Skip to content

Commit

Permalink
fixing par deps with units
Browse files Browse the repository at this point in the history
  • Loading branch information
andreatramacere committed Apr 18, 2024
1 parent b60e76e commit d869587
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 55 deletions.
26 changes: 19 additions & 7 deletions jetset/base_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -258,15 +258,22 @@ def load_model(cls, file_name):
raise RuntimeError(e)

def _fix_par_dep_on_load(self,):
#print("\n \n ========> fix dep on load START")
for p in self.parameters.par_array:
if p._is_dependent is True and p._linked is False:
#print('==> _master_par_list',p._master_par_list)
#print('==> _depending_par_expr',p._depending_par_expr)
_master_par_list=copy.deepcopy(p._master_par_list)
#print('==> _master_par_list',p._master_par_list, " for", p.name)
#print('==> _depending_par_expr',p._depending_par_expr, " for", p.name)
_master_par_list=[p for p in p._master_par_list]
_depending_par_expr=copy.deepcopy(p._depending_par_expr)
p.reset_dependencies()
self.make_dependent_par(p.name, _master_par_list, _depending_par_expr,set_par_expr_source_code=True)


#def _set_pars_dep(self):
# for p in self.parameters.par_array:
# if


#def _set_pars_dep(self):
# for p in self.parameters.par_array:
# if
Expand Down Expand Up @@ -356,7 +363,8 @@ def _test_par_expr(self,master_par_list,par_expr):
except:
raise RuntimeError('the parameter expression is not valid')

def make_dependent_par(self, par, depends_on, par_expr,verbose=True,set_par_expr_source_code=True):
def make_dependent_par(self, par, depends_on, par_expr,verbose=True,set_par_expr_source_code=True,master_pars=None):
#print("\n ===> make par: ",par, "depending on : ",depends_on, " START")
master_par_list = depends_on

dep_par=self.parameters.get_par_by_name(par)
Expand All @@ -380,16 +388,20 @@ def make_dependent_par(self, par, depends_on, par_expr,verbose=True,set_par_expr
message='problem with parameter name: %s'%p
message+='\nexception:%s'%str(e)
raise RuntimeError(message)

for p in master_par_list:

m = self.parameters.get_par_by_name(p)
if m._is_dependent is False:
m.val=m.val
if m._is_dependent is False and m.par_type == 'user_defined':
try:
m.val=m.val
except:
pass
if set_par_expr_source_code is True:
dep_par._set_par_expr_source_code()
if verbose is True:
dep_par.par_expression_source_code
#print(" ===> make par: ",par, "depending on : ",depends_on, " END\n")

def add_user_par(self,name,val,units='',val_min=None,val_max=None):
self.parameters.add_par(ModelParameter(name=name,units=units,val=val,val_min=val_min,val_max=val_max,par_type='user_defined'))
Expand Down
1 change: 0 additions & 1 deletion jetset/jet_paramters.py
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,6 @@ def set(self,**keywords):
super(JetParameter,self).set(**keywords )

if 'val' in keywords.keys():
#print('val',keywords['val'],'to',(self.name))
self.assign_val_to_jetkernel(self.name,keywords['val'])

#This instruction is not requested since it is performed in
Expand Down
6 changes: 3 additions & 3 deletions jetset/model_manager.py
Original file line number Diff line number Diff line change
Expand Up @@ -364,9 +364,9 @@ def _build_model(c):
#print(p.name,p._root_par,[p.model],p._linked_root_model,p.immutable)
c.parameters.link_par(p._root_par.name,[p.model.name],p._linked_root_model.name)

for m in c.components.components_list:
if isinstance(m,Jet):
m._fix_par_dep_on_load()
#for m in c.components.components_list:
# if isinstance(m,Jet):
# m._fix_par_dep_on_load()
if isinstance(c, Model):
c.eval()
return c
Expand Down
7 changes: 4 additions & 3 deletions jetset/model_parameters.py
Original file line number Diff line number Diff line change
Expand Up @@ -326,6 +326,7 @@ def reset_dependencies(self):
self._is_dependent = False
self._master_par_list=[]
self._depending_pars=[]
self._master_pars=[]
self.par_expr=None


Expand All @@ -336,6 +337,7 @@ def _add_depending_par(self,par):

def _add_master_par(self,par):
if par not in self._master_pars :
print("adding par:",par.name,"to ",self.name)
self._master_pars.append(par)

#def make_dependent_par(self, master_par, func, root_model=None):
Expand Down Expand Up @@ -408,7 +410,7 @@ def _eval_par_func(self):
if _user_par_.adimensional:
_par_values[ID] = _user_par_.val_lin
else:
_par_values[ID] = _user_par_.val_lin*_user_par_.units
_par_values[ID] = _user_par_.val_lin*u.Unit(str(_user_par_.units))
exec(_user_par_.name + '=_par_values[ID]')
res = eval(self._depending_par_expr)
elif callable(self._depending_par_expr) is True:
Expand All @@ -417,7 +419,7 @@ def _eval_par_func(self):
if _user_par_.adimensional:
_par_values[_user_par_.name] = _user_par_.val_lin
else:
_par_values[_user_par_.name] = _user_par_.val_lin*_user_par_.units
_par_values[_user_par_.name] = _user_par_.val_lin*u.Unit(str(_user_par_.units))
res=self._depending_par_expr(**_par_values)

_unit = None
Expand Down Expand Up @@ -469,7 +471,6 @@ def set(self, *args, skip_dep_par_warning=False, **keywords):
self._val.val = keywords[kw]
if self._depending_pars is not []:
for p in self._depending_pars:
#print("==> setting dep par",p.name, 'to',p._func(),'p=',p,'master',self,'name',self.name)
p.set(val=p._func(),skip_dep_par_warning=True)

elif kw == 'log':
Expand Down
56 changes: 55 additions & 1 deletion jetset/tests/test_depending_parameters.py
Original file line number Diff line number Diff line change
Expand Up @@ -178,4 +178,58 @@ def par_func(R0,B0,R_H,m_B):

np.testing.assert_allclose(new_fit_model.jet_leptonic.parameters.B.val, par_func(B0,R0,R_H,m_B))


def test_dep_par_fit_model(self,plot=False):
import numpy as np
from jetset.model_manager import FitModel
from jetset.jet_model import Jet
j=Jet()

j.add_EC_component(['EC_BLR','EC_Disk','EC_DT'],disk_type='MultiBB')

#kaspi+ 2007:https://iopscience.iop.org/article/10.1086/512094/pdf
j.make_dependent_par(par='R_BLR_in', depends_on=['L_Disk'], par_expr='1E17*(L_Disk/1E45)**0.5')

j.make_dependent_par(par='R_BLR_out', depends_on=['R_BLR_in'], par_expr='R_BLR_in*1.1')

#Cleary+ 2007:https://iopscience.iop.org/article/10.1086/511969/pdf
j.make_dependent_par(par='R_DT', depends_on=['L_Disk'], par_expr='2.5E18*(L_Disk/1E45)**0.5')

j.add_user_par(name='theta_open',val=5, units='deg')
def f_par(R_H,theta_open):
return np.tan( theta_open)*R_H

j.make_dependent_par(par='R', depends_on=['R_H','theta_open'],
par_expr=f_par)

j.eval()

R_H_val=1E21
j.parameters.R_H.val=R_H_val
R_val=j.parameters.R.val

j.save_model('test_jet_EC.pkl')

fit_model = FitModel(jet=j, name='EC-best-fit-lsb', template=None)
fit_model.save_model('test_fit_model_EC.pkl')

new_jet=Jet.load_model('test_jet_EC.pkl')

new_jet.eval()
print("units of theta_open",new_jet.parameters.theta_open.units)


new_jet.parameters.R_H.val=R_H_val

assert(new_jet.parameters.R.val==R_val)
new_jet.show_model()

new_fit_model=FitModel.load_model('test_fit_model_EC.pkl')


new_fit_model.eval()
print("units of theta_open",new_fit_model.jet_leptonic.parameters.theta_open.units)

new_fit_model.jet_leptonic.parameters.R_H.val=R_H_val

assert(new_fit_model.jet_leptonic.parameters.R.val==R_val)
new_fit_model.show_model()
41 changes: 1 addition & 40 deletions jetset/tests/test_jet_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ def integration_suite(self,plot=False):
self.test_build_bessel(plot=plot)

self.test_set_N_from_nuFnu(plot=plot)
self.test_EC(plot=plot)
#self.test_EC(plot=plot)

def test_build_bessel(self,plot=False):
print('--------> test_build_bessel',plot)
Expand Down Expand Up @@ -55,46 +55,7 @@ def test_set_N_from_nuFnu(self,plot=False):
y = j.eval(nu=[nu], get_model=True)
np.testing.assert_allclose(y, nuFnu, rtol=1E-2)


def test_EC(self,plot=False):
print('--------> test_EC', plot)
j=Jet()
j.add_EC_component(['EC_BLR','EC_Disk','EC_DT'],disk_type='MultiBB')

#kaspi+ 2007:https://iopscience.iop.org/article/10.1086/512094/pdf
j.make_dependent_par(par='R_BLR_in', depends_on=['L_Disk'], par_expr='1E17*(L_Disk/1E45)**0.5')

j.make_dependent_par(par='R_BLR_out', depends_on=['R_BLR_in'], par_expr='R_BLR_in*1.1')

#Cleary+ 2007:https://iopscience.iop.org/article/10.1086/511969/pdf
j.make_dependent_par(par='R_DT', depends_on=['L_Disk'], par_expr='2.5E18*(L_Disk/1E45)**0.5')

j.add_user_par(name='theta_open',val=5, units='deg')
def f_par(R_H,theta_open):
return np.tan( theta_open)*R_H

j.make_dependent_par(par='R', depends_on=['R_H','theta_open'],
par_expr=f_par)

j.eval()

j.parameters.R_H.val=1E20
R_val=j.parameters.R.val

j.save_model('test_jet_EC.pkl')

new_jet=Jet.load_model('test_jet_EC.pkl')

#new_jet.make_dependent_par(par='R', depends_on=['R_H','theta_open'],
# par_expr=f_par)
new_jet.eval()
print("units of theta_open",new_jet.parameters.theta_open.units)

new_jet.parameters.R_H.val=1E20

assert(new_jet.parameters.R.val==R_val)
new_jet.show_model()

class TestJetHadronic(TestBase):

def integration_suite(self,plot=False):
Expand Down

0 comments on commit d869587

Please sign in to comment.