Index: /branches/rel1/psLib/src/dataManip/Makefile =================================================================== --- /branches/rel1/psLib/src/dataManip/Makefile (revision 1039) +++ /branches/rel1/psLib/src/dataManip/Makefile (revision 1040) @@ -9,5 +9,4 @@ SRC_OBJS = psStats.o \ psFunctions.o \ - psMinimize.o \ psImageStats.o \ psMatrix.o \ Index: anches/rel1/psLib/src/dataManip/psMinimize.c =================================================================== --- /branches/rel1/psLib/src/dataManip/psMinimize.c (revision 1039) +++ (revision ) @@ -1,304 +1,0 @@ -#include -#include -#include -#include -#include -#include - -#include "psMemory.h" -#include "psVector.h" -#include "psTrace.h" -#include "psError.h" -#include "psAbort.h" -#include "psFunctions.h" -#include "psSort.h" -#include "float.h" -#include -// NOTE: rewrite so there is no maximum order for the polynomials. -#define MAX_POLY_ORDER 10 -#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) -int MyInfoLevel = 0; -/****************************************************************************** - This routine must minimize an arbitrary function. - *****************************************************************************/ -psVector * -psMinimize(float (*myFunction)(const psVector *restrict), - psVector *restrict initialGuess, - psVector *restrict paramMask) -{ - return(NULL); -} - -/****************************************************************************** - This routine must minimize an arbitrary function. - *****************************************************************************/ -psVector * -psMinimizeChi2(float (*evalModel)(const psVector *restrict, - const psVector *restrict), - const psVector *restrict domain, - const psVector *restrict data, - const psVector *restrict errors, - psVector *restrict initialGuess, - const psVector *restrict paramMask) -{ - return(NULL); -} - -/** @brief This procedure calculates various combinations of powers of x and y - * and stores them in the data structure sums[][]. After it completes: - * sums[i][j] == x^i * y^j - */ -void buildSums(double x, - double y, - /*@out@*/double sums[MAX_POLY_ORDER+1][MAX_POLY_ORDER+1], - int polyOrder) -{ - int i = 0; // loop index variable - int j = 0; // loop index variable - double xSum = 0.0; // The running sum of X terms - double ySum = 0.0; // The running sum of Y terms - - xSum = 1.0; - ySum = 1.0; - for(i=0;i<=polyOrder;i++) { - ySum = xSum; - for(j=0;j<=polyOrder;j++) { - sums[i][j] = ySum; - ySum*= y; - } - xSum*= x; - } -} - -/** @brief The coefficients of the matrix in equation (7) from the ADD will - * be very large if the x and y values are in the 0-511 range (ie: the sum y^7 - * for all 0 - *
  • height - *
  • width - * - * - * Output: - *
      - *
    • imageScalingFactors - *
    - * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ -void buildImageScalingFactors(int height, - int width, - float **imageScalingFactors) -{ - int maxDim = 0; // The largest dimension of the image. - int i = 0; // loop index variable. - - // Calculate the maximum dimensional extent of the image. - if (height > width) { - maxDim = height; - } else { - maxDim = width; - } - - - // Allocate memory for the output array. - *imageScalingFactors = (float *) psAlloc((maxDim+10) * sizeof(float)); - - // This code is somewhat arbitrary. For an image with a height/width - // of 512x512, the scaling factors will be between 0.0-1.0. - for (i=0;i. - * polyTerms[poly][i][1] = the power to which Y is raised in the - * i-th term of in an poly-order sky - * background polynomial

    . - * - * NOTE: the C_0 term defined in the ADD begins at i=2 in our data - * structures (ie. the x/y powers of the i-th term in the sky model - * polynomial are actually stored at polyTerms[][i+2][]. There are two - * reasons for this. First, there is a term prior to C_0 in equation - * (7) of the ADD. Second, our linear algebra codes assume data is - * stored offset from index 1. - * - * Input: - *
      - *
    • polyTerms[][][] - *
    - * - * Output: - *
      - *
    • polyTerms[][][] - *
    - * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ -void buildPolyTerms(/*@out@*/ int polyTerms[MAX_POLY_ORDER+1][(MAX_POLYNOMIAL_TERMS+2)][2]) -{ - int polyOrder=0; // loop index variable. - int i=0; // loop index variable. - int term = 0; // loop index variable. - int num=0; // loop index variable. - - for(polyOrder=0;polyOrder<=MAX_POLY_ORDER;polyOrder++) { - // The following 4 terms should not be used in any of the subsequent - // computation. We initialize them to zero in order to produce stable - // results for debugging purposes should they mistakenly be used. - polyTerms[polyOrder][0][0] = 0; - polyTerms[polyOrder][0][1] = 0; - polyTerms[polyOrder][1][0] = 0; - polyTerms[polyOrder][1][1] = 0; - - // This code segment loops through each term i in the polynomial and - // calculates the power to which x/y are raised in that i-th term. - i=2; - for (term=0;term<=polyOrder;term++) { - for (num=0;num<=term;num++) { - polyTerms[polyOrder][i][0] = term-num; - polyTerms[polyOrder][i][1] = num; - if (MyInfoLevel > 2) { - printf("%d-th order Sky polynomial term %d is x^%d y^%d\n", - polyOrder, i, - polyTerms[polyOrder][i][0], polyTerms[polyOrder][i][1]); - } - i++; - } - } - } -} - - -/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th - * order sky background polynomial defined by the coefficients in the array B[] - * are consistent with zero. If true, then *flag is set to 1. Otherwise, - * *flag is set to 0. The matrix inversion code in the middle of this - * procedure draws from Numerical Recipes in C page 48. - * - * Input: - *
      - *
    • A This is the LUD decomposition of the original matrix A. - *
    • N The size of the matrix (plus 1, actually, since offset 1). - *
    • indx misc Numerical Recipes data structure. - *
    • B The coefficients of the sky polynomial. - *
    • polyOrder The degree of the sky polynomial. - *
    - * Output: - *
      - *
    • *flag Set this to 1 if we must recalculate the coefficients. - *
    - * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ -void polyOrderCheck(float **A, - int N, - int *indx, - float *B, - int polyOrder, - int *flag) -{ - float **y = NULL; // This 2-D matrix will hold A^-1 - float *col = NULL; // misc NumerRecipes data structure - float *error=NULL; // will hold the sqrt() of the - // diagonal of y[][]. - int i=0; // loop-index variable - int j=0; // loop-index variable - int numPolyTerms = 0; // The number of terms in the - // polynomial. - int lastTerm = 0; // The index location of the first - // n-th order term in array B[]. - int firstTerm = 0; // Index location of last such term. - - // Allocate the necessary data structures for this procedure... - error = (float *) psAlloc((N + 1) * sizeof(float)); - col = (float *) psAlloc((N + 1) * sizeof(float)); - y = (float **) psAlloc((N + 1) * sizeof(float *)); - for(i=1;i<=N;i++) { - y[i] = (float *) psAlloc((N + 1) * sizeof(float)); - } - - // Invert the matrix A and put the result in y[][]. This code is taken - // from Numerical Recipes in C page 48. - for(j=1;j<=N;j++) { - for(i=1;i<=N;i++) { - col[i] = 0.0; - } - col[j] = 1.0; - // NOTE: substitue the LUD rotine - // lubksb(A, N, indx, col); - for(i=1;i<=N;i++) { - y[i][j] = col[i]; - } - } - - // Determine where the first n-th order (in this comment, n equals - // polyOrder) polynomial term is stored in the matrix B[], and also were - // the last n-order term is stored. Then we loop over all the n-order - // terms and check if they are consistent with zero. - - numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); - lastTerm = numPolyTerms + 1; - firstTerm = lastTerm - polyOrder; - *flag = 1; - for (i=firstTerm; i<=lastTerm; i++) { - #ifdef DARWIN - error[i] = (float)sqrt(y[i][i]); - #else - - error[i] = sqrtf(y[i][i]); - #endif - - if (!((B[i] <= (2.0f * error[i])) && - ((-2.0f * error[i]) <= B[i]))) { - *flag = 0; - } - } - - // Free all memory allocated in this routine. - psFree(error); - psFree(col); - for(j=1;j<=N;j++) { - psFree(y[j]); - } - psFree(y); -} - - - - - - - - - -/****************************************************************************** - This routine must fit a polynomial of degree myPoly to the data points - (x, y) and return the coefficients of that polynomial, as well as the - error for each data poiny (yErr). - *****************************************************************************/ -psPolynomial1D * -psGetArrayPolynomial(psPolynomial1D myPoly, - const psVector *restrict x, - const psVector *restrict y, - const psVector *restrict yErr) -{ - return(NULL); -} Index: anches/rel1/psLib/src/dataManip/psMinimize.d =================================================================== --- /branches/rel1/psLib/src/dataManip/psMinimize.d (revision 1039) +++ (revision ) @@ -1,4 +1,0 @@ -psMinimize.o psMinimize.d : psMinimize.c ../sysUtils/psMemory.h \ - ../collections/psVector.h ../collections/psType.h ../sysUtils/psTrace.h \ - ../sysUtils/psError.h ../sysUtils/psAbort.h psFunctions.h \ - ../collections/psSort.h Index: anches/rel1/psLib/src/dataManip/psMinimize.h =================================================================== --- /branches/rel1/psLib/src/dataManip/psMinimize.h (revision 1039) +++ (revision ) @@ -1,43 +1,0 @@ -#if !defined(PS_MINIMIZE_H) -#define PS_MINIMIZE_H - -/** \file psMinimize.h - * \brief minimization operations - * \ingroup Stats - */ - -/** Functions **************************************************************/ -/** \addtogroup Stats - * \{ - */ - -/** This routine must ninimize a particular non-linear function */ -psFloatArray * -psMinimize(float (*myFunction)(const psFloatArray *restrict), ///< Function to minimize - psFloatArray *restrict initialGuess, ///< Initial guess - psIntArray *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant - ); - - -/** Minimize chi^2 for input data */ -psFloatArray * -psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict, - const psFloatArray *restrict), ///< Model to fit; (domain and params) - const psFloatArray *restrict domain, ///< The domain values for the corresponding measurements - const psFloatArray *restrict data, ///< Data to fit - const psFloatArray *restrict errors, ///< Errors in the data - psFloatArray *restrict initialGuess, ///< Initial guess - const psIntArray *restrict paramMask ///< 1 = fit for parameter, 0 = hold parameter constant - ); - -/** Derive a polynomial fit by chi^2 minimisation --- can be done analytically */ -psPolynomial1D * -psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit - const psFloatArray *restrict x, ///< Ordinates (or NULL to just use the indices) - const psFloatArray *restrict y, ///< Coordinates - const psFloatArray *restrict yErr ///< Errors in coordinates, or NULL - ); - -/* \} */ // End of MathGroup Functions - -#endif Index: anches/rel1/psLib/src/math/psMinimize.c =================================================================== --- /branches/rel1/psLib/src/math/psMinimize.c (revision 1039) +++ (revision ) @@ -1,304 +1,0 @@ -#include -#include -#include -#include -#include -#include - -#include "psMemory.h" -#include "psVector.h" -#include "psTrace.h" -#include "psError.h" -#include "psAbort.h" -#include "psFunctions.h" -#include "psSort.h" -#include "float.h" -#include -// NOTE: rewrite so there is no maximum order for the polynomials. -#define MAX_POLY_ORDER 10 -#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) -int MyInfoLevel = 0; -/****************************************************************************** - This routine must minimize an arbitrary function. - *****************************************************************************/ -psVector * -psMinimize(float (*myFunction)(const psVector *restrict), - psVector *restrict initialGuess, - psVector *restrict paramMask) -{ - return(NULL); -} - -/****************************************************************************** - This routine must minimize an arbitrary function. - *****************************************************************************/ -psVector * -psMinimizeChi2(float (*evalModel)(const psVector *restrict, - const psVector *restrict), - const psVector *restrict domain, - const psVector *restrict data, - const psVector *restrict errors, - psVector *restrict initialGuess, - const psVector *restrict paramMask) -{ - return(NULL); -} - -/** @brief This procedure calculates various combinations of powers of x and y - * and stores them in the data structure sums[][]. After it completes: - * sums[i][j] == x^i * y^j - */ -void buildSums(double x, - double y, - /*@out@*/double sums[MAX_POLY_ORDER+1][MAX_POLY_ORDER+1], - int polyOrder) -{ - int i = 0; // loop index variable - int j = 0; // loop index variable - double xSum = 0.0; // The running sum of X terms - double ySum = 0.0; // The running sum of Y terms - - xSum = 1.0; - ySum = 1.0; - for(i=0;i<=polyOrder;i++) { - ySum = xSum; - for(j=0;j<=polyOrder;j++) { - sums[i][j] = ySum; - ySum*= y; - } - xSum*= x; - } -} - -/** @brief The coefficients of the matrix in equation (7) from the ADD will - * be very large if the x and y values are in the 0-511 range (ie: the sum y^7 - * for all 0 - *
  • height - *
  • width - * - * - * Output: - *
      - *
    • imageScalingFactors - *
    - * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ -void buildImageScalingFactors(int height, - int width, - float **imageScalingFactors) -{ - int maxDim = 0; // The largest dimension of the image. - int i = 0; // loop index variable. - - // Calculate the maximum dimensional extent of the image. - if (height > width) { - maxDim = height; - } else { - maxDim = width; - } - - - // Allocate memory for the output array. - *imageScalingFactors = (float *) psAlloc((maxDim+10) * sizeof(float)); - - // This code is somewhat arbitrary. For an image with a height/width - // of 512x512, the scaling factors will be between 0.0-1.0. - for (i=0;i. - * polyTerms[poly][i][1] = the power to which Y is raised in the - * i-th term of in an poly-order sky - * background polynomial

    . - * - * NOTE: the C_0 term defined in the ADD begins at i=2 in our data - * structures (ie. the x/y powers of the i-th term in the sky model - * polynomial are actually stored at polyTerms[][i+2][]. There are two - * reasons for this. First, there is a term prior to C_0 in equation - * (7) of the ADD. Second, our linear algebra codes assume data is - * stored offset from index 1. - * - * Input: - *
      - *
    • polyTerms[][][] - *
    - * - * Output: - *
      - *
    • polyTerms[][][] - *
    - * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ -void buildPolyTerms(/*@out@*/ int polyTerms[MAX_POLY_ORDER+1][(MAX_POLYNOMIAL_TERMS+2)][2]) -{ - int polyOrder=0; // loop index variable. - int i=0; // loop index variable. - int term = 0; // loop index variable. - int num=0; // loop index variable. - - for(polyOrder=0;polyOrder<=MAX_POLY_ORDER;polyOrder++) { - // The following 4 terms should not be used in any of the subsequent - // computation. We initialize them to zero in order to produce stable - // results for debugging purposes should they mistakenly be used. - polyTerms[polyOrder][0][0] = 0; - polyTerms[polyOrder][0][1] = 0; - polyTerms[polyOrder][1][0] = 0; - polyTerms[polyOrder][1][1] = 0; - - // This code segment loops through each term i in the polynomial and - // calculates the power to which x/y are raised in that i-th term. - i=2; - for (term=0;term<=polyOrder;term++) { - for (num=0;num<=term;num++) { - polyTerms[polyOrder][i][0] = term-num; - polyTerms[polyOrder][i][1] = num; - if (MyInfoLevel > 2) { - printf("%d-th order Sky polynomial term %d is x^%d y^%d\n", - polyOrder, i, - polyTerms[polyOrder][i][0], polyTerms[polyOrder][i][1]); - } - i++; - } - } - } -} - - -/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th - * order sky background polynomial defined by the coefficients in the array B[] - * are consistent with zero. If true, then *flag is set to 1. Otherwise, - * *flag is set to 0. The matrix inversion code in the middle of this - * procedure draws from Numerical Recipes in C page 48. - * - * Input: - *
      - *
    • A This is the LUD decomposition of the original matrix A. - *
    • N The size of the matrix (plus 1, actually, since offset 1). - *
    • indx misc Numerical Recipes data structure. - *
    • B The coefficients of the sky polynomial. - *
    • polyOrder The degree of the sky polynomial. - *
    - * Output: - *
      - *
    • *flag Set this to 1 if we must recalculate the coefficients. - *
    - * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ -void polyOrderCheck(float **A, - int N, - int *indx, - float *B, - int polyOrder, - int *flag) -{ - float **y = NULL; // This 2-D matrix will hold A^-1 - float *col = NULL; // misc NumerRecipes data structure - float *error=NULL; // will hold the sqrt() of the - // diagonal of y[][]. - int i=0; // loop-index variable - int j=0; // loop-index variable - int numPolyTerms = 0; // The number of terms in the - // polynomial. - int lastTerm = 0; // The index location of the first - // n-th order term in array B[]. - int firstTerm = 0; // Index location of last such term. - - // Allocate the necessary data structures for this procedure... - error = (float *) psAlloc((N + 1) * sizeof(float)); - col = (float *) psAlloc((N + 1) * sizeof(float)); - y = (float **) psAlloc((N + 1) * sizeof(float *)); - for(i=1;i<=N;i++) { - y[i] = (float *) psAlloc((N + 1) * sizeof(float)); - } - - // Invert the matrix A and put the result in y[][]. This code is taken - // from Numerical Recipes in C page 48. - for(j=1;j<=N;j++) { - for(i=1;i<=N;i++) { - col[i] = 0.0; - } - col[j] = 1.0; - // NOTE: substitue the LUD rotine - // lubksb(A, N, indx, col); - for(i=1;i<=N;i++) { - y[i][j] = col[i]; - } - } - - // Determine where the first n-th order (in this comment, n equals - // polyOrder) polynomial term is stored in the matrix B[], and also were - // the last n-order term is stored. Then we loop over all the n-order - // terms and check if they are consistent with zero. - - numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); - lastTerm = numPolyTerms + 1; - firstTerm = lastTerm - polyOrder; - *flag = 1; - for (i=firstTerm; i<=lastTerm; i++) { - #ifdef DARWIN - error[i] = (float)sqrt(y[i][i]); - #else - - error[i] = sqrtf(y[i][i]); - #endif - - if (!((B[i] <= (2.0f * error[i])) && - ((-2.0f * error[i]) <= B[i]))) { - *flag = 0; - } - } - - // Free all memory allocated in this routine. - psFree(error); - psFree(col); - for(j=1;j<=N;j++) { - psFree(y[j]); - } - psFree(y); -} - - - - - - - - - -/****************************************************************************** - This routine must fit a polynomial of degree myPoly to the data points - (x, y) and return the coefficients of that polynomial, as well as the - error for each data poiny (yErr). - *****************************************************************************/ -psPolynomial1D * -psGetArrayPolynomial(psPolynomial1D myPoly, - const psVector *restrict x, - const psVector *restrict y, - const psVector *restrict yErr) -{ - return(NULL); -} Index: anches/rel1/psLib/src/math/psMinimize.h =================================================================== --- /branches/rel1/psLib/src/math/psMinimize.h (revision 1039) +++ (revision ) @@ -1,43 +1,0 @@ -#if !defined(PS_MINIMIZE_H) -#define PS_MINIMIZE_H - -/** \file psMinimize.h - * \brief minimization operations - * \ingroup Stats - */ - -/** Functions **************************************************************/ -/** \addtogroup Stats - * \{ - */ - -/** This routine must ninimize a particular non-linear function */ -psFloatArray * -psMinimize(float (*myFunction)(const psFloatArray *restrict), ///< Function to minimize - psFloatArray *restrict initialGuess, ///< Initial guess - psIntArray *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant - ); - - -/** Minimize chi^2 for input data */ -psFloatArray * -psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict, - const psFloatArray *restrict), ///< Model to fit; (domain and params) - const psFloatArray *restrict domain, ///< The domain values for the corresponding measurements - const psFloatArray *restrict data, ///< Data to fit - const psFloatArray *restrict errors, ///< Errors in the data - psFloatArray *restrict initialGuess, ///< Initial guess - const psIntArray *restrict paramMask ///< 1 = fit for parameter, 0 = hold parameter constant - ); - -/** Derive a polynomial fit by chi^2 minimisation --- can be done analytically */ -psPolynomial1D * -psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit - const psFloatArray *restrict x, ///< Ordinates (or NULL to just use the indices) - const psFloatArray *restrict y, ///< Coordinates - const psFloatArray *restrict yErr ///< Errors in coordinates, or NULL - ); - -/* \} */ // End of MathGroup Functions - -#endif