IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psModules

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1212pmSourceIO_CMF_PS1_V1.v1.c
        1313pmSourceIO_CMF_PS1_V4.c
         14pmSourceIO_CMF_PS1_V5.c
        1415pmSourceIO_CMF_PS1_SV1.c
        1516pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
        1618pmSourceIO_CMF_PS1_DV1.c
        1719pmSourceIO_CMF_PS1_DV2.c
        1820pmSourceIO_CMF_PS1_DV3.c
        19 
         21pmSourceIO_CMF_PS1_DV4.c

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects/pmPCM_MinimizeChisq.c

    r35768 r37066  
    3131#include "pmMoments.h"
    3232#include "pmModelFuncs.h"
     33#include "pmModelClass.h"
    3334#include "pmModel.h"
    3435#include "pmModelUtils.h"
    35 #include "pmModelClass.h"
    3636#include "pmSourceMasks.h"
    3737#include "pmSourceExtendedPars.h"
    3838#include "pmSourceDiffStats.h"
    3939#include "pmSourceSatstar.h"
     40#include "pmSourceLensing.h"
    4041#include "pmSource.h"
    4142#include "pmSourceFitModel.h"
    4243#include "pmPCMdata.h"
     44
     45# define SAVE_IMAGES 0
     46# if (SAVE_IMAGES)
     47int psphotSaveImage (psMetadata *header, psImage *image, char *filename);
     48# endif
    4349
    4450# define FACILITY "psModules.objects"
    … …  
    9197    psF32 lambda = 0.001;
    9298    psF32 dLinear = 0.0;
    93     psF32 nu = 2.0;
     99    psF32 nu = 3.0;
    94100
    95101# if (USE_FFT && PRE_CONVOLVE)
    … …  
    130136        }
    131137
     138        if (min->isInteractive) {
     139            fprintf (stderr, "%d : ", min->iter);
     140            for (int ti = 0; ti < params->n; ti++) {
     141                fprintf (stderr, "%f  ", params->data.F32[ti]);
     142            }
     143            fprintf (stderr, " : %f\n", min->value);
     144        }
     145
     146        char key[10]; // used for interactive responses
     147        bool testValue = false;
     148
    132149        // set a new guess for Alpha, Beta, Params
    133150        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
     151            if (false && min->isInteractive) {
     152                fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] ");
     153                if (!fgets(key, 8, stdin)) {
     154                    psWarning("Unable to read option");
     155                }
     156                switch (key[0]) {
     157                  case 'n':
     158                  case 'N':
     159                    done = true;
     160                    break;
     161                  case 'y':
     162                  case 'Y':
     163                  case '\n':
     164                    lambda *= 10.0;
     165                    continue;
     166                  default:
     167                    lambda *= 10.0;
     168                    continue;
     169                }
     170                if (done) break;
     171            }
    134172            min->iter ++;
    135173            if (min->iter >=  min->maxIter) break;
    … …  
    138176        }
    139177
     178        if (false && min->isInteractive) {
     179            p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: ");
     180            fprintf (stdout, "last chisq : %f\n", min->value);
     181            bool getOptions = true;
     182            while (getOptions) {
     183                fprintf (stdout, "options: (m)odify, (g)o, (q)uit: ");
     184                if (!fgets(key, 8, stdin)) {
     185                    psWarning("Unable to read option");
     186                }
     187                switch (key[0]) {
     188                  case 'm':
     189                  case 'M':
     190                    testValue = TRUE;
     191                    fprintf (stdout, "enter (Npar) (value): ");
     192                    int Npar = 0;
     193                    float value= 0;
     194                    int Nscan = fscanf (stdin, "%d %f", &Npar, &value);
     195                    if (Nscan != 2) {
     196                      fprintf (stderr, "scan failure\n");
     197                    }
     198                    Params->data.F32[Npar] = value;
     199                    break;
     200                  case 'g':
     201                  case 'G':
     202                  case '\n':
     203                    getOptions = false;
     204                    break;
     205                  default:
     206                    done = true;
     207                    break;
     208                }
     209                fprintf (stderr, "foo\n");
     210            }
     211            if (done) break;
     212        }
     213           
    140214        // dump some useful info if trace is defined
    141215        if (psTraceGetLevel(FACILITY) >= 6) {
    … …  
    202276        // XXX : Madsen gives suggestion for better use of rho
    203277        // rho is positive if the new chisq is smaller
    204         if (rho >= -1e-6) {
     278        if (testValue || (rho >= -1e-6)) {
    205279            min->value = Chisq;
    206280            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
    … …  
    215289          case 0:
    216290            if (rho >= -1e-6) {
    217                 lambda *= 0.25;
     291                lambda *= 0.1;
    218292            } else {
    219293                lambda *= 10.0;
    … …  
    234308            if (rho > 0.0) {
    235309                lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
    236                 nu = 2.0;
     310                nu = 3.0;
    237311            } else {
    238312                lambda *= nu;
    239                 nu *= 2.0;
     313                nu *= 3.0;
    240314            }
    241315            break;
    … …  
    381455            coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0);
    382456
    383             pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord);
     457            pcm->modelFlux->data.F32[i][j] = pcm->modelConv->class->modelFunc (deriv, params, coord);
    384458
    385459            for (int n = 0; n < params->n; n++) {
    … …  
    408482# else
    409483    if (pcm->use1Dgauss) {
    410         // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    411         // * the model flux is not masked
    412         // * threading takes place above this level
    413         pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
    414         psImageSmooth_PreAlloc_F32 (pcm->modelConvFlux, pcm->smdata);
    415         // psImageSmooth (pcm->modelConvFlux, pcm->sigma, pcm->nsigma);
     484
     485        if (USE_1D_CACHE) {
     486            // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     487            // * the model flux is not masked
     488            // * threading takes place above this level
     489            pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
     490            psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata);
     491        } else {
     492            pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
     493            psImageSmooth2dCache_F32 (pcm->modelConvFlux, pcm->smdata2d);
     494        }
    416495    } else {
    417496        psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
    … …  
    428507# else
    429508        if (pcm->use1Dgauss) {
    430             // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    431             // * the model flux is not masked
    432             // * threading takes place above this level
    433             dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
    434             psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata);
    435             // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma);
     509            if (USE_1D_CACHE) {
     510                // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     511                // * the model flux is not masked
     512                // * threading takes place above this level
     513                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     514                psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
     515            } else {
     516                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     517                psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d);
     518            }
    436519        } else {
    437520            psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
    … …  
    449532
    450533        if (pcm->use1Dgauss) {
    451             // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    452             // * the model flux is not masked
    453             // * threading takes place above this level
    454             dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
    455             psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata);
    456             // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma);
     534            if (USE_1D_CACHE) {
     535                // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     536                // * the model flux is not masked
     537                // * threading takes place above this level
     538                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     539                psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
     540            } else {
     541                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     542                psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d);
     543            }
    457544        } else {
    458545            psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf);
    … …  
    474561    // XXX TEST : SAVE IMAGES
    475562# if (SAVE_IMAGES)
    476     psphotSaveImage (NULL, pcm->psf->image, "psf.fits");
    477     psphotSaveImage (NULL, pcm->modelFlux, "model.fits");
    478     psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits");
    479     psphotSaveImage (NULL, source->pixels, "obj.fits");
    480     psphotSaveImage (NULL, source->maskObj, "mask.fits");
    481     psphotSaveImage (NULL, source->variance, "variance.fits");
     563    static int Npass = 0;
     564    char name[128];
     565    if (!pcm->use1Dgauss) {
     566      snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name);
     567    }
     568    snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name);
     569    snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name);
     570    snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name);
     571    snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name);
     572    snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name);
     573    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
     574        psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
     575        if (!dmodelConv) continue;
     576        snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name);
     577    }
     578    Npass ++;
    482579# endif
    483580
    … …  
    506603
    507604            float ymodel  = pcm->modelConvFlux->data.F32[i][j];
    508             float yweight = 1.0 / source->variance->data.F32[i][j];
     605
     606            // XXXX note this point here:::
     607            float yweight = pcm->poissonErrors ? 1.0 / source->variance->data.F32[i][j] : 1.0;
    509608            float delta = ymodel - source->pixels->data.F32[i][j];
    510609
Note: See TracChangeset for help on using the changeset viewer.