/home66/gary/public_html/cloudy/c08_branch/source/optimize_func.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*optimize_func actual function called during evaluation of optimization run */
#include "cddefines.h"
#include "init.h"
#include "lines.h"
#include "prt.h"
#include "called.h"
#include "radius.h"
#include "input.h"
#include "cloudy.h"
#include "cddrive.h"
#include "optimize.h"
#include "grid.h"
/* used below to derive chi2 */
STATIC double chi2_func(double,double,double);

double optimize_func(realnum param[])
{

#define MAXCAT 4

        char /* chCAPLB[5], */
          chFind[5];

        bool lgBAD,
                lgHIT, 
          lgLimOK;

        long int cat,
          i, 
          j, 
          nfound, 
          nobs_cat[MAXCAT],
          np;

        static long int ipobs[NOBSLM];

        double chi1, 
          chi2_cat[MAXCAT],
          chisq, 
          func_v, 
          help,
          predin,
          scld,
          snorm,
          theocl,
          temp_theory;

        static char name_cat[MAXCAT][13] = 
        {
          "rel flux    ",
          "column dens ",
          "abs flux    ",
          "mean Temp   "
        };

        static bool lgLinSet = false;

        DEBUG_ENTRY( "optimize_func()" );

        /* This routine is called by optimizer with values of the
         * variable parameters for CLOUDY in the array p(i). It returns
         * the value FUNC = SUM (obs-model)**2/sig**2 for the lines
         * specified in the observational data file, values held in the
         * common blocks /OBSLIN/ & /OBSINT/
         * replacement input strings for CLOUDY READR held in /chCardSav/
         * parameter information for setting chCardSav held in /parmv/
         * additional variables
         * Gary's variables
         */

        /* zero out lots of variables */
        zero();

        if( optimize.lgOptimFlow )
        {
                fprintf( ioQQQ, " trace, optimize_func variables" );
                for( i=0; i < optimize.nvary; i++ )
                {
                        fprintf( ioQQQ, "%.2e", param[i] );
                }
                fprintf( ioQQQ, "\n" );
        }

        for( i=0; i < optimize.nvary; i++ )
        {
                optimize.vparm[0][i] = param[i];
        }

        /* call routine to pack variables with appropriate
         * CLOUDY input lines given the array of variable parameters p(i) */
        vary_input(&lgLimOK);

        if( !lgLimOK )
        {
                /* these parameters are not within limits of parameter search
                 * >>chng 96 apr 26, as per Peter van Hoof comment */
                fprintf( ioQQQ, " Iteration %ld not within range.\n", 
                  optimize.nOptimiz );

                /* this is error; very bad since not within range of parameters */
                func_v = (double)FLT_MAX;

                /* always increment nOptimiz, even if parameters are out of bounds,
                 * this prevents optimizer to get stuck in infinite loop */
                ++optimize.nOptimiz;
                return( func_v );
        }

        for( i=0; i < optimize.nvary; i++ )
        {
                optimize.varmax[i] = (realnum)MAX2(optimize.varmax[i],optimize.vpused[i]);
                optimize.varmin[i] = (realnum)MIN2(optimize.varmin[i],optimize.vpused[i]);
        }

        lgBAD = cloudy();
        if( lgBAD )
        {
                fprintf( ioQQQ, " PROBLEM Cloudy returned error condition - what happened?\n" );
        }

        if( grid.lgGrid )
        {
                /* this is the function's return value */
                chisq = 0.;
        }
        else
        {
                /* this branch optimizing, not grid 
                / * extract line fluxes and compare with observations */
                chisq = 0.0;
                for( i=0; i < MAXCAT; i++ )
                {
                        nobs_cat[i] = 0;
                        chi2_cat[i] = 0.0;
                }

                if( LineSave.ipNormWavL < 0 )
                {
                        fprintf( ioQQQ, 
                                " Normalization line array index is bad.  What has gone wrong?\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                if( (snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) == 0. )
                {
                        fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity.  What has gone wrong?\n" );
                        fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* first print all warnings */
                cdWarnings(ioQQQ);

                /* cycle through the observational values */
                nfound = 0;

                /* first is to optimize relative emission line spectrum */
                if( optimize.lgOptLin )
                {
                        /* set pointers to all optimized lines if first call */
                        if( !lgLinSet )
                        {
                                lgHIT = true;
                                lgLinSet = true;
                                for( i=0; i < optimize.nlobs; i++ )
                                {
                                        double temp1, temp2;
                                        cap4(chFind , (char*)optimize.chLineLabel[i]);
                                        /* j = 0; */

                                        /* >> chng 06 may 04, use cdLine instead of ad hoc treatment.
                                         * no need to complain, cdLine will do it automatically.  */
                                        /* this is an intensity, get the line, returns false if could not find it */
                                        j=cdLine( chFind , optimize.wavelength[i] , &temp1, &temp2 );
                                        if( j<=0 )
                                        {
                                                fprintf( ioQQQ, "\n" );
                                                lgHIT = false;
                                        }
                                        else
                                        {
                                                ipobs[i] = j;
                                        }
                                }

                                /* we did not find the line */
                                if( !lgHIT )
                                {
                                        fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
                                        fprintf( ioQQQ, " Sorry.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }

                        for( i=0; i < 10; i++ )
                        {
                                optimize.SavGenericData[i] = 0.;
                        }

                        for( i=0; i < optimize.nlobs; i++ )
                        {
                                /* and find corresponding model value by straight search */
                                nfound += 1;
                                scld = (double)LineSv[ipobs[i]].sumlin[LineSave.lgLineEmergent]/(double)snorm*LineSave.ScaleNormLine;
                                chi1 = chi2_func(scld,(double)optimize.xLineInt_Obs[i],(double)optimize.xLineInt_error[i]);
                                cat = 0;
                                nobs_cat[cat]++;
                                chi2_cat[cat] += chi1;

                                fprintf( ioQQQ, " %4.4s ", 
                                  LineSv[ipobs[i]].chALab);

                                prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength);

                                fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity", 
                                  scld, 
                                  optimize.xLineInt_Obs[i], 
                                  optimize.xLineInt_error[i], 
                                  chi1 );

                                fprintf( ioQQQ, "\n" );

                                if( i<10 )
                                {
                                        optimize.SavGenericData[i] = chi1;
                                }
                        }
                }

                /* >>chng 06 may 04, moved this from above so that it only
                 * gets printed if all lines are found */
                /*if( optimize.lgOptLum || optimize.lgOptCol || optimize.lgOptTemp || optimize.lgOptLin )*/
                if( optimize.lgOptimize )
                {
                        fprintf( ioQQQ, "  ID               Model    Observed       error      chi**2     Type\n" );
                }
                else
                {
                        ASSERT( grid.lgGrid );
                }

                /* this is to optimize a mean temperature */
                if( optimize.lgOptTemp )
                {
                        for( i=0; i < optimize.nTempObs; i++ )
                        {
                                if( cdTemp(/*(char*)*/optimize.chTempLab[i],optimize.ionTemp[i], &temp_theory, optimize.chTempWeight[i]) )
                                {
                                        /* did not find column density */
                                        fprintf(ioQQQ," optimizer did not find column density %s %li \n",
                                                optimize.chTempLab[i],optimize.ionTemp[i] );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                nfound += 1;
                                chi1 = chi2_func(temp_theory,(double)optimize.temp_obs[i],(double)optimize.temp_error[i]);
                                cat = 3;
                                nobs_cat[cat]++;
                                chi2_cat[cat] += chi1;

                                fprintf( ioQQQ, " %4.4s %2ld        ",
                                  optimize.chTempLab[i], 
                                  optimize.ionTemp[i] );
                                PrintE82( ioQQQ, temp_theory );
                                fprintf( ioQQQ, "    ");
                                PrintE82( ioQQQ, optimize.temp_obs[i] );
                                fprintf( ioQQQ, "    %.5f   %.2e", 
                                        optimize.temp_error[i],  chi1 );
                                fprintf( ioQQQ, " Temperature\n");
                        }
                }

                /* option to optimize column densities */
                if( optimize.lgOptCol )
                {
                        for( i=0; i < optimize.ncobs; i++ )
                        {
                                if( cdColm((char*)optimize.chColDen_label[i],optimize.ion_ColDen[i], &theocl) )
                                {
                                        /* did not find column density */
                                        fprintf(ioQQQ," optimizer did not find column density %s %li \n",
                                                optimize.chColDen_label[i],optimize.ion_ColDen[i] );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                nfound += 1;
                                chi1 = chi2_func(theocl,(double)optimize.ColDen_Obs[i],(double)optimize.chColDen_error[i]);
                                cat = 1;
                                nobs_cat[cat]++;
                                chi2_cat[cat] += chi1;

                                fprintf( ioQQQ, " %4.4s%6ld%12.4e%12.4e%12.5f%12.2e Column density\n", 
                                  optimize.chColDen_label[i], optimize.ion_ColDen[i], theocl, 
                                  optimize.ColDen_Obs[i], optimize.chColDen_error[i], chi1 );
                        }
                }

                /* option to optimize line flux */
                if( optimize.lgOptLum )
                {
                        nfound += 1;
                        if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0.f )
                        {
                                predin = log10(LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten;
                                help = pow(10.,predin-(double)optimize.optint);
                                chi1 = chi2_func(help,1.,(double)optimize.optier);
                        }
                        else
                        {
                                predin = -999.99999;
                                chi1 = (double)FLT_MAX;
                        }
                        cat = 2;
                        nobs_cat[cat]++;
                        chi2_cat[cat] += chi1;

                        fprintf( ioQQQ, " %4.4s%6f%12.5f%12.5f%12.5f%12.2e Line intensity\n", 
                          LineSv[LineSave.ipNormWavL].chALab, 
                          LineSv[LineSave.ipNormWavL].wavelength, 
                          predin, 
                          optimize.optint, 
                          optimize.optier, 
                          chi1 );
                }

                /* >>chng 06 apr 26, do not have to have line matches if doing grid. */
                if( nfound <= 0 && !grid.lgGrid )
                {
                        fprintf( ioQQQ, " WARNING; no line matches found\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* write out chisquared for this iteration */
                fprintf( ioQQQ, "\n" );
                for( i=0; i < MAXCAT; i++ )
                {
                        if( nobs_cat[i] > 0 )
                        {
                                chisq += chi2_cat[i]/nobs_cat[i];
                                fprintf( ioQQQ, " Category %s #obs.%3ld  Total Chi**2%11.3e  Average Chi**2%11.3e\n",
                                  name_cat[i],nobs_cat[i],chi2_cat[i],chi2_cat[i]/nobs_cat[i] );
                        }
                }
                if( nfound )
                {
                        fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq );
                }
        }

        /* increment nOptimiz, the grid / optimizer counter */
        ++optimize.nOptimiz;

        /* only print this if output has been turned on */
        if( called.lgTalk )
        {
                fprintf( ioQQQ, "\n" );
                for( i=0; i < optimize.nvary; i++ )
                {
                        optimize.vparm[0][i] = (realnum)MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
                        optimize.vparm[0][i] = (realnum)MAX2(optimize.vparm[0][i],optimize.varang[i][0]);
                        param[i] = optimize.vparm[0][i];
                        np = optimize.nvfpnt[i];

                        /* now generate the actual command with parameter,
                         * there will be from 1 to 3 numbers on the line */
                        if( optimize.nvarxt[i] == 1 )
                        {
                                /* case with 1 parameter */
                                sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] );
                        }

                        else if( optimize.nvarxt[i] == 2 )
                        {
                                /* case with 2 parameter */
                                sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]);
                        }

                        else if( optimize.nvarxt[i] == 3 )
                        {
                                /* case with 3 parameter */
                                sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 
                                        optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] );
                        }

                        else if( optimize.nvarxt[i] == 4 )
                        {
                                /* case with 4 parameter */
                                sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 
                                        optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] );
                        }

                        else if( optimize.nvarxt[i] == 5 )
                        {
                                /* case with 5 parameter */
                                sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 
                                        optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i],
                                        optimize.vparm[3][i] , optimize.vparm[4][i]);
                        }

                        else
                        {
                                fprintf(ioQQQ,"The number of variable options on this line makes no sense to me4\n");
                                cdEXIT(EXIT_FAILURE);
                        }

                        fprintf( ioQQQ, " Varied command: %s\n", 
                          input.chCardSav[np] );
                }
        }

        func_v = MIN2(chisq,(double)FLT_MAX);
        return( func_v );
}

/* ============================================================================== */
STATIC double chi2_func(double ymodl,
        double ymeas,
        double yerr)
{
        double chi2_func_v,
                temp;

        DEBUG_ENTRY( "chi2_func()" );

        /* compute chi**2 by comparing model quantity ymodl with a measured
         * quantity ymeas with relative error yerr (negative means upper limit)
         */

        if( ymeas <= 0. )
        {
                fprintf( ioQQQ, "chi2_func: non-positive observed quantity, this should not happen\n" );
                cdEXIT(EXIT_FAILURE);
        }

        if( yerr > 0. )
        {
                if( ymodl > 0. )
                {
                        temp = POW2((ymodl-ymeas)/(MIN2(ymodl,ymeas)*yerr));
                        chi2_func_v = MIN2( temp , (double)FLT_MAX );
                }
                else
                        chi2_func_v = (double)FLT_MAX;
        }
        else if( yerr < 0. )
        {
                /* value quoted is an upper limit, so add to chisq
                 * only if limit exceeded, otherwise return zero.
                 */
                if( ymodl > ymeas )
                {
                        temp = POW2((ymodl-ymeas)/(ymeas*yerr));
                        chi2_func_v = MIN2(temp,(double)FLT_MAX);
                }
                else
                        chi2_func_v = 0.;
        }
        else
        {
                fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
                cdEXIT(EXIT_FAILURE);
        }
        return chi2_func_v;
}

---