public software.psfex

[/] [branches/] [rhl/] [src/] [makeit.c] - Diff between revs 188 and 191

Go to most recent revision | Show entire file | Details | Blame | View Log

Rev 188 Rev 191
Line 38... Line 38...
 
 
#include        <math.h>
#include        <math.h>
#include        <stdio.h>
#include        <stdio.h>
#include        <stdlib.h>
#include        <stdlib.h>
#include        <string.h>
#include        <string.h>
#include        <time.h>
 
 
 
#include        "define.h"
#include        "define.h"
#include        "types.h"
#include        "types.h"
#include        "globals.h"
#include        "globals.h"
#include        "fits/fitscat.h"
#include        "fits/fitscat.h"
Line 56... Line 55...
#include        "prefs.h"
#include        "prefs.h"
#include        "psf.h"
#include        "psf.h"
#include        "sample.h"
#include        "sample.h"
#include        "xml.h"
#include        "xml.h"
 
 
psfstruct       *make_psf(setstruct *set, float psfstep,
 
                        float *basis, int nbasis, contextstruct *context);
 
void            write_error(char *msg1, char *msg2);
 
time_t          thetime, thetime2;
time_t          thetime, thetime2;
 
void            write_error(char *msg1, char *msg2);
 
 
/********************************** makeit ***********************************/
/********************************** makeit ***********************************/
/*
/*
*/
*/
void    makeit(void)
void    makeit(void)
Line 145... Line 142...
        fields[c]->next, fields[c]->next>1 ? "s":"",
        fields[c]->next, fields[c]->next>1 ? "s":"",
        fields[c]->ndet, fields[c]->ndet>1 ? "s":"");
        fields[c]->ndet, fields[c]->ndet>1 ? "s":"");
    }
    }
  QPRINTF(OUTPUT, "\n");
  QPRINTF(OUTPUT, "\n");
 
 
  next = fields[0]->next;
 
 
 
  if (prefs.xml_flag)
  if (prefs.xml_flag)
    init_xml(ncat);
    init_xml(ncat);
 
 
  psfstep = prefs.psf_step;
  makeit_body(fields, &context, &fullcontext, 1);
  psfsteps = NULL;
  next = fields[0]->next;
  nbasis = 0;
 
  psfbasis = NULL;
 
  psfbasiss = NULL;
 
 
 
/* Initialize context */
 
  NFPRINTF(OUTPUT, "Initializing contexts...");
 
  context = context_init(prefs.context_name, prefs.context_group,
 
                prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
 
                CONTEXT_REMOVEHIDDEN);
 
  fullcontext = context->npc?
 
                  context_init(prefs.context_name, prefs.context_group,
 
                    prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
 
                    CONTEXT_KEEPHIDDEN)
 
                : context;
 
 
 
  if (context->npc && ncat<2)
 
    warning("Hidden dependencies cannot be derived from",
 
        " a single catalog");
 
  else if (context->npc && prefs.stability_type == STABILITY_EXPOSURE)
 
    warning("Hidden dependencies have no effect",
 
        " in STABILITY_TYPE EXPOSURE mode");
 
 
 
/* Compute PSF steps */
 
  if (!prefs.psf_step)
 
    {
 
    NFPRINTF(OUTPUT, "Computing optimum PSF sampling steps...");
 
    if (prefs.newbasis_type==NEWBASIS_PCACOMMON
 
        || (prefs.stability_type == STABILITY_SEQUENCE
 
                && prefs.psf_mef_type == PSF_MEF_COMMON))
 
      {
 
      set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
 
      psfstep = (float)((set->fwhm/2.35)*0.5);
 
      end_set(set);
 
      }
 
/*-- Need to derive a common pixel step for each ext */
 
    else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT
 
        || context->npc
 
        || (prefs.stability_type == STABILITY_SEQUENCE
 
                && prefs.psf_mef_type == PSF_MEF_INDEPENDENT))
 
      {
 
/*-- Run through all samples to derive a different pixel step for each extension */
 
      QMALLOC(psfsteps, float, next);
 
      for (ext=0 ; ext<next; ext++)
 
        {
 
        set = load_samples(incatnames, 0, ncat, ext, next, context);
 
        psfsteps[ext] = (float)(psfstep? psfstep : (set->fwhm/2.35)*0.5);
 
        end_set(set);
 
        }
 
      }
 
    }
 
 
 
/* Derive a new common PCA basis for all extensions */
 
  if (prefs.newbasis_type==NEWBASIS_PCACOMMON)
 
    {
 
    QMALLOC(cpsf, psfstruct *, ncat*next);
 
    for (ext=0 ; ext<next; ext++)
 
      for (c=0; c<ncat; c++)
 
        {
 
        sprintf(str, "Computing new PCA image basis from %s...",
 
                fields[c]->rtcatname);
 
        NFPRINTF(OUTPUT, str);
 
        set = load_samples(incatnames, c, 1, ext, next, context);
 
        step = psfstep;
 
        cpsf[c+ext*ncat] = make_psf(set, psfstep, NULL, 0, context);
 
        end_set(set);
 
        }
 
    nbasis = prefs.newbasis_number;
 
    psfbasis = pca_onsnaps(cpsf, ncat*next, nbasis);
 
    for (i=0 ; i<ncat*next; i++)
 
      psf_end(cpsf[i]);
 
    free(cpsf);
 
    }
 
/* Derive a new PCA basis for each extension */
 
  else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT)
 
    {
 
    nbasis = prefs.newbasis_number;
 
    QMALLOC(psfbasiss, float *, next);
 
    for (ext=0; ext<next; ext++)
 
      {
 
      if (psfsteps)
 
        step = psfsteps[ext];
 
      else
 
        step = psfstep;
 
      QMALLOC(cpsf, psfstruct *, ncat);
 
      for (c=0; c<ncat; c++)
 
        {
 
        sprintf(str, "Computing new PCA image basis from %s...",
 
                fields[c]->rtcatname);
 
        NFPRINTF(OUTPUT, str);
 
        set = load_samples(incatnames, c, 1, ext, next, context);
 
        cpsf[c] = make_psf(set, step, NULL, 0, context);
 
        end_set(set);
 
        }
 
      psfbasiss[ext] = pca_onsnaps(cpsf, ncat, nbasis);
 
      for (c=0 ; c<ncat; c++)
 
        psf_end(cpsf[c]);
 
      free(cpsf);
 
      }
 
    }
 
 
 
  if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_COMMON)
 
/*-- Derive principal components of PSF variation from the whole mosaic */
 
    {
 
    p = 0;
 
    QMALLOC(cpsf, psfstruct *, ncat*next);
 
    for (c=0; c<ncat; c++)
 
      {
 
      sprintf(str, "Computing hidden dependency parameter(s) from %s...",
 
                fields[c]->rtcatname);
 
      NFPRINTF(OUTPUT, str);
 
      for (ext=0 ; ext<next; ext++)
 
        {
 
        set = load_samples(incatnames, c, 1, ext, next, context);
 
        if (psfsteps)
 
          step = psfsteps[ext];
 
        else
 
          step = psfstep;
 
        basis = psfbasiss? psfbasiss[ext] : psfbasis;
 
        cpsf[p++] = make_psf(set, step, basis, nbasis, context);
 
        end_set(set);
 
        }
 
      }
 
    free(fullcontext->pc);
 
    fullcontext->pc = pca_oncomps(cpsf, next, ncat, context->npc);
 
    for (c=0 ; c<ncat*next; c++)
 
      psf_end(cpsf[c]);
 
    free(cpsf);
 
    }
 
 
 
/* Compute "final" PSF models */
 
  if (prefs.psf_mef_type == PSF_MEF_COMMON)
 
    {
 
    if (prefs.stability_type == STABILITY_SEQUENCE)
 
      {
 
/*---- Load all the samples at once */
 
      set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
 
      step = psfstep;
 
      basis = psfbasis;
 
      field_count(fields, set, COUNT_LOADED);
 
      psf = make_psf(set, step, basis, nbasis, context);
 
      field_count(fields, set, COUNT_ACCEPTED);
 
      end_set(set);
 
      NFPRINTF(OUTPUT, "Computing final PSF model...");
 
      context_apply(context, psf, fields, ALL_EXTENSIONS, 0, ncat);
 
      psf_end(psf);
 
      }
 
    else
 
      for (c=0; c<ncat; c++)
 
        {
 
/*------ Load the samples for current exposure */
 
        sprintf(str, "Computing final PSF model from %s...",
 
                fields[c]->rtcatname);
 
        NFPRINTF(OUTPUT, str);
 
        set = load_samples(incatnames, c, 1, ALL_EXTENSIONS, next, context);
 
        if (psfstep)
 
          step = psfstep;
 
        else
 
          step = (float)((set->fwhm/2.35)*0.5);
 
        basis = psfbasis;
 
        field_count(fields, set, COUNT_LOADED);
 
        psf = make_psf(set, step, basis, nbasis, fullcontext);
 
        field_count(fields, set, COUNT_ACCEPTED);
 
        end_set(set);
 
        context_apply(fullcontext, psf, fields, ALL_EXTENSIONS, c, 1);
 
        psf_end(psf);
 
        }
 
    }
 
  else
 
    for (ext=0 ; ext<next; ext++)
 
      {
 
      basis = psfbasiss? psfbasiss[ext] : psfbasis;
 
      if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_INDEPENDENT)
 
/*------ Derive principal components of PSF components */
 
        {
 
        QMALLOC(cpsf, psfstruct *, ncat);
 
        if (psfsteps)
 
          step = psfsteps[ext];
 
        else
 
          step = psfstep;
 
        for (c=0; c<ncat; c++)
 
          {
 
          if (next>1)
 
            sprintf(str,
 
                "Computing hidden dependency parameter(s) from %s[%d/%d]...",
 
                fields[c]->rtcatname, ext+1, next);
 
          else
 
            sprintf(str, "Computing hidden dependency parameter(s) from %s...",
 
                fields[c]->rtcatname);
 
          NFPRINTF(OUTPUT, str);
 
          set = load_samples(incatnames, c, 1, ext, next, context);
 
          field_count(fields, set, COUNT_LOADED);
 
          cpsf[c] = make_psf(set, step, basis, nbasis, context);
 
          field_count(fields, set, COUNT_ACCEPTED);
 
          end_set(set);
 
          }
 
        free(fullcontext->pc);
 
        fullcontext->pc = pca_oncomps(cpsf, 1, ncat, context->npc);
 
        for (c=0 ; c<ncat; c++)
 
          psf_end(cpsf[c]);
 
        free(cpsf);
 
        }
 
 
 
      if (prefs.stability_type == STABILITY_SEQUENCE)
 
        {
 
/*------ Load all the samples at once */
 
        if (next>1)
 
          {
 
          sprintf(str, "Computing final PSF model for extension [%d/%d]...",
 
                ext+1, next);
 
          NFPRINTF(OUTPUT, str);
 
          }
 
        else
 
          NFPRINTF(OUTPUT, "Computing final PSF model...");
 
        set = load_samples(incatnames, 0, ncat, ext, next, fullcontext);
 
        if (psfstep)
 
          step = psfstep;
 
        else if (psfsteps)
 
          step = psfsteps[ext];
 
        else
 
          step = (float)((set->fwhm/2.35)*0.5);
 
        field_count(fields, set, COUNT_LOADED);
 
        psf = make_psf(set, step, basis, nbasis, fullcontext);
 
        field_count(fields, set, COUNT_ACCEPTED);
 
        end_set(set);
 
        context_apply(fullcontext, psf, fields, ext, 0, ncat);
 
        psf_end(psf);
 
        }
 
      else
 
        for (c=0; c<ncat; c++)
 
          {
 
/*-------- Load the samples for current exposure */
 
          if (next>1)
 
            sprintf(str, "Reading data from %s[%d/%d]...",
 
                fields[c]->rtcatname, ext+1, next);
 
          else
 
            sprintf(str, "Reading data from %s...",
 
                fields[c]->rtcatname);
 
          NFPRINTF(OUTPUT, str);
 
          set = load_samples(incatnames, c, 1, ext, next, context);
 
          if (psfstep)
 
            step = psfstep;
 
          else if (psfsteps)
 
            step = psfsteps[ext];
 
          else
 
            step = (float)((set->fwhm/2.35)*0.5);
 
          if (next>1)
 
            sprintf(str, "Computing final PSF model for %s[%d/%d]...",
 
                fields[c]->rtcatname, ext+1, next);
 
          else
 
            sprintf(str, "Computing final PSF model for %s...",
 
                fields[c]->rtcatname);
 
          NFPRINTF(OUTPUT, str);
 
          field_count(fields, set, COUNT_LOADED);
 
          psf = make_psf(set, step, basis, nbasis, context);
 
          field_count(fields, set, COUNT_ACCEPTED);
 
          end_set(set);
 
          context_apply(context, psf, fields, ext, c, 1);
 
          psf_end(psf);
 
          }
 
      }
 
 
 
  free(psfsteps);
 
  if (psfbasiss)
 
    {
 
    for (ext=0; ext<next; ext++)
 
      free(psfbasiss[ext]);
 
    free(psfbasiss);
 
    }
 
  else if (psfbasis)
 
    free(psfbasis);
 
 
 
/* Compute diagnostics and check-images */
/* Write XML */
#ifdef USE_THREADS
  if (prefs.xml_flag)
/* Force MKL using a single thread as diagnostic code is already multithreaded*/
 
 #ifdef HAVE_MKL
 
  if (prefs.context_nsnap>2)
 
    mkl_set_num_threads(1);
 
 #endif
 
#endif
 
  QIPRINTF(OUTPUT,
 
        "   filename      [ext] accepted/total samp. chi2/dof FWHM ellip."
 
        " resi. asym.");
 
  for (c=0; c<ncat; c++)
 
    {
 
    field = fields[c];
 
    for (ext=0 ; ext<next; ext++)
 
      {
 
      psf = field->psf[ext];
 
      wcs = field->wcs[ext];
 
      if (next>1)
 
        sprintf(str, "Computing diagnostics for %s[%d/%d]...",
 
                field->rtcatname, ext+1, next);
 
      else
 
        sprintf(str, "Computing diagnostics for %s...",
 
                field->rtcatname);
 
      NFPRINTF(OUTPUT, str);
 
/*---- Check PSF with individual datasets */
 
      set2 = load_samples(incatnames, c, 1, ext, next, context);
 
      psf->samples_loaded = set2->nsample;
 
      if (set2->nsample>1)
 
        {
 
/*------ Remove bad PSF candidates */
 
        psf_clean(psf, set2, prefs.prof_accuracy);
 
        psf->chi2 = set2->nsample? psf_chi2(psf, set2) : 0.0;
 
        }
 
      psf->samples_accepted = set2->nsample;
 
/*---- Compute diagnostics and field statistics */
 
      psf_diagnostic(psf);
 
      psf_wcsdiagnostic(psf, wcs);
 
      nmed = psf->nmed;
 
      field_stats(fields, set2);
 
/*---- Display stats for current catalog/extension */
 
      if (next>1)
 
        sprintf(str, "[%d/%d]", ext+1, next);
 
      else
 
        str[0] = '\0';
 
      QPRINTF(OUTPUT, "%-17.17s%-7.7s %5d/%-5d %6.2f %6.2f %6.2f  %4.2f"
 
        " %5.2f %5.2f\n",
 
        ext==0? field->rtcatname : "",
 
        str,
 
        psf->samples_accepted, psf->samples_loaded,
 
        psf->pixstep,
 
        psf->chi2,
 
        psf->moffat_fwhm,
 
        psf->moffat_ellipticity,
 
        psf->pfmoffat_residuals,
 
        psf->sym_residuals);
 
/*---- Save "Check-images" */
 
      for (i=0; i<prefs.ncheck_type; i++)
 
        if (prefs.check_type[i])
 
          {
          {
          sprintf(str, "Saving CHECK-image #%d...", i+1);
    NFPRINTF(OUTPUT, "Writing XML file...");
          NFPRINTF(OUTPUT, str);
    write_xml(prefs.xml_name);
          check_write(field, set2, prefs.check_name[i], prefs.check_type[i],
    end_xml();
                ext, next, prefs.check_cubeflag);
 
          }
 
/*---- Free memory */
 
      end_set(set2);
 
      }
 
    }
    }
 
 
#ifdef USE_THREADS
 
/* Back to multithreaded MKL */
 
 #ifdef HAVE_MKL
 
   mkl_set_num_threads(prefs.nthreads);
 
 #endif
 
#endif
 
 
 
/* Save result */
/* Save result */
  for (c=0; c<ncat; c++)
  for (c=0; c<ncat; c++)
    {
    {
    sprintf(str, "Saving PSF model and metadata for %s...",
    sprintf(str, "Saving PSF model and metadata for %s...",
        fields[c]->rtcatname);
        fields[c]->rtcatname);
Line 561... Line 217...
    cplot_counts(fields[c]);
    cplot_counts(fields[c]);
    cplot_countfrac(fields[c]);
    cplot_countfrac(fields[c]);
    cplot_modchi2(fields[c]);
    cplot_modchi2(fields[c]);
    cplot_modresi(fields[c]);
    cplot_modresi(fields[c]);
#endif
#endif
 
 
/*-- Update XML */
/*-- Update XML */
    if (prefs.xml_flag)
    if (prefs.xml_flag)
      update_xml(fields[c]);
      update_xml(fields[c]);
    }
    }
 
 
Line 575... Line 232...
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
  sprintf(prefs.stime_end,"%02d:%02d:%02d",
  sprintf(prefs.stime_end,"%02d:%02d:%02d",
        tm->tm_hour, tm->tm_min, tm->tm_sec);
        tm->tm_hour, tm->tm_min, tm->tm_sec);
  prefs.time_diff = difftime(thetime2, thetime);
  prefs.time_diff = difftime(thetime2, thetime);
 
 
/* Write XML */
 
  if (prefs.xml_flag)
 
    {
 
    NFPRINTF(OUTPUT, "Writing XML file...");
 
    write_xml(prefs.xml_name);
 
    end_xml();
 
    }
 
 
 
/* Free memory */
/* Free memory */
  for (c=0; c<ncat; c++)
  for (c=0; c<ncat; c++)
    field_end(fields[c]);
    field_end(fields[c]);
  free(fields);
  free(fields);
 
 
Line 596... Line 245...
 
 
  return;
  return;
  }
  }
 
 
 
 
/****** make_psf *************************************************************
 
PROTO   psfstruct *make_psf(setstruct *set, float psfstep,
 
                        float *basis, int nbasis, contextstruct *context)
 
PURPOSE Make PSFs from a set of FITS binary catalogs.
 
INPUT   Pointer to a sample set,
 
        PSF sampling step,
 
        Pointer to basis image vectors,
 
        Number of basis vectors,
 
        Pointer to context structure.
 
OUTPUT  Pointer to the PSF structure.
 
NOTES   Diagnostics are computed only if diagflag != 0.
 
AUTHOR  E. Bertin (IAP)
 
VERSION 03/09/2009
 
 ***/
 
psfstruct       *make_psf(setstruct *set, float psfstep,
 
                        float *basis, int nbasis, contextstruct *context)
 
  {
 
   psfstruct            *psf;
 
   basistypenum         basistype;
 
   float                pixsize[2];
 
 
 
  pixsize[0] = (float)prefs.psf_pixsize[0];
 
  pixsize[1] = (float)prefs.psf_pixsize[1];
 
//  NFPRINTF(OUTPUT,"Initializing PSF modules...");
 
  psf = psf_init(context, prefs.psf_size, psfstep, pixsize, set->nsample);
 
 
 
  psf->samples_loaded = set->nsample;
 
  psf->fwhm = set->fwhm;
 
 
 
/* Make the basic PSF-model (1st pass) */
 
//  NFPRINTF(OUTPUT,"Modeling the PSF (1/3)...");
 
  psf_make(psf, set, 0.2);
 
  if (basis && nbasis)
 
    {
 
    QMEMCPY(basis, psf->basis, float, nbasis*psf->size[0]*psf->size[1]);
 
    psf->nbasis = nbasis;
 
    }
 
  else
 
    {
 
//    NFPRINTF(OUTPUT,"Generating the PSF model...");
 
    basistype = prefs.basis_type;
 
    if (basistype==BASIS_PIXEL_AUTO)
 
      basistype = (psf->fwhm < PSF_AUTO_FWHM)? BASIS_PIXEL : BASIS_NONE;
 
    psf_makebasis(psf, set, basistype, prefs.basis_number);
 
    }
 
  psf_refine(psf, set);
 
 
 
/* Remove bad PSF candidates */
 
  if (set->nsample>1)
 
    {
 
    psf_clean(psf, set, 0.2);
 
 
 
/*-- Make the basic PSF-model (2nd pass) */
 
//    NFPRINTF(OUTPUT,"Modeling the PSF (2/3)...");
 
    psf_make(psf, set, 0.1);
 
    psf_refine(psf, set);
 
    }
 
 
 
/* Remove bad PSF candidates */
 
  if (set->nsample>1)
 
    {
 
    psf_clean(psf, set, 0.1);
 
 
 
/*-- Make the basic PSF-model (3rd pass) */
 
//    NFPRINTF(OUTPUT,"Modeling the PSF (3/3)...");
 
    psf_make(psf, set, 0.05);
 
    psf_refine(psf, set);
 
    }
 
 
 
/* Remove bad PSF candidates */
 
  if (set->nsample>1)
 
    psf_clean(psf, set, 0.05);
 
 
 
  psf->samples_accepted = set->nsample;
 
 
 
/* Refine the PSF-model */
 
  psf_make(psf, set, prefs.prof_accuracy);
 
  psf_refine(psf, set);
 
 
 
/* Clip the PSF-model */
 
  psf_clip(psf);
 
 
 
/*-- Just check the Chi2 */
 
  psf->chi2 = set->nsample? psf_chi2(psf, set) : 0.0;
 
 
 
  return psf;
 
  }
 
 
 
 
 
/****** write_error ********************************************************
/****** write_error ********************************************************
PROTO   void    write_error(char *msg1, char *msg2)
PROTO   void    write_error(char *msg1, char *msg2)
PURPOSE Manage files in case of a catched error
PURPOSE Manage files in case of a catched error
Line 707... Line 269...
  end_xml();
  end_xml();
 
 
  return;
  return;
  }
  }
 
 
 
 
 
 
 No newline at end of file
 No newline at end of file