/home66/gary/public_html/cloudy/c08_branch/source/optimize_do.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 */
/*lgOptimize_do main driver for optimization runs*/
#include "cddefines.h"
#define NPLXMX  (LIMPAR*(LIMPAR+6)+1)
#include "input.h"
#include "called.h"
#include "prt.h"
#include "punch.h"
#include "optimize.h"

/* called by cdDrive, this returns false if things went ok, true for disaster */
bool lgOptimize_do(void)
{
        long int i, 
          iflag, 
          ii, 
          iworke[NPLXMX], 
          j, 
          mode, 
          need, 
          nfe, 
          np;
        realnum fret, 
          fx, 
          param[LIMPAR], 
          ptem[LIMPAR], 
          delta[LIMPAR], 
          toler,
          worke[NPLXMX];
/*      double _e0[210]; */
/*      realnum (*const p)[21] = (realnum(*)[21])_e0; */
/*      realnum (*const xi)[20] = (realnum(*)[20])_e0; */

        DEBUG_ENTRY( "lgOptimize_do()" );

        /* main driver for optimization runs
         * Drives cloudy to optimize variables;*/

        /* code originally written by R.F. Carswell, IOA Cambridge */

        toler = (realnum)log10(1. + optimize.OptGlobalErr);

        /*>>chng 07 feb 19, rm powell's method
         * this is phymir */
        if( strcmp(optimize.chOptRtn,"PHYM") == 0 )
        {
                /* Phymir method */
                for( j=0; j < optimize.nvary; j++ )
                {
                        ptem[j] = optimize.vparm[0][j];
                        delta[j] = optimize.vincr[j];
                }
                /* >>chng 06 jan 02, fix uninitialized var problem detected by valgrind/purify, PvH */
                for( j=optimize.nvary; j < LIMPAR; j++ )
                {
                        ptem[j] = -FLT_MAX;
                        delta[j] = -FLT_MAX;
                }
                optimize_phymir(ptem,delta,optimize.nvary,&fret,toler);
                for( j=0; j < optimize.nvary; j++ )
                {
                        optimize.vparm[0][j] = ptem[j];
                }

        }

        else if( strcmp(optimize.chOptRtn,"SUBP") == 0 )
        {
                fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" );
                need = 2*optimize.nvary + optimize.nvary*(optimize.nvary + 4) + 1;
                if( need > NPLXMX )
                {
                        fprintf( ioQQQ, " Increase size of NPLXMX in parameter statements to handle this many variables.\n" );
                        fprintf( ioQQQ, " I need at least %5ld\n", need );
                        cdEXIT(EXIT_FAILURE);
                }
                for( j=0; j < optimize.nvary; j++ )
                {
                        ptem[j] = optimize.vparm[0][j];
                }

                /* The subroutine SUBPLX input into cloudy 8/4/94.
                 * The program itself is very well commented.
                 * The mode must set to 0 for the default values.
                 * The switch iflag tells if the program terminated normally.   */
                mode = 0;

                /*  >>chng 97 dec 08, remove first arg, optimize_func, since not used in routines */
                optimize_subplex(
                        /* the number of parameters to vary */
                        optimize.nvary,
                        /* the relative error, single number */
                        toler,
                        /* maximum number of function evaluations before giving up */
                        optimize.nIterOptim,
                        /* mode of operation, we simply set to zero */
                        mode,
                        /* the initial changes in the guessed best coefficients, typically 0.2 to 1  */
                        optimize.vincr,
                        /* a vector of nvary initial parameters that are the starting guesses for the parameters */
                        ptem,
                        /* a realnum, this is simply ignored */
                        &fx,
                        /* another parameter that is simply ignored, a long int */
                        &nfe,
                        /* a realnum that is NPLXMX long, used for working space by the routine */
                        /* an array that is 20*26 + 1 elements long, used for working space */
                        worke,
                        /* a long int that is NPLXMX long, used for working space by the routine */
                        /* an array that is 20*26 + 1 elements long, used for working space */
                        iworke,
                        /* a long int - says what happened, if -1 then exceeded nIterOptim iteration */
                        &iflag);

                if( iflag == -1 )
                {
                        fprintf( ioQQQ, " SUBPLEX exceeding maximum iterations.\n This can be reset with the OPTIMZE ITERATIONS command.\n" );
                }

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

                if( optimize.lgOptimFlow )
                {
                        fprintf( ioQQQ, " trace return optimize_subplex:\n" );
                        for( j=0; j < optimize.nvary; j++ )
                        {
                                fprintf( ioQQQ, " Values:" );
                                for( ii=1; ii <= optimize.nvarxt[j]; ii++ )
                                {
                                        fprintf( ioQQQ, "%10.2e", optimize.vparm[ii-1][j] );
                                }
                                fprintf( ioQQQ, "\n" );
                        }
                }
        }
        else
                TotalInsanity();
        /*>>chng 07 feb 08m rm optimize_amoeba due to Robin & Peter license
         * concerns - R837*/

        fprintf( ioQQQ, " **************************************************\n" );
        fprintf( ioQQQ, " **************************************************\n" );
        fprintf( ioQQQ, " **************************************************\n" );
        fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", optimize.nOptimiz );

        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 me.\n");
                        cdEXIT(EXIT_FAILURE);
                }


                fprintf( ioQQQ, " Optimal command: %s\n", input.chCardSav[np]);
                fprintf( ioQQQ, "  Smallest value:%10.2e Largest value:%10.2e Allowed range %10.2e to %10.2e\n", 
                  optimize.varmin[i], optimize.varmax[i], optimize.varang[i][0], 
                  optimize.varang[i][1] );
        }

        called.lgTalk = true;
        called.lgTalkIsOK = true;
        prt.lgFaintOn = true;

        /* though a page eject */
        fprintf( ioQQQ, "\f" );

        /* punch optimal parameters on unit ioOptim */
        if( optimize.ioOptim == NULL )
        {
                /* open default file name, optimal.in */
                optimize.ioOptim = open_data( chOptimFileName, "w", AS_LOCAL_ONLY );
        }

        for( i=0; i <= input.nSave; i++ )
        {
                fprintf( optimize.ioOptim, "%s\n", input.chCardSav[i]);
        }
        fclose( optimize.ioOptim );

        /* recalculate values in cloudy for the best fit, and print out
         * all the information */
        fret = (realnum)optimize_func(param);


        if( lgAbort )
        {
                /* busted set means there were serious problems somewhere */
                return true;
        }
        else
        {
                /* return 0 if everything is ok */
                return false;
        }
}

---