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