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