Optimization & Simulation functions¶
Fitting spectra¶
Evaluate the model and convolve with response¶
The smallest step in a fitting process is calculating the model, fold the model spectrum through the response matrix, and calculate the likelihood compared to a loaded dataset.
The eval procedure calculates the model and folds it through the response matrix:
call spex%eval()
The interface looks like this:
module subroutine spex_fit_eval(this)
class(sapi), intent(inout) :: this
end subroutine
If you need the fit statistics and degrees of freedom, they
can be returned using the eval_stat procedure:
call spex%eval_stat(stat,dof)
The interface looks like this:
module subroutine spex_fit_eval_stat(this, stat, dof)
class(sapi), intent(inout) :: this
real, intent(out) :: stat !! Statistics value
integer, intent(out) :: dof !! Degrees of freedom
end subroutine
Fitting¶
Spectral fitting is done with the fit procedure:
call spex%fit(ier)
or
call spex%fit(ier,iter=10)
The interface looks like this:
module subroutine spex_fit_fit(this, ier, iter)
class(sapi), intent(inout) :: this
integer, intent(inout) :: ier !!Error status
integer, intent(inout), optional :: iter !!Number of iterations
end subroutine
One can tweak the fits with the commands below. The fit method can be set with:
call spex%fit_method('classical',ier)
Interface:
module subroutine spex_fit_set_method(this, imethod, ier)
class(sapi), intent(inout) :: this
character*(*), intent(in) :: imethod !! Method fitting
integer, intent(out) :: ier !! Error status
end subroutine
The terminal output for each iteration can be set with:
call spex%fit_print(1,ier)
Interface:
module subroutine spex_fit_set_print(this, i, ier)
class(sapi), intent(inout) :: this
integer, intent(in) :: i !!printing the intermediate step (yes if i=1)
integer, intent(out) :: ier !!error status
end subroutine
The fit statistics can be set using:
call spex%fit_stat('cstat',ier)
Interface:
module subroutine spex_fit_set_statistic(this, stat, ier)
class(sapi), intent(inout) :: this
character*(*), intent(in) :: stat !! Fit statistics (chi2, cstat, or wstat)
integer, intent(out) :: ier !! error status
end subroutine
And setting the statistics for each instrument separately:
call spex%fit_stat_inst('cstat',1,1,ier)
Interface:
module subroutine spex_fit_set_statistic_inst(this, stat, inst, ireg, ier)
class(sapi), intent(inout) :: this
character*(*), intent(in) :: stat !! Fit statistics (chi2, cstat, or wstat)
integer, intent(in) :: inst !! Instrument number
integer, intent(in) :: ireg !! Region number
integer, intent(out) :: ier !! error status
end subroutine
There is also a deprecated option to set the weigths for the fit:
call spex%fit_weight('model',ier)
Interface:
module subroutine spex_fit_set_weight(this, fweight, ier)
class(sapi), intent(inout) :: this
character*(*), intent(in) :: fweight !!fit weight: 'data' or 'model'
integer, intent(out) :: ier !!error status
end subroutine
Error calculation¶
The error on a parameter can be calculated with:
call spex%error(isect,icomp,pnam,pvalue,perr_low,perr_upp,perr_lc,perr_cmin,perr_pmin,ier)
Interface:
module subroutine spex_fit_error(this, isect, icomp, pnam, pvalue, perr_low, perr_upp, &
perr_lc, perr_cmin, perr_pmin, ier)
class(sapi), intent(inout) :: this
integer, intent(in) :: isect !! The instrument number
integer, intent(in) :: icomp !! The component number
character*(*), intent(inout) :: pnam !! Parameter name
real(dp), intent(out) :: pvalue !! Parameter value
real(dp), intent(out) :: perr_low !! Lower error
real(dp), intent(out) :: perr_upp !! Upper error
logical, intent(out) :: perr_lc !! Lower chi/c statistics found (True/False)
real(dp), intent(out) :: perr_cmin!! Lowest chi/c statistics found
real(dp), intent(out) :: perr_pmin!! Parameter value for new minimum
integer, intent(inout) :: ier !! status error
end subroutine
Example:
call spex%error(1,1,'norm',pvalue,perr_low,perr_upp,perr_lc,perr_cmin,perr_pmin,ier)
One can also change the target delta-chi:
call spex%error_dchi(dchi)
Interface:
module subroutine spex_fit_set_error_dchi(this,dchi)
class(sapi), intent(inout) :: this
real, intent(in) :: dchi !! Delta chisq of the error to search for
end subroutine
Example:
call spex%error_dchi(1.0)
Setting the error start value can be done with:
call spex%error_start(start)
Interface:
module subroutine spex_fit_set_error_start(this,start)
class(sapi), intent(inout) :: this
real, intent(in) :: start !! Starting value of error calculation
end subroutine
Example:
call spex%error_start(1E+5)
Step¶
Similar to the step command in SPEX, the step process is set up with multiple procedures:
call spex%step_axis(iaxis,isect,icomp,npar,lim1,lim2,nstep,log)
The above sets up one stepping axis for a parameter, including the limits and stepsize.
The interface:
module subroutine spex_fit_step_axis(this, iaxis, isect, icomp, npar, lim1, lim2, nstep, log)
class(sapi), intent(inout) :: this
integer, intent(in) :: iaxis !! Axis number to set
integer, intent(in) :: isect !! Sector number
integer, intent(in) :: icomp !! Component number
character*4, intent(inout) :: npar !! Parameter name
real, intent(in) :: lim1 !! Lower limit of axis
real, intent(in) :: lim2 !! Upper limit of axis
integer, intent(in) :: nstep !! Number of steps
logical, intent(in) :: log !! Axis logarithmic (yes=.true., no=.false.)
end subroutine
Check if the axes are properly defined:
call spex%step_check_axis(isect,icomp,npar,lim1,lim2,nstep,log,ier)
If ier is non-zero, the axes were not set up properly.
The interface:
module subroutine spex_fit_step_check_axis(this, isect, icomp, npar, lim1, lim2, nstep, log, ier)
class(sapi), intent(inout) :: this
integer, intent(in) :: isect !! Sector number
integer, intent(in) :: icomp !! Component number
character*4, intent(inout) :: npar !! Parameter name
real, intent(in) :: lim1 !! Lower limit of axis
real, intent(in) :: lim2 !! Upper limit of axis
integer, intent(in) :: nstep !! Number of steps
logical, intent(in) :: log !! Axis logarithmic (yes=.true., no=.false.)
integer, intent(inout) :: ier
end subroutine
Example:
call spex%step_check_axis(1,1,'norm',)
Set the number of axes for the step session:
call spex%step_dim(ndim,ier)
Interface:
module subroutine spex_fit_step_dim(this, ndim, ier)
class(sapi), intent(inout) :: this
integer, intent(in) :: ndim !! Step dimension
integer, intent(out) :: ier
end subroutine
Setting an output file for step is done with:
call spex%step_file('stepname')
Interface:
module subroutine spex_fit_step_file(this, stpname)
class(sapi), intent(inout) :: this
character*(*), intent(in) :: stpname
end subroutine
The final step is started with:
call spex%step(ier)
Interfaces:
module subroutine spex_fit_step(this,ier,pstep)
class(sapi), intent(inout) :: this
integer, intent(inout) :: ier !error status
type(step), optional, pointer :: pstep
end subroutine
module subroutine spex_fit_step_check(this, ier, pstep)
class(sapi), intent(inout) :: this
integer, intent(inout) :: ier
type(step), pointer, optional :: pstep
end subroutine
Simulating spectra¶
Following the simulate command in SPEX, the Fortran API has procedures to set the different simulation properties if needed.
Simulating noise on the background or source can be switched on and off:
call spex%sim_noise(.false.)
call spex%sim_bnoise(.false.)
The interfaces for these procedures:
module subroutine spex_fit_sim_noise(this,noise)
class(sapi), intent(inout) :: this
logical, intent(in) :: noise !! Set source noise for simulation (true/false)
end subroutine
module subroutine spex_fit_sim_bnoise(this,bnoise)
class(sapi), intent(inout) :: this
logical, intent(in) :: bnoise !! Set background noise for simulation (true/false)
end subroutine
The random seed can be set to a value or randomized with:
call spex%sim_seed(.true.,0) ! Set random seed to a random number based on computer clock
call spex%sim_seed(.false.,12345678) ! Set random seed to a particular number
The interface:
module subroutine spex_fit_sim_seed(this,rand,seed)
class(sapi), intent(inout) :: this
logical, intent(in) :: rand !! Should the seed be randomly set? True/False
integer, intent(in) :: seed !! Set seed manually if rand is False
end subroutine
Additional systematic errors can be added (not recommended):
call spex%sim_syserr(sys_src,sys_bkg)
Interface:
module subroutine spex_fit_sim_syserr(this,sys_src,sys_bkg)
class(sapi), intent(inout) :: this
real, intent(in) :: sys_src !! Systematic error on the source spectrum
real, intent(in) :: sys_bkg !! Systematic error on the background spectrum
end subroutine
Select the instruments that need to be simulated:
call spex%sim_inst(in1,in2,ic1,ic2,ier)
The interface:
module subroutine spex_fit_sim_inst(this,in1,in2,ic1,ic2,ier)
class(sapi), intent(inout) :: this
integer, intent(in) :: in1
integer, intent(in) :: in2
integer, intent(in) :: ic1
integer, intent(in) :: ic2
integer, intent(inout) :: ier
end subroutine
And finally, the main simulate command:
call spex%sim(exposure)
The interface:
module subroutine spex_fit_simulate(this, exposure)
class(sapi), intent(inout) :: this
real(dp), intent(in) :: exposure !! Exposure time
end subroutine
Example:
call spex%sim(1E+6)