cloudy
trunk
Main Page
Related Pages
Data Structures
Files
File List
Globals
All
Data Structures
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
source
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
{
113
optimize
.
varmax
[i] = (
realnum
)
MAX2
(
optimize
.
varmax
[i],
optimize
.
vpused
[i]);
114
optimize
.
varmin
[i] = (
realnum
)
MIN2
(
optimize
.
varmin
[i],
optimize
.
vpused
[i]);
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,
219
optimize
.
xLineInt_Obs
[i],
220
optimize
.
xLineInt_error
[i],
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"
,
253
optimize
.
chTempLab
[i],
optimize
.
ionTemp
[i] );
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
,
" "
);
267
PrintE82
(
ioQQQ
,
optimize
.
temp_obs
[i] );
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"
,
283
optimize
.
chColDen_label
[i],
optimize
.
ion_ColDen
[i] );
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,
294
optimize
.
ColDen_Obs
[i],
optimize
.
chColDen_error
[i], chi1 );
295
}
296
}
297
298
/* option to optimize line flux */
299
if
(
optimize
.
lgOptLum
)
300
{
301
nfound += 1;
302
if
(
LineSv
[
LineSave
.
ipNormWavL
].
sumlin
[
LineSave
.
lgLineEmergent
] > 0.f )
303
{
304
predin = log10(
LineSv
[
LineSave
.
ipNormWavL
].
sumlin
[
LineSave
.
lgLineEmergent
]) +
radius
.
Conv2PrtInten
;
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"
,
318
LineSv
[
LineSave
.
ipNormWavL
].
chALab
,
319
LineSv
[
LineSave
.
ipNormWavL
].
wavelength
,
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
1.8.1.1