[PYSPEX] how to get sector and component index of a parameter

[drogantini]

Hi Jelle, is there a way I can get the sector/component number of a specific parameter?

s.com('hot')
s.par(1, 1, 'nh', 4E-4, thawn=True)
hot_nh = s.par_get(1,2,'norm')

hot_nh has these variables

index – Parameter number
name  – Name of parameter
type – Type of parameter (norm, abun, fitp, cons or text)
desc  – Description of parameter
value  – Value
low  – Lower boundary of parameter range
upp  – Upper boundary of parameter range
err_low  – Lower error boundary
err_upp  – Upper error boundary
free  – True if parameter is free
linked – True if parameter is linked
link_sector – Sector to which parameter is linked
link_comp – Component to which parameter is linked
link_par  – Parameter to which parameter is linked
coupling  – Coupling factor

Many thanks, Daniele

[drogantini]

Follow-up question: can I get which statistic method I am using? I know that there isn't a proper command in SPEX. Wondering if it is easier with PySPEX. Thanks

[jdeplaa]

Hi Jelle, is there a way I can get the sector/component number of a specific parameter?

s.com('hot')
s.par(1, 1, 'nh', 4E-4, thawn=True)
hot_nh = s.par_get(1,2,'norm')

hot_nh has these variables

index – Parameter number
name  – Name of parameter
type – Type of parameter (norm, abun, fitp, cons or text)
desc  – Description of parameter
value  – Value
low  – Lower boundary of parameter range
upp  – Upper boundary of parameter range
err_low  – Lower error boundary
err_upp  – Upper error boundary
free  – True if parameter is free
linked – True if parameter is linked
link_sector – Sector to which parameter is linked
link_comp – Component to which parameter is linked
link_par  – Parameter to which parameter is linked
coupling  – Coupling factor

Many thanks, Daniele

Hi Daniele,

Currently, that is not possible. Could you explain what you would expect from such a feature exactly? Do you want to be able to search the parameter list for certain properties? If you have multiple components, it is very likely to find parameters multiple times (multiple 'nh' parameters, for example). How do you suggest to find the right one?

In the current version of pyspex, the parameters are organised in a nested tree with first a sector number, then the component number, and last the list of parameters. You could loop over them to find the parameter that you need. This tree is linked in the s.mod object:

>>> dir(s.mod.sect[0].comp[0].par[0])

The indices differ from the true sector and component numbers by i-1, because Python starts counting at 0. Otherwise, this contains the entire list of parameters in order. You can use this list to create a function to do what you want.

NOTE: If you just made a number of changes to the parameters, make sure that you run an update command before you start the search to make sure that the pyspex tree is the same as the tree in SPEX. This command copies the settings in SPEX to the pyspex variables that you want to search:

>>> s.mod.update()

Thanks for your feedback. If this results in a useful function, I am interested to add it to the interface!

[jdeplaa]

Follow-up question: can I get which statistic method I am using? I know that there isn't a proper command in SPEX. Wondering if it is easier with PySPEX. Thanks

There isn't a dedicated function for that, but it is easy to get. A copy of the fit settings is kept in the s.opt_fit object. You just need to make sure that the copy is still up-to-date. You can get the current fit statistic with these two commands:

>>> s.opt_fit.update()   # updates the python copy of the SPEX settings
>>> print(s.opt_fit.stat)  # prints the statistic setting
csta

This gives you the general setting: 'csta', 'chi2', or 'wsta'. Recently, it became possible in SPEX to assign a fit statistic per spectrum. This setting does not have a copy in pyspex yet, unfortunately.

[drogantini]

Hi Daniele,

Currently, that is not possible. Could you explain what you would expect from such a feature exactly? Do you want to be able to search the parameter list for certain properties? If you have multiple components, it is very likely to find parameters multiple times (multiple 'nh' parameters, for example). How do you suggest to find the right one?

Hi Jelle, thanks for your quick reply. I am finally writing the Bayesian tool for SPEX using the python interface. It would be great if the parameter class saves the section and the component of the specific parameter. To be clear: I want to transfer the parameter info to a function (based on the piece of code I posted before):

hot_nh = s.par_get(1,2,'norm')
solver = prior_function(model, hot_nh)

where prior_function is something like:

def prior_function(model, par):
      low = log10(par.low)
      spread = log10(par.upp) - log10(par.low)
      [...]
      # get new values from an external algorithm and update the value in SPEX
     s.par(par.sec, par.comp, par.name, new_par_value)

Since I cannot access to the sector and component of the parameter (hypothetically par.sec and par.comp), I provide them manually and the function looks like this:

solver = prior_function(model, hot_nh)
def prior_function(model, sec, comp, par):

It is fine, but somehow it increases the chances of typos/mistakes, especially if the number of parameters is large. Is it complicated to include the sector and component number in the parameter class? I see that you do for the linked parameter (link_sector, link_comp) so I thought it is should be doable. Many thanks

[jdeplaa]

Thanks for the more detailed explanation. This should be relatively easy to do. I will make a note to add it to SPEX 3.08.00.

[drogantini]

s.opt_fit.update() # updates the python copy of the SPEX settings print(s.opt_fit.stat) # prints the statistic setting

This works well, many thanks!

[drogantini]

Hi Jelle,

I am afraid I will ask you a pretty large number of questions in the coming days regarding pyspex. Let me know if you prefer that I send you private emails or if I should keep using this help page (in this case, shall I open a new issue for different questions? It would be easier to search for them). I am trying to get the total flux for a model. I understand that pyspex only gets the flux of single components. Therefore, is the following small script the easy/quick way to get the total flux?

flux = numpy.array(0.0)
# loop over all the components 
for i in range(s.mod.sect[0].ncomp):
    # check if it is additive component
    if s.mod.sect[0].comp[i].add == True:
        # sum the flux
        flux += numpy.array(s.flux_get(1,i+1).photlum)

[jdeplaa]

Hi Daniele,

Yes, it is better to open an issue for each question that you have. That makes it easier to find for others and we keep better track of questions that have been answered. Your last question has issue number #30 now. I will answer it there.

[jdeplaa]

Hi Daniele,

I added the sector and region number to the parameter object in pyspex. It is quite simple, so if you want, you can change it in your local version as well. This is the diff of the file python/pyspex/model.py with the modifications:

--- a/python/pyspex/model.py
+++ b/python/pyspex/model.py
@@ -626,7 +626,11 @@ class Component:

 class Parameter:
     """Class for model parameters.
-    
+
+    :ivar isect: Sector number
+    :vartype index: int
+    :ivar icomp: Component number
+    :vartype index: int
     :ivar index: Parameter number
     :vartype index: int
     :ivar name: Name of parameter
@@ -660,6 +664,8 @@ class Parameter:
     """

     def __init__(self):
+        self.isect = 0        # Sector number
+        self.icomp = 0        # Component number
         self.index = 0        # Parameter number
         self.name = ''        # Name of parameter
         self.type = ''        # Type of parameter (norm, abun, fitp, cons or text)
@@ -689,6 +695,8 @@ class Parameter:
         """
         pyspex_f2py.api_model_parameter.api_par_update(float(isect), float(icomp), float(ipar))

+        self.isect = isect
+        self.icomp = icomp
         self.index = int(pyspex_f2py.api_model_parameter.par_index)
         self.name = pyspex_f2py.api_model_parameter.par_acronym.view('S4').tobytes().decode('utf-8')
         self.name = self.name.strip()

Since this is just a change in the python interface, there is no need to recompile the Fortran part. Let me know whether this works.