cloudy trunk
Loading...
Searching...
No Matches
optimize_do.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/*lgOptimize_do main driver for optimization runs*/
4#include "cddefines.h"
5#define NPLXMX (LIMPAR*(LIMPAR+6)+1)
6#include "input.h"
7#include "called.h"
8#include "prt.h"
9#include "save.h"
10#include "optimize.h"
11
12/* called by cdDrive, this returns false if things went ok, true for disaster */
13bool lgOptimize_do(void)
14{
15 long int i,
16 iflag,
17 ii,
18 iworke[NPLXMX],
19 j,
20 mode,
21 need,
22 nfe,
23 np;
24 chi2_type fret;
25 realnum fx,
26 ptem[LIMPAR],
27 delta[LIMPAR],
28 toler,
29 worke[NPLXMX];
30
31 DEBUG_ENTRY( "lgOptimize_do()" );
32
33 /* main driver for optimization runs
34 * Drives cloudy to optimize variables;*/
35
36 /* code originally written by R.F. Carswell, IOA Cambridge */
37
38 toler = (realnum)log10(1. + optimize.OptGlobalErr);
39
40 if( strcmp(optimize.chOptRtn,"PHYM") == 0 )
41 {
42 /* Phymir method */
43 for( j=0; j < optimize.nvary; j++ )
44 {
45 ptem[j] = optimize.vparm[0][j];
46 delta[j] = optimize.vincr[j];
47 }
48 /* >>chng 06 jan 02, fix uninitialized var problem detected by valgrind/purify, PvH */
49 for( j=optimize.nvary; j < LIMPAR; j++ )
50 {
51 ptem[j] = -FLT_MAX;
52 delta[j] = -FLT_MAX;
53 }
54 optimize_phymir(ptem,delta,optimize.nvary,&fret,toler);
55 for( j=0; j < optimize.nvary; j++ )
56 optimize.vparm[0][j] = ptem[j];
57 }
58
59 else if( strcmp(optimize.chOptRtn,"SUBP") == 0 )
60 {
61 fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" );
62 need = 2*optimize.nvary + optimize.nvary*(optimize.nvary + 4) + 1;
63 if( need > NPLXMX )
64 {
65 fprintf( ioQQQ, " Increase size of NPLXMX in parameter statements to handle this many variables.\n" );
66 fprintf( ioQQQ, " I need at least %5ld\n", need );
68 }
69 for( j=0; j < optimize.nvary; j++ )
70 ptem[j] = optimize.vparm[0][j];
71
72 /* The subroutine SUBPLX input into cloudy 8/4/94.
73 * The program itself is very well commented.
74 * The mode must set to 0 for the default values.
75 * The switch iflag tells if the program terminated normally. */
76 mode = 0;
77
78 /* >>chng 97 dec 08, remove first arg, optimize_func, since not used in routines */
80 /* the number of parameters to vary */
81 optimize.nvary,
82 /* the relative error, single number */
83 toler,
84 /* maximum number of function evaluations before giving up */
85 optimize.nIterOptim,
86 /* mode of operation, we simply set to zero */
87 mode,
88 /* the initial changes in the guessed best coefficients, typically 0.2 to 1 */
89 optimize.vincr,
90 /* a vector of nvary initial parameters that are the starting guesses for the parameters */
91 ptem,
92 /* a realnum, this is simply ignored */
93 &fx,
94 /* another parameter that is simply ignored, a long int */
95 &nfe,
96 /* a realnum that is NPLXMX long, used for working space by the routine */
97 /* an array that is 20*26 + 1 elements long, used for working space */
98 worke,
99 /* a long int that is NPLXMX long, used for working space by the routine */
100 /* an array that is 20*26 + 1 elements long, used for working space */
101 iworke,
102 /* a long int - says what happened, if -1 then exceeded nIterOptim iteration */
103 &iflag);
104
105 if( iflag == -1 )
106 {
107 fprintf( ioQQQ, " SUBPLEX exceeding maximum iterations.\n"
108 " This can be reset with the OPTIMZE ITERATIONS command.\n" );
109 }
110
111 for( j=0; j < optimize.nvary; j++ )
112 optimize.vparm[0][j] = ptem[j];
113
114 if( optimize.lgOptimFlow )
115 {
116 fprintf( ioQQQ, " trace return optimize_subplex:\n" );
117 for( j=0; j < optimize.nvary; j++ )
118 {
119 fprintf( ioQQQ, " Values:" );
120 for( ii=1; ii <= optimize.nvarxt[j]; ii++ )
121 fprintf( ioQQQ, " %.2e", optimize.vparm[ii-1][j] );
122 fprintf( ioQQQ, "\n" );
123 }
124 }
125 }
126 else
128
129 // turn output back on for final model
130 called.lgTalk = cpu.i().lgMPI_talk();
131 called.lgTalkIsOK = cpu.i().lgMPI_talk();
132 prt.lgFaintOn = true;
133
134 if( called.lgTalk )
135 {
136 fprintf( ioQQQ, " **************************************************\n" );
137 fprintf( ioQQQ, " **************************************************\n" );
138 fprintf( ioQQQ, " **************************************************\n" );
139 fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", optimize.nOptimiz );
140
141 for( i=0; i < optimize.nvary; i++ )
142 {
143 np = optimize.nvfpnt[i];
144
145 /* now generate the actual command with parameter */
146 if( optimize.nvarxt[i] == 1 )
147 {
148 sprintf( input.chCardSav[np], optimize.chVarFmt[i], optimize.vparm[0][i] );
149 }
150
151 else if( optimize.nvarxt[i] == 2 )
152 {
153 sprintf( input.chCardSav[np], optimize.chVarFmt[i], optimize.vparm[0][i],
154 optimize.vparm[1][i]);
155 }
156
157 else if( optimize.nvarxt[i] == 3 )
158 {
159 sprintf( input.chCardSav[np], optimize.chVarFmt[i],
160 optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i] );
161 }
162
163 else if( optimize.nvarxt[i] == 4 )
164 {
165 sprintf( input.chCardSav[np], optimize.chVarFmt[i],
166 optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i],
167 optimize.vparm[3][i] );
168 }
169
170 else if( optimize.nvarxt[i] == 5 )
171 {
172 sprintf( input.chCardSav[np], optimize.chVarFmt[i],
173 optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i],
174 optimize.vparm[3][i], optimize.vparm[4][i]);
175 }
176
177 else
178 {
179 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me.\n");
181 }
182
183 fprintf( ioQQQ, " Optimal command: %s\n", input.chCardSav[np]);
184 fprintf( ioQQQ, " Smallest value:%10.2e Largest value:%10.2e Allowed range %10.2e to %10.2e\n",
185 optimize.varmin[i], optimize.varmax[i], optimize.varang[i][0],
186 optimize.varang[i][1] );
187 }
188 }
189
190 if( cpu.i().lgMaster() )
191 {
192 /* save optimal parameters */
193 FILE *ioOptim = open_data( chOptimFileName, "w", AS_LOCAL_ONLY );
194 for( i=0; i <= input.nSave; i++ )
195 fprintf( ioOptim, "%s\n", input.chCardSav[i]);
196 fclose( ioOptim );
197
198 /* recalculate values in cloudy for the best fit, and print out
199 * all the information */
200 fprintf( ioQQQ, "\f" );
201
202 for( i=0; i < optimize.nvary; i++ )
203 ptem[i] = optimize.vparm[0][i];
204
205 (void)optimize_func(ptem);
206 }
207
208 return lgAbort;
209}
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
static t_cpu cpu
Definition cpu.h:355
@ AS_LOCAL_ONLY
Definition cpu.h:208
t_input input
Definition input.cpp:12
char chOptimFileName[INPUT_LINE_LENGTH]
Definition optimize.cpp:6
t_optimize optimize
Definition optimize.cpp:5
void optimize_phymir(realnum[], const realnum[], long, chi2_type *, realnum)
double chi2_type
Definition optimize.h:49
void optimize_subplex(long int n, double tol, long int maxnfe, long int mode, realnum scale[], realnum x[], realnum *fx, long int *nfe, realnum work[], long int iwork[], long int *iflag)
chi2_type optimize_func(const realnum param[], int grid_index=-1)
const long LIMPAR
Definition optimize.h:61
#define NPLXMX
bool lgOptimize_do(void)
t_prt prt
Definition prt.cpp:10