/home66/gary/public_html/cloudy/c08_branch/source/parse_optimize.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 */
/*ParseOptimize parse the optimize and grid command lines */
/*GetOptColDen read observed column densities & errors for optimizer */
/*GetOptLineInt parse observed line intensities for optimization routines */
/*GetOptTemp read observed temperatures & errors for optimizer */
#include "cddefines.h"
#include "trace.h"
#include "optimize.h"
#include "grid.h"
#include "input.h"
#include "prt.h"
#include "parse.h"
#include "lines_service.h"

static const realnum DEFERR = 0.05f;

/*GetOptLineInt parse observed line intensities for optimization routines */
STATIC void GetOptLineInt(char *chCard );

/*GetOptColDen read observed column densities & errors for optimizer */
STATIC void GetOptColDen(char *chCard );

/*GetOptTemp read observed temperatures & errors for optimizer */
STATIC void GetOptTemp(char *chCard );

/* ParseOptimize - called from ParseCommands if GRID or OPTIMIZE command found */
void ParseOptimize(
        /* command line, which was changed to all caps in main parsing routine */
        char *chCard)
{
        bool lgEOL;
        long int i;

        DEBUG_ENTRY( "ParseOptimize()" );

        /* this must come first so that contents of filename do not trigger wrong command */
        if( nMatch("FILE",chCard) )
        {
                /* option to send final set of parameters to an input file 
                 * get name within double quotes, 
                 * and set to blanks in chCard and OrgCard */
                /* file handle set NULL in optimize_do when first set up,
                 * this routine called many times during optimization process,
                 * only want to open file one time */
                if( optimize.ioOptim == NULL )
                {
                        /* chCard is all caps at this point.  GetQuote will work with 
                         * original version of command line, which preserves case of
                         * characters.  Also removes string between quotes */
                        if( GetQuote( chOptimFileName , chCard , true ) )
                                /* this can't happen since GetQuote exits if quote not
                                 * found, since true sets abort mode */
                                TotalInsanity();

                        /* open the file */
                        optimize.ioOptim = open_data( chOptimFileName, "w", AS_LOCAL_ONLY );
                }
        }

        else if( nMatch("COLU",chCard) )
        {
                /* optimize column density */
                optimize.lgOptCol = true;
                optimize.lgOptimize = true;

                /* read column densities to match */
                GetOptColDen(chCard);
        }

        else if( nMatch("CONT",chCard) )
        {
                /* set flag saying that optimization should start from continue file */
                optimize.lgOptCont = true;
                optimize.lgOptimize = true;
        }

        /* >>chng 06 sep 4, exclude lines with "GRID", handled below. */
        else if( nMatch("INCR",chCard) && !nMatch("GRID",chCard) )
        {
                /* scan off increments for the previously selected parameter */
                if( optimize.nparm > 0 )
                {
                        /* also called during optimization process, ignore then */
                        i = 5;
                        optimize.OptIncrm[optimize.nparm-1] = 
                                (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                }
        }

        else if( nMatch("LUMI",chCard) || nMatch("INTE",chCard) )
        {
                /* scan off intensity or luminosity of normalization line */
                i = 5;
                optimize.optint = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                optimize.optier = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        optimize.optier = DEFERR;

                /* set flag to say that intensity or luminosity of line set */
                optimize.lgOptLum = true;
                optimize.lgOptimize = true;
        }

        else if( nMatch("ITER",chCard) )
        {
                /* scan off number of iterations */
                i = 5;
                optimize.nIterOptim = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
        }

        else if( nMatch("LINE",chCard) )
        {
                /* read lines to match */
                GetOptLineInt(chCard);

                /* set flag saying that something has been set */
                optimize.lgOptLin = true;
                optimize.lgOptimize = true;
        }

        else if( nMatch("PHYM",chCard) )
        {
                /* use PvH's PHYMIR to optimize parameters */
                strcpy( optimize.chOptRtn, "PHYM" );
#               ifdef __unix
                optimize.lgParallel = ! nMatch("SEQU",chCard);
#               else
                optimize.lgParallel = false;
#               endif
                if( optimize.lgParallel ) 
                {
                        i = 5;
                        long dum = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        /* default is the total number of available CPUs */
                        optimize.useCPU = lgEOL ? cpu.nCPU() : dum;
                }
                else 
                {
                        optimize.useCPU = 1;
                }
        }

        /* >>chng 06 aug 22, grid now done differently, requires two range parameters
         * and the number of points, parsed below as "GRID". */
        else if( nMatch("RANG",chCard) && !nMatch("GRID",chCard)  )
        {
                /* scan off range for the previously selected variable */
                if( optimize.nparm > 0 )
                {
                        bool lgFirstOneReal = false;

                        i = 5;
                        optimize.varang[optimize.nparm-1][0] = 
                                (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                optimize.varang[optimize.nparm-1][0] = -1e38f;
                        }
                        else
                        {
                                lgFirstOneReal = true;
                        }

                        optimize.varang[optimize.nparm-1][1] = 
                                (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                optimize.varang[optimize.nparm-1][1] = 1e38f;
                        }
                        else if( lgFirstOneReal )
                        {
                                /* >>chng 06 aug 22, swap if second range parameter is less than the first,
                                 * and always increment optimize.nRangeSet */
                                ++optimize.nRangeSet;
                                if( optimize.varang[optimize.nparm-1][1] < optimize.varang[optimize.nparm-1][0] )
                                {
                                        realnum temp = optimize.varang[optimize.nparm-1][0];
                                        optimize.varang[optimize.nparm-1][0] = optimize.varang[optimize.nparm-1][1];
                                        optimize.varang[optimize.nparm-1][1] = temp;
                                }
                        }
                }
        }

        else if( nMatch("SUBP",chCard) )
        {
                /* use subplex to optimize parameters */
                strcpy( optimize.chOptRtn, "SUBP" );
        }

        /* match a temperature */
        else if( nMatch("TEMP",chCard) )
        {
                /* read temperatures to match */
                GetOptTemp(chCard);

                /* set flag saying that optimize temps has been set */
                optimize.lgOptTemp = true;
                optimize.lgOptimize = true;
        }

        else if( nMatch("TOLE",chCard) )
        {
                /* scan off tolerance of fit, sum of residuals must be smaller than this
                 * default is 0.10 */
                i = 5;
                optimize.OptGlobalErr = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
        }

        else if( nMatch("TRAC",chCard) )
        {
                if( nMatch("STAR",chCard) )
                {
                        /* trace start iteration number
                         * turn on trace printout starting on nth call to cloudy */
                        i = 5;
                        optimize.nTrOpt = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " optimize trace start command:\n" );
                                fprintf( ioQQQ, " The iteration number must appear.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        optimize.lgTrOpt = true;
                }
                else if( nMatch("FLOW",chCard) )
                {
                        /* trace flow
                         * follow logical flow within code */
                        optimize.lgOptimFlow = true;
                }
                else
                {
                        fprintf( ioQQQ, " optimize trace flow command:\n" );
                        fprintf( ioQQQ, " One of the sub keys START or FLOW must appear.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else
        {
                fprintf( ioQQQ, " ==%s== is unrecognized keyword, consult HAZY.\n", chCard  );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

/*GetOptColDen read observed column densities & errors for optimizer */
STATIC void GetOptColDen(char *chCard )
{
        char chCap[INPUT_LINE_LENGTH];
        bool lgEOF, 
          lgEOL;
        long int i;

        DEBUG_ENTRY( "GetOptColDen()" );

        /* read observed column densities & errors */
        optimize.ncobs = 0;

        /* get first line */
        input_readarray(chCard,&lgEOF);
        if( lgEOF )
        {
                fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* now copy this line to chCap, then convert to caps */
        strcpy( chCap, chCard );
        caps(chCap);

        while( !lgEOF )
        {
                if( optimize.ncobs > NCOLLM )
                {
                        fprintf( ioQQQ, " Too many column densities have been entered; the limit is%4ld.  Increase variable NCOLLM.\n", 
                          NCOLLM );
                        cdEXIT(EXIT_FAILURE);
                }

                /* order on line is element label (col 1-4), ionization stage, column density, err */
                /* copy cap'd version of first 4 char of chCard to chColDen_label */
                cap4((char*)optimize.chColDen_label[optimize.ncobs] , chCard );

                /* now get the ion stage, this should be 1 for atom, up to element
                 * number plus one */
                i = 5;
                optimize.ion_ColDen[optimize.ncobs] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " %s\n", chCard );
                        fprintf( ioQQQ, " The ionization stage MUST appear on this line.  Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* the ion must be 1 or greater unless requesting a special,
                 * like a molecule or excited state population, in which
                 * case ion = 0
                 * can't check on upper limit yet since have not resolved element name */
                if( optimize.ion_ColDen[optimize.ncobs] < 0 )
                {
                        fprintf( ioQQQ, " %s\n", chCard );
                        fprintf( ioQQQ, " An ionization stage of%4ld does not make sense.  Sorry.\n", 
                          optimize.ion_ColDen[optimize.ncobs] );
                        cdEXIT(EXIT_FAILURE);
                }

                optimize.ColDen_Obs[optimize.ncobs] = (realnum)pow(10.,FFmtRead(chCard,&i,
                  INPUT_LINE_LENGTH,&lgEOL));
                if( lgEOL )
                {
                        fprintf( ioQQQ, " %80.80s\n", chCard );
                        fprintf( ioQQQ, " An observed column density MUST be entered.  Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                optimize.chColDen_error[optimize.ncobs] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( optimize.chColDen_error[optimize.ncobs] <= 0.0 )
                {
                        /* this is the relative error allowed */
                        optimize.chColDen_error[optimize.ncobs] = (realnum)DEFERR;
                }

                /* check if number is a limit - if '<' appears on the line then it is */
                if( strchr( chCard , '<' ) != NULL )
                {
                        /* value is an upper limit only, use negative error to flag */
                        optimize.chColDen_error[optimize.ncobs] = -optimize.chColDen_error[optimize.ncobs];
                }

                input_readarray(chCard,&lgEOF);
                strcpy( chCap, chCard );
                caps(chCap);
                if( lgEOF )
                {
                        fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                if( strncmp( chCap , "END" , 3) == 0 )
                {
                        lgEOF = true;
                }

                /* now increment the number of columns we have entered */
                optimize.ncobs += 1;
        }

        if( trace.lgTrace && optimize.lgTrOpt )
        {
                fprintf( ioQQQ, "%4ld columns were entered, they were;\n", 
                  optimize.ncobs );
                for( i=0; i < optimize.ncobs; i++ )
                {
                        fprintf( ioQQQ, " %4.4s ion=%5ld%10.2e%10.2e\n", 
                          optimize.chColDen_label[i], optimize.ion_ColDen[i], optimize.ColDen_Obs[i], 
                          optimize.chColDen_error[i] );
                }
        }
        return;
}

/*GetOptLineInt parse observed line intensities for optimization routines */
STATIC void GetOptLineInt(char *chCard )
{
        char chCap[INPUT_LINE_LENGTH];

        bool lgEOF, 
          lgEOL;
        long int i;

        DEBUG_ENTRY( "GetOptLineInt()" );

        /* read observed line fluxes & errors */
        optimize.nlobs = 0;

        input_readarray(chCard,&lgEOF);
        if( lgEOF )
        {
                fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        strcpy( chCap, chCard );
        caps(chCap);

        while( !lgEOF )
        {
                if( optimize.nlobs >= NOBSLM )
                {
                        fprintf( ioQQQ, 
                                " Too many lines have been entered; the limit is %ld.  Increase variable NOBSLM.\n", 
                          NOBSLM );
                        cdEXIT(EXIT_FAILURE);
                }

                /* order on line is label (col 1-4), wavelength, flux, error */
                strncpy( optimize.chLineLabel[optimize.nlobs], chCard , 4 );
                /* null terminate the label*/
                optimize.chLineLabel[optimize.nlobs][4] = 0;

                i = 5;
                /* next get the wavelength */
                optimize.wavelength[optimize.nlobs] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                /* now convert wavelength to angstroms */
                /* line was entered, look for possible micron or cm label */
                if( input.chCARDCAPS[i-1] == 'M' )
                {
                        /* microns */
                        optimize.wavelength[optimize.nlobs] *= 1e4f;
                }
                else if( input.chCARDCAPS[i-1] == 'C' )
                {
                        /* microns */
                        optimize.wavelength[optimize.nlobs] *= 1e8f;
                }

                /* get the error associated with 4 significant figures */
                optimize.errorwave[optimize.nlobs] = 
                        WavlenErrorGet( optimize.wavelength[optimize.nlobs] );

                /* next get the observed intensity */
                optimize.xLineInt_Obs[optimize.nlobs] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " %s\n", chCard );
                        fprintf( ioQQQ, " The wavelength and relative intensity MUST be entered on this line.  Sorry.\n" );
                        fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
                        cdEXIT(EXIT_FAILURE);
                }

                if( optimize.xLineInt_Obs[optimize.nlobs] <= 0. )
                {
                        fprintf( ioQQQ, " An observed intensity of %.2e is not allowed.  Sorry.\n", 
                          optimize.xLineInt_Obs[optimize.nlobs] );
                        fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
                        cdEXIT(EXIT_FAILURE);
                }

                /* finally the optional error */
                optimize.xLineInt_error[optimize.nlobs] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                /* most often will use the default */
                if( optimize.xLineInt_error[optimize.nlobs] <= 0.0 )
                {
                        /* this is the relative error allowed */
                        optimize.xLineInt_error[optimize.nlobs] = (realnum)DEFERR;
                }

                /* check if number is a limit - if '<' appears on the line then it is */
                if( strchr( chCard , '<' ) != NULL )
                {
                        /* value is an upper limit only, use negative error to flag */
                        optimize.xLineInt_error[optimize.nlobs] = -optimize.xLineInt_error[optimize.nlobs];
                }

                /* get next line */
                input_readarray(chCard,&lgEOF);
                if( lgEOF )
                {
                        fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                strcpy( chCap, chCard );
                caps(chCap);
                if( strncmp( chCap ,"END" , 3 ) == 0 )
                        lgEOF = true;

                /* finally increment the number of observed lines */
                ++optimize.nlobs;
        }

        if( trace.lgTrace && trace.lgTrOptm )
        {
                fprintf( ioQQQ, "%4ld lines were entered, they were;\n", 
                  optimize.nlobs );

                for( i=0; i < optimize.nlobs; i++ )
                {
                        fprintf( ioQQQ, " %4.4s ", optimize.chLineLabel[i] );
                        prt_wl( ioQQQ, optimize.wavelength[i] );

                        fprintf( ioQQQ, " %10.2e%10.2e\n", 
                                optimize.xLineInt_Obs[i], 
                                optimize.xLineInt_error[i] );
                }
        }
        return;
}

/*GetOptTemp parse observed line intensities for optimization routines */
STATIC void GetOptTemp(char *chCard )
{
        char chCap[INPUT_LINE_LENGTH];

        bool lgEOF, 
          lgEOL;
        long int i;

        DEBUG_ENTRY( "GetOptTemp()" );

        /* read observed line fluxes & errors - first set total number of observe temps */
        optimize.nTempObs = 0;

        input_readarray(chCard,&lgEOF);
        if( lgEOF )
        {
                fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* make line caps so we can parse it */
        strcpy( chCap, chCard );
        caps(chCap);

        while( !lgEOF )
        {
                if( optimize.nTempObs >= NOBSLM )
                {
                        fprintf( ioQQQ, 
                                " Too many temperatures have been entered; the limit is %ld.  Increase variable NOBSLM.\n", 
                          NOBSLM );
                        cdEXIT(EXIT_FAILURE);
                }

                /* order on line is label (col 1-4), ion, temperature, error */
                strncpy( optimize.chTempLab[optimize.nTempObs], chCard , 4 );
                /* null terminate the label*/
                optimize.chTempLab[optimize.nTempObs][4] = 0;

                i = 5;
                /* next get the ion stage */
                optimize.ionTemp[optimize.nTempObs] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                /* next get the observed temperature */
                optimize.temp_obs[optimize.nTempObs] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " %s\n", chCard );
                        fprintf( ioQQQ, " The ion stage and temperature MUST be entered on this line.  Sorry.\n" );
                        fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
                        cdEXIT(EXIT_FAILURE);
                }

                /* temperatures less than or equal to 10 are logs */
                if( optimize.temp_obs[optimize.nTempObs] <= 10. )
                {
                        optimize.temp_obs[optimize.nTempObs] = (realnum)pow( 10. , (double)optimize.temp_obs[optimize.nTempObs] );
                }

                /* finally the optional error */
                optimize.temp_error[optimize.nTempObs] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                /* most often will use the default */
                if( optimize.temp_error[optimize.nTempObs] <= 0.0 )
                {
                        /* this is the relative error allowed */
                        optimize.temp_error[optimize.nTempObs] = (realnum)DEFERR;
                }

                /* check if number is a limit - if '<' appears on the line then it is */
                if( strchr( chCard , '<' ) != NULL )
                {
                        /* value is an upper limit only, use negative error to flag */
                        optimize.temp_error[optimize.nTempObs] = -optimize.temp_error[optimize.nTempObs];
                }

                /* check for radius or volume for how to weight the mean temp 
                 * this will be the default */
                strcpy( optimize.chTempWeight[optimize.nTempObs] , "radius" );
                /* >>chng 05 dec 29, from chCard to chCap, unlike much of code in this file,
                 * chCard has been read in by this routine and contains the original form of 
                 * the command line.  It was converted to caps and stored in chCap above.
                 * As it was written only VOLUME was matched, would not match volume.
                 * Bug caught and corrected by Bohdan Melekh */
                /*if( nMatch( "VOLUME" , chCard ) )*/
                if( nMatch( "VOLUME" , chCap ) )
                {
                        strcpy( optimize.chTempWeight[optimize.nTempObs] , "volume" );
                }

                /* get next line */
                input_readarray(chCard,&lgEOF);
                if( lgEOF )
                {
                        fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                strcpy( chCap, chCard );
                caps(chCap);
                if( strncmp( chCap ,"END" , 3 ) == 0 )
                        lgEOF = true;

                /* finally increment the number of observed lines */
                ++optimize.nTempObs;
        }

        if( trace.lgTrace && trace.lgTrOptm )
        {
                fprintf( ioQQQ, "%4ld temperatures were entered, they were;\n", 
                  optimize.nTempObs );

                for( i=0; i < optimize.nTempObs; i++ )
                {
                        fprintf( ioQQQ, " %4.4s ", optimize.chTempLab[i] );
                        fprintf( ioQQQ, " %li " , optimize.ionTemp[i] );

                        fprintf( ioQQQ, " %.2e %.2e\n", 
                                optimize.temp_obs[i], 
                                optimize.temp_error[i] );
                }
        }
        return;
}

---