00001
00002
00003
00004 #include "cddefines.h"
00005 #include "zero.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 #include "iso.h"
00016 #include "taulines.h"
00017 #include "mole.h"
00018
00019 static double chi2_func(double,double,double);
00020
00021 double optimize_func(float param[])
00022 {
00023
00024 #define MAXCAT 4
00025
00026 char
00027 chFind[5];
00028
00029 bool lgBAD,
00030 lgHIT,
00031 lgLimOK;
00032
00033 long int cat,
00034 i,
00035 j,
00036 nfound,
00037 nobs_cat[MAXCAT],
00038 np;
00039
00040 static long int ipobs[NOBSLM];
00041
00042 double chi1,
00043 chi2_cat[MAXCAT],
00044 chisq,
00045 func_v,
00046 help,
00047 predin,
00048 scld,
00049 snorm,
00050 theocl,
00051 temp_theory;
00052
00053 static char name_cat[MAXCAT][13] =
00054 {
00055 "rel flux ",
00056 "column dens ",
00057 "abs flux ",
00058 "mean Temp "
00059 };
00060
00061 static bool lgLinSet = false;
00062
00063 DEBUG_ENTRY( "optimize_func()" );
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 zero();
00081
00082
00083
00084
00085
00086 if( optimize.lgOptimFlow )
00087 {
00088 fprintf( ioQQQ, " trace, optimize_func variables" );
00089 for( i=0; i < optimize.nvary; i++ )
00090 {
00091 fprintf( ioQQQ, "%10.2e", param[i] );
00092 }
00093 fprintf( ioQQQ, "\n" );
00094 }
00095
00096 for( i=0; i < optimize.nvary; i++ )
00097 {
00098 optimize.vparm[0][i] = param[i];
00099 }
00100
00101
00102
00103 ++optimize.nOptimiz;
00104
00105
00106
00107 vary_input(&lgLimOK);
00108
00109 if( !lgLimOK )
00110 {
00111
00112
00113 fprintf( ioQQQ, " Iteration %ld not within range.\n",
00114 optimize.nOptimiz );
00115
00116
00117 func_v = (double)FLT_MAX;
00118
00119 DEBUG_EXIT( "optimize_func()" );
00120 return( func_v );
00121 }
00122
00123 for( i=0; i < optimize.nvary; i++ )
00124 {
00125 optimize.varmax[i] = (float)MAX2(optimize.varmax[i],optimize.vpused[i]);
00126 optimize.varmin[i] = (float)MIN2(optimize.varmin[i],optimize.vpused[i]);
00127 }
00128
00129 lgBAD = cloudy();
00130 if( lgBAD )
00131 {
00132 fprintf( ioQQQ, " Cloudy returned error condition - what happened?\n" );
00133 }
00134
00135 if( grid.lgGrid )
00136 {
00137 GridInitialize();
00138 }
00139
00140
00141
00142 chisq = 0.0;
00143 for( i=0; i < MAXCAT; i++ )
00144 {
00145 nobs_cat[i] = 0;
00146 chi2_cat[i] = 0.0;
00147 }
00148
00149 if( LineSave.ipNormWavL < 0 )
00150 {
00151 fprintf( ioQQQ,
00152 " Normalization line array index is bad. What has gone wrong?\n" );
00153 puts( "[Stop in optimize_func]" );
00154 cdEXIT(EXIT_FAILURE);
00155 }
00156 snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent];
00157
00158 if( snorm == 0. )
00159 {
00160 fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity. What has gone wrong?\n" );
00161 fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
00162 puts( "[Stop in optimize_func]" );
00163 cdEXIT(EXIT_FAILURE);
00164 }
00165
00166
00167 cdWarnings(ioQQQ);
00168
00169
00170 nfound = 0;
00171
00172
00173 if( optimize.lgOptLin )
00174 {
00175
00176 if( !lgLinSet )
00177 {
00178 lgHIT = true;
00179 lgLinSet = true;
00180 for( i=0; i < optimize.nlobs; i++ )
00181 {
00182 double temp1, temp2;
00183 cap4(chFind , (char*)optimize.chLineLabel[i]);
00184
00185
00186
00187
00188
00189 j=cdLine( chFind , optimize.wavelength[i] , &temp1, &temp2 );
00190 if( j<=0 )
00191 {
00192 fprintf( ioQQQ, "\n" );
00193 lgHIT = false;
00194 }
00195 else
00196 {
00197 ipobs[i] = j;
00198 }
00199 }
00200
00201
00202 if( !lgHIT )
00203 {
00204 fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
00205 fprintf( ioQQQ, " Sorry.\n");
00206 puts( "[Stop in optimize_func]" );
00207 cdEXIT(EXIT_FAILURE);
00208 }
00209 }
00210
00211 for( i=0; i < 10; i++ )
00212 {
00213 optimize.SavGenericData[i] = 0.;
00214 }
00215
00216 for( i=0; i < optimize.nlobs; i++ )
00217 {
00218
00219 nfound += 1;
00220 scld = (double)LineSv[ipobs[i]].sumlin[LineSave.lgLineEmergent]/(double)snorm*LineSave.ScaleNormLine;
00221 chi1 = chi2_func(scld,(double)optimize.xLineInt_Obs[i],(double)optimize.xLineInt_error[i]);
00222 cat = 0;
00223 nobs_cat[cat]++;
00224 chi2_cat[cat] += chi1;
00225
00226 fprintf( ioQQQ, " %4.4s ",
00227 LineSv[ipobs[i]].chALab);
00228
00229 prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength);
00230
00231 fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity",
00232 scld,
00233 optimize.xLineInt_Obs[i],
00234 optimize.xLineInt_error[i],
00235 chi1 );
00236
00237 {
00238
00239 enum {DEBUG_LOC=false};
00240
00241
00242 long ipHi, ipLo;
00243
00244 for( ipLo = ipHe1s1S; ipLo <=20; ipLo ++ )
00245 {
00246 for( ipHi = ipLo+1; ipHi < iso.numLevels_local[ipHE_LIKE][ipHELIUM] -
00247 iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipHi++ )
00248 {
00249 if( fabs(EmisLines[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng - optimize.wavelength[i])
00250 < optimize.errorwave[i] && ipLo!=ipHe2p3P0 && ipLo!=ipHe2p3P1 )
00251 {
00252 fprintf( ioQQQ, "\t%li %li %li %li %li %li",
00253 iso.quant_desig[ipHE_LIKE][ipHELIUM][ipHi].n,
00254 iso.quant_desig[ipHE_LIKE][ipHELIUM][ipHi].l,
00255 iso.quant_desig[ipHE_LIKE][ipHELIUM][ipHi].s,
00256 iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].n,
00257 iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].l,
00258 iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].s );
00259 break;
00260 }
00261 }
00262 }
00263 }
00264
00265 fprintf( ioQQQ, "\n" );
00266
00267 if( i<10 )
00268 {
00269 optimize.SavGenericData[i] = chi1;
00270 }
00271 }
00272 }
00273
00274
00275
00276 if( optimize.lgOptLum || optimize.lgOptCol || optimize.lgOptTemp || optimize.lgOptLin )
00277 {
00278 fprintf( ioQQQ, " ID Model Observed error chi**2 Type\n" );
00279 }
00280 else
00281 {
00282 ASSERT( grid.lgGrid );
00283 }
00284
00285
00286 if( optimize.lgOptTemp )
00287 {
00288 for( i=0; i < optimize.nTempObs; i++ )
00289 {
00290 if( cdTemp(optimize.chTempLab[i],optimize.ionTemp[i], &temp_theory, optimize.chTempWeight[i]) )
00291 {
00292
00293 fprintf(ioQQQ," optimizer did not find column density %s %li \n",
00294 optimize.chTempLab[i],optimize.ionTemp[i] );
00295 puts( "[Stop in optimize_func]" );
00296 cdEXIT(EXIT_FAILURE);
00297 }
00298 nfound += 1;
00299 chi1 = chi2_func(temp_theory,(double)optimize.temp_obs[i],(double)optimize.temp_error[i]);
00300 cat = 3;
00301 nobs_cat[cat]++;
00302 chi2_cat[cat] += chi1;
00303
00304 fprintf( ioQQQ, " %4.4s %2ld ",
00305 optimize.chTempLab[i],
00306 optimize.ionTemp[i] );
00307 PrintE82( ioQQQ, temp_theory );
00308 fprintf( ioQQQ, " ");
00309 PrintE82( ioQQQ, optimize.temp_obs[i] );
00310 fprintf( ioQQQ, " %.5f %.2e",
00311 optimize.temp_error[i], chi1 );
00312 fprintf( ioQQQ, " Temperature\n");
00313 }
00314 }
00315
00316
00317 if( optimize.lgOptCol )
00318 {
00319 for( i=0; i < optimize.ncobs; i++ )
00320 {
00321 if( cdColm((char*)optimize.chColDen_label[i],optimize.ion_ColDen[i], &theocl) )
00322 {
00323
00324 fprintf(ioQQQ," optimizer did not find column density %s %li \n",
00325 optimize.chColDen_label[i],optimize.ion_ColDen[i] );
00326 puts( "[Stop in optimize_func]" );
00327 cdEXIT(EXIT_FAILURE);
00328 }
00329 nfound += 1;
00330 chi1 = chi2_func(theocl,(double)optimize.ColDen_Obs[i],(double)optimize.chColDen_error[i]);
00331 cat = 1;
00332 nobs_cat[cat]++;
00333 chi2_cat[cat] += chi1;
00334
00335 fprintf( ioQQQ, " %4.4s%6ld%12.4e%12.4e%12.5f%12.2e Column density\n",
00336 optimize.chColDen_label[i], optimize.ion_ColDen[i], theocl,
00337 optimize.ColDen_Obs[i], optimize.chColDen_error[i], chi1 );
00338 }
00339 }
00340
00341
00342 if( optimize.lgOptLum )
00343 {
00344 nfound += 1;
00345 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0.f )
00346 {
00347 predin = log10(LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten;
00348 help = pow(10.,predin-(double)optimize.optint);
00349 chi1 = chi2_func(help,1.,(double)optimize.optier);
00350 }
00351 else
00352 {
00353 predin = -999.99999;
00354 chi1 = (double)FLT_MAX;
00355 }
00356 cat = 2;
00357 nobs_cat[cat]++;
00358 chi2_cat[cat] += chi1;
00359
00360 fprintf( ioQQQ, " %4.4s%6f%12.5f%12.5f%12.5f%12.2e Line intensity\n",
00361 LineSv[LineSave.ipNormWavL].chALab,
00362 LineSv[LineSave.ipNormWavL].wavelength,
00363 predin,
00364 optimize.optint,
00365 optimize.optier,
00366 chi1 );
00367 }
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377 if( nfound <= 0 && !grid.lgGrid )
00378 {
00379 fprintf( ioQQQ, " WARNING; no line matches found\n" );
00380 puts( "[Stop in optimize_func]" );
00381 cdEXIT(EXIT_FAILURE);
00382 }
00383
00384
00385 fprintf( ioQQQ, "\n" );
00386 for( i=0; i < MAXCAT; i++ )
00387 {
00388 if( nobs_cat[i] > 0 )
00389 {
00390 chisq += chi2_cat[i]/nobs_cat[i];
00391 fprintf( ioQQQ, " Category %s #obs.%3ld Total Chi**2%11.3e Average Chi**2%11.3e\n",
00392 name_cat[i],nobs_cat[i],chi2_cat[i],chi2_cat[i]/nobs_cat[i] );
00393 }
00394 }
00395
00396 if( nfound )
00397 {
00398 fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq );
00399 }
00400
00401
00402 if( called.lgTalk )
00403 {
00404 fprintf( ioQQQ, "\n" );
00405 for( i=0; i < optimize.nvary; i++ )
00406 {
00407 optimize.vparm[0][i] = (float)MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
00408 optimize.vparm[0][i] = (float)MAX2(optimize.vparm[0][i],optimize.varang[i][0]);
00409 param[i] = optimize.vparm[0][i];
00410 np = optimize.nvfpnt[i];
00411
00412
00413
00414 if( optimize.nvarxt[i] == 1 )
00415 {
00416
00417 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] );
00418 }
00419
00420 else if( optimize.nvarxt[i] == 2 )
00421 {
00422
00423 sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]);
00424 }
00425
00426 else if( optimize.nvarxt[i] == 3 )
00427 {
00428
00429 sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
00430 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] );
00431 }
00432
00433 else if( optimize.nvarxt[i] == 4 )
00434 {
00435
00436 sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
00437 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] );
00438 }
00439
00440 else if( optimize.nvarxt[i] == 5 )
00441 {
00442
00443 sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
00444 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i],
00445 optimize.vparm[3][i] , optimize.vparm[4][i]);
00446 }
00447
00448 else
00449 {
00450 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me4\n");
00451 puts( "[Stop in optimize_func]" );
00452 cdEXIT(EXIT_FAILURE);
00453 }
00454
00455 fprintf( ioQQQ, " Optimal command: %s\n",
00456 input.chCardSav[np] );
00457 }
00458 }
00459
00460 func_v = MIN2(chisq,(double)FLT_MAX);
00461
00462 DEBUG_EXIT( "optimize_func()" );
00463 return( func_v );
00464 }
00465
00466
00467 static double chi2_func(double ymodl,
00468 double ymeas,
00469 double yerr)
00470 {
00471 double chi2_func_v,
00472 temp;
00473
00474 DEBUG_ENTRY( "chi2_func()" );
00475
00476
00477
00478
00479
00480 if( ymeas <= 0. )
00481 {
00482 fprintf( ioQQQ, "chi2_func: non-positive observed quantity, this should not happen\n" );
00483 puts( "[Stop]" );
00484 cdEXIT(EXIT_FAILURE);
00485 }
00486
00487 if( yerr > 0. )
00488 {
00489 if( ymodl > 0. )
00490 {
00491 temp = POW2((ymodl-ymeas)/(MIN2(ymodl,ymeas)*yerr));
00492 chi2_func_v = MIN2( temp , (double)FLT_MAX );
00493 }
00494 else
00495 chi2_func_v = (double)FLT_MAX;
00496 }
00497 else if( yerr < 0. )
00498 {
00499
00500
00501
00502 if( ymodl > ymeas )
00503 {
00504 temp = POW2((ymodl-ymeas)/(ymeas*yerr));
00505 chi2_func_v = MIN2(temp,(double)FLT_MAX);
00506 }
00507 else
00508 chi2_func_v = 0.;
00509 }
00510 else
00511 {
00512 fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
00513 puts( "[Stop]" );
00514 cdEXIT(EXIT_FAILURE);
00515 }
00516
00517 DEBUG_EXIT( "chi2_func()" );
00518 return chi2_func_v;
00519 }