Commit abb4c515047cfe0f1788ef1cb3905bd6a6d1554b

Authored by Annie Hughes
1 parent 08e3260e
Exists in master

fixed logical bugs in initialising params

src/idl/dustem_init_fixed_params.pro
1 -PRO dustem_init_fixed_params,fpd,fiv,help=help 1 +PRO dustem_init_fixed_params,fpd,fiv,help=help,clear=clear
2 2
3 ;+ 3 ;+
4 ; NAME: 4 ; NAME:
@@ -28,6 +28,7 @@ PRO dustem_init_fixed_params,fpd,fiv,help=help @@ -28,6 +28,7 @@ PRO dustem_init_fixed_params,fpd,fiv,help=help
28 ; 28 ;
29 ; ACCEPTED KEY-WORDS: 29 ; ACCEPTED KEY-WORDS:
30 ; help = if set, print this help 30 ; help = if set, print this help
  31 +; clear = if set, clear all existing fixed parameters
31 ; 32 ;
32 ; COMMON BLOCKS: 33 ; COMMON BLOCKS:
33 ; None 34 ; None
@@ -57,9 +58,16 @@ IF keyword_set(help) THEN BEGIN @@ -57,9 +58,16 @@ IF keyword_set(help) THEN BEGIN
57 goto,the_end 58 goto,the_end
58 END 59 END
59 60
  61 +if keyword_set(clear) then begin
  62 +; message,'Clearing the fixed parameters',/info
  63 +; stop
  64 + (*!dustem_fit).fixed_param_descs=ptr_new()
  65 + (*!dustem_fit).fixed_param_init_values=ptr_new()
  66 + goto,the_end
  67 +end
  68 +
60 ;Save the values of fit parameter descriptions and values 69 ;Save the values of fit parameter descriptions and values
61 need_save=0 70 need_save=0
62 -  
63 ;This routine will use (*!dustem_fit).param_descs) 71 ;This routine will use (*!dustem_fit).param_descs)
64 ;Therefore, if the variable already exists, its content is saved, and will be put back in place later 72 ;Therefore, if the variable already exists, its content is saved, and will be put back in place later
65 IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN 73 IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN
@@ -69,7 +77,8 @@ IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN @@ -69,7 +77,8 @@ IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN
69 need_save=1 77 need_save=1
70 ENDIF 78 ENDIF
71 79
72 -;temporarily use the regular parameter system variables to update 80 +;TEMPORARILY put FPD and FIV into the regular parameter system
  81 +;variables. This allows us to the use the dustem_set_params routine
73 (*!dustem_fit).param_descs=ptr_new(fpd) 82 (*!dustem_fit).param_descs=ptr_new(fpd)
74 (*!dustem_fit).param_init_values=ptr_new(fiv) 83 (*!dustem_fit).param_init_values=ptr_new(fiv)
75 ;Trimmed JPB on 08 March 2024 84 ;Trimmed JPB on 08 March 2024
@@ -77,9 +86,8 @@ ENDIF @@ -77,9 +86,8 @@ ENDIF
77 ind=where(fpd NE '',count) 86 ind=where(fpd NE '',count)
78 87
79 ;fixing the parameters 88 ;fixing the parameters
80 -IF count NE 0 THEN BEGIN 89 +IF count NE 0 THEN $
81 dustem_set_params,[fiv] 90 dustem_set_params,[fiv]
82 -ENDIF  
83 91
84 ;set the fixed system variables, used to keep memory of the fixed 92 ;set the fixed system variables, used to keep memory of the fixed
85 ;parameters description and values 93 ;parameters description and values
@@ -92,6 +100,11 @@ IF need_save EQ 1 THEN BEGIN @@ -92,6 +100,11 @@ IF need_save EQ 1 THEN BEGIN
92 (*!dustem_fit).param_init_values=ptr_new(iv) 100 (*!dustem_fit).param_init_values=ptr_new(iv)
93 ;(*!dustem_fit).param_func=ptr_new(find) 101 ;(*!dustem_fit).param_func=ptr_new(find)
94 ENDIF 102 ENDIF
  103 +IF need_save EQ 0 THEN BEGIN
  104 + (*!dustem_fit).param_descs=ptr_new()
  105 + (*!dustem_fit).param_init_values=ptr_new()
  106 + ;(*!dustem_fit).param_func=ptr_new(find)
  107 +ENDIF
95 108
96 the_end: 109 the_end:
97 END 110 END
src/idl/dustem_init_params.pro
@@ -4,7 +4,9 @@ PRO dustem_init_params,model,pd,iv, $ @@ -4,7 +4,9 @@ PRO dustem_init_params,model,pd,iv, $
4 ,tied=tied,fixed=fixed $ 4 ,tied=tied,fixed=fixed $
5 ,polarization=polarization $ 5 ,polarization=polarization $
6 ,no_reset_plugin_structure=no_reset_plugin_structure $ 6 ,no_reset_plugin_structure=no_reset_plugin_structure $
7 - ,force_reset=force_reset 7 + ,force_reset=force_reset $
  8 + ,help=help,debug=debug
  9 +
8 ;+ 10 ;+
9 ; NAME: 11 ; NAME:
10 ; dustem_init_params 12 ; dustem_init_params
@@ -184,148 +186,217 @@ if keyword_set(polarization) and polexists eq 0 then $ @@ -184,148 +186,217 @@ if keyword_set(polarization) and polexists eq 0 then $
184 ; basic sanity checking of PLUGINS 186 ; basic sanity checking of PLUGINS
185 ; to be written 187 ; to be written
186 188
187 -;== PLAY NICE -- check what already exists in the parameter structure  
188 -defsysv,'!dustem_parinfo',exists=par_status 189 +; first decide what needs to be done
  190 +fpar_status=ptr_valid((*!dustem_fit).fixed_param_descs)
  191 +par_status=ptr_valid((*!dustem_fit).param_descs)
  192 +fpd_set = keyword_set(fpd)
  193 +pd_set = keyword_set(pd)
189 194
190 -if keyword_set(force_reset) and !dustem_verbose NE 0 then $  
191 - message, 'Reinitialising all previous parameter options (FORCE_RESET=1)',/info 195 +old_pd=[] & old_iv=[] & old_parinfo=[]
  196 +old_ulimed=[] & old_llimed=[] & old_ulims=[] & old_llims=[]
  197 +old_fixed=[] & old_tied=[]
  198 +nold=0.
192 199
193 -if par_status eq 0 and !dustem_verbose NE 0 then $  
194 - message, 'No existing parameters found',/info 200 +new_pd=[] & new_iv=[] & new_parinfo=[]
  201 +new_ulimed=[] & new_llimed=[] & new_ulims=[] & new_llims=[]
  202 +new_fixed=[] & new_tied=[]
  203 +nnew=0.
195 204
196 -if par_status eq 1 and not keyword_set(force_reset) then begin  
197 - if !dustem_verbose NE 0 then message,'Updating/adding model parameters',/info  
198 - new_pd=pd & new_iv=iv & nnewpars=n_elements(new_pd)  
199 - new_ulimed=intarr(nnewpars)  
200 - new_llimed=intarr(nnewpars)  
201 - new_ulims=intarr(nnewpars)  
202 - new_llims=intarr(nnewpars)  
203 - new_fixed=intarr(nnewpars)  
204 - new_tied=intarr(nnewpars)  
205 - if keyword_set(ulimed) then new_ulimed=ulimed  
206 - if keyword_set(llimed) then new_llimed=llimed  
207 - if keyword_set(ulims) then new_ulims=ulims  
208 - if keyword_set(llims) then new_llims=llims  
209 - if keyword_set(fixed) then new_fixed=fixed  
210 - if keyword_set(tied) then new_tied=tied 205 +old_fpd=[] & old_fiv=[]
  206 +nold_fixed=0.
211 207
  208 +new_fpd=[] & new_fiv=[]
  209 +nnew_fixed=0.
  210 +
  211 +if fpar_status eq 1 then begin
  212 + old_fpd=(*(*!dustem_fit).fixed_param_descs)
  213 + old_fiv=(*(*!dustem_fit).fixed_param_init_values)
  214 + nold_fixed=n_elements(old_fpd)
  215 +end
  216 +
  217 +if par_status eq 1 then begin
212 old_pd=(*(*!dustem_fit).param_descs) 218 old_pd=(*(*!dustem_fit).param_descs)
213 old_iv=(*(*!dustem_fit).param_init_values) 219 old_iv=(*(*!dustem_fit).param_init_values)
  220 + nold=n_elements(old_pd)
214 old_parinfo=*!dustem_parinfo 221 old_parinfo=*!dustem_parinfo
215 old_fixed=old_parinfo.fixed 222 old_fixed=old_parinfo.fixed
216 old_tied=old_parinfo.tied 223 old_tied=old_parinfo.tied
217 noldpars=n_elements(old_parinfo) 224 noldpars=n_elements(old_parinfo)
218 - old_ulimed=intarr(noldpars)  
219 - old_llimed=intarr(noldpars)  
220 - old_ulims=fltarr(noldpars)  
221 - old_llims=fltarr(noldpars)  
222 - for i=0,noldpars-1 do begin 225 + old_ulimed=intarr(nold)
  226 + old_llimed=intarr(nold)
  227 + old_ulims=fltarr(nold)
  228 + old_llims=fltarr(nold)
  229 + for i=0,nold-1 do begin
223 old_ulimed[i]=old_parinfo[i].limited[1] 230 old_ulimed[i]=old_parinfo[i].limited[1]
224 old_llimed[i]=old_parinfo[i].limited[0] 231 old_llimed[i]=old_parinfo[i].limited[0]
225 old_ulims[i]=(old_parinfo[i].limits[1])*old_iv[i] 232 old_ulims[i]=(old_parinfo[i].limits[1])*old_iv[i]
226 old_llims[i]=(old_parinfo[i].limits[0])*old_iv[i] 233 old_llims[i]=(old_parinfo[i].limits[0])*old_iv[i]
227 endfor 234 endfor
228 - ; check for duplicates  
229 - noldpars=n_elements(old_pd)  
230 - match,strupcase(new_pd),strupcase(old_pd),a,b,count=ndups  
231 - if ndups gt noldpars then begin  
232 - message, 'Number of duplicates exceeds number of existing parameters. This should not happen ?'/info  
233 - stop  
234 - endif  
235 - if ndups gt 0 and ndups le noldpars then begin  
236 - if !dustem_verbose NE 0 then begin 235 +end
  236 +
  237 +if pd_set eq 1 then begin
  238 + new_pd=pd & new_iv=iv & nnew=n_elements(new_pd)
  239 + new_ulimed=intarr(nnew)
  240 + new_llimed=intarr(nnew)
  241 + new_ulims=intarr(nnew)
  242 + new_llims=intarr(nnew)
  243 + new_fixed=intarr(nnew)
  244 + new_tied=intarr(nnew)
  245 + if keyword_set(ulimed) then new_ulimed=ulimed
  246 + if keyword_set(llimed) then new_llimed=llimed
  247 + if keyword_set(ulims) then new_ulims=ulims
  248 + if keyword_set(llims) then new_llims=llims
  249 + if keyword_set(fixed) then new_fixed=fixed
  250 + if keyword_set(tied) then new_tied=tied
  251 +end
  252 +
  253 +if fpd_set eq 1 then begin
  254 + new_fpd=fpd & new_fiv=fiv & nnew_fixed=n_elements(new_fpd)
  255 +end
  256 +
  257 +if keyword_set(force_reset) then goto, initialise
  258 +
  259 +if pd_set eq 1 then begin
  260 + if par_status eq 1 then begin
  261 + match,strupcase(new_pd),strupcase(old_pd),a,b,count=ndups
  262 + if ndups gt nold then begin
  263 + message, 'Number of duplicates exceeds number of existing free parameters. This should not happen ?'/info
  264 + stop
  265 + end else if ndups eq 0 then begin
  266 + message,'No duplicate free parameters found...',/info
  267 + end else if ndups eq nold then begin
  268 + message,'Updating all old values of free parameters...',/info
  269 + old_pd=[]
  270 + old_iv=[]
  271 + old_parinfo=[]
  272 + old_ulimed=[]
  273 + old_llimed=[]
  274 + old_ulims=[]
  275 + old_llims=[]
  276 + old_fixed=[]
  277 + old_tied=[]
  278 + end else if ndups gt 0 and ndups lt nold then begin
237 message,string(ndups)+' matching parameters found. Updating values...',/info 279 message,string(ndups)+' matching parameters found. Updating values...',/info
238 for i=0,ndups-1 do print,old_pd[b[i]] 280 for i=0,ndups-1 do print,old_pd[b[i]]
239 - end  
240 - remove,b,old_pd  
241 - remove,b,old_iv  
242 - remove,b,old_parinfo  
243 - remove,b,old_ulimed  
244 - remove,b,old_llimed  
245 - remove,b,old_ulims  
246 - remove,b,old_llims  
247 - remove,b,old_fixed  
248 - remove,b,old_tied 281 + remove,b,old_pd
  282 + remove,b,old_iv
  283 + remove,b,old_parinfo
  284 + remove,b,old_ulimed
  285 + remove,b,old_llimed
  286 + remove,b,old_ulims
  287 + remove,b,old_llims
  288 + remove,b,old_fixed
  289 + remove,b,old_tied
  290 + endif
249 endif 291 endif
250 - ;merge the structures  
251 - pd=[new_pd,old_pd]  
252 - iv=[new_iv,old_iv]  
253 - ulimed=[new_ulimed,old_ulimed]  
254 - llimed=[new_llimed,old_llimed]  
255 - ulims=[new_ulims,old_ulims]  
256 - llims=[new_llims,old_llims]  
257 - fixed=[new_fixed,old_fixed]  
258 - tied=[new_tied,old_tied]  
259 -endif  
260 292
  293 + if fpar_status eq 1 then begin
  294 + match,strupcase(new_pd),strupcase(old_fpd),a,b,count=ndups
  295 + if ndups gt nold_fixed then begin
  296 + message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
  297 + stop
  298 + end else if ndups eq 0 then begin
  299 + message,'No duplicate free parameters found in old fixed params...',/info
  300 + end else if ndups eq nold_fixed then begin
  301 + message,'Updating all old values of fixed parameters found in new free params...',/info
  302 + old_fpd=[]
  303 + old_fiv=[]
  304 + end else if ndups gt 0 and ndups lt nold_fixed then begin
  305 + message,string(ndups)+' matching fixed parameters found. Updating values...',/info
  306 + for i=0,ndups-1 do print,old_fpd[b[i]]
  307 + remove,b,old_fpd
  308 + remove,b,old_fiv
  309 + endif
  310 + endif
  311 +endif
261 312
262 -;; stop  
263 -;; message,'line 262',/info  
264 -;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)  
265 -;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)  
266 -;; stop 313 +if fpd_set eq 1 then begin
  314 + if par_status eq 1 then begin
  315 + match,strupcase(new_fpd),strupcase(old_pd),a,b,count=ndups
  316 + if ndups gt nold then begin
  317 + message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
  318 + stop
  319 + end else if ndups eq 0 then begin
  320 + message,'No duplicate fixed parameters found in old free params...',/info
  321 + end else if ndups eq nold then begin
  322 + message,'Updating all old values of free parameters found in fixed params...',/info
  323 + old_pd=[]
  324 + old_iv=[]
  325 + old_parinfo=[]
  326 + old_ulimed=[]
  327 + old_llimed=[]
  328 + old_ulims=[]
  329 + old_llims=[]
  330 + old_fixed=[]
  331 + old_tied=[]
  332 + end else if ndups gt 0 and ndups lt nold then begin
  333 + message,string(ndups)+' matching parameters found. Updating values...',/info
  334 + for i=0,ndups-1 do print,old_pd[b[i]]
  335 + remove,b,old_pd
  336 + remove,b,old_iv
  337 + remove,b,old_parinfo
  338 + remove,b,old_ulimed
  339 + remove,b,old_llimed
  340 + remove,b,old_ulims
  341 + remove,b,old_llims
  342 + remove,b,old_fixed
  343 + remove,b,old_tied
  344 + endif
  345 + endif
267 346
268 -;== INITIALIZE DUST MODEL PARAMETERS and PLUGIN FREE PARAMETERS  
269 -dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims,fixed=fixed,tied=tied  
270 -  
271 -;; stop  
272 -;; message,'line 271',/info  
273 -;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)  
274 -;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)  
275 -;; stop  
276 -  
277 -;== Initialize ANY FIXED PARAMETERS FOR DUST MODEL AND PLUGINS  
278 -IF keyword_set(fpd) THEN BEGIN  
279 -;== PLAY NICE -- check what already exists in the fixed parameter structure  
280 - fpar_status=ptr_valid((*!dustem_fit).fixed_param_descs)  
281 -  
282 - if keyword_set(force_reset) and !dustem_verbose NE 0 then $  
283 - message, 'Reinitialising all previous fixed parameter options (FORCE_RESET=1)',/info  
284 -  
285 - if fpar_status eq 0 and !dustem_verbose NE 0 then $  
286 - message, 'No existing fixed parameters found',/info  
287 -  
288 - if fpar_status eq 1 and not keyword_set(force_reset) then begin  
289 - if !dustem_verbose NE 0 then message,'Updating/adding fixed parameters',/info  
290 - new_fpd=fpd & new_fiv=fiv & nnewfixed=n_elements(new_fpd)  
291 -  
292 - old_fpd=(*(*!dustem_fit).fixed_param_descs)  
293 - old_fiv=(*(*!dustem_fit).fixed_param_init_values)  
294 - noldfixed=n_elements(old_fpd)  
295 - 347 + if fpar_status eq 1 then begin
296 match,strupcase(new_fpd),strupcase(old_fpd),a,b,count=ndups 348 match,strupcase(new_fpd),strupcase(old_fpd),a,b,count=ndups
297 - if ndups gt noldfixed then begin  
298 - message, 'Number of duplicates exceeds number of fixed parameters. This should not happen ?'/info 349 + if ndups gt nold_fixed then begin
  350 + message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
299 stop 351 stop
300 - endif  
301 - if ndups gt 0 and ndups le noldfixed then begin  
302 - if !dustem_verbose NE 0 then message,string(ndups)+' matching fixed parameters found. Updating values for '+old_fpd[b],/info 352 + end else if ndups eq 0 then begin
  353 + message,'No duplicate fixed parameters found in old fixed params...',/info
  354 + end else if ndups eq nold_fixed then begin
  355 + message,'Updating all old values of fixed parameters...',/info
  356 + old_fpd=[]
  357 + old_fiv=[]
  358 + end else if ndups gt 0 and ndups lt nold_fixed then begin
  359 + message,string(ndups)+' matching fixed parameters found. Updating values...',/info
  360 + for i=0,ndups-1 do print,old_fpd[b[i]]
303 remove,b,old_fpd 361 remove,b,old_fpd
304 remove,b,old_fiv 362 remove,b,old_fiv
305 endif 363 endif
306 - fpd=[new_fpd,old_fpd]  
307 - fiv=[new_fiv,old_fiv]  
308 endif 364 endif
309 - dustem_init_fixed_params,fpd,fiv  
310 -ENDIF  
311 -;; stop  
312 -;; message,'line 312',/info  
313 -;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)  
314 -;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)  
315 -;; stop 365 +endif
  366 +
  367 +;merge the pd structures
  368 +pd=[new_pd,old_pd]
  369 +iv=[new_iv,old_iv]
  370 +ulimed=[new_ulimed,old_ulimed]
  371 +llimed=[new_llimed,old_llimed]
  372 +ulims=[new_ulims,old_ulims]
  373 +llims=[new_llims,old_llims]
  374 +fixed=[new_fixed,old_fixed]
  375 +tied=[new_tied,old_tied]
  376 +
  377 +fpd=[new_fpd,old_fpd]
  378 +fiv=[new_fiv,old_fiv]
  379 +
  380 +initialise:
  381 +;== INITIALIZE DUST MODEL PARAMETERS and PLUGIN FREE PARAMETERS
  382 +; currently this works without checking if PD is non-NULL because we
  383 +; require there alway to be a PD vector... one day we should change this
  384 +npd=n_elements(pd) & nfpd=n_elements(fpd)
  385 +
  386 +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
  387 +if npd eq 0 and nold gt 0 then $
  388 + dustem_init_parinfo,pd,iv,/clear
  389 +
  390 +if nfpd ge 1 then dustem_init_fixed_params,fpd,fiv
  391 +if nfpd eq 0 and nold_fixed gt 0 then $
  392 + dustem_init_fixed_params,fpd,fiv,/clear
316 393
317 ;== INITIALIZE ANY PLUGINS 394 ;== INITIALIZE ANY PLUGINS
318 IF not keyword_set(no_reset_plugin_structure) THEN BEGIN 395 IF not keyword_set(no_reset_plugin_structure) THEN BEGIN
319 dustem_init_plugins,pd,fpd=fpd 396 dustem_init_plugins,pd,fpd=fpd
320 ENDIF 397 ENDIF
321 -;; stop  
322 -;; message,'line 322',/info  
323 -;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)  
324 -;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)  
325 -;; stop  
326 398
327 if !dustem_verbose NE 0 then message,'Finished initializing free and fixed parameters of dust model and plugins',/info 399 if !dustem_verbose NE 0 then message,'Finished initializing free and fixed parameters of dust model and plugins',/info
328 -;;stop  
329 400
330 the_end: 401 the_end:
331 402
src/idl/dustem_init_parinfo.pro
1 PRO dustem_init_parinfo,pd,iv, $ 1 PRO dustem_init_parinfo,pd,iv, $
2 fixed=fixed,tied=tied, $ 2 fixed=fixed,tied=tied, $
3 up_limited=up_limited,lo_limited=lo_limited,up_limits=up_limits,lo_limits=lo_limits, $ 3 up_limited=up_limited,lo_limited=lo_limited,up_limits=up_limits,lo_limits=lo_limits, $
4 - relstep=relstep,step=step 4 + relstep=relstep,step=step, $
  5 + clear=clear
5 6
6 ;+ 7 ;+
7 ; NAME: 8 ; NAME:
@@ -48,6 +49,16 @@ IF keyword_set(help) THEN BEGIN @@ -48,6 +49,16 @@ IF keyword_set(help) THEN BEGIN
48 goto,the_end 49 goto,the_end
49 END 50 END
50 51
  52 +if keyword_set(clear) then begin
  53 + message,'Normally you should never get here before we overhaul the parameter description structure',/info
  54 + stop
  55 + (*!dustem_fit).param_descs=ptr_new()
  56 + (*!dustem_fit).param_init_values=ptr_new()
  57 + (*!dustem_fit).current_param_values = ptr_new()
  58 + defsysv,'!dustem_parinfo',ptr_new()
  59 + goto, the_end
  60 +end
  61 +
51 62
52 (*!dustem_fit).param_descs=ptr_new(pd) 63 (*!dustem_fit).param_descs=ptr_new(pd)
53 64