diff --git a/jetset/base_model.py b/jetset/base_model.py index 7395db95..f10bcf05 100644 --- a/jetset/base_model.py +++ b/jetset/base_model.py @@ -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 @@ -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) @@ -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')) diff --git a/jetset/jet_paramters.py b/jetset/jet_paramters.py index ba07d47a..a43a4513 100644 --- a/jetset/jet_paramters.py +++ b/jetset/jet_paramters.py @@ -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 diff --git a/jetset/model_manager.py b/jetset/model_manager.py index e941a33e..61b22f54 100644 --- a/jetset/model_manager.py +++ b/jetset/model_manager.py @@ -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 diff --git a/jetset/model_parameters.py b/jetset/model_parameters.py index 418e5fb5..eb0c9f27 100644 --- a/jetset/model_parameters.py +++ b/jetset/model_parameters.py @@ -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 @@ -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): @@ -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: @@ -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 @@ -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': diff --git a/jetset/tests/test_depending_parameters.py b/jetset/tests/test_depending_parameters.py index c6df86c5..4beef9cf 100644 --- a/jetset/tests/test_depending_parameters.py +++ b/jetset/tests/test_depending_parameters.py @@ -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)) - \ No newline at end of file + 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() \ No newline at end of file diff --git a/jetset/tests/test_jet_model.py b/jetset/tests/test_jet_model.py index 9ef5a414..03eceeb6 100644 --- a/jetset/tests/test_jet_model.py +++ b/jetset/tests/test_jet_model.py @@ -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) @@ -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):