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 /*optimize_func actual function called during evaluation of optimization run */ 00004 #include "cddefines.h" 00005 #include "init.h" 00006 #include "lines.h" 00007 #include "prt.h" 00008 #include "called.h" 00009 #include "radius.h" 00010 #include "input.h" 00011 #include "cloudy.h" 00012 #include "cddrive.h" 00013 #include "optimize.h" 00014 #include "grid.h" 00015 /* used below to derive chi2 */ 00016 STATIC double chi2_func(double,double,double); 00017 00018 double optimize_func(realnum param[]) 00019 { 00020 00021 #define MAXCAT 4 00022 00023 char /* chCAPLB[5], */ 00024 chFind[5]; 00025 00026 bool lgBAD, 00027 lgHIT, 00028 lgLimOK; 00029 00030 long int cat, 00031 i, 00032 j, 00033 nfound, 00034 nobs_cat[MAXCAT], 00035 np; 00036 00037 static long int ipobs[NOBSLM]; 00038 00039 double chi1, 00040 chi2_cat[MAXCAT], 00041 chisq, 00042 func_v, 00043 help, 00044 predin, 00045 scld, 00046 snorm, 00047 theocl, 00048 temp_theory; 00049 00050 static char name_cat[MAXCAT][13] = 00051 { 00052 "rel flux ", 00053 "column dens ", 00054 "abs flux ", 00055 "mean Temp " 00056 }; 00057 00058 static bool lgLinSet = false; 00059 00060 DEBUG_ENTRY( "optimize_func()" ); 00061 00062 /* This routine is called by optimizer with values of the 00063 * variable parameters for CLOUDY in the array p(i). It returns 00064 * the value FUNC = SUM (obs-model)**2/sig**2 for the lines 00065 * specified in the observational data file, values held in the 00066 * common blocks /OBSLIN/ & /OBSINT/ 00067 * replacement input strings for CLOUDY READR held in /chCardSav/ 00068 * parameter information for setting chCardSav held in /parmv/ 00069 * additional variables 00070 * Gary's variables 00071 */ 00072 00073 /* zero out lots of variables */ 00074 zero(); 00075 00076 if( optimize.lgOptimFlow ) 00077 { 00078 fprintf( ioQQQ, " trace, optimize_func variables" ); 00079 for( i=0; i < optimize.nvary; i++ ) 00080 { 00081 fprintf( ioQQQ, "%.2e", param[i] ); 00082 } 00083 fprintf( ioQQQ, "\n" ); 00084 } 00085 00086 for( i=0; i < optimize.nvary; i++ ) 00087 { 00088 optimize.vparm[0][i] = param[i]; 00089 } 00090 00091 /* call routine to pack variables with appropriate 00092 * CLOUDY input lines given the array of variable parameters p(i) */ 00093 vary_input(&lgLimOK); 00094 00095 if( !lgLimOK ) 00096 { 00097 /* these parameters are not within limits of parameter search 00098 * >>chng 96 apr 26, as per Peter van Hoof comment */ 00099 fprintf( ioQQQ, " Iteration %ld not within range.\n", 00100 optimize.nOptimiz ); 00101 00102 /* this is error; very bad since not within range of parameters */ 00103 func_v = (double)FLT_MAX; 00104 00105 /* always increment nOptimiz, even if parameters are out of bounds, 00106 * this prevents optimizer to get stuck in infinite loop */ 00107 ++optimize.nOptimiz; 00108 return( func_v ); 00109 } 00110 00111 for( i=0; i < optimize.nvary; i++ ) 00112 { 00113 optimize.varmax[i] = (realnum)MAX2(optimize.varmax[i],optimize.vpused[i]); 00114 optimize.varmin[i] = (realnum)MIN2(optimize.varmin[i],optimize.vpused[i]); 00115 } 00116 00117 lgBAD = cloudy(); 00118 if( lgBAD ) 00119 { 00120 fprintf( ioQQQ, " PROBLEM Cloudy returned error condition - what happened?\n" ); 00121 } 00122 00123 if( grid.lgGrid ) 00124 { 00125 /* this is the function's return value */ 00126 chisq = 0.; 00127 } 00128 else 00129 { 00130 /* this branch optimizing, not grid 00131 / * extract line fluxes and compare with observations */ 00132 chisq = 0.0; 00133 for( i=0; i < MAXCAT; i++ ) 00134 { 00135 nobs_cat[i] = 0; 00136 chi2_cat[i] = 0.0; 00137 } 00138 00139 if( LineSave.ipNormWavL < 0 ) 00140 { 00141 fprintf( ioQQQ, 00142 " Normalization line array index is bad. What has gone wrong?\n" ); 00143 cdEXIT(EXIT_FAILURE); 00144 } 00145 00146 if( (snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) == 0. ) 00147 { 00148 fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity. What has gone wrong?\n" ); 00149 fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" ); 00150 cdEXIT(EXIT_FAILURE); 00151 } 00152 00153 /* first print all warnings */ 00154 cdWarnings(ioQQQ); 00155 00156 /* cycle through the observational values */ 00157 nfound = 0; 00158 00159 /* first is to optimize relative emission line spectrum */ 00160 if( optimize.lgOptLin ) 00161 { 00162 /* set pointers to all optimized lines if first call */ 00163 if( !lgLinSet ) 00164 { 00165 lgHIT = true; 00166 lgLinSet = true; 00167 for( i=0; i < optimize.nlobs; i++ ) 00168 { 00169 double temp1, temp2; 00170 cap4(chFind , (char*)optimize.chLineLabel[i]); 00171 /* j = 0; */ 00172 00173 /* >> chng 06 may 04, use cdLine instead of ad hoc treatment. 00174 * no need to complain, cdLine will do it automatically. */ 00175 /* this is an intensity, get the line, returns false if could not find it */ 00176 j=cdLine( chFind , optimize.wavelength[i] , &temp1, &temp2 ); 00177 if( j<=0 ) 00178 { 00179 fprintf( ioQQQ, "\n" ); 00180 lgHIT = false; 00181 } 00182 else 00183 { 00184 ipobs[i] = j; 00185 } 00186 } 00187 00188 /* we did not find the line */ 00189 if( !lgHIT ) 00190 { 00191 fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" ); 00192 fprintf( ioQQQ, " Sorry.\n"); 00193 cdEXIT(EXIT_FAILURE); 00194 } 00195 } 00196 00197 for( i=0; i < 10; i++ ) 00198 { 00199 optimize.SavGenericData[i] = 0.; 00200 } 00201 00202 for( i=0; i < optimize.nlobs; i++ ) 00203 { 00204 /* and find corresponding model value by straight search */ 00205 nfound += 1; 00206 scld = (double)LineSv[ipobs[i]].sumlin[LineSave.lgLineEmergent]/(double)snorm*LineSave.ScaleNormLine; 00207 chi1 = chi2_func(scld,(double)optimize.xLineInt_Obs[i],(double)optimize.xLineInt_error[i]); 00208 cat = 0; 00209 nobs_cat[cat]++; 00210 chi2_cat[cat] += chi1; 00211 00212 fprintf( ioQQQ, " %4.4s ", 00213 LineSv[ipobs[i]].chALab); 00214 00215 prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength); 00216 00217 fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity", 00218 scld, 00219 optimize.xLineInt_Obs[i], 00220 optimize.xLineInt_error[i], 00221 chi1 ); 00222 00223 fprintf( ioQQQ, "\n" ); 00224 00225 if( i<10 ) 00226 { 00227 optimize.SavGenericData[i] = chi1; 00228 } 00229 } 00230 } 00231 00232 /* >>chng 06 may 04, moved this from above so that it only 00233 * gets printed if all lines are found */ 00234 /*if( optimize.lgOptLum || optimize.lgOptCol || optimize.lgOptTemp || optimize.lgOptLin )*/ 00235 if( optimize.lgOptimize ) 00236 { 00237 fprintf( ioQQQ, " ID Model Observed error chi**2 Type\n" ); 00238 } 00239 else 00240 { 00241 ASSERT( grid.lgGrid ); 00242 } 00243 00244 /* this is to optimize a mean temperature */ 00245 if( optimize.lgOptTemp ) 00246 { 00247 for( i=0; i < optimize.nTempObs; i++ ) 00248 { 00249 if( cdTemp(/*(char*)*/optimize.chTempLab[i],optimize.ionTemp[i], &temp_theory, optimize.chTempWeight[i]) ) 00250 { 00251 /* did not find column density */ 00252 fprintf(ioQQQ," optimizer did not find column density %s %li \n", 00253 optimize.chTempLab[i],optimize.ionTemp[i] ); 00254 cdEXIT(EXIT_FAILURE); 00255 } 00256 nfound += 1; 00257 chi1 = chi2_func(temp_theory,(double)optimize.temp_obs[i],(double)optimize.temp_error[i]); 00258 cat = 3; 00259 nobs_cat[cat]++; 00260 chi2_cat[cat] += chi1; 00261 00262 fprintf( ioQQQ, " %4.4s %2ld ", 00263 optimize.chTempLab[i], 00264 optimize.ionTemp[i] ); 00265 PrintE82( ioQQQ, temp_theory ); 00266 fprintf( ioQQQ, " "); 00267 PrintE82( ioQQQ, optimize.temp_obs[i] ); 00268 fprintf( ioQQQ, " %.5f %.2e", 00269 optimize.temp_error[i], chi1 ); 00270 fprintf( ioQQQ, " Temperature\n"); 00271 } 00272 } 00273 00274 /* option to optimize column densities */ 00275 if( optimize.lgOptCol ) 00276 { 00277 for( i=0; i < optimize.ncobs; i++ ) 00278 { 00279 if( cdColm((char*)optimize.chColDen_label[i],optimize.ion_ColDen[i], &theocl) ) 00280 { 00281 /* did not find column density */ 00282 fprintf(ioQQQ," optimizer did not find column density %s %li \n", 00283 optimize.chColDen_label[i],optimize.ion_ColDen[i] ); 00284 cdEXIT(EXIT_FAILURE); 00285 } 00286 nfound += 1; 00287 chi1 = chi2_func(theocl,(double)optimize.ColDen_Obs[i],(double)optimize.chColDen_error[i]); 00288 cat = 1; 00289 nobs_cat[cat]++; 00290 chi2_cat[cat] += chi1; 00291 00292 fprintf( ioQQQ, " %4.4s%6ld%12.4e%12.4e%12.5f%12.2e Column density\n", 00293 optimize.chColDen_label[i], optimize.ion_ColDen[i], theocl, 00294 optimize.ColDen_Obs[i], optimize.chColDen_error[i], chi1 ); 00295 } 00296 } 00297 00298 /* option to optimize line flux */ 00299 if( optimize.lgOptLum ) 00300 { 00301 nfound += 1; 00302 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0.f ) 00303 { 00304 predin = log10(LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten; 00305 help = pow(10.,predin-(double)optimize.optint); 00306 chi1 = chi2_func(help,1.,(double)optimize.optier); 00307 } 00308 else 00309 { 00310 predin = -999.99999; 00311 chi1 = (double)FLT_MAX; 00312 } 00313 cat = 2; 00314 nobs_cat[cat]++; 00315 chi2_cat[cat] += chi1; 00316 00317 fprintf( ioQQQ, " %4.4s%6f%12.5f%12.5f%12.5f%12.2e Line intensity\n", 00318 LineSv[LineSave.ipNormWavL].chALab, 00319 LineSv[LineSave.ipNormWavL].wavelength, 00320 predin, 00321 optimize.optint, 00322 optimize.optier, 00323 chi1 ); 00324 } 00325 00326 /* >>chng 06 apr 26, do not have to have line matches if doing grid. */ 00327 if( nfound <= 0 && !grid.lgGrid ) 00328 { 00329 fprintf( ioQQQ, " WARNING; no line matches found\n" ); 00330 cdEXIT(EXIT_FAILURE); 00331 } 00332 00333 /* write out chisquared for this iteration */ 00334 fprintf( ioQQQ, "\n" ); 00335 for( i=0; i < MAXCAT; i++ ) 00336 { 00337 if( nobs_cat[i] > 0 ) 00338 { 00339 chisq += chi2_cat[i]/nobs_cat[i]; 00340 fprintf( ioQQQ, " Category %s #obs.%3ld Total Chi**2%11.3e Average Chi**2%11.3e\n", 00341 name_cat[i],nobs_cat[i],chi2_cat[i],chi2_cat[i]/nobs_cat[i] ); 00342 } 00343 } 00344 if( nfound ) 00345 { 00346 fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq ); 00347 } 00348 } 00349 00350 /* increment nOptimiz, the grid / optimizer counter */ 00351 ++optimize.nOptimiz; 00352 00353 /* only print this if output has been turned on */ 00354 if( called.lgTalk ) 00355 { 00356 fprintf( ioQQQ, "\n" ); 00357 for( i=0; i < optimize.nvary; i++ ) 00358 { 00359 optimize.vparm[0][i] = (realnum)MIN2(optimize.vparm[0][i],optimize.varang[i][1]); 00360 optimize.vparm[0][i] = (realnum)MAX2(optimize.vparm[0][i],optimize.varang[i][0]); 00361 param[i] = optimize.vparm[0][i]; 00362 np = optimize.nvfpnt[i]; 00363 00364 /* now generate the actual command with parameter, 00365 * there will be from 1 to 3 numbers on the line */ 00366 if( optimize.nvarxt[i] == 1 ) 00367 { 00368 /* case with 1 parameter */ 00369 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] ); 00370 } 00371 00372 else if( optimize.nvarxt[i] == 2 ) 00373 { 00374 /* case with 2 parameter */ 00375 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]); 00376 } 00377 00378 else if( optimize.nvarxt[i] == 3 ) 00379 { 00380 /* case with 3 parameter */ 00381 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 00382 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] ); 00383 } 00384 00385 else if( optimize.nvarxt[i] == 4 ) 00386 { 00387 /* case with 4 parameter */ 00388 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 00389 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] ); 00390 } 00391 00392 else if( optimize.nvarxt[i] == 5 ) 00393 { 00394 /* case with 5 parameter */ 00395 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 00396 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], 00397 optimize.vparm[3][i] , optimize.vparm[4][i]); 00398 } 00399 00400 else 00401 { 00402 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me4\n"); 00403 cdEXIT(EXIT_FAILURE); 00404 } 00405 00406 fprintf( ioQQQ, " Varied command: %s\n", 00407 input.chCardSav[np] ); 00408 } 00409 } 00410 00411 func_v = MIN2(chisq,(double)FLT_MAX); 00412 return( func_v ); 00413 } 00414 00415 /* ============================================================================== */ 00416 STATIC double chi2_func(double ymodl, 00417 double ymeas, 00418 double yerr) 00419 { 00420 double chi2_func_v, 00421 temp; 00422 00423 DEBUG_ENTRY( "chi2_func()" ); 00424 00425 /* compute chi**2 by comparing model quantity ymodl with a measured 00426 * quantity ymeas with relative error yerr (negative means upper limit) 00427 */ 00428 00429 if( ymeas <= 0. ) 00430 { 00431 fprintf( ioQQQ, "chi2_func: non-positive observed quantity, this should not happen\n" ); 00432 cdEXIT(EXIT_FAILURE); 00433 } 00434 00435 if( yerr > 0. ) 00436 { 00437 if( ymodl > 0. ) 00438 { 00439 temp = POW2((ymodl-ymeas)/(MIN2(ymodl,ymeas)*yerr)); 00440 chi2_func_v = MIN2( temp , (double)FLT_MAX ); 00441 } 00442 else 00443 chi2_func_v = (double)FLT_MAX; 00444 } 00445 else if( yerr < 0. ) 00446 { 00447 /* value quoted is an upper limit, so add to chisq 00448 * only if limit exceeded, otherwise return zero. 00449 */ 00450 if( ymodl > ymeas ) 00451 { 00452 temp = POW2((ymodl-ymeas)/(ymeas*yerr)); 00453 chi2_func_v = MIN2(temp,(double)FLT_MAX); 00454 } 00455 else 00456 chi2_func_v = 0.; 00457 } 00458 else 00459 { 00460 fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" ); 00461 cdEXIT(EXIT_FAILURE); 00462 } 00463 return chi2_func_v; 00464 }