public software.psfex

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

Show entire file | Details | Blame | View Log

Rev 188 Rev 191
Line 1... Line 1...
/*
/*
*                               makeit.c
*                               makeit2.c
*
*
* Main loop.
* Main loop, with the I/O parts removed
*
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*
*       This file part of:      PSFEx
*       This file part of:      PSFEx
*
*
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 58... Line 57...
#include        "sample.h"
#include        "sample.h"
#include        "xml.h"
#include        "xml.h"
 
 
psfstruct       *make_psf(setstruct *set, float psfstep,
psfstruct       *make_psf(setstruct *set, float psfstep,
                        float *basis, int nbasis, contextstruct *context);
                        float *basis, int nbasis, contextstruct *context);
void            write_error(char *msg1, char *msg2);
 
time_t          thetime, thetime2;
 
 
 
/********************************** makeit ***********************************/
/********************************** makeit_body ******************************/
/*
/*
*/
*/
void    makeit(void)
void
 
makeit_body(
 
   fieldstruct          **fields,
 
   contextstruct        **context_,
 
   contextstruct        **fullcontext_,
 
   int                  free_sets       /* if false, don't free any sets as we don't own them */
 
   )
  {
  {
   wcsstruct            *wcs;
   wcsstruct            *wcs;
   fieldstruct          **fields,
   fieldstruct          *field;
                        *field;
 
   psfstruct            **cpsf,
   psfstruct            **cpsf,
                        *psf;
                        *psf;
   setstruct            *set, *set2;
   setstruct            *set, *set2;
   contextstruct        *context, *fullcontext;
   contextstruct        *context, *fullcontext;
   struct tm            *tm;
 
   char                 str[MAXCHAR];
   char                 str[MAXCHAR];
   char                 **incatnames,
   char                 **incatnames;
                        *pstr;
 
   float                **psfbasiss,
   float                **psfbasiss,
                        *psfsteps, *psfbasis, *basis,
                        *psfsteps, *psfbasis, *basis,
                        psfstep, step;
                        psfstep, step;
   int                  c,i,p, ncat, ext, next, nmed, nbasis;
   int                  c,i,p, ncat, ext, next, nmed, nbasis;
 
 
/* Install error logging */
 
  error_installfunc(write_error);
 
 
 
  incatnames = prefs.incat_name;
  incatnames = prefs.incat_name;
  ncat = prefs.ncat;
  ncat = prefs.ncat;
 
 
/* Processing start date and time */
 
  thetime = time(NULL);
 
  tm = localtime(&thetime);
 
  sprintf(prefs.sdate_start,"%04d-%02d-%02d",
 
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
 
  sprintf(prefs.stime_start,"%02d:%02d:%02d",
 
        tm->tm_hour, tm->tm_min, tm->tm_sec);
 
 
 
  NFPRINTF(OUTPUT, "");
 
  QPRINTF(OUTPUT,
 
        "----- %s %s started on %s at %s with %d thread%s\n\n",
 
                BANNER,
 
                MYVERSION,
 
                prefs.sdate_start,
 
                prefs.stime_start,
 
                prefs.nthreads,
 
                prefs.nthreads>1? "s":"");
 
 
 
 
 
/* End here if no filename has been provided */
 
  if (!ncat)
 
    {
 
/*-- Processing end date and time */
 
    thetime2 = time(NULL);
 
    tm = localtime(&thetime2);
 
    sprintf(prefs.sdate_end,"%04d-%02d-%02d",
 
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
 
    sprintf(prefs.stime_end,"%02d:%02d:%02d",
 
        tm->tm_hour, tm->tm_min, tm->tm_sec);
 
    prefs.time_diff = difftime(thetime2, thetime);
 
 
 
/*-- Write XML */
 
    if (prefs.xml_flag)
 
      {
 
      init_xml(0);
 
      write_xml(prefs.xml_name);
 
      end_xml();
 
      }
 
    return;
 
    }
 
 
 
/* Create an array of PSFs (one PSF for each extension) */
 
  QMALLOC(fields, fieldstruct *, ncat);
 
 
 
  NFPRINTF(OUTPUT, "");
 
  QPRINTF(OUTPUT, "----- %d input catalogues:\n", ncat);
 
  for (c=0; c<ncat; c++)
 
    {
 
    fields[c] = field_init(incatnames[c]);
 
    QPRINTF(OUTPUT, "%-20.20s:  \"%-16.16s\"  %3d extension%s %7d detection%s\n",
 
        fields[c]->rcatname, fields[c]->ident,
 
        fields[c]->next, fields[c]->next>1 ? "s":"",
 
        fields[c]->ndet, fields[c]->ndet>1 ? "s":"");
 
    }
 
  QPRINTF(OUTPUT, "\n");
 
 
 
  next = fields[0]->next;
  next = fields[0]->next;
 
 
  if (prefs.xml_flag)
 
    init_xml(ncat);
 
 
 
  psfstep = prefs.psf_step;
  psfstep = prefs.psf_step;
  psfsteps = NULL;
  psfsteps = NULL;
  nbasis = 0;
  nbasis = 0;
  psfbasis = NULL;
  psfbasis = NULL;
  psfbasiss = NULL;
  psfbasiss = NULL;
Line 166... Line 102...
  fullcontext = context->npc?
  fullcontext = context->npc?
                  context_init(prefs.context_name, prefs.context_group,
                  context_init(prefs.context_name, prefs.context_group,
                    prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
                    prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
                    CONTEXT_KEEPHIDDEN)
                    CONTEXT_KEEPHIDDEN)
                : context;
                : context;
 
  *context_ = context;
 
  *fullcontext_ = fullcontext;
 
 
  if (context->npc && ncat<2)
  if (context->npc && ncat<2)
    warning("Hidden dependencies cannot be derived from",
    warning("Hidden dependencies cannot be derived from",
        " a single catalog");
        " a single catalog");
  else if (context->npc && prefs.stability_type == STABILITY_EXPOSURE)
  else if (context->npc && prefs.stability_type == STABILITY_EXPOSURE)
Line 184... Line 122...
        || (prefs.stability_type == STABILITY_SEQUENCE
        || (prefs.stability_type == STABILITY_SEQUENCE
                && prefs.psf_mef_type == PSF_MEF_COMMON))
                && prefs.psf_mef_type == PSF_MEF_COMMON))
      {
      {
      set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
      set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
      psfstep = (float)((set->fwhm/2.35)*0.5);
      psfstep = (float)((set->fwhm/2.35)*0.5);
      end_set(set);
      if (free_sets) end_set(set);
      }
      }
/*-- Need to derive a common pixel step for each ext */
/*-- Need to derive a common pixel step for each ext */
    else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT
    else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT
        || context->npc
        || context->npc
        || (prefs.stability_type == STABILITY_SEQUENCE
        || (prefs.stability_type == STABILITY_SEQUENCE
Line 198... Line 136...
      QMALLOC(psfsteps, float, next);
      QMALLOC(psfsteps, float, next);
      for (ext=0 ; ext<next; ext++)
      for (ext=0 ; ext<next; ext++)
        {
        {
        set = load_samples(incatnames, 0, ncat, ext, next, context);
        set = load_samples(incatnames, 0, ncat, ext, next, context);
        psfsteps[ext] = (float)(psfstep? psfstep : (set->fwhm/2.35)*0.5);
        psfsteps[ext] = (float)(psfstep? psfstep : (set->fwhm/2.35)*0.5);
        end_set(set);
        if (free_sets) end_set(set);
        }
        }
      }
      }
    }
    }
 
 
/* Derive a new common PCA basis for all extensions */
/* Derive a new common PCA basis for all extensions */
Line 216... Line 154...
                fields[c]->rtcatname);
                fields[c]->rtcatname);
        NFPRINTF(OUTPUT, str);
        NFPRINTF(OUTPUT, str);
        set = load_samples(incatnames, c, 1, ext, next, context);
        set = load_samples(incatnames, c, 1, ext, next, context);
        step = psfstep;
        step = psfstep;
        cpsf[c+ext*ncat] = make_psf(set, psfstep, NULL, 0, context);
        cpsf[c+ext*ncat] = make_psf(set, psfstep, NULL, 0, context);
        end_set(set);
        if (free_sets) end_set(set);
        }
        }
    nbasis = prefs.newbasis_number;
    nbasis = prefs.newbasis_number;
    psfbasis = pca_onsnaps(cpsf, ncat*next, nbasis);
    psfbasis = pca_onsnaps(cpsf, ncat*next, nbasis);
    for (i=0 ; i<ncat*next; i++)
    for (i=0 ; i<ncat*next; i++)
      psf_end(cpsf[i]);
      psf_end(cpsf[i]);
Line 243... Line 181...
        sprintf(str, "Computing new PCA image basis from %s...",
        sprintf(str, "Computing new PCA image basis from %s...",
                fields[c]->rtcatname);
                fields[c]->rtcatname);
        NFPRINTF(OUTPUT, str);
        NFPRINTF(OUTPUT, str);
        set = load_samples(incatnames, c, 1, ext, next, context);
        set = load_samples(incatnames, c, 1, ext, next, context);
        cpsf[c] = make_psf(set, step, NULL, 0, context);
        cpsf[c] = make_psf(set, step, NULL, 0, context);
        end_set(set);
        if (free_sets) end_set(set);
        }
        }
      psfbasiss[ext] = pca_onsnaps(cpsf, ncat, nbasis);
      psfbasiss[ext] = pca_onsnaps(cpsf, ncat, nbasis);
      for (c=0 ; c<ncat; c++)
      for (c=0 ; c<ncat; c++)
        psf_end(cpsf[c]);
        psf_end(cpsf[c]);
      free(cpsf);
      free(cpsf);
Line 271... Line 209...
          step = psfsteps[ext];
          step = psfsteps[ext];
        else
        else
          step = psfstep;
          step = psfstep;
        basis = psfbasiss? psfbasiss[ext] : psfbasis;
        basis = psfbasiss? psfbasiss[ext] : psfbasis;
        cpsf[p++] = make_psf(set, step, basis, nbasis, context);
        cpsf[p++] = make_psf(set, step, basis, nbasis, context);
        end_set(set);
        if (free_sets) end_set(set);
        }
        }
      }
      }
    free(fullcontext->pc);
    free(fullcontext->pc);
    fullcontext->pc = pca_oncomps(cpsf, next, ncat, context->npc);
    fullcontext->pc = pca_oncomps(cpsf, next, ncat, context->npc);
    for (c=0 ; c<ncat*next; c++)
    for (c=0 ; c<ncat*next; c++)
Line 293... Line 231...
      step = psfstep;
      step = psfstep;
      basis = psfbasis;
      basis = psfbasis;
      field_count(fields, set, COUNT_LOADED);
      field_count(fields, set, COUNT_LOADED);
      psf = make_psf(set, step, basis, nbasis, context);
      psf = make_psf(set, step, basis, nbasis, context);
      field_count(fields, set, COUNT_ACCEPTED);
      field_count(fields, set, COUNT_ACCEPTED);
      end_set(set);
      if (free_sets) end_set(set);
      NFPRINTF(OUTPUT, "Computing final PSF model...");
      NFPRINTF(OUTPUT, "Computing final PSF model...");
      context_apply(context, psf, fields, ALL_EXTENSIONS, 0, ncat);
      context_apply(context, psf, fields, ALL_EXTENSIONS, 0, ncat);
      psf_end(psf);
      psf_end(psf);
      }
      }
    else
    else
Line 314... Line 252...
          step = (float)((set->fwhm/2.35)*0.5);
          step = (float)((set->fwhm/2.35)*0.5);
        basis = psfbasis;
        basis = psfbasis;
        field_count(fields, set, COUNT_LOADED);
        field_count(fields, set, COUNT_LOADED);
        psf = make_psf(set, step, basis, nbasis, fullcontext);
        psf = make_psf(set, step, basis, nbasis, fullcontext);
        field_count(fields, set, COUNT_ACCEPTED);
        field_count(fields, set, COUNT_ACCEPTED);
        end_set(set);
        if (free_sets) end_set(set);
        context_apply(fullcontext, psf, fields, ALL_EXTENSIONS, c, 1);
        context_apply(fullcontext, psf, fields, ALL_EXTENSIONS, c, 1);
        psf_end(psf);
        psf_end(psf);
        }
        }
    }
    }
  else
  else
Line 345... Line 283...
          NFPRINTF(OUTPUT, str);
          NFPRINTF(OUTPUT, str);
          set = load_samples(incatnames, c, 1, ext, next, context);
          set = load_samples(incatnames, c, 1, ext, next, context);
          field_count(fields, set, COUNT_LOADED);
          field_count(fields, set, COUNT_LOADED);
          cpsf[c] = make_psf(set, step, basis, nbasis, context);
          cpsf[c] = make_psf(set, step, basis, nbasis, context);
          field_count(fields, set, COUNT_ACCEPTED);
          field_count(fields, set, COUNT_ACCEPTED);
          end_set(set);
          if (free_sets) end_set(set);
          }
          }
        free(fullcontext->pc);
        free(fullcontext->pc);
        fullcontext->pc = pca_oncomps(cpsf, 1, ncat, context->npc);
        fullcontext->pc = pca_oncomps(cpsf, 1, ncat, context->npc);
        for (c=0 ; c<ncat; c++)
        for (c=0 ; c<ncat; c++)
          psf_end(cpsf[c]);
          psf_end(cpsf[c]);
Line 375... Line 313...
        else
        else
          step = (float)((set->fwhm/2.35)*0.5);
          step = (float)((set->fwhm/2.35)*0.5);
        field_count(fields, set, COUNT_LOADED);
        field_count(fields, set, COUNT_LOADED);
        psf = make_psf(set, step, basis, nbasis, fullcontext);
        psf = make_psf(set, step, basis, nbasis, fullcontext);
        field_count(fields, set, COUNT_ACCEPTED);
        field_count(fields, set, COUNT_ACCEPTED);
        end_set(set);
        if (free_sets) end_set(set);
        context_apply(fullcontext, psf, fields, ext, 0, ncat);
        context_apply(fullcontext, psf, fields, ext, 0, ncat);
        psf_end(psf);
        psf_end(psf);
        }
        }
      else
      else
        for (c=0; c<ncat; c++)
        for (c=0; c<ncat; c++)
Line 407... Line 345...
                fields[c]->rtcatname);
                fields[c]->rtcatname);
          NFPRINTF(OUTPUT, str);
          NFPRINTF(OUTPUT, str);
          field_count(fields, set, COUNT_LOADED);
          field_count(fields, set, COUNT_LOADED);
          psf = make_psf(set, step, basis, nbasis, context);
          psf = make_psf(set, step, basis, nbasis, context);
          field_count(fields, set, COUNT_ACCEPTED);
          field_count(fields, set, COUNT_ACCEPTED);
          end_set(set);
          if (free_sets) end_set(set);
          context_apply(context, psf, fields, ext, c, 1);
          context_apply(context, psf, fields, ext, c, 1);
          psf_end(psf);
          psf_end(psf);
          }
          }
      }
      }
 
 
Line 489... Line 427...
          NFPRINTF(OUTPUT, str);
          NFPRINTF(OUTPUT, str);
          check_write(field, set2, prefs.check_name[i], prefs.check_type[i],
          check_write(field, set2, prefs.check_name[i], prefs.check_type[i],
                ext, next, prefs.check_cubeflag);
                ext, next, prefs.check_cubeflag);
          }
          }
/*---- Free memory */
/*---- Free memory */
      end_set(set2);
      if (free_sets) end_set(set2);
      }
      }
    }
    }
 
 
#ifdef USE_THREADS
#ifdef USE_THREADS
/* Back to multithreaded MKL */
/* Back to multithreaded MKL */
 #ifdef HAVE_MKL
 #ifdef HAVE_MKL
   mkl_set_num_threads(prefs.nthreads);
   mkl_set_num_threads(prefs.nthreads);
 #endif
 #endif
#endif
#endif
 
 
/* Save result */
 
  for (c=0; c<ncat; c++)
 
    {
 
    sprintf(str, "Saving PSF model and metadata for %s...",
 
        fields[c]->rtcatname);
 
    NFPRINTF(OUTPUT, str);
 
/*-- Create a file name with a "PSF" extension */
 
    if (*prefs.psf_dir)
 
      {
 
      if ((pstr = strrchr(incatnames[c], '/')))
 
        pstr++;
 
      else
 
        pstr = incatnames[c];
 
      sprintf(str, "%s/%s", prefs.psf_dir, pstr);
 
      }
 
    else
 
      strcpy(str, incatnames[c]);
 
    if (!(pstr = strrchr(str, '.')))
 
      pstr = str+strlen(str);
 
    sprintf(pstr, "%s", prefs.psf_suffix);
 
    field_psfsave(fields[c], str);
 
/* Create homogenisation kernels */
 
    if (prefs.homobasis_type != HOMOBASIS_NONE)
 
      {
 
      for (ext=0; ext<next; ext++)
 
        {
 
        if (next>1)
 
          sprintf(str, "Computing homogenisation kernel for %s[%d/%d]...",
 
                fields[c]->rtcatname, ext+1, next);
 
        else
 
          sprintf(str, "Computing homogenisation kernel for %s...",
 
                fields[c]->rtcatname);
 
        NFPRINTF(OUTPUT, str);
 
        if (*prefs.homokernel_dir)
 
          {
 
          if ((pstr = strrchr(incatnames[c], '/')))
 
            pstr++;
 
          else
 
            pstr = incatnames[c];
 
          sprintf(str, "%s/%s", prefs.homokernel_dir, pstr);
 
          }
 
        else
 
          strcpy(str, incatnames[c]);
 
        if (!(pstr = strrchr(str, '.')))
 
          pstr = str+strlen(str);
 
        sprintf(pstr, "%s", prefs.homokernel_suffix);
 
        psf_homo(fields[c]->psf[ext], str, prefs.homopsf_params,
 
                prefs.homobasis_number, prefs.homobasis_scale, ext, next);
 
        }
 
      }
 
#ifdef HAVE_PLPLOT
 
/* Plot diagnostic maps for all catalogs */
 
    cplot_ellipticity(fields[c]);
 
    cplot_fwhm(fields[c]);
 
    cplot_moffatresi(fields[c]);
 
    cplot_asymresi(fields[c]);
 
    cplot_counts(fields[c]);
 
    cplot_countfrac(fields[c]);
 
    cplot_modchi2(fields[c]);
 
    cplot_modresi(fields[c]);
 
#endif
 
/*-- Update XML */
 
    if (prefs.xml_flag)
 
      update_xml(fields[c]);
 
    }
 
 
 
/* Processing end date and time */
 
  thetime2 = time(NULL);
 
  tm = localtime(&thetime2);
 
  sprintf(prefs.sdate_end,"%04d-%02d-%02d",
 
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
 
  sprintf(prefs.stime_end,"%02d:%02d:%02d",
 
        tm->tm_hour, tm->tm_min, tm->tm_sec);
 
  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 */
 
  for (c=0; c<ncat; c++)
 
    field_end(fields[c]);
 
  free(fields);
 
 
 
  if (context->npc)
 
    context_end(fullcontext);
 
  context_end(context);
 
 
 
  return;
  return;
  }
  }
 
 
 
 
/****** make_psf *************************************************************
/****** make_psf *************************************************************
PROTO   psfstruct *make_psf(setstruct *set, float psfstep,
PROTO   psfstruct *make_psf(setstruct *set, float psfstep,
                        float *basis, int nbasis, contextstruct *context)
                        float *basis, int nbasis, contextstruct *context)
PURPOSE Make PSFs from a set of FITS binary catalogs.
PURPOSE Make PSFs from a set of FITS binary catalogs.
INPUT   Pointer to a sample set,
INPUT   Pointer to a sample set,
Line 684... Line 529...
  psf->chi2 = set->nsample? psf_chi2(psf, set) : 0.0;
  psf->chi2 = set->nsample? psf_chi2(psf, set) : 0.0;
 
 
  return psf;
  return psf;
  }
  }
 
 
 
 
/****** write_error ********************************************************
 
PROTO   void    write_error(char *msg1, char *msg2)
 
PURPOSE Manage files in case of a catched error
 
INPUT   a character string,
 
        another character string
 
OUTPUT  RETURN_OK if everything went fine, RETURN_ERROR otherwise.
 
NOTES   -.
 
AUTHOR  E. Bertin (IAP)
 
VERSION 23/02/2007
 
 ***/
 
void    write_error(char *msg1, char *msg2)
 
  {
 
   char error[MAXCHAR];
 
 
 
  sprintf(error, "%s%s", msg1,msg2);
 
  if (prefs.xml_flag)
 
    write_xmlerror(prefs.xml_name, error);
 
  end_xml();
 
 
 
  return;
 
  }
 
 
 
 
 
 
 
 No newline at end of file
 No newline at end of file