cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*lgOptimize_do main driver for optimization runs*/ 00004 #include "cddefines.h" 00005 #define NPLXMX (LIMPAR*(LIMPAR+6)+1) 00006 #include "input.h" 00007 #include "called.h" 00008 #include "prt.h" 00009 #include "punch.h" 00010 #include "optimize.h" 00011 00012 /* called by cdDrive, this returns false if things went ok, true for disaster */ 00013 bool lgOptimize_do(void) 00014 { 00015 long int i, 00016 iflag, 00017 ii, 00018 iworke[NPLXMX], 00019 j, 00020 mode, 00021 need, 00022 nfe, 00023 np; 00024 realnum fret, 00025 fx, 00026 param[LIMPAR], 00027 ptem[LIMPAR], 00028 delta[LIMPAR], 00029 toler, 00030 worke[NPLXMX]; 00031 /* double _e0[210]; */ 00032 /* realnum (*const p)[21] = (realnum(*)[21])_e0; */ 00033 /* realnum (*const xi)[20] = (realnum(*)[20])_e0; */ 00034 00035 DEBUG_ENTRY( "lgOptimize_do()" ); 00036 00037 /* main driver for optimization runs 00038 * Drives cloudy to optimize variables;*/ 00039 00040 /* code originally written by R.F. Carswell, IOA Cambridge */ 00041 00042 toler = (realnum)log10(1. + optimize.OptGlobalErr); 00043 00044 /*>>chng 07 feb 19, rm powell's method 00045 * this is phymir */ 00046 if( strcmp(optimize.chOptRtn,"PHYM") == 0 ) 00047 { 00048 /* Phymir method */ 00049 for( j=0; j < optimize.nvary; j++ ) 00050 { 00051 ptem[j] = optimize.vparm[0][j]; 00052 delta[j] = optimize.vincr[j]; 00053 } 00054 /* >>chng 06 jan 02, fix uninitialized var problem detected by valgrind/purify, PvH */ 00055 for( j=optimize.nvary; j < LIMPAR; j++ ) 00056 { 00057 ptem[j] = -FLT_MAX; 00058 delta[j] = -FLT_MAX; 00059 } 00060 optimize_phymir(ptem,delta,optimize.nvary,&fret,toler); 00061 for( j=0; j < optimize.nvary; j++ ) 00062 { 00063 optimize.vparm[0][j] = ptem[j]; 00064 } 00065 00066 } 00067 00068 else if( strcmp(optimize.chOptRtn,"SUBP") == 0 ) 00069 { 00070 fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" ); 00071 need = 2*optimize.nvary + optimize.nvary*(optimize.nvary + 4) + 1; 00072 if( need > NPLXMX ) 00073 { 00074 fprintf( ioQQQ, " Increase size of NPLXMX in parameter statements to handle this many variables.\n" ); 00075 fprintf( ioQQQ, " I need at least %5ld\n", need ); 00076 cdEXIT(EXIT_FAILURE); 00077 } 00078 for( j=0; j < optimize.nvary; j++ ) 00079 { 00080 ptem[j] = optimize.vparm[0][j]; 00081 } 00082 00083 /* The subroutine SUBPLX input into cloudy 8/4/94. 00084 * The program itself is very well commented. 00085 * The mode must set to 0 for the default values. 00086 * The switch iflag tells if the program terminated normally. */ 00087 mode = 0; 00088 00089 /* >>chng 97 dec 08, remove first arg, optimize_func, since not used in routines */ 00090 optimize_subplex( 00091 /* the number of parameters to vary */ 00092 optimize.nvary, 00093 /* the relative error, single number */ 00094 toler, 00095 /* maximum number of function evaluations before giving up */ 00096 optimize.nIterOptim, 00097 /* mode of operation, we simply set to zero */ 00098 mode, 00099 /* the initial changes in the guessed best coefficients, typically 0.2 to 1 */ 00100 optimize.vincr, 00101 /* a vector of nvary initial parameters that are the starting guesses for the parameters */ 00102 ptem, 00103 /* a realnum, this is simply ignored */ 00104 &fx, 00105 /* another parameter that is simply ignored, a long int */ 00106 &nfe, 00107 /* a realnum that is NPLXMX long, used for working space by the routine */ 00108 /* an array that is 20*26 + 1 elements long, used for working space */ 00109 worke, 00110 /* a long int that is NPLXMX long, used for working space by the routine */ 00111 /* an array that is 20*26 + 1 elements long, used for working space */ 00112 iworke, 00113 /* a long int - says what happened, if -1 then exceeded nIterOptim iteration */ 00114 &iflag); 00115 00116 if( iflag == -1 ) 00117 { 00118 fprintf( ioQQQ, " SUBPLEX exceeding maximum iterations.\n This can be reset with the OPTIMZE ITERATIONS command.\n" ); 00119 } 00120 00121 for( j=0; j < optimize.nvary; j++ ) 00122 { 00123 optimize.vparm[0][j] = ptem[j]; 00124 } 00125 00126 if( optimize.lgOptimFlow ) 00127 { 00128 fprintf( ioQQQ, " trace return optimize_subplex:\n" ); 00129 for( j=0; j < optimize.nvary; j++ ) 00130 { 00131 fprintf( ioQQQ, " Values:" ); 00132 for( ii=1; ii <= optimize.nvarxt[j]; ii++ ) 00133 { 00134 fprintf( ioQQQ, "%10.2e", optimize.vparm[ii-1][j] ); 00135 } 00136 fprintf( ioQQQ, "\n" ); 00137 } 00138 } 00139 } 00140 else 00141 TotalInsanity(); 00142 /*>>chng 07 feb 08m rm optimize_amoeba due to Robin & Peter license 00143 * concerns - R837*/ 00144 00145 fprintf( ioQQQ, " **************************************************\n" ); 00146 fprintf( ioQQQ, " **************************************************\n" ); 00147 fprintf( ioQQQ, " **************************************************\n" ); 00148 fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", optimize.nOptimiz ); 00149 00150 for( i=0; i < optimize.nvary; i++ ) 00151 { 00152 optimize.vparm[0][i] = (realnum)MIN2(optimize.vparm[0][i],optimize.varang[i][1]); 00153 optimize.vparm[0][i] = (realnum)MAX2(optimize.vparm[0][i],optimize.varang[i][0]); 00154 param[i] = optimize.vparm[0][i]; 00155 np = optimize.nvfpnt[i]; 00156 00157 /* now generate the actual command with parameter, 00158 * there will be from 1 to 3 numbers on the line */ 00159 if( optimize.nvarxt[i] == 1 ) 00160 { 00161 /* case with 1 parameter */ 00162 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] ); 00163 } 00164 00165 else if( optimize.nvarxt[i] == 2 ) 00166 { 00167 /* case with 2 parameter */ 00168 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]); 00169 } 00170 00171 else if( optimize.nvarxt[i] == 3 ) 00172 { 00173 /* case with 3 parameter */ 00174 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 00175 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] ); 00176 } 00177 00178 else if( optimize.nvarxt[i] == 4 ) 00179 { 00180 /* case with 4 parameter */ 00181 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 00182 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] ); 00183 } 00184 00185 else if( optimize.nvarxt[i] == 5 ) 00186 { 00187 /* case with 5 parameter */ 00188 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 00189 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] , 00190 optimize.vparm[3][i] , optimize.vparm[4][i]); 00191 } 00192 00193 else 00194 { 00195 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me.\n"); 00196 cdEXIT(EXIT_FAILURE); 00197 } 00198 00199 00200 fprintf( ioQQQ, " Optimal command: %s\n", input.chCardSav[np]); 00201 fprintf( ioQQQ, " Smallest value:%10.2e Largest value:%10.2e Allowed range %10.2e to %10.2e\n", 00202 optimize.varmin[i], optimize.varmax[i], optimize.varang[i][0], 00203 optimize.varang[i][1] ); 00204 } 00205 00206 called.lgTalk = true; 00207 called.lgTalkIsOK = true; 00208 prt.lgFaintOn = true; 00209 00210 /* though a page eject */ 00211 fprintf( ioQQQ, "\f" ); 00212 00213 /* punch optimal parameters on unit ioOptim */ 00214 if( optimize.ioOptim == NULL ) 00215 { 00216 /* open default file name, optimal.in */ 00217 optimize.ioOptim = open_data( chOptimFileName, "w", AS_LOCAL_ONLY ); 00218 } 00219 00220 for( i=0; i <= input.nSave; i++ ) 00221 { 00222 fprintf( optimize.ioOptim, "%s\n", input.chCardSav[i]); 00223 } 00224 fclose( optimize.ioOptim ); 00225 00226 /* recalculate values in cloudy for the best fit, and print out 00227 * all the information */ 00228 fret = (realnum)optimize_func(param); 00229 00230 00231 if( lgAbort ) 00232 { 00233 /* busted set means there were serious problems somewhere */ 00234 return true; 00235 } 00236 else 00237 { 00238 /* return 0 if everything is ok */ 00239 return false; 00240 } 00241 }