Skip to content

GitLab

Commit 99ad09ea authored by Matthieu Tristram's avatar Matthieu Tristram
Browse files

Correct bug

parent 2eabc76a
Hide whitespace changes
Inline Side-by-side
Showing with 53 additions and 40 deletions
work/tools/python/camel.py
View file @ 99ad09ea
......@@ -110,6 +110,9 @@ def Parnames( parfile, par_type=None):
def Interval( chain, par, CL=68, symmetrical=False):
return( getci( np.percentile( chain[par], [(100-CL)/2,50,100-(100-CL)/2]), symmetrical=symmetrical))
def Upper( chain, par, CL=95):
return( np.percentile( chain[par], CL))
def scan(filename):
f = open( filename, 'r')
......@@ -439,12 +442,15 @@ def posterior1d( chains, params, nbin=50, smooth=1,
#legend
if names != []:
# if len(params) >= 2:
if len(params) == subplot[0]*subplot[1]:
ax.legend(names, loc='upper center', frameon=False, ncol=len(names), bbox_to_anchor=(0.,1.2))
plt.subplots_adjust( top=0.85)
if len(params) ==1:
ax.legend(names, loc='upper center')
else:
ax.legend(names, bbox_to_anchor=(1.05, 0), loc='lower left', borderaxespad=0.)
plt.subplots_adjust( left=0.1, right=0.85)
if len(params) == subplot[0]*subplot[1]:
ax.legend(names, loc='upper center', frameon=False, ncol=len(names), bbox_to_anchor=(0.,1.2))
plt.subplots_adjust( top=0.85)
else:
ax.legend(names, bbox_to_anchor=(1.05, 0), loc='lower left', borderaxespad=0.)
plt.subplots_adjust( left=0.1, right=0.85)
# ax.legend(names, loc='upper right')
......@@ -594,7 +600,8 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
tmp = []
for chain in chains:
tmp = tmp+list(chain[par])
rge[par] = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4]))
rge[par] = list([min(tmp),max(tmp)])
# rge[par] = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4]))
for par in params:
......@@ -605,7 +612,7 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
chain = chains[c]
weight = weights[c]
h,X=np.histogram( chain[par], bins=nbin, range=rge[par], normed=True, weights=weight)
Y = nd.gaussian_filter( h, smooth)
Y = nd.gaussian_filter1d( h, smooth,mode='nearest')
Y = Y/np.max(Y) if Norm else Y
plt.plot( (X[1:]+X[:-1])/2., Y, color=colors[c], linestyle=linestyles[c])
......@@ -762,19 +769,17 @@ def rectangle( chains, par0, params, nbin=50, parnames=None, names=[], colors=[]
def hist1d(x, bins=50, weights=None, smooth=0., normmax=True, **kwargs):
# ax = kwargs.get("ax", plt.gca())
extents = kwargs.get("extents",[x.min(),x.max()])
n, b = np.histogram(x, weights=weights, bins=bins,normed=True)
ns = nd.gaussian_filter1d(n,smooth,mode='nearest')
if normmax:
ns /= ns.max()
xx = np.insert( (b[:-1]+b[1:])/2., 0, b[0])
yy = np.insert( ns, 0, 0.)
xx = (b[:-1]+b[1:])/2.
yy = ns
plt.plot( xx, yy, **kwargs)
plt.xlim( extents)
# plt.show()
......@@ -819,7 +824,7 @@ def ctr_level(histo2d, lvl):
alvl = np.searchsorted(cum_h, lvl)
clist = h[-alvl]
return clist
return clist[::-1]
def FeldmanCousins( xmin, sigma, CL=95):
......@@ -1392,8 +1397,11 @@ class prof:
self.dbdir = dbdir
self.parname = parname
self._nameforclass = nameforclass if nameforclass else self.parname
self._result = False
self._result = False
self._str_result = ""
self._is_fitted = False
self._offset = offset
self._CL = 68
path, dirs, files = os.walk("%s" % (self.dbdir)).next()
if sum([self.parname+"_" in d for d in dirs])==0:
......@@ -1409,25 +1417,24 @@ class prof:
if writeBestfit:
self._writeBestFit(idir,f)
# for f2 in [f1 for f1 in f if ".o" in f1]:
#check there are more than nestimate bestfit in directory idir
if sum([s.startswith("best_fit") for s in f]) >= nestimate:
# gets the value indicated in the subfolder's name
# gets the value indicated in the subfolder's name
keep=True
for ii in remove:
if ii == idir.replace(self.parname+'_',''):
keep=False
if keep:
# values in string
self._s.append( idir.replace(self.parname+'_',''))
# same thing but casts in double
# values in double + offset
self._v.append( np.double(idir.replace(self.parname+'_',''))+self._offset)
if len(self._s) <= 2:
raise ValueError( "not enough profiled-value available (need at least two)")
raise ValueError( "not enough profiled-value available:\n %s/%s" % (self.dbdir,idir))
sortidx = np.argsort(self._v)
self._s = [self._s[s] for s in sortidx]
......@@ -1609,9 +1616,9 @@ class prof:
from scipy.interpolate import interp1d,UnivariateSpline,InterpolatedUnivariateSpline
if CL:
self.CL = CL
self._CL = CL
else:
self.CL = 95 if upper else 68
self._CL = 95 if upper else 68
datav = np.array(self._v)[(self._p-min(self._p)) < deltaMax]
datap = np.array(self._p)[(self._p-min(self._p)) < deltaMax]
......@@ -1637,28 +1644,34 @@ class prof:
if upper:
xlow,xhigh = self._deltaChi2( CL=68)
res = FeldmanCousins( xmin, xhigh-xmin, self.CL)
self._result="%s < %f (%d%% C.L.)" % (parname.get(self._nameforclass,self._nameforclass),res, self.CL)
print(self._result)
self._res = FeldmanCousins( xmin, xhigh-xmin, self._CL)
self._str_res="%s < %f (%d%% C.L.)" % (parname.get(self._nameforclass,self._nameforclass),self._res, self._CL)
print(self._str_res)
else:
xlow,xhigh = self._deltaChi2( CL=self.CL)
res = (xmin,(xlow,xhigh))
self._result="$%f_{%+f}^{%+f}$" % (xmin,xlow-xmin,xhigh-xmin)
print(self._result)
xlow,xhigh = self._deltaChi2( CL=self._CL)
self._res = (xmin,(xlow,xhigh))
self._str_res="$%f_{%+f}^{%+f}$" % (xmin,xlow-xmin,xhigh-xmin)
print(self._str_res)
return( res)
self._is_fitted = True
return( self._res)
def get_fit(self):
"""
Return the fitted function
"""
return(self.f)
if self._is_fitted:
return(self.f)
else:
print( "Please fit before...")
return()
def get_result(self):
"""
Return the result from the fit in string format
Return the result from the fit
"""
return(self._result)
print( self._str_res)
return(self._res)
def _minimum( self, bounds=None):
......@@ -1725,7 +1738,7 @@ class prof:
ax = plt.gca()
from matplotlib import __version__
if __version__ >= '1.5':
colordef = next(ax._get_lines.prop_cycler)
colordef = next(ax._get_lines.prop_cycler)['color']
else:
colordef = next(ax._get_lines.color_cycle)
......@@ -1739,7 +1752,7 @@ class prof:
#plot profile
refmin = self.ymin
if self._result:
if self._is_fitted:
#use min(ftopot) to set the y shift
vartoplot=np.linspace(*extent, num=1000)
ftoplot=self.f(vartoplot)
......@@ -1752,17 +1765,17 @@ class prof:
ax.set_ylim((-0.1,5))
color = baseplt.get_color()
if draw_CL and self.CL==68:
if draw_CL and self._CL==68:
plt.axhline(y=1,color='k',linestyle='--')
#plot fit
if self._result:
if self._is_fitted:
ax.plot(vartoplot,ftoplot-min(ftoplot), color=color, linestyle=linestyle, label=label, **kwargs)
CI = (self._xlow,self._xhigh)
if draw_CL:
if np.isinf(CI[0]):
res = FeldmanCousins( self._xmin, self._xhigh-self._xmin, self.CL)
res = FeldmanCousins( self._xmin, self._xhigh-self._xmin, self._CL)
ax.plot([res,res],[-0.1,1], color=color, **kwargs)
else:
ax.plot([CI[0],CI[0]],[-0.1,1], color=color, **kwargs)
......@@ -1770,8 +1783,8 @@ class prof:
if result:
if np.isinf(CI[0]):
ax.text(max([self._xmin,0.]), 4, ' '+self._result, horizontalalignment='left' )
ax.text(max([self._xmin,0.]), 4, ' '+self._str_res, horizontalalignment='left' )
else:
ax.text(self._xmin, 4, self._result, horizontalalignment='center')
ax.text(self._xmin, 4, self._str_res, horizontalalignment='center')
########################################################################
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment