cloudy  trunk
 All Data Structures 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-2008 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 "radius.h"
10 #include "input.h"
11 #include "cloudy.h"
12 #include "cddrive.h"
13 #include "optimize.h"
14 #include "grid.h"
15 /* used below to derive chi2 */
16 STATIC double chi2_func(double,double,double);
17 
18 double optimize_func(realnum param[])
19 {
20 
21 #define MAXCAT 4
22 
23  char /* chCAPLB[5], */
24  chFind[5];
25 
26  bool lgBAD,
27  lgHIT,
28  lgLimOK;
29 
30  long int cat,
31  i,
32  j,
33  nfound,
34  nobs_cat[MAXCAT],
35  np;
36 
37  static long int ipobs[NOBSLM];
38 
39  double chi1,
40  chi2_cat[MAXCAT],
41  chisq,
42  func_v,
43  help,
44  predin,
45  scld,
46  snorm,
47  theocl,
48  temp_theory;
49 
50  static char name_cat[MAXCAT][13] =
51  {
52  "rel flux ",
53  "column dens ",
54  "abs flux ",
55  "mean Temp "
56  };
57 
58  static bool lgLinSet = false;
59 
60  DEBUG_ENTRY( "optimize_func()" );
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  /* zero out lots of variables */
74  zero();
75 
76  if( optimize.lgOptimFlow )
77  {
78  fprintf( ioQQQ, " trace, optimize_func variables" );
79  for( i=0; i < optimize.nvary; i++ )
80  {
81  fprintf( ioQQQ, "%.2e", param[i] );
82  }
83  fprintf( ioQQQ, "\n" );
84  }
85 
86  for( i=0; i < optimize.nvary; i++ )
87  {
88  optimize.vparm[0][i] = param[i];
89  }
90 
91  /* call routine to pack variables with appropriate
92  * CLOUDY input lines given the array of variable parameters p(i) */
93  vary_input(&lgLimOK);
94 
95  if( !lgLimOK )
96  {
97  /* these parameters are not within limits of parameter search
98  * >>chng 96 apr 26, as per Peter van Hoof comment */
99  fprintf( ioQQQ, " Iteration %ld not within range.\n",
100  optimize.nOptimiz );
101 
102  /* this is error; very bad since not within range of parameters */
103  func_v = (double)FLT_MAX;
104 
105  /* always increment nOptimiz, even if parameters are out of bounds,
106  * this prevents optimizer to get stuck in infinite loop */
107  ++optimize.nOptimiz;
108  return( func_v );
109  }
110 
111  for( i=0; i < optimize.nvary; i++ )
112  {
115  }
116 
117  lgBAD = cloudy();
118  if( lgBAD )
119  {
120  fprintf( ioQQQ, " PROBLEM Cloudy returned error condition - what happened?\n" );
121  }
122 
123  if( grid.lgGrid )
124  {
125  /* this is the function's return value */
126  chisq = 0.;
127  }
128  else
129  {
130  /* this branch optimizing, not grid
131  / * extract line fluxes and compare with observations */
132  chisq = 0.0;
133  for( i=0; i < MAXCAT; i++ )
134  {
135  nobs_cat[i] = 0;
136  chi2_cat[i] = 0.0;
137  }
138 
139  if( LineSave.ipNormWavL < 0 )
140  {
141  fprintf( ioQQQ,
142  " Normalization line array index is bad. What has gone wrong?\n" );
143  cdEXIT(EXIT_FAILURE);
144  }
145 
146  if( (snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) == 0. )
147  {
148  fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity. What has gone wrong?\n" );
149  fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
150  cdEXIT(EXIT_FAILURE);
151  }
152 
153  /* first print all warnings */
154  cdWarnings(ioQQQ);
155 
156  /* cycle through the observational values */
157  nfound = 0;
158 
159  /* first is to optimize relative emission line spectrum */
160  if( optimize.lgOptLin )
161  {
162  /* set pointers to all optimized lines if first call */
163  if( !lgLinSet )
164  {
165  lgHIT = true;
166  lgLinSet = true;
167  for( i=0; i < optimize.nlobs; i++ )
168  {
169  double temp1, temp2;
170  cap4(chFind , (char*)optimize.chLineLabel[i]);
171  /* j = 0; */
172 
173  /* >> chng 06 may 04, use cdLine instead of ad hoc treatment.
174  * no need to complain, cdLine will do it automatically. */
175  /* this is an intensity, get the line, returns false if could not find it */
176  j=cdLine( chFind , optimize.wavelength[i] , &temp1, &temp2 );
177  if( j<=0 )
178  {
179  fprintf( ioQQQ, "\n" );
180  lgHIT = false;
181  }
182  else
183  {
184  ipobs[i] = j;
185  }
186  }
187 
188  /* we did not find the line */
189  if( !lgHIT )
190  {
191  fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
192  fprintf( ioQQQ, " Sorry.\n");
193  cdEXIT(EXIT_FAILURE);
194  }
195  }
196 
197  for( i=0; i < 10; i++ )
198  {
199  optimize.SavGenericData[i] = 0.;
200  }
201 
202  for( i=0; i < optimize.nlobs; i++ )
203  {
204  /* and find corresponding model value by straight search */
205  nfound += 1;
206  scld = (double)LineSv[ipobs[i]].sumlin[LineSave.lgLineEmergent]/(double)snorm*LineSave.ScaleNormLine;
207  chi1 = chi2_func(scld,(double)optimize.xLineInt_Obs[i],(double)optimize.xLineInt_error[i]);
208  cat = 0;
209  nobs_cat[cat]++;
210  chi2_cat[cat] += chi1;
211 
212  fprintf( ioQQQ, " %4.4s ",
213  LineSv[ipobs[i]].chALab);
214 
215  prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength);
216 
217  fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity",
218  scld,
221  chi1 );
222 
223  fprintf( ioQQQ, "\n" );
224 
225  if( i<10 )
226  {
227  optimize.SavGenericData[i] = chi1;
228  }
229  }
230  }
231 
232  /* >>chng 06 may 04, moved this from above so that it only
233  * gets printed if all lines are found */
234  /*if( optimize.lgOptLum || optimize.lgOptCol || optimize.lgOptTemp || optimize.lgOptLin )*/
235  if( optimize.lgOptimize )
236  {
237  fprintf( ioQQQ, " ID Model Observed error chi**2 Type\n" );
238  }
239  else
240  {
241  ASSERT( grid.lgGrid );
242  }
243 
244  /* this is to optimize a mean temperature */
245  if( optimize.lgOptTemp )
246  {
247  for( i=0; i < optimize.nTempObs; i++ )
248  {
249  if( cdTemp(/*(char*)*/optimize.chTempLab[i],optimize.ionTemp[i], &temp_theory, optimize.chTempWeight[i]) )
250  {
251  /* did not find column density */
252  fprintf(ioQQQ," optimizer did not find column density %s %li \n",
254  cdEXIT(EXIT_FAILURE);
255  }
256  nfound += 1;
257  chi1 = chi2_func(temp_theory,(double)optimize.temp_obs[i],(double)optimize.temp_error[i]);
258  cat = 3;
259  nobs_cat[cat]++;
260  chi2_cat[cat] += chi1;
261 
262  fprintf( ioQQQ, " %4.4s %2ld ",
263  optimize.chTempLab[i],
264  optimize.ionTemp[i] );
265  PrintE82( ioQQQ, temp_theory );
266  fprintf( ioQQQ, " ");
268  fprintf( ioQQQ, " %.5f %.2e",
269  optimize.temp_error[i], chi1 );
270  fprintf( ioQQQ, " Temperature\n");
271  }
272  }
273 
274  /* option to optimize column densities */
275  if( optimize.lgOptCol )
276  {
277  for( i=0; i < optimize.ncobs; i++ )
278  {
279  if( cdColm((char*)optimize.chColDen_label[i],optimize.ion_ColDen[i], &theocl) )
280  {
281  /* did not find column density */
282  fprintf(ioQQQ," optimizer did not find column density %s %li \n",
284  cdEXIT(EXIT_FAILURE);
285  }
286  nfound += 1;
287  chi1 = chi2_func(theocl,(double)optimize.ColDen_Obs[i],(double)optimize.chColDen_error[i]);
288  cat = 1;
289  nobs_cat[cat]++;
290  chi2_cat[cat] += chi1;
291 
292  fprintf( ioQQQ, " %4.4s%6ld%12.4e%12.4e%12.5f%12.2e Column density\n",
293  optimize.chColDen_label[i], optimize.ion_ColDen[i], theocl,
295  }
296  }
297 
298  /* option to optimize line flux */
299  if( optimize.lgOptLum )
300  {
301  nfound += 1;
303  {
305  help = pow(10.,predin-(double)optimize.optint);
306  chi1 = chi2_func(help,1.,(double)optimize.optier);
307  }
308  else
309  {
310  predin = -999.99999;
311  chi1 = (double)FLT_MAX;
312  }
313  cat = 2;
314  nobs_cat[cat]++;
315  chi2_cat[cat] += chi1;
316 
317  fprintf( ioQQQ, " %4.4s%6f%12.5f%12.5f%12.5f%12.2e Line intensity\n",
320  predin,
321  optimize.optint,
322  optimize.optier,
323  chi1 );
324  }
325 
326  /* >>chng 06 apr 26, do not have to have line matches if doing grid. */
327  if( nfound <= 0 && !grid.lgGrid )
328  {
329  fprintf( ioQQQ, " WARNING; no line matches found\n" );
330  cdEXIT(EXIT_FAILURE);
331  }
332 
333  /* write out chisquared for this iteration */
334  fprintf( ioQQQ, "\n" );
335  for( i=0; i < MAXCAT; i++ )
336  {
337  if( nobs_cat[i] > 0 )
338  {
339  chisq += chi2_cat[i]/nobs_cat[i];
340  fprintf( ioQQQ, " Category %s #obs.%3ld Total Chi**2%11.3e Average Chi**2%11.3e\n",
341  name_cat[i],nobs_cat[i],chi2_cat[i],chi2_cat[i]/nobs_cat[i] );
342  }
343  }
344  if( nfound )
345  {
346  fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq );
347  }
348  }
349 
350  /* increment nOptimiz, the grid / optimizer counter */
351  ++optimize.nOptimiz;
352 
353  /* only print this if output has been turned on */
354  if( called.lgTalk )
355  {
356  fprintf( ioQQQ, "\n" );
357  for( i=0; i < optimize.nvary; i++ )
358  {
359  optimize.vparm[0][i] = (realnum)MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
360  optimize.vparm[0][i] = (realnum)MAX2(optimize.vparm[0][i],optimize.varang[i][0]);
361  param[i] = optimize.vparm[0][i];
362  np = optimize.nvfpnt[i];
363 
364  /* now generate the actual command with parameter,
365  * there will be from 1 to 3 numbers on the line */
366  if( optimize.nvarxt[i] == 1 )
367  {
368  /* case with 1 parameter */
369  sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] );
370  }
371 
372  else if( optimize.nvarxt[i] == 2 )
373  {
374  /* case with 2 parameter */
375  sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]);
376  }
377 
378  else if( optimize.nvarxt[i] == 3 )
379  {
380  /* case with 3 parameter */
381  sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
382  optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] );
383  }
384 
385  else if( optimize.nvarxt[i] == 4 )
386  {
387  /* case with 4 parameter */
388  sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
389  optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] );
390  }
391 
392  else if( optimize.nvarxt[i] == 5 )
393  {
394  /* case with 5 parameter */
395  sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
396  optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i],
397  optimize.vparm[3][i] , optimize.vparm[4][i]);
398  }
399 
400  else
401  {
402  fprintf(ioQQQ,"The number of variable options on this line makes no sense to me4\n");
403  cdEXIT(EXIT_FAILURE);
404  }
405 
406  fprintf( ioQQQ, " Varied command: %s\n",
407  input.chCardSav[np] );
408  }
409  }
410 
411  func_v = MIN2(chisq,(double)FLT_MAX);
412  return( func_v );
413 }
414 
415 /* ============================================================================== */
416 STATIC double chi2_func(double ymodl,
417  double ymeas,
418  double yerr)
419 {
420  double chi2_func_v,
421  temp;
422 
423  DEBUG_ENTRY( "chi2_func()" );
424 
425  /* compute chi**2 by comparing model quantity ymodl with a measured
426  * quantity ymeas with relative error yerr (negative means upper limit)
427  */
428 
429  if( ymeas <= 0. )
430  {
431  fprintf( ioQQQ, "chi2_func: non-positive observed quantity, this should not happen\n" );
432  cdEXIT(EXIT_FAILURE);
433  }
434 
435  if( yerr > 0. )
436  {
437  if( ymodl > 0. )
438  {
439  temp = POW2((ymodl-ymeas)/(MIN2(ymodl,ymeas)*yerr));
440  chi2_func_v = MIN2( temp , (double)FLT_MAX );
441  }
442  else
443  chi2_func_v = (double)FLT_MAX;
444  }
445  else if( yerr < 0. )
446  {
447  /* value quoted is an upper limit, so add to chisq
448  * only if limit exceeded, otherwise return zero.
449  */
450  if( ymodl > ymeas )
451  {
452  temp = POW2((ymodl-ymeas)/(ymeas*yerr));
453  chi2_func_v = MIN2(temp,(double)FLT_MAX);
454  }
455  else
456  chi2_func_v = 0.;
457  }
458  else
459  {
460  fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
461  cdEXIT(EXIT_FAILURE);
462  }
463  return chi2_func_v;
464 }

Generated for cloudy by doxygen 1.8.1.1