Fixed CONCENTRATION_MODEL issue.

Removed singular matrix warning.
Added MATRIX_MODELERR measurement parameters.
Added DURATION_ANALYSIS measurement parameter.

config.h.in

@@ -62,6 +62,9 @@
 /* Define to 1 if you have the `getpagesize' function. */
 #undef HAVE_GETPAGESIZE
 
+/* Define to 1 if you have the `gettimeofday' function. */
+#undef HAVE_GETTIMEOFDAY
+
 /* Define to 1 if you have the <inttypes.h> header file. */
 #undef HAVE_INTTYPES_H
 

configure

@@ -22427,8 +22427,9 @@ done
+
 for ac_func in atexit getenv memcpy memmove memset mkdir munmap strstr \
-		sincosf logf
+		sincosf logf gettimeofday
 do
 as_ac_var=`$as_echo "ac_cv_func_$ac_func" | $as_tr_sh`
 { $as_echo "$as_me:$LINENO: checking for $ac_func" >&5

configure.ac

@@ -93,7 +93,7 @@ AC_TYPE_SIGNAL
 AC_FUNC_STAT
 AC_FUNC_STRFTIME
 AC_CHECK_FUNCS([atexit getenv memcpy memmove memset mkdir munmap strstr \
-		sincosf logf])
+		sincosf logf gettimeofday])
 
 # Check support for large files
 AC_SYS_LARGEFILE

src/analyse.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	analyse(), endobject()...: measurements on detections.
 *
-*	Last modify:	11/09/2009
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -385,10 +385,14 @@ Final processing of object data, just before saving it to the catalog.
 void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
 		picstruct *dwfield, int n, objliststruct *objlist)
   {
-   objstruct	*obj;
-   checkstruct	*check;
-   double	rawpos[NAXIS];
-   int		i,j, ix,iy,selecflag, newnumber,nsub;
+   objstruct		*obj;
+   checkstruct		*check;
+   double		rawpos[NAXIS],
+			analtime1;
+   int			i,j, ix,iy,selecflag, newnumber,nsub;
+
+  if (FLAG(obj2.analtime))
+    analtime1 = counter_seconds();
 
   obj = &objlist->obj[n];
 
@@ -737,6 +741,9 @@ void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
           obj->number = ++thecat.ntotal;
         }
 
+      if (FLAG(obj2.analtime) && !j)
+        obj2->analtime = (float)(counter_seconds() - analtime1);
+
       FPRINTF(OUTPUT, "%8d %6.1f %6.1f %5.1f %5.1f %12g "
 			"%c%c%c%c%c%c%c%c\n",
 	obj->number, obj->mx+1.0, obj->my+1.0,

src/catout.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	functions for output of catalog data.
 *
-*	Last modify:	16/09/2009
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -378,6 +378,7 @@ void	updateparamflags()
 			| FLAG(obj2.prof_spheroid_peak);
   prefs.prof_flag |= FLAG(obj2.prof_chi2) | FLAG(obj2.prof_niter)
 			| FLAG(obj2.prof_vector) | FLAG(obj2.prof_errvector)
+			| FLAG(obj2.prof_errmatrix)
 			| FLAG(obj2.x_prof) | FLAG(obj2.y_prof)
 			| FLAG(obj2.prof_mx2)
 			| FLAG(obj2.peak_prof)

src/globals.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	global declarations.
 *
-*	Last modify:	11/05/2008
+*	Last modify:	14/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -73,6 +73,8 @@ extern void	alloccatparams(void),
 		write_error(char *msg1, char *msg2),
 		write_vo_fields(FILE *file);
 
+extern double	counter_seconds(void);
+
 extern float	hmedian(float *, int);
 
 extern int	addobj(int, objliststruct *, objliststruct *),

src/levmar/Axb_core.c

@@ -923,7 +923,7 @@ LM_REAL *a, *b, *work;
 			exit(1);
 		}
 		else{
-      fprintf(stderr, RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
+//      fprintf(stderr, RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
 #ifndef LINSOLVERS_RETAIN_MEMORY
       free(buf);
 #endif

src/makeit.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	main program.
 *
-*	Last modify:	20/11/2008
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -62,7 +62,7 @@ void	makeit()
    patternstruct	*pattern;
    static time_t        thetime1, thetime2;
    struct tm		*tm;
-   int			nflag[MAXFLAG],
+   int			nflag[MAXFLAG], nparam2[2],
 			i, nok, ntab, next, ntabmax, forcextflag,
 			nima0,nima1, nweight0,nweight1, npat;
 
@@ -110,8 +110,10 @@ void	makeit()
 /* Create profiles at full resolution */
     NFPRINTF(OUTPUT, "Preparing profile models");
     theprofit = profit_init(thepsf);
-    changecatparamarrays("VECTOR_PROF", &theprofit->nparam, 1);
-    changecatparamarrays("VECTOR_PROFERR", &theprofit->nparam, 1);
+    changecatparamarrays("VECTOR_MODEL", &theprofit->nparam, 1);
+    changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1);
+    nparam2[0] = nparam2[1] = theprofit->nparam;
+    changecatparamarrays("MATRIX_MODELERR", &nparam2, 2);
     if (prefs.pattern_flag)
       {
       npat = prefs.prof_disk_patternvectorsize;

src/misc.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	miscellaneous functions.
 *
-*	Last modify:	13/12/2002
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -18,6 +18,9 @@
 #include        "config.h"
 #endif
 
+#include	<time.h>
+#include	<sys/time.h>
+
 #include	"define.h"
 #include	"globals.h"
 
@@ -70,3 +73,23 @@ float	hmedian(float *ra, int n)
   }
 
 
+/****** counter_seconds *******************************************************
+PROTO	double counter_seconds(void)
+PURPOSE	Count the number of seconds (with an arbitrary offset).
+INPUT	-.
+OUTPUT	Returns a number of seconds.
+NOTES	Results are meaningful only for tasks that take one microsec or more.
+AUTHOR	E. Bertin (IAP)
+VERSION	24/09/2009
+ ***/
+double	counter_seconds(void)
+  {
+   struct timeval	tp;
+   struct timezone	tzp;
+   int			dummy;
+
+  dummy = gettimeofday(&tp,&tzp);
+  return (double) tp.tv_sec + (double) tp.tv_usec * 1.0e-6;
+  }
+
+

src/param.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	parameter list for catalog data.
 *
-*	Last modify:	20/07/2009
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -851,6 +851,10 @@ keystruct	objkey[] = {
 	&outobj2.poserrtheta1950_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg",
 	"pos.posAng;pos.errorEllipse", "deg"},
 
+  {"DURATION_ANALYSIS", "Duration of analysis for this source",
+	&outobj2.analtime, H_FLOAT, T_FLOAT, "%9.4g", "s",
+	"time.duration;time.processing", "s"},
+
 /* Temporarily (at least) deactivated not to confuse users with model-fitting
   {"X2PC_IMAGE", "PC variance along x",
 	&outobj2.mx2_pc, H_EXPO, T_DOUBLE, "%18.10e", "pixel**2",

src/paramprofit.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	Model-fitting parameter list for catalog data.
 *
-*	Last modify:	16/09/2009
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -21,6 +21,10 @@
 	&outobj2.prof_errvector, H_FLOAT, T_FLOAT, "%12.4g", "",
 	"stat.stdev;stat.fit;src.morph.param", "", 1,
 	&prefs.prof_errvectorsize},
+  {"MATRIX_MODELERR", "Model-fitting covariance matrix",
+	&outobj2.prof_errmatrix, H_FLOAT, T_FLOAT, "%12.4g", "",
+	"stat.covariance;stat.fit;src.morph.param", "", 2,
+	prefs.prof_errmatrixsize},
   {"CHI2_MODEL", "Reduced Chi2 of the fit",
 	&outobj2.prof_chi2, H_FLOAT, T_FLOAT, "%12.7g", "",
 	"stat.fit.chi2;src.morph", ""},

src/prefs.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	Keywords for the configuration file.
 *
-*	Last modify:	20/11/2008
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -215,6 +215,7 @@ typedef struct
 /*----- Profile-fitting */
   int		prof_vectorsize;			/* nb of params */
   int		prof_errvectorsize;			/* nb of params */
+  int		prof_errmatrixsize[2];			/* nb of params */
   int		prof_disk_patternvectorsize;		/* nb of params */
   int		prof_disk_patternncomp;			/* nb of params */
   int		prof_disk_patternmodvectorsize;		/* nb of params */

src/profit.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	Fit an arbitrary profile combination to a detection.
 *
-*	Last modify:	20/09/2009
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -163,7 +163,7 @@ OUTPUT	Pointer to an allocated fit structure (containing details about the
 	fit).
 NOTES	It is a modified version of the lm_minimize() of lmfit.
 AUTHOR	E. Bertin (IAP)
-VERSION	21/09/2009
+VERSION	24/09/2009
  ***/
 void	profit_fit(profitstruct *profit,
 		picstruct *field, picstruct *wfield,
@@ -176,9 +176,10 @@ void	profit_fit(profitstruct *profit,
     checkstruct		*check;
     double		emx2,emy2,emxy, a , cp,sp, cn, bn, n;
     float		psf_fwhm, dchi2, err;
-    int			i,j,p, nparam, ncomp, nprof;
+    int			i,j,p, nparam, nparam2, ncomp, nprof;
 
   nparam = profit->nparam;
+  nparam2 = nparam*nparam;
   nprof = profit->nprof;
   if (profit->psfdft)
     {
@@ -218,6 +219,12 @@ void	profit_fit(profitstruct *profit,
     if (FLAG(obj2.prof_vector))
       for (p=0; p<nparam; p++)
         obj2->prof_vector[p] = 0.0;
+    if (FLAG(obj2.prof_errvector))
+      for (p=0; p<nparam; p++)
+        obj2->prof_errvector[p] = 0.0;
+    if (FLAG(obj2.prof_errmatrix))
+      for (p=0; p<nparam2; p++)
+        obj2->prof_errmatrix[p] = 0.0;
     obj2->prof_niter = 0;
     return;
     }
@@ -255,26 +262,6 @@ void	profit_fit(profitstruct *profit,
   profit->niter = profit_minimize(profit, PROFIT_MAXITER);
   profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi);
 
-/*
-  QMEMCPY(profit->paraminit, oldparaminit, float, nparam);
-  if (profit_setparam(profit, PARAM_ARMS_PITCH, 160.0, 130.0, 175.0)==RETURN_OK)
-    {
-    oldchi2 = profit->chi2;
-    oldniter = profit->niter;
-    profit_resetparams(profit);
-    profit_setparam(profit, PARAM_ARMS_PITCH, 160.0, 130.0, 175.0);
-    profit->niter = profit_minimize(profit, PROFIT_MAXITER);
-    if (profit->chi2 > oldchi2)
-      {
-      memcpy(profit->paraminit, oldparaminit, nparam*sizeof(float));
-      profit->chi2 = oldchi2;
-      profit->niter = oldniter;
-      }
-    else
-      obj2->prof_flag |= PROFIT_FLIPPED;
-    }  
-*/
-
 /* Convert covariance matrix to bound space */
   profit_covarunboundtobound(profit);
   for (p=0; p<nparam; p++)
@@ -309,6 +296,11 @@ void	profit_fit(profitstruct *profit,
     for (p=0; p<nparam; p++)
       obj2->prof_errvector[p]= profit->paramerr[p];
     }
+  if (FLAG(obj2.prof_errmatrix))
+    {
+    for (p=0; p<nparam2; p++)
+      obj2->prof_errmatrix[p]= profit->covar[p];
+    }
 
   obj2->prof_niter = profit->niter;
   obj2->flux_prof = profit->flux;
@@ -601,7 +593,7 @@ void	profit_fit(profitstruct *profit,
     }
 
 /* Star/galaxy classification */
-  if (FLAG(obj2.prof_class_star))
+  if (FLAG(obj2.prof_class_star) || FLAG(obj2.prof_concentration))
     {
     pprofit = *profit;
     memset(pprofit.paramindex, 0, PARAM_NPARAM*sizeof(int));
@@ -620,15 +612,21 @@ void	profit_fit(profitstruct *profit,
     pprofit.paraminit[pprofit.paramindex[PARAM_DISK_FLUX]] = profit->flux;
     pprofit.niter = profit_minimize(&pprofit, PROFIT_MAXITER);
     profit_residuals(&pprofit,field,wfield, 10.0, pprofit.param,pprofit.resi);
-    dchi2 = 0.5*(pprofit.chi2 - profit->chi2);
-    obj2->prof_class_star = dchi2 < 50.0?
+    if (FLAG(obj2.prof_class_star))
+      {
+      dchi2 = 0.5*(pprofit.chi2 - profit->chi2);
+      obj2->prof_class_star = dchi2 < 50.0?
 	(dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0;
-    if (profit->flux > 0.0 && pprofit.flux > 0.0)
-      obj2->prof_concentration = -2.5*log10(pprofit.flux / profit->flux);
-    else  if (profit->flux > 0.0)
-      obj2->prof_concentration = 99.0;
-    else  if (pprofit.flux > 0.0)
-      obj2->prof_concentration = -99.0;
+      }
+    if (FLAG(obj2.prof_concentration))
+      {
+      if (profit->flux > 0.0 && pprofit.flux > 0.0)
+        obj2->prof_concentration = -2.5*log10(pprofit.flux / profit->flux);
+      else  if (profit->flux > 0.0)
+        obj2->prof_concentration = 99.0;
+      else  if (pprofit.flux > 0.0)
+        obj2->prof_concentration = -99.0;
+      }
     prof_end(pprofit.prof[0]);
     free(pprofit.prof);
     free(pprofit.covar);
@@ -642,9 +640,7 @@ void	profit_fit(profitstruct *profit,
   free(profit->objpix);
   free(profit->objweight);
   free(profit->resi);
-/*
-  free(oldparaminit);
-*/
+
   return;
   }
 
@@ -970,7 +966,7 @@ INPUT	Profile-fitting structure,
 OUTPUT	Vector of residuals.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	20/03/2009
+VERSION	24/03/2009
  ***/
 float	*profit_residuals(profitstruct *profit, picstruct *field,
 		picstruct *wfield, float dynparam, float *param, float *resi)
@@ -983,8 +979,11 @@ float	*profit_residuals(profitstruct *profit, picstruct *field,
     profit->param[p] = param[p];
 /* Simple PSF shortcut */
   if (profit->nprof == 1 && profit->prof[0]->code == PROF_DIRAC)
+    {
     profit_resample(profit, profit->psfpix, profit->lmodpix,
 		*profit->prof[0]->flux);
+    profit->flux = *profit->prof[0]->flux;
+    }
   else
     {
     profit->flux = 0.0;

src/types.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	global type definitions.
 *
-*	Last modify:	16/09/2009
+*	Last modify:	24/09/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -307,6 +307,7 @@ typedef struct
 /* ---- Profile-fitting */
   float		*prof_vector;			/* Model parameters */
   float		*prof_errvector;		/* Model parameter errors */
+  float		*prof_errmatrix;		/* Model parameter covariances*/
   float		prof_chi2;			/* Reduced chi2 */
   BYTE		prof_flag;			/* Model-fitting flags */
   BYTE		prof_flagw;			/* Model-fitting WORLD flag */
@@ -465,8 +466,10 @@ typedef struct
   float		prof_arms_starterrw;		/* RMS error */
   float		prof_arms_quadfrac;		/* Arms quadrature fraction */
   float		prof_arms_quadfracerr;		/* RMS error */
-/* ---- MEF */
+/* ---- MEF ----*/
   short		ext_number;			/* FITS extension number */
+/* ---- Time ---- */
+  float		analtime;			/* Analysis time (s) */
   }	obj2struct;
 
 /*----------------------------- lists of objects ----------------------------*/