IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1103


Ignore:
Timestamp:
Jun 25, 2004, 2:24:30 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/collections/psVector.c

    r1073 r1103  
    88 *  @author Ross Harman, MHPCC
    99 *
    10  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-06-23 23:00:15 $
     10 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-06-26 00:24:30 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    7474    psVec->type.type = elemType;
    7575    psVec->nalloc = nalloc;
    76     psVec->n = 0;
     76    psVec->n = nalloc;
    7777
    7878    // Create vector data array

  • trunk/psLib/src/dataManip/psMinimize.c

    r1093 r1103  
    3737}
    3838psModelData;
     39
    3940// The first argument to evalModel() and
    4041// d_evalModel() specifies the data point.
    … …  
    7778                   gsl_vector *f)
    7879{
     80    return 0;
    7981    int i;    // Loop index variable.
    8082    int j;    // Loop index variable.
    … …  
    148150    // have those parameters removed.  Here will create a new parameter list
    149151    // with the masked parameters added (we expand initialGuess).
     152
     153    psMemCheckCorruption(1);
     154    printf("Calling psVectorAlloc(%d)\n", initialGuess->n);
     155
    150156    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     157    printf("Called psVectorAlloc(%d)\n", initialGuess->n);
     158
    151159    if (mask != NULL) {
    152160        j = 0;
    … …  
    208216
    209217    int numData = domain->numCols; // Number of data points
    210     gsl_multifit_fdfsolver *s = NULL; // GSL data structure.
    211218    int status;    // Return status for the GSL solver.
    212219    int i = 0;    // Loop index variable.
    213220    int j = 0;    // Loop index variable.
    214221    int iter = 0;   // Iteration counter.
    215     psModelData *inputData = NULL; // Contains data that the user-supplied
    216     // function/derivate need.
    217222    gsl_multifit_function_fdf f; // GSL structure that contains the
    218223    // functions/derivative to be solved.
    219     double *xInit = NULL;  // The initial guess at the parameters
     224    double *xInit = NULL;     // The initial guess at the parameters
    220225    // with masked parameters removed.
    221226    const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
    … …  
    223228    // Marquardt algorithm for chi2
    224229    // minimization.
    225 
    226     inputData = (psModelData *) psAlloc(sizeof(psModelData));
    227     inputData->n = numData;
    228     inputData->count = 0;
    229     inputData->initialGuess = initialGuess;
    230     inputData->domain = domain;
    231     //    inputData->data = data;
    232     //    inputData->errors = errors;
    233     //    inputData->paramMask = paramMask;
    234     //    inputData->evalModel = evalModel;
    235     //    inputData->d_evalModel = DevalModel;
    236 
    237     inputData->count = 0;
     230    gsl_multifit_fdfsolver *s; // GSL data structure.
     231
     232    //    inputData = (psModelData *) psAlloc(sizeof(psModelData));
     233    psModelData inputData;
     234    inputData.n = numData;
     235    inputData.count = 0;
     236    inputData.initialGuess = initialGuess;
     237    inputData.domain = domain;
     238    inputData.count = 0;
    238239    if (paramMask != NULL) {
    239240        for (i=0;i<paramMask->n;i++) {
    240241            if (paramMask->data.U8[i] != 0) {
    241                 inputData->count++;
     242                inputData.count++;
    242243            }
    243244        }
    244245    } else {
    245         inputData->count= initialGuess->n;
     246        inputData.count= initialGuess->n;
    246247    }
    247248
    … …  
    250251    // parameters are masked out.
    251252
    252     xInit = (double *) psAlloc(inputData->count * sizeof(double));
     253    xInit = (double *) psAlloc(inputData.count * sizeof(double));
    253254    if (paramMask != NULL) {
    254255        j = 0;
    255256        for (i=0;i<initialGuess->n;i++) {
    256257            if (paramMask->data.U8[i] == 0) {
    257                 xInit[j++] = initialGuess->data.F32[i];
     258                xInit[j++] = (double) initialGuess->data.F32[i];
    258259            }
    259260        }
    260261    } else {
    261262        for (i=0;i<initialGuess->n;i++) {
    262             xInit[i] = initialGuess->data.F32[i];
    263         }
    264     }
    265 
    266     // Creates the vector for x which GSL uses.  Must deallocate.
    267     gsl_vector_view x = gsl_vector_view_array(xInit, inputData->count);
     263            xInit[i] = (double) initialGuess->data.F32[i];
     264        }
     265    }
    268266
    269267    const gsl_rng_type *type;
    … …  
    283281    f.fdf = &gsl_my_function_fdf;
    284282    f.n = numData;
    285     f.p = inputData->count;
     283    f.p = inputData.count;
    286284    f.params = &inputData;
    287285
     286    printf("inputData.count is %d\n", inputData.count);
     287    printf("numData is %d\n", numData);
     288    // Creates the vector for x which GSL uses.  Must deallocate.
     289    gsl_vector_view x = gsl_vector_view_array(xInit, inputData.count);
     290
    288291    // Create an instance of the GSL solver that we will be iterating on.
    289     // It will have numData data points and inputData->count parameters.
    290     s = gsl_multifit_fdfsolver_alloc(T, numData, inputData->count);
    291 
     292    // It will have numData data points and inputData.count parameters.
     293    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.count);
     294
     295    printf("HERE 04\n");
    292296    // Initialize the GSL minimizer to use function defined by the data
    293297    // structure "f" and x.vector as an initial guess for the parameters.
    … …  
    295299
    296300
     301    printf("HERE 05\n");
    297302    // Each iteration of the following loop will perform one step in an
    298303    // attempt to minimized chi-squared for the function.  The loop exits
    … …  
    302307        iter++;
    303308        // Perform an iteration of the GSL solver.
     309        printf("HERE 06\n");
    304310        status = gsl_multifit_fdfsolver_iterate(s);
    305 
     311        printf("HERE 07\n");
    306312        // If there was a problem, abort.
    307313        if (status) {
     314            exit(1);
    308315            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
    309316        }
     317        printf("HERE 08\n");
    310318
    311319        // Test if the parameters changed by a small enough amount.
    312320        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     321        printf("HERE 09\n");
    313322    } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
     323    printf("HERE 12\n");
    314324
    315325
    … …  
    339349    gsl_multifit_fdfsolver_free(s);
    340350    psFree(xInit);
    341     psFree(inputData);
    342351
    343352    // Bye bye.
    … …  
    621630    return(NULL);
    622631}
     632
     633
     634
     635
     636
     637

  • trunk/psLib/src/dataManip/psMinimize.d

    r986 r1103  
    11psMinimize.o psMinimize.d : psMinimize.c ../sysUtils/psMemory.h \
    2   ../collections/psVector.h ../collections/psType.h ../sysUtils/psTrace.h \
    3   ../sysUtils/psError.h ../sysUtils/psAbort.h psFunctions.h \
    4   ../collections/psSort.h
     2  ../collections/psVector.h ../collections/psType.h \
     3  ../collections/psImage.h ../sysUtils/psTrace.h ../sysUtils/psError.h \
     4  ../sysUtils/psAbort.h psFunctions.h ../collections/psSort.h \
     5  psMinimize.h

  • trunk/psLib/src/math/psMinimize.c

    r1093 r1103  
    3737}
    3838psModelData;
     39
    3940// The first argument to evalModel() and
    4041// d_evalModel() specifies the data point.
    … …  
    7778                   gsl_vector *f)
    7879{
     80    return 0;
    7981    int i;    // Loop index variable.
    8082    int j;    // Loop index variable.
    … …  
    148150    // have those parameters removed.  Here will create a new parameter list
    149151    // with the masked parameters added (we expand initialGuess).
     152
     153    psMemCheckCorruption(1);
     154    printf("Calling psVectorAlloc(%d)\n", initialGuess->n);
     155
    150156    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     157    printf("Called psVectorAlloc(%d)\n", initialGuess->n);
     158
    151159    if (mask != NULL) {
    152160        j = 0;
    … …  
    208216
    209217    int numData = domain->numCols; // Number of data points
    210     gsl_multifit_fdfsolver *s = NULL; // GSL data structure.
    211218    int status;    // Return status for the GSL solver.
    212219    int i = 0;    // Loop index variable.
    213220    int j = 0;    // Loop index variable.
    214221    int iter = 0;   // Iteration counter.
    215     psModelData *inputData = NULL; // Contains data that the user-supplied
    216     // function/derivate need.
    217222    gsl_multifit_function_fdf f; // GSL structure that contains the
    218223    // functions/derivative to be solved.
    219     double *xInit = NULL;  // The initial guess at the parameters
     224    double *xInit = NULL;     // The initial guess at the parameters
    220225    // with masked parameters removed.
    221226    const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
    … …  
    223228    // Marquardt algorithm for chi2
    224229    // minimization.
    225 
    226     inputData = (psModelData *) psAlloc(sizeof(psModelData));
    227     inputData->n = numData;
    228     inputData->count = 0;
    229     inputData->initialGuess = initialGuess;
    230     inputData->domain = domain;
    231     //    inputData->data = data;
    232     //    inputData->errors = errors;
    233     //    inputData->paramMask = paramMask;
    234     //    inputData->evalModel = evalModel;
    235     //    inputData->d_evalModel = DevalModel;
    236 
    237     inputData->count = 0;
     230    gsl_multifit_fdfsolver *s; // GSL data structure.
     231
     232    //    inputData = (psModelData *) psAlloc(sizeof(psModelData));
     233    psModelData inputData;
     234    inputData.n = numData;
     235    inputData.count = 0;
     236    inputData.initialGuess = initialGuess;
     237    inputData.domain = domain;
     238    inputData.count = 0;
    238239    if (paramMask != NULL) {
    239240        for (i=0;i<paramMask->n;i++) {
    240241            if (paramMask->data.U8[i] != 0) {
    241                 inputData->count++;
     242                inputData.count++;
    242243            }
    243244        }
    244245    } else {
    245         inputData->count= initialGuess->n;
     246        inputData.count= initialGuess->n;
    246247    }
    247248
    … …  
    250251    // parameters are masked out.
    251252
    252     xInit = (double *) psAlloc(inputData->count * sizeof(double));
     253    xInit = (double *) psAlloc(inputData.count * sizeof(double));
    253254    if (paramMask != NULL) {
    254255        j = 0;
    255256        for (i=0;i<initialGuess->n;i++) {
    256257            if (paramMask->data.U8[i] == 0) {
    257                 xInit[j++] = initialGuess->data.F32[i];
     258                xInit[j++] = (double) initialGuess->data.F32[i];
    258259            }
    259260        }
    260261    } else {
    261262        for (i=0;i<initialGuess->n;i++) {
    262             xInit[i] = initialGuess->data.F32[i];
    263         }
    264     }
    265 
    266     // Creates the vector for x which GSL uses.  Must deallocate.
    267     gsl_vector_view x = gsl_vector_view_array(xInit, inputData->count);
     263            xInit[i] = (double) initialGuess->data.F32[i];
     264        }
     265    }
    268266
    269267    const gsl_rng_type *type;
    … …  
    283281    f.fdf = &gsl_my_function_fdf;
    284282    f.n = numData;
    285     f.p = inputData->count;
     283    f.p = inputData.count;
    286284    f.params = &inputData;
    287285
     286    printf("inputData.count is %d\n", inputData.count);
     287    printf("numData is %d\n", numData);
     288    // Creates the vector for x which GSL uses.  Must deallocate.
     289    gsl_vector_view x = gsl_vector_view_array(xInit, inputData.count);
     290
    288291    // Create an instance of the GSL solver that we will be iterating on.
    289     // It will have numData data points and inputData->count parameters.
    290     s = gsl_multifit_fdfsolver_alloc(T, numData, inputData->count);
    291 
     292    // It will have numData data points and inputData.count parameters.
     293    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.count);
     294
     295    printf("HERE 04\n");
    292296    // Initialize the GSL minimizer to use function defined by the data
    293297    // structure "f" and x.vector as an initial guess for the parameters.
    … …  
    295299
    296300
     301    printf("HERE 05\n");
    297302    // Each iteration of the following loop will perform one step in an
    298303    // attempt to minimized chi-squared for the function.  The loop exits
    … …  
    302307        iter++;
    303308        // Perform an iteration of the GSL solver.
     309        printf("HERE 06\n");
    304310        status = gsl_multifit_fdfsolver_iterate(s);
    305 
     311        printf("HERE 07\n");
    306312        // If there was a problem, abort.
    307313        if (status) {
     314            exit(1);
    308315            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
    309316        }
     317        printf("HERE 08\n");
    310318
    311319        // Test if the parameters changed by a small enough amount.
    312320        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     321        printf("HERE 09\n");
    313322    } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
     323    printf("HERE 12\n");
    314324
    315325
    … …  
    339349    gsl_multifit_fdfsolver_free(s);
    340350    psFree(xInit);
    341     psFree(inputData);
    342351
    343352    // Bye bye.
    … …  
    621630    return(NULL);
    622631}
     632
     633
     634
     635
     636
     637

  • trunk/psLib/src/mathtypes/psVector.c

    r1073 r1103  
    88 *  @author Ross Harman, MHPCC
    99 *
    10  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-06-23 23:00:15 $
     10 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-06-26 00:24:30 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    7474    psVec->type.type = elemType;
    7575    psVec->nalloc = nalloc;
    76     psVec->n = 0;
     76    psVec->n = nalloc;
    7777
    7878    // Create vector data array

  • trunk/psLib/test/dataManip/tst_psMinimize00.c

    r1095 r1103  
    1515#define NUM_PARAMS 4
    1616
    17 float *myFunc(psVector *myData,
    18               psVector *myParams)
     17float myFunc(const psVector *restrict myData,
     18             const psVector *restrict myParams)
    1919{
    2020    float x = myData->data.F32[0];
    … …  
    2424    float C = myParams->data.F32[2];
    2525    float D = myParams->data.F32[3];
    26     float *tmp = (float *) psAlloc(sizeof(float));
    2726
    28     *tmp = A + (B*x) + (C*y) + (D*x*y);
    29     return(tmp);
     27    return(A + (B*x) + (C*y) + (D*x*y));
    3028}
    3129
    32 float *myFuncDeriv(psVector *myData,
    33                    psVector *myParams,
    34                    int whichParamDeriv)
     30float myFuncDeriv(const psVector *restrict myData,
     31                  const psVector *restrict myParams,
     32                  int whichParamDeriv)
    3533{
    3634    float x = myData->data.F32[0];
    3735    float y = myData->data.F32[1];
    38     float *tmp = (float *) psAlloc(sizeof(float));
     36    float tmp = 0.0;
    3937
    4038    if (whichParamDeriv == 0) {
    41         *tmp = 1.0;
     39        tmp = 1.0;
    4240    } else if (whichParamDeriv == 1) {
    43         *tmp = x;
     41        tmp = x;
    4442    } else if (whichParamDeriv == 2) {
    45         *tmp = y;
     43        tmp = y;
    4644    } else if (whichParamDeriv == 3) {
    47         *tmp = x * y;
     45        tmp = x * y;
    4846    }
    4947    return(tmp);
    … …  
    6664    int memLeaks = 0;
    6765
    68     domain = psImageAlloc(NUM_DATA, DATA_WIDTH, PS_TYPE_F32);
     66    domain = psImageAlloc(DATA_WIDTH, NUM_DATA, PS_TYPE_F32);
    6967    data = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    7068    errors = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    … …  
    7371    tmpVecPtr = psVectorAlloc(DATA_WIDTH, PS_TYPE_F32);
    7472
     73    printf("FUCK: initialGuess->n is %d\n", initialGuess->n);
    7574    // Build the data.
    7675    initialGuess->data.F32[0] = 2.0;
    … …  
    9190            tmpVecPtr->data.F32[i] = domain->data.F32[i][j];
    9291        }
    93         data->data.F32[i] = *myFunc(tmpVecPtr, initialGuess);
     92        data->data.F32[i] = myFunc(tmpVecPtr, initialGuess);
    9493    }
    9594
    … …  
    104103                   errors,
    105104                   initialGuess,
    106                    paramMask,
     105                   NULL,
    107106                   &chiSq);
    108107
Note: See TracChangeset for help on using the changeset viewer.