Blame view

src/idl/dustem_init_params.pro 12.5 KB
a4c9699a   Annie Hughes   first commit
1
2
3
4
PRO dustem_init_params,model,pd,iv, $
                       fpd=fpd,fiv=fiv $
                       ,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims $
                       ,tied=tied,fixed=fixed $
4a2d1650   Jean-Philippe Bernard   improved
5
                       ,polarization=polarization $
b57d58a1   Annie Hughes   addition of force...
6
                       ,no_reset_plugin_structure=no_reset_plugin_structure $
abb4c515   Annie Hughes   fixed logical bug...
7
8
9
                       ,force_reset=force_reset $
                       ,help=help,debug=debug

a4c9699a   Annie Hughes   first commit
10
11
12
13
14
;+
; NAME:
;    dustem_init_params
;
; PURPOSE:
b57d58a1   Annie Hughes   addition of force...
15
16
17
18
19
; This routine initialises free and fixed parameters of both the
; dust models and plugins. It does some basic sanity checking. By
; default, it preserves parameters and related information that are
; already stored in memory. To re-initialise, run with the
; /force_reset keyword
a4c9699a   Annie Hughes   first commit
20
21
22
23
24
;  
; CATEGORY:
;    DustEMWrap, Distributed, Mid-Level, Initialization
;
; CALLING SEQUENCE:
b57d58a1   Annie Hughes   addition of force...
25
26
27
;    dustem_init_params, model, pd, iv, $
;          [,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,lolims=llims,tied=tied,fixed=fixed $
;          ,pol=pol,no_reset_plugin_structure=no_reset_plugin_structure,force_reset=force_reset]
a4c9699a   Annie Hughes   first commit
28
29
30
31
32
33
34
;  
; INPUTS:
;    model -- the dust model
;    pd -- the parameter description vector
;    iv -- the initial values of the free parameters 
;
; OPTIONAL INPUT PARAMETERS:
041ca6b6   Annie Hughes   updated help
35
36
37
38
39
40
;    fpd -- the fixed parameter description vector
;    fiv -- the initial values of the fixed parameters 
;    ulimed -- flag 1/0 if upper limits ON/OFF
;    llimed -- flag 1/0 if lower limits ON/OFF
;    ulims -- upper limit values
;    llims -- lower limit values
b57d58a1   Annie Hughes   addition of force...
41
;
a4c9699a   Annie Hughes   first commit
42
43
44
45
46
47
48
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    None
;
; ACCEPTED KEY-WORDS:
b57d58a1   Annie Hughes   addition of force...
49
;    polarization = check that the model makes a polarisation prediction
4a2d1650   Jean-Philippe Bernard   improved
50
;    no_reset_plugin_structure = if set, does not reset the plugin structure
b57d58a1   Annie Hughes   addition of force...
51
52
53
54
55
56
;    force_reset = if set, resets existing parameters in memory to the
;                  default values of the model (removes all trace of
;                  limits, no-default values, free/fixed from previous
;                  instructions)
;    help      = if set, print this help
;
a4c9699a   Annie Hughes   first commit
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURES AND SUBROUTINES USED:
;   dustem_init_parinfo
;   dustem_init_plugins
;   dustem_init_fixed_params
;  
; EXAMPLES
;
; MODIFICATION HISTORY:
;    Written by AH October-2022
b57d58a1   Annie Hughes   addition of force...
76
77
78
;    Significant changes by AH March-2024:
;       addition of force_reset keyword and nicer management of
;       existing parinfo in memory
a4c9699a   Annie Hughes   first commit
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_init_params'
  goto,the_end
END

Npar=n_elements(pd)
Niv=n_elements(iv)

if Npar ne Niv then begin
   message, 'Number of parameters does not equal number of initial values provided',/info
   stop
end

if keyword_set(ulimed) then begin
   Nulimed=n_elements(ulimed)
   if Npar ne Nulimed then begin
      message, 'Number of parameters does not equal number of upper limit flags',/info
      stop
   end
end

if keyword_set(ulims) then begin
   Nulims=n_elements(ulims)
   if Npar ne Nulims then begin
      message, 'Number of upper limits must equal number of parameters (use ulimed to specify ON/OFF for each parameter)',/info
      stop
   end
4e556e9f   Annie Hughes   Some upper and lo...
110
111
112
113
114
115
116
117
118
119
   ; compare ulims and ivs
   on=where(ulimed eq 1,ct)
   if ct gt 0 then begin
      test=total(ulims[on] lt iv[on])
      if test ne 0 then begin
         message, 'Upper limits are less than requested initial values',/info
         stop
      end
      
   end
a4c9699a   Annie Hughes   first commit
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
end

if keyword_set(llimed) then begin
   Nllimed=n_elements(llimed)
   if Npar ne Nllimed then begin
      message, 'Number of parameters does not equal number of lower limit flags',/info
      stop
   end
end

if keyword_set(llims) then begin
   Nllims=n_elements(llims)
   if Npar ne Nllims then begin
      message, 'Number of lower limits must equal number of parameters (use llimed to specify ON/OFF for each parameter)',/info
      stop
   end
4e556e9f   Annie Hughes   Some upper and lo...
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
   ; compare llims and ivs
   on=where(llimed eq 1,ct)
   if ct gt 0 then begin
      test=total(llims[on] gt iv[on])
      if test ne 0 then begin
         message, 'Lower limits are greater than requested initial values',/info
         stop
      end
      
   end
end

if keyword_set(llims) and keyword_set(ulims) then begin
   on=where(llimed eq 1 and ulimed eq 1,ct)
   if ct gt 0 then begin
      test=total(llims[on] gt ulims[on])
      if test ne 0 then begin
         message, 'Requested lower limits are greater than upper limits',/info
         stop
      end
   end
a4c9699a   Annie Hughes   first commit
157
158
159
160
161
162
163
164
165
166
167
end


if keyword_set(fpd) then begin
   Nfpar=n_elements(fpd)
   Nfiv=n_elements(fiv)
   if Nfpar ne Nfiv then begin
      message, 'Number of fixed parameters does not equal number of fixed initial values',/info
      stop
   end
;== CHECK IF USER HAS MADE SOMETHING FIXED AND FREE
87fc4ae7   Annie Hughes   removed worst bugs
168
   match,strupcase(pd),strupcase(fpd),a,b,count=ct
a4c9699a   Annie Hughes   first commit
169
   if ct ne 0 then begin
b57d58a1   Annie Hughes   addition of force...
170
      message, 'Parameter(s) that are both free and fixed?',/info
a4c9699a   Annie Hughes   first commit
171
172
173
174
175
176
      print,'Matching free parameters: '+pd[a]
      print,'Matching fixed parameters: '+fpd[b]
      stop
   end
end

041ca6b6   Annie Hughes   updated help
177

eef37ef9   Annie Hughes   do not remember, ...
178
exists=dustem_test_model_exists(model,/silent)
041ca6b6   Annie Hughes   updated help
179
180
181
if exists ne 1 then $
   message,'Unknown dust model '+model

eef37ef9   Annie Hughes   do not remember, ...
182
polexists=dustem_test_model_exists(model,/pol,/silent)
041ca6b6   Annie Hughes   updated help
183
184
if keyword_set(polarization) and polexists eq 0 then $
   message,'You have requested polarization mode, but the '+model+' model does not predict polarization'
a4c9699a   Annie Hughes   first commit
185

b57d58a1   Annie Hughes   addition of force...
186
; basic sanity checking of PLUGINS
041ca6b6   Annie Hughes   updated help
187
; to be written
a4c9699a   Annie Hughes   first commit
188

abb4c515   Annie Hughes   fixed logical bug...
189
190
191
192
193
; first decide what needs to be done
fpar_status=ptr_valid((*!dustem_fit).fixed_param_descs)
par_status=ptr_valid((*!dustem_fit).param_descs)
fpd_set = keyword_set(fpd)
pd_set = keyword_set(pd)
b57d58a1   Annie Hughes   addition of force...
194

abb4c515   Annie Hughes   fixed logical bug...
195
196
197
198
old_pd=[] & old_iv=[] & old_parinfo=[]
old_ulimed=[] & old_llimed=[] & old_ulims=[] & old_llims=[]
old_fixed=[] & old_tied=[]
nold=0.
b57d58a1   Annie Hughes   addition of force...
199

abb4c515   Annie Hughes   fixed logical bug...
200
201
202
203
new_pd=[] & new_iv=[] & new_parinfo=[]
new_ulimed=[] & new_llimed=[] & new_ulims=[] & new_llims=[]
new_fixed=[] & new_tied=[]
nnew=0.
b57d58a1   Annie Hughes   addition of force...
204

abb4c515   Annie Hughes   fixed logical bug...
205
206
old_fpd=[] & old_fiv=[]
nold_fixed=0.
b57d58a1   Annie Hughes   addition of force...
207

abb4c515   Annie Hughes   fixed logical bug...
208
209
210
211
212
213
214
215
216
217
new_fpd=[] & new_fiv=[]
nnew_fixed=0.

if fpar_status eq 1 then begin
   old_fpd=(*(*!dustem_fit).fixed_param_descs)
   old_fiv=(*(*!dustem_fit).fixed_param_init_values)
   nold_fixed=n_elements(old_fpd)
end

if par_status eq 1 then begin
b57d58a1   Annie Hughes   addition of force...
218
219
   old_pd=(*(*!dustem_fit).param_descs)
   old_iv=(*(*!dustem_fit).param_init_values)
abb4c515   Annie Hughes   fixed logical bug...
220
   nold=n_elements(old_pd)
b57d58a1   Annie Hughes   addition of force...
221
222
223
224
   old_parinfo=*!dustem_parinfo
   old_fixed=old_parinfo.fixed
   old_tied=old_parinfo.tied
   noldpars=n_elements(old_parinfo)
abb4c515   Annie Hughes   fixed logical bug...
225
226
227
228
229
   old_ulimed=intarr(nold)
   old_llimed=intarr(nold)
   old_ulims=fltarr(nold)
   old_llims=fltarr(nold)
   for i=0,nold-1 do begin
b57d58a1   Annie Hughes   addition of force...
230
231
232
233
234
      old_ulimed[i]=old_parinfo[i].limited[1]
      old_llimed[i]=old_parinfo[i].limited[0]
      old_ulims[i]=(old_parinfo[i].limits[1])*old_iv[i]
      old_llims[i]=(old_parinfo[i].limits[0])*old_iv[i]
   endfor
abb4c515   Annie Hughes   fixed logical bug...
235
236
237
238
239
240
241
242
243
end

if pd_set eq 1 then begin
   new_pd=pd & new_iv=iv & nnew=n_elements(new_pd)
   new_ulimed=intarr(nnew)
   new_llimed=intarr(nnew)
   new_ulims=intarr(nnew)
   new_llims=intarr(nnew)
   new_fixed=intarr(nnew)
05ebb391   Annie Hughes   fixed handling of...
244
   new_tied=strarr(nnew)
abb4c515   Annie Hughes   fixed logical bug...
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
   if keyword_set(ulimed) then new_ulimed=ulimed
   if keyword_set(llimed) then new_llimed=llimed
   if keyword_set(ulims) then new_ulims=ulims
   if keyword_set(llims) then new_llims=llims
   if keyword_set(fixed) then new_fixed=fixed
   if keyword_set(tied) then new_tied=tied
end

if fpd_set eq 1 then begin
   new_fpd=fpd & new_fiv=fiv & nnew_fixed=n_elements(new_fpd)
end

if keyword_set(force_reset) then goto, initialise

if pd_set eq 1 then begin
   if par_status eq 1 then begin
      match,strupcase(new_pd),strupcase(old_pd),a,b,count=ndups
      if ndups gt nold then begin
         message, 'Number of duplicates exceeds number of existing free parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate free parameters found...',/info
      end else if ndups eq nold then begin
         message,'Updating all old values of free parameters...',/info
         old_pd=[]
         old_iv=[]
         old_parinfo=[]
         old_ulimed=[]
         old_llimed=[]
         old_ulims=[]
         old_llims=[]
         old_fixed=[]
         old_tied=[]
      end else if ndups gt 0 and ndups lt nold then begin
87fc4ae7   Annie Hughes   removed worst bugs
279
280
         message,string(ndups)+' matching parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_pd[b[i]]
abb4c515   Annie Hughes   fixed logical bug...
281
282
283
284
285
286
287
288
289
290
         remove,b,old_pd
         remove,b,old_iv
         remove,b,old_parinfo
         remove,b,old_ulimed
         remove,b,old_llimed
         remove,b,old_ulims
         remove,b,old_llims
         remove,b,old_fixed
         remove,b,old_tied
      endif
b57d58a1   Annie Hughes   addition of force...
291
   endif
b57d58a1   Annie Hughes   addition of force...
292

abb4c515   Annie Hughes   fixed logical bug...
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
   if fpar_status eq 1 then begin
      match,strupcase(new_pd),strupcase(old_fpd),a,b,count=ndups
      if ndups gt nold_fixed then begin
         message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate free parameters found in old fixed params...',/info
      end else if ndups eq nold_fixed then begin
         message,'Updating all old values of fixed parameters found in new free params...',/info
         old_fpd=[]
         old_fiv=[]
      end else if ndups gt 0 and ndups lt nold_fixed then begin
         message,string(ndups)+' matching fixed parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_fpd[b[i]]
         remove,b,old_fpd
         remove,b,old_fiv
      endif
   endif
endif
87fc4ae7   Annie Hughes   removed worst bugs
312

abb4c515   Annie Hughes   fixed logical bug...
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
if fpd_set eq 1 then begin
   if par_status eq 1 then begin
      match,strupcase(new_fpd),strupcase(old_pd),a,b,count=ndups
      if ndups gt nold then begin
         message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate fixed parameters found  in old free params...',/info
      end else if ndups eq nold then begin
         message,'Updating all old values of free parameters found in fixed params...',/info
         old_pd=[]
         old_iv=[]
         old_parinfo=[]
         old_ulimed=[]
         old_llimed=[]
         old_ulims=[]
         old_llims=[]
         old_fixed=[]
         old_tied=[]
      end else if ndups gt 0 and ndups lt nold then begin
         message,string(ndups)+' matching parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_pd[b[i]]
         remove,b,old_pd
         remove,b,old_iv
         remove,b,old_parinfo
         remove,b,old_ulimed
         remove,b,old_llimed
         remove,b,old_ulims
         remove,b,old_llims
         remove,b,old_fixed
         remove,b,old_tied
      endif
   endif
87fc4ae7   Annie Hughes   removed worst bugs
346

abb4c515   Annie Hughes   fixed logical bug...
347
   if fpar_status eq 1 then begin
87fc4ae7   Annie Hughes   removed worst bugs
348
      match,strupcase(new_fpd),strupcase(old_fpd),a,b,count=ndups
abb4c515   Annie Hughes   fixed logical bug...
349
350
      if ndups gt nold_fixed then begin
         message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
b57d58a1   Annie Hughes   addition of force...
351
         stop
abb4c515   Annie Hughes   fixed logical bug...
352
353
354
355
356
357
358
359
360
      end else if ndups eq 0 then begin
         message,'No duplicate fixed parameters found in old fixed params...',/info
      end else if ndups eq nold_fixed then begin
         message,'Updating all old values of fixed parameters...',/info
         old_fpd=[]
         old_fiv=[]
      end else if ndups gt 0 and ndups lt nold_fixed then begin
         message,string(ndups)+' matching fixed parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_fpd[b[i]]
b57d58a1   Annie Hughes   addition of force...
361
362
         remove,b,old_fpd
         remove,b,old_fiv
b57d58a1   Annie Hughes   addition of force...
363
      endif
b57d58a1   Annie Hughes   addition of force...
364
   endif
abb4c515   Annie Hughes   fixed logical bug...
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
endif

;merge the pd structures
pd=[new_pd,old_pd]
iv=[new_iv,old_iv]
ulimed=[new_ulimed,old_ulimed]
llimed=[new_llimed,old_llimed]
ulims=[new_ulims,old_ulims]
llims=[new_llims,old_llims]
fixed=[new_fixed,old_fixed]
tied=[new_tied,old_tied]

fpd=[new_fpd,old_fpd]
fiv=[new_fiv,old_fiv]

initialise:
;== INITIALIZE DUST MODEL PARAMETERS and PLUGIN FREE PARAMETERS
; currently this works without checking if PD is non-NULL because we
; require there alway to be a PD vector... one day we should change this
npd=n_elements(pd) & nfpd=n_elements(fpd)

if npd ge 1 then dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims,fixed=fixed,tied=tied
05ebb391   Annie Hughes   fixed handling of...
387

abb4c515   Annie Hughes   fixed logical bug...
388
389
390
391
if npd eq 0 and nold gt 0 then $
   dustem_init_parinfo,pd,iv,/clear

if nfpd ge 1 then dustem_init_fixed_params,fpd,fiv
05ebb391   Annie Hughes   fixed handling of...
392

abb4c515   Annie Hughes   fixed logical bug...
393
394
if nfpd eq 0 and nold_fixed gt 0 then $
   dustem_init_fixed_params,fpd,fiv,/clear
87fc4ae7   Annie Hughes   removed worst bugs
395
396
397
398
399

;== INITIALIZE ANY PLUGINS
IF not keyword_set(no_reset_plugin_structure) THEN BEGIN
   dustem_init_plugins,pd,fpd=fpd
ENDIF
a4c9699a   Annie Hughes   first commit
400

b57d58a1   Annie Hughes   addition of force...
401
if !dustem_verbose NE 0 then message,'Finished initializing free and fixed parameters of dust model and plugins',/info
a4c9699a   Annie Hughes   first commit
402
403
404
405

the_end:

END