cloudy trunk
Loading...
Searching...
No Matches
optimize_func.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*optimize_func actual function called during evaluation of optimization run */
4#include "cddefines.h"
5#include "init.h"
6#include "lines.h"
7#include "prt.h"
8#include "called.h"
9#include "predcont.h"
10#include "radius.h"
11#include "rfield.h"
12#include "input.h"
13#include "cloudy.h"
14#include "cddrive.h"
15#include "optimize.h"
16#include "grid.h"
17/* used below to derive chi2 */
18STATIC double chi2_func(double,double,double);
19
21 int grid_index)
22{
23 const int MAXCAT = 6;
24
25 static const char name_cat[MAXCAT][13] =
26 {
27 "rel flux ",
28 "column dens ",
29 "abs flux ",
30 "mean temp ",
31 "ang diameter",
32 "photometry "
33 };
34
35 char chFind[5];
36
37 bool lgBAD,
38 lgLimOK;
39
40 long int cat,
41 i,
42 j,
43 nfound,
44 nobs_cat[MAXCAT],
45 np;
46
47 chi2_type chi1,
48 chi2_cat[MAXCAT],
49 chisq,
50 help,
51 predin,
52 scld,
53 snorm,
54 theocl,
55 temp_theory;
56
57 DEBUG_ENTRY( "optimize_func()" );
58
59 if( grid_index >= 0 )
60 optimize.nOptimiz = grid_index;
61
62 /* This routine is called by optimizer with values of the
63 * variable parameters for CLOUDY in the array p(i). It returns
64 * the value FUNC = SUM (obs-model)**2/sig**2 for the lines
65 * specified in the observational data file, values held in the
66 * common blocks /OBSLIN/ & /OBSINT/
67 * replacement input strings for CLOUDY READR held in /chCardSav/
68 * parameter information for setting chCardSav held in /parmv/
69 * additional variables
70 * Gary's variables
71 */
72
73 if( optimize.lgOptimFlow )
74 {
75 fprintf( ioQQQ, " trace, optimize_func variables" );
76 for( i=0; i < optimize.nvary; i++ )
77 {
78 fprintf( ioQQQ, "%.2e", param[i] );
79 }
80 fprintf( ioQQQ, "\n" );
81 }
82
83 for( i=0; i < optimize.nvary; i++ )
84 {
85 optimize.vparm[0][i] = param[i];
86 }
87
88 /* call routine to pack variables with appropriate
89 * CLOUDY input lines given the array of variable parameters p(i) */
90 vary_input( &lgLimOK, grid_index );
91
92 // nothing more to be done...
93 if( strcmp(optimize.chOptRtn,"XSPE") == 0 )
94 return 0.;
95
96 /* zero out lots of variables */
97 zero();
98
99 for( i=0; i < optimize.nvary; i++ )
100 {
101 optimize.varmax[i] = max(optimize.varmax[i],min(param[i],optimize.varang[i][1]));
102 optimize.varmin[i] = min(optimize.varmin[i],max(param[i],optimize.varang[i][0]));
103 }
104
105 if( !lgLimOK )
106 {
107 /* these parameters are not within limits of parameter search
108 * >>chng 96 apr 26, as per Peter van Hoof comment */
109 fprintf( ioQQQ, " Iteration %ld not within range.\n",
110 optimize.nOptimiz );
111
112 /* always increment nOptimiz, even if parameters are out of bounds,
113 * this prevents optimizer to get stuck in infinite loop */
114 ++optimize.nOptimiz;
115
116 /* this is error; very bad since not within range of parameters */
117 return BIG_CHI2;
118 }
119
120 lgBAD = cloudy();
121 if( lgBAD )
122 {
123 fprintf( ioQQQ, " PROBLEM Cloudy returned error condition - what happened?\n" );
124 }
125
126 if( grid.lgGrid )
127 {
128 /* this is the function's return value */
129 chisq = 0.;
130 }
131 else
132 {
133 /* this branch optimizing, not grid
134 / * extract line fluxes and compare with observations */
135 chisq = 0.0;
136 for( i=0; i < MAXCAT; i++ )
137 {
138 nobs_cat[i] = 0;
139 chi2_cat[i] = 0.0;
140 }
141
142 if( LineSave.ipNormWavL < 0 )
143 {
144 fprintf( ioQQQ,
145 " Normalization line array index is bad. What has gone wrong?\n" );
147 }
148
149 if( (snorm = LineSv[LineSave.ipNormWavL].SumLine[optimize.nEmergent]) == 0. )
150 {
151 fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity. What has gone wrong?\n" );
152 fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
154 }
155
156 /* first print all warnings */
158
159 /* print header before any of the individual chi2 values are printed */
160 if( optimize.lgOptimize )
161 fprintf( ioQQQ, " ID Model Observed error chi**2 Type\n" );
162 else
163 ASSERT( grid.lgGrid );
164
165 /* cycle through the observational values */
166 nfound = 0;
167
168 /* first is to optimize relative emission line spectrum */
169 if( optimize.xLineInt_Obs.size() > 0 )
170 {
171 /* set pointers to all optimized lines if first call */
172 if( optimize.ipobs.size() == 0 )
173 {
174 optimize.ipobs.resize( optimize.xLineInt_Obs.size() );
175 bool lgHIT = true;
176 for( i=0; i < long(optimize.xLineInt_Obs.size()); i++ )
177 {
178 chi2_type temp1, temp2;
179 cap4( chFind, optimize.chLineLabel[i].c_str() );
180
181 /* >> chng 06 may 04, use cdLine instead of ad hoc treatment.
182 * no need to complain, cdLine will do it automatically. */
183 /* this is an intensity, get the line, returns false if could not find it */
184 j = cdLine( chFind, optimize.wavelength[i], &temp1, &temp2 );
185 if( j <= 0 )
186 {
187 fprintf( ioQQQ, "\n" );
188 lgHIT = false;
189 }
190 else
191 {
192 optimize.ipobs[i] = j;
193 }
194 }
195
196 /* we did not find the line */
197 if( !lgHIT )
198 {
199 fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
200 fprintf( ioQQQ, " Sorry.\n");
202 }
203 }
204
205 for( i=0; i < 10; i++ )
206 {
207 optimize.SavGenericData[i] = 0.;
208 }
209
210 for( i=0; i < long(optimize.xLineInt_Obs.size()); i++ )
211 {
212 /* and find corresponding model value by straight search */
213 nfound += 1;
214 scld = (chi2_type)LineSv[optimize.ipobs[i]].SumLine[optimize.nEmergent]/
215 (chi2_type)snorm*LineSave.ScaleNormLine;
216 chi1 = chi2_func(scld,(chi2_type)optimize.xLineInt_Obs[i],
217 (chi2_type)optimize.xLineInt_error[i]);
218 cat = 0;
219 nobs_cat[cat]++;
220 chi2_cat[cat] += chi1;
221
222 fprintf( ioQQQ, " %4.4s ", LineSv[optimize.ipobs[i]].chALab);
223
224 prt_wl( ioQQQ,LineSv[optimize.ipobs[i]].wavelength);
225
226 fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity",
227 scld,
228 optimize.xLineInt_Obs[i],
229 optimize.xLineInt_error[i],
230 chi1 );
231
232 fprintf( ioQQQ, "\n" );
233
234 if( i<10 )
235 {
236 optimize.SavGenericData[i] = chi1;
237 }
238 }
239 }
240
241 /* this is to optimize a mean temperature */
242 for( i=0; i < long(optimize.temp_obs.size()); i++ )
243 {
244 if( cdTemp( optimize.chTempLab[i].c_str(), optimize.ionTemp[i],
245 &temp_theory, optimize.chTempWeight[i].c_str()) )
246 {
247 /* did not find column density */
248 fprintf(ioQQQ," optimizer did not find column density %s %li \n",
249 optimize.chTempLab[i].c_str(),optimize.ionTemp[i] );
251 }
252 nfound += 1;
253 chi1 = chi2_func(temp_theory,(chi2_type)optimize.temp_obs[i],
254 (chi2_type)optimize.temp_error[i]);
255 cat = 3;
256 nobs_cat[cat]++;
257 chi2_cat[cat] += chi1;
258
259 fprintf( ioQQQ, " %4.4s%7ld%12.4e%12.4e%12.5f%12.2e Temperature\n",
260 optimize.chTempLab[i].c_str(), optimize.ionTemp[i], temp_theory,
261 optimize.temp_obs[i], optimize.temp_error[i], chi1 );
262 }
263
264 /* option to optimize column densities */
265 for( i=0; i < long(optimize.ColDen_Obs.size()); i++ )
266 {
267 if( cdColm(optimize.chColDen_label[i].c_str(),optimize.ion_ColDen[i], &theocl) )
268 {
269 /* did not find column density */
270 fprintf(ioQQQ," optimizer did not find column density %s %li \n",
271 optimize.chColDen_label[i].c_str(), optimize.ion_ColDen[i] );
273 }
274 nfound++;
275 chi1 = chi2_func(theocl,(chi2_type)optimize.ColDen_Obs[i],
276 (chi2_type)optimize.ColDen_error[i]);
277 cat = 1;
278 nobs_cat[cat]++;
279 chi2_cat[cat] += chi1;
280
281 fprintf( ioQQQ, " %4.4s%7ld%12.4e%12.4e%12.5f%12.2e Column density\n",
282 optimize.chColDen_label[i].c_str(), optimize.ion_ColDen[i], theocl,
283 optimize.ColDen_Obs[i], optimize.ColDen_error[i], chi1 );
284 }
285
286 /* option to optimize line flux */
287 if( optimize.lgOptLum )
288 {
289 ++nfound;
290 if( LineSv[LineSave.ipNormWavL].SumLine[optimize.nOptLum] > 0.f )
291 {
292 predin = log10(LineSv[LineSave.ipNormWavL].SumLine[optimize.nOptLum]) +
293 radius.Conv2PrtInten;
294 help = pow(10.,predin-(chi2_type)optimize.optint);
295 chi1 = chi2_func(help,1.,(chi2_type)optimize.optier);
296 }
297 else
298 {
299 predin = -999.99999;
300 chi1 = BIG_CHI2;
301 }
302 cat = 2;
303 nobs_cat[cat]++;
304 chi2_cat[cat] += chi1;
305
306 fprintf( ioQQQ, " %4.4s ",
307 LineSv[LineSave.ipNormWavL].chALab );
308
309 prt_wl( ioQQQ, LineSv[LineSave.ipNormWavL].wavelength );
310
311 fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Line intensity\n",
312 predin,
313 optimize.optint,
314 optimize.optier,
315 chi1 );
316 }
317
318 /* option to optimize the absolute continuum flux */
319 for( size_t k=0; k < optimize.ContIndex.size(); k++ )
320 {
321 nfound++;
322 // there are 4 entries for each wavelength: nFnu, nInu, InwT, InwC
323 long ind = t_PredCont::Inst().offset() + 4*optimize.ContIndex[k];
324 chi2_type nFnu_model = 0.;
325 if( LineSv[ind].SumLine[0] > SMALLFLOAT )
326 {
327 nFnu_model = chi2_type( log10(LineSv[ind].SumLine[0]) + radius.Conv2PrtInten );
328 nFnu_model = pow(10.,nFnu_model);
329 }
330 Flux F_model( optimize.ContNFnu[k].E(), nFnu_model );
331
332 chi1 = chi2_func(nFnu_model,optimize.ContNFnu[k].get("erg/s/cm2"),optimize.ContNFnuErr[k]);
333 const char* catstr;
334 // treat radio continuum flux as absolute flux so that it can be used
335 // as a more accurate replacement of the normalization line intensity
336 if( optimize.ContEner[k].mm() <= 1. )
337 {
338 cat = 5;
339 catstr = "Photometry";
340 }
341 else
342 {
343 cat = 2;
344 catstr = "Radio intensity";
345 }
346 nobs_cat[cat]++;
347 chi2_cat[cat] += chi1;
348
349 fprintf( ioQQQ, " %4.4s ", LineSv[ind].chALab);
351 string unit = optimize.ContNFnu[k].uu();
352 fprintf( ioQQQ, "%12.4g%12.4g%12.5f%12.2e %s [%s]\n",
353 F_model.get(unit),
354 optimize.ContNFnu[k].get(unit),
355 optimize.ContNFnuErr[k], chi1,
356 catstr, unit.c_str() );
357 }
358
359 /* option to optimize angular diamater */
360 if( optimize.lgOptDiam )
361 {
362 nfound++;
363 chi2_type diam_model;
364 // get diameter in cm
365 if( rfield.lgUSphON )
366 diam_model = 2.*rfield.rstrom; // ionization bounded -> use Stroemgren radius
367 else
368 diam_model = 2.*radius.Radius; // density bounded -> use outer radius
369 // now convert to arcsec if necessary
370 if( !optimize.lgDiamInCM && radius.distance > 0. )
371 diam_model *= AS1RAD/radius.distance;
372
373 chi1 = chi2_func(diam_model,optimize.optDiam,optimize.optDiamErr);
374 cat = 4;
375 nobs_cat[cat]++;
376 chi2_cat[cat] += chi1;
377
378 fprintf( ioQQQ, " %12.4g%12.4g%12.5f%12.2e Angular diameter\n",
379 diam_model, optimize.optDiam, optimize.optDiamErr, chi1 );
380 }
381
382 /* do not have to have line matches if doing grid. */
383 if( nfound <= 0 && !grid.lgGrid )
384 {
385 fprintf( ioQQQ, " WARNING; no line matches found\n" );
387 }
388
389 /* write out chisquared for this iteration */
390 fprintf( ioQQQ, "\n" );
391 for( i=0; i < MAXCAT; i++ )
392 {
393 if( nobs_cat[i] > 0 )
394 {
395 chisq += chi2_cat[i]/nobs_cat[i];
396 fprintf( ioQQQ, " Category %s #obs.%3ld Total Chi**2%11.3e Average Chi**2%11.3e\n",
397 name_cat[i],nobs_cat[i],chi2_cat[i],chi2_cat[i]/nobs_cat[i] );
398 }
399 }
400 if( nfound )
401 {
402 fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq );
403 }
404 }
405
406 /* increment nOptimiz, the grid / optimizer counter */
407 ++optimize.nOptimiz;
408
409 /* only print this if output has been turned on */
410 if( called.lgTalk )
411 {
412 fprintf( ioQQQ, "\n" );
413 for( i=0; i < optimize.nvary; i++ )
414 {
415 np = optimize.nvfpnt[i];
416
417 /* now generate the actual command with parameter,
418 * there will be from 1 to 3 numbers on the line */
419 if( optimize.nvarxt[i] == 1 )
420 {
421 /* case with 1 parameter */
422 sprintf( input.chCardSav[np], optimize.chVarFmt[i], optimize.vparm[0][i] );
423 }
424
425 else if( optimize.nvarxt[i] == 2 )
426 {
427 /* case with 2 parameter */
428 sprintf( input.chCardSav[np], optimize.chVarFmt[i], optimize.vparm[0][i],
429 optimize.vparm[1][i]);
430 }
431
432 else if( optimize.nvarxt[i] == 3 )
433 {
434 /* case with 3 parameter */
435 sprintf( input.chCardSav[np], optimize.chVarFmt[i],
436 optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i] );
437 }
438
439 else if( optimize.nvarxt[i] == 4 )
440 {
441 /* case with 4 parameter */
442 sprintf( input.chCardSav[np], optimize.chVarFmt[i],
443 optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i],
444 optimize.vparm[3][i] );
445 }
446
447 else if( optimize.nvarxt[i] == 5 )
448 {
449 /* case with 5 parameter */
450 sprintf( input.chCardSav[np], optimize.chVarFmt[i],
451 optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i],
452 optimize.vparm[3][i], optimize.vparm[4][i]);
453 }
454
455 else
456 {
457 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me4\n");
459 }
460
461 fprintf( ioQQQ, " Varied command: %s\n",
462 input.chCardSav[np] );
463 }
464 }
465
466 return min(chisq,BIG_CHI2);
467}
468
469/* ============================================================================== */
471 chi2_type ymeas,
472 chi2_type yerr)
473{
474 chi2_type chi2_func_v,
475 temp;
476
477 DEBUG_ENTRY( "chi2_func()" );
478
479 /* compute chi**2 by comparing model quantity ymodl with a measured
480 * quantity ymeas with relative error yerr (negative means upper limit)
481 */
482
483 if( ymeas <= 0. )
484 {
485 fprintf( ioQQQ, "chi2_func: non-positive observed quantity, this should not happen\n" );
487 }
488
489 if( yerr > 0. )
490 {
491 if( ymodl > 0. )
492 {
493 temp = pow2((ymodl-ymeas)/(min(ymodl,ymeas)*yerr));
494 chi2_func_v = min(temp,BIG_CHI2);
495 }
496 else
497 chi2_func_v = BIG_CHI2;
498 }
499 else if( yerr < 0. )
500 {
501 /* value quoted is an upper limit, so add to chisq
502 * only if limit exceeded, otherwise return zero.
503 */
504 if( ymodl > ymeas )
505 {
506 temp = pow2((ymodl-ymeas)/(ymeas*yerr));
507 chi2_func_v = min(temp,BIG_CHI2);
508 }
509 else
510 chi2_func_v = 0.;
511 }
512 else
513 {
514 fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
516 }
517 return chi2_func_v;
518}
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
T pow2(T a)
Definition cddefines.h:931
#define EXIT_FAILURE
Definition cddefines.h:140
void cap4(char *chCAP, const char *chLab)
Definition service.cpp:240
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
long max(int a, long b)
Definition cddefines.h:775
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition cddrive.cpp:636
void cdWarnings(FILE *ioPNT)
Definition cddrive.cpp:198
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition cddrive.cpp:1602
LinSv * LineSv
Definition cdinit.cpp:70
Definition flux.h:10
double get() const
Definition flux.h:51
static t_PredCont & Inst()
Definition cddefines.h:175
long offset() const
Definition predcont.h:36
bool cloudy()
Definition cloudy.cpp:38
const realnum SMALLFLOAT
Definition cpu.h:191
t_grid grid
Definition grid.cpp:5
void zero(void)
Definition zero.cpp:73
t_input input
Definition input.cpp:12
t_LineSave LineSave
Definition lines.cpp:5
static realnum * wavelength
t_optimize optimize
Definition optimize.cpp:5
double chi2_type
Definition optimize.h:49
const chi2_type BIG_CHI2
Definition optimize.h:51
void vary_input(bool *lgLimOK, int grid_index)
STATIC double chi2_func(double, double, double)
chi2_type optimize_func(const realnum param[], int grid_index)
UNUSED const double AS1RAD
Definition physconst.h:129
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8