Page MenuHomeHEPForge
Files F7879480

No OneTemporary
Actions

View Options

diff --git a/bin/prof2-burst b/bin/prof2-burst
new file mode 100755
--- /dev/null
+++ b/bin/prof2-burst
@@ -0,0 +1,564 @@
+#! /usr/bin/env python
+
+"""\
+%prog rundir ipolfile datadir [opts]
+
+TODO:
+ * Update this description
+"""
+
+import optparse, os, sys
+op = optparse.OptionParser(usage=__doc__)
+op.add_option("--ierr", dest="IERR", default="symm", help="Whether to interpolate MC errors: none, mean, median, symm (default: %default)") #< add rel, asymm
+op.add_option("--pname", "--pfile", dest="PNAME", default="params.dat", help="Name of the params file to be found in each run directory (default: %default)")
+op.add_option("-n", dest="NSAMPLES", type=int, default=1, help="Number of samples")
+op.add_option("-j", dest="MULTI", type=int, default=1, help="Number of threads to use")
+op.add_option("--minos", dest="MINOS", default=False, action="store_true", help="Run Minos algorithm after minimisation")
+op.add_option("--wfile", dest="WFILE", default=None, help="Path to a weight file, used to restrict ipol building to a subset of bins (default: %default)")
+op.add_option("--limits", dest="LIMITS", default=None, help="Simple text file with parameter limits and fixed parameters")
+op.add_option("-g", "--gradient", dest="GRADIENT", action="store_true", default=False, help="Run minimisation with analytic gradient (EXPERIMENTAL!)")
+op.add_option("-s", "--strategy", dest="STRATEGY", default=1, type=int, help="Set Minuit strategy [0 fast, 1 default, 2 slow]")
+op.add_option("-r", "--result", dest="RESULT", default=None, help="Minimisation result to use for eigentunes calculation")
+op.add_option("--filter", dest="FILTER", action="store_true", default=False, help="Filter out data bins that have 0 error")
+op.add_option("--perc", dest="PERCENTILE", type=float, default=95, help="Percentile for eigentunes")
+op.add_option("-v", "--debug", dest="DEBUG", action="store_true", default=False, help="Turn on some debug messages")
+op.add_option("-q", "--quiet", dest="QUIET", action="store_true", default=False, help="turn off messages")
+op.add_option("-o", "--outfile", dest="OUTFILE", default="gofs.txt", help="Output file name for the gof measures")
+op.add_option("-p", "--prefix", dest="PREFIX", default="BURST", help="Output file name prefix")
+op.add_option("-O", "--outdir", dest="OUTDIR", default="MYBURST", help="Output directory")
+op.add_option("-S", "--subrun", dest="SUBRUN", default=0, type=int, help="Subrun in case of distributed computing")
+opts, args = op.parse_args()
+
+
+
+## Get mandatory arguments
+if len(args) < 1:
+ print "Argument missing... exiting\n\n"
+ op.print_usage()
+ sys.exit(1)
+RUNSDIR = args[0]
+IFILE = args[1]
+
+## Load the Professor machinery
+import professor2 as prof
+if not opts.QUIET:
+ print prof.logo
+
+
+IHISTOS, IMETA = prof.read_ipoldata(IFILE)
+PARAMS = prof.read_params(RUNSDIR, opts.PNAME)
+# Throw away those points not used in the original ipol
+USEDRUNS = IMETA["Runs"].split()
+from collections import OrderedDict
+USEDPARAMS = OrderedDict()
+for k, v in PARAMS.iteritems():
+ if k in USEDRUNS:
+ USEDPARAMS[k]=v
+
+## Weight file parsing to select a histos subset
+if opts.WFILE:
+ matchers = prof.read_pointmatchers(opts.WFILE)
+ for hn in IHISTOS.keys():
+ if not any(m.match_path(hn) for m in matchers.keys()):
+ del IHISTOS[hn]
+ elif opts.DEBUG:
+ print "Observable %s passed weight file path filter" % hn
+ print "Filtered observables by path, %d remaining" % len(IHISTOS)
+
+
+REFDIR = args[2]
+DHISTOS_raw = prof.read_all_histos(REFDIR)
+DHISTOS ={}
+
+# Throw away all data histos not needed
+for k in IHISTOS.keys():
+ DHISTOS[k]=DHISTOS_raw[k]
+del DHISTOS_raw
+
+matchers = prof.read_pointmatchers(opts.WFILE) if opts.WFILE else None
+
+def mkSmearedIpolSet(inputHistos, inputParams, seed=1):
+
+
+ myHISTOS ={} #{'HISTONAME': {'0000': <Histo with 0 bins>}}
+
+ # Smearing and preparation of ipol input
+ for hname, ihist in inputHistos.iteritems():
+ temp = {}
+ for run, params in inputParams.iteritems():
+ temp[run] = ihist.toDataHisto(params).mkSmearedCopy(seed)
+ myHISTOS[hname] = temp
+
+ myUSEDRUNS, myPARAMNAMES, myPARAMSLIST = prof.mk_ipolinputs(inputParams)
+ myHNAMES = sorted(inputHistos.keys())
+
+
+ valorder = inputHistos[inputHistos.keys()[0]].bins[0].ival.order
+ errorder = inputHistos[inputHistos.keys()[0]].bins[0].ierrs.order
+
+ print opts.MULTI
+ CONFIG = {"MULTI":opts.MULTI, "ORDER":valorder, "IERR":opts.IERR, "ERR_ORDER":valorder}
+ tempDict = prof.mkStandardIpols(myHISTOS, myHNAMES, myUSEDRUNS, myPARAMSLIST, CONFIG)
+
+ return tempDict, [myPARAMNAMES, myPARAMSLIST], myUSEDRUNS
+
+
+def smearDataHistos(histDict, SEED):
+ rdict = {}
+ for k, v in histDict.iteritems():
+ DH = prof.DataHisto(v.bins, v.path).mkSmearedCopy(SEED)
+ rdict[k] = DH
+ return rdict
+
+def mkFitFunc(dhistos, ipolfname, MAXERRDICT, MATCHERS, doFilter):
+ import professor2 as prof
+ ## Take parameter names directly from ifile, or fallback
+ PNAMES = METADATA["ParamNames"].split()
+ if not PNAMES:
+ PNAMES = ["A%03i" % i for i in xrange(int(METADATA["Dimension"]))]
+
+
+ ## Function definition wrapper
+ funcdef = prof.mk_fitfunc("prof.simpleGoF", PNAMES, "profGoF", ["DBINS", "IBINS", "MAXERRS"])
+ return funcdef
+
+def mkSingleTune(DHISTOS, ipolfname, MAXERRDICT, MATCHERS, doFilter):
+ import professor2 as prof
+
+ IHISTOS, METADATA = prof.read_ipoldata(ipolfname)
+ DBINS, IBINS, MAXERRS = prof.prepareBins(DHISTOS, IHISTOS, MAXERRDICT, MATCHERS, doFilter)
+
+
+
+ ## Take parameter names directly from ifile, or fallback
+ PNAMES = METADATA["ParamNames"].split()
+ if not PNAMES:
+ PNAMES = ["A%03i" % i for i in xrange(int(METADATA["Dimension"]))]
+
+
+ ## Function definition wrapper
+ funcdef = prof.mk_fitfunc("prof.simpleGoF", PNAMES, "profGoF", ["DBINS", "IBINS", "MAXERRS"])
+ exec funcdef in locals()
+ if opts.DEBUG:
+ print "Built GoF wrapper from:\n '%s'" % funcdef
+
+
+
+ try:
+ from iminuit import Minuit
+ except ImportError, e:
+ print "Unable to import iminuit, exiting", e
+ print "Try installing iminuit with pip: pip install iminuit --user"
+ import sys
+ sys.exit(1)
+
+ if not opts.QUIET:
+ print "\n"
+ print 66*"*"
+ print "* Using iminuit, please visit https://github.com/iminuit/iminuit *"
+ print 66*"*"
+ print "\n"
+
+
+ # IHISTOS, METADATA = prof.read_ipoldata(ipolfname)
+ PMINS = [float(x) for x in METADATA["MinParamVals"].split()]
+ PMAXS = [float(x) for x in METADATA["MaxParamVals"].split()]
+ FARG=prof.setupMinuitFitarg(PNAMES,PMINS,PMAXS,opts.LIMITS)
+
+
+
+ # TODO: errordef as CL params?
+ PRINTLEVEL = 0 if opts.QUIET else 1
+
+ if opts.GRADIENT:
+ graddef = prof.mk_fitfunc("simpleGoFGradient", PNAMES, "myGrad")
+ exec graddef in locals()
+ minuit = Minuit(profGoF, grad_fcn=myGrad, errordef=1, print_level=PRINTLEVEL, forced_parameters=PNAMES, **FARG)
+ else:
+ minuit = Minuit(profGoF, errordef=1, print_level=PRINTLEVEL, forced_parameters=PNAMES, **FARG)
+ minuit.strategy = opts.STRATEGY
+
+ import time
+ start_time = time.time()
+ ## Lift off
+ minuit.migrad()
+ if opts.MINOS:
+ minuit.minos()
+ if not opts.QUIET:
+ print("Minimisation finished after %s seconds" % (time.time() - start_time))
+
+ # from IPython import embed
+ # embed()
+ # import sys
+ # sys.exit(1)
+ MIN = [minuit.values[p] for p in PNAMES]
+
+ ## Goodness of fit etc
+ chi2 = minuit.fval
+
+ return chi2, MIN, minuit.covariance, len(IBINS)
+
+
+def mkGoFPlot(X, dof=1295):
+ import numpy as np
+ g_68 = np.percentile(X, 68.8)
+ g_95 = np.percentile(X, 95)
+ g_99 = np.percentile(X, 99.9)
+ import pylab
+ pylab.axvspan(min(X),g_68, label="68.8 pct", facecolor='b', alpha=0.3)
+ pylab.axvspan(g_68, g_95, label="95 pct", facecolor='r', alpha=0.3)
+ pylab.axvspan(g_95, g_99, label="99.9 pct", facecolor='g', alpha=0.3)
+ pylab.hist(X, 50, histtype="step")
+
+ # from scipy.stats import chi2
+ # x = np.linspace(chi2.ppf(0.01, dof), chi2.ppf(0.99, dof), 100)
+ # pylab.plot(x, chi2.pdf(x, dof), 'r-', lw=5, alpha=0.6, label='chi2 pdf dof=%i'%dof)
+
+ pylab.legend()
+ pylab.savefig("%s-myprofgof.pdf"%opts.PREFIX)
+
+ pylab.clf()
+ pylab.hist(X, 50, cumulative=True, histtype="step")
+ pylab.savefig("%s-myprofgofcum.pdf"%opts.PREFIX)
+
+
+def mySolve(center_t, direction_t, TRAFO, GOFdef, target):
+ exec GOFdef in globals() # Note globals!
+
+ def getVal(a):
+ temp_t = center_t + a*direction_t
+ temp = TRAFO * temp_t.transpose()
+ temp_r = temp.transpose().tolist()[0]
+ return profGoF(*temp_r) - target
+
+ def getP(a):
+ temp_t = center_t + a*direction_t
+ temp = TRAFO * temp_t.transpose()
+ return temp.transpose().tolist()[0]
+
+ from scipy.optimize import fsolve
+ x=fsolve(getVal,1)
+ print "Scaling is", x
+ return getP(x)
+
+def mkNewParamCorrelation(T_trans, T, point, GOFdef, target):# center_t, TRAFO, GOFdef, target, oldCOV, TRAFO_trans):
+ exec GOFdef in globals() # Note globals!
+
+ from numpy import matrix
+ rv = matrix(point.values())
+ center_t = (T_trans * rv.transpose()).transpose()
+
+
+ DIM=len(point.values())
+ from scipy.optimize import fsolve
+ from numpy import zeros,diag
+
+ def getVal(a, direction):
+ temp_t = center_t + a*direction
+ temp = T * temp_t.transpose()
+ temp_r = temp.transpose().tolist()[0]
+ return profGoF(*temp_r) - target
+
+ def getX(a, direction):
+ return center_t + a*direction
+
+ newdiag=[]
+ for i in xrange(DIM):
+ ev = zeros(DIM)
+ ev[i] = 1
+ temp = fsolve(lambda x:getVal(x, ev), 1)
+ temp2 = fsolve(lambda x:getVal(x, -ev), 1)
+ if temp > temp2:
+ newdiag.append(temp)
+ else:
+ newdiag.append(temp2)
+
+ N = diag([x[0] for x in newdiag])
+ return T_trans*N*T
+
+
+
+
+
+
+def mkEigentunes(T_trans, T, point, GOFdef, target, plus=True):
+ """
+ COV ... real symmetric covariance matrix
+ point ... could be any point in the true parameter space but should be
+ the minimisation result i.e. the center of COV
+ """
+ from numpy import sqrt, zeros, matrix
+ # Trsnform the minimisation result into the other coordinate system
+ rv = matrix(point.values())
+ rv_trans = (T_trans * rv.transpose()).transpose()
+
+ ret = []
+
+
+ # Construct all base vectors (in rotated system) with pos and neg directions
+ dim = len(point.values())
+ EVS = []
+ for i in xrange(dim):
+ ev = zeros(dim) # A zero vector in len(S) dimensions
+ # Set one of the coordinates to 1 or -1
+ ev[i] = 1 if plus else -1
+ EVS.append(ev)
+ print ev
+
+ # Get the eigentunes
+ for num, ev in enumerate(EVS):
+ thisEigentune = mySolve(rv_trans, ev, T, GOFdef, target)
+ ret.append([int(ev[num]*(num+1)), thisEigentune])
+
+ return ret
+
+def calcHistoCov(h, COV_P, result):
+ """
+ Propagate the parameter covariance onto the histogram covariance
+ using the ipol gradients.
+ """
+ IBINS=h.bins
+ from numpy import zeros
+ COV_H = zeros((h.nbins, h.nbins))
+ from numpy import array
+ for i in xrange(len(IBINS)):
+ GRD_i = array(IBINS[i].grad(result))
+ for j in xrange(len(IBINS)):
+ GRD_j = array(IBINS[j].grad(result))
+ pc =GRD_i.dot(COV_P).dot(GRD_j)
+ COV_H[i][j] = pc
+ return COV_H
+
+def mkErrorPropagationHistos(IHISTOS, point, COV, combine=False):
+ covs={}
+ properrs = {}
+ ipolerrs = {}
+ ipolvals = {}
+ from numpy import sqrt
+ for k, v in IHISTOS.iteritems():
+ covs[k] = calcHistoCov(v, COV, point)
+ properrs[k] = sqrt(covs[k].diagonal())
+ ipolerrs[k] = [b.err for b in v.toDataHisto(point).bins]
+ ipolvals[k] = [b.val for b in v.toDataHisto(point).bins]
+
+ scatters=[]
+ for k, v in IHISTOS.iteritems():
+ T=v.toDataHisto(point)
+ for i in xrange(T.nbins):
+ T.bins[i].errs=properrs[k][i] if combine is False else sqrt(properrs[k][i]**2 + ipolerrs[k][i]**2)
+ scatters.append(T.toScatter2D())
+ return scatters
+
+
+def mkScatters(ipolH, ppoint):
+ scatters_e =[]
+ for k in sorted(ipolH.keys()):
+ v = ipolH[k]
+ T=v.toDataHisto(ppoint)
+ scatters_e.append(T.toScatter2D())
+ return scatters_e
+
+def covarianceToList(cov, pnames):
+ """
+ Reformat minuit covariance matrix dict to single line, e.g. 3 dimensions
+ 0 1 2
+ 1 3 4 --> [0 1 2 3 4 5]
+ 2 4 5
+ """
+ covline = []
+ for i, px in enumerate(pnames):
+ for j, py in enumerate(pnames):
+ if i<=j:
+ covline.append(cov[(px, py)])
+ return covline
+
+
+def mkEnvelopes(central, etunes):
+ ret = {}
+ for i in xrange(1, len(etunes.keys())-2):
+ ret[i] = []
+ Hplus = etunes[i]
+ Hminus = etunes[-i]
+ for num_h, h in enumerate(central):
+ temp = h.clone()
+ for num_p, p in enumerate(temp.points):
+ yplus = Hplus[num_h].points[num_p].y
+ yminus = Hminus[num_h].points[num_p].y
+ if yplus > p.y:
+ eplus = yplus - p.y
+ eminus = p.y - yminus
+ else:
+ eplus = yminus - p.y
+ eminus = p.y - yplus
+ p.yErrs = (eminus, eplus)
+ ret[i].append(temp)
+ return ret
+
+def writeToFile(items, nbins, params, outdir, fname, suffix):
+ """
+ Write out gof values to file as well as some meta info in the first lines
+ """
+ head = "# Nbins %i\n# Parameters"%nbins
+ for p in params:
+ head += " %s"%p
+ head += "\n"
+
+ import os
+ if not os.path.exists(outdir):
+ os.makedirs(outdir)
+ fname += "-%i.txt"%suffix
+ outname = os.path.join(outdir, fname)
+ with open(outname, "w") as f:
+ f.write(head)
+
+ for g in items:
+ if type(g)==list:
+ temp =""
+ for i in g:
+ temp+=" %e"%i
+ f.write(temp.strip()+"\n")
+ else:
+ f.write("%e\n"% g)
+
+def writeMINS(mins, outdir, params, suffix=None):
+ """
+ Write out gof values to file as well as some meta info in the first lines
+ """
+ head = "# Parameters"
+ for p in params:
+ head += " %s"%p
+ head += "\n"
+
+ import os
+ if not os.path.exists(outdir):
+ os.makedirs(outdir)
+ fname = "results.txt" if suffix is None else "results-%s.txt"%suffix
+ outname = os.path.join(outdir, fname)
+ with open(outname, "w") as f:
+ f.write(head)
+
+ for m in mins:
+ temp =""
+ for i in m:
+ temp+=" %e"%i
+ f.write(temp.strip()+"\n")
+
+def writeCOVS(covs, outdir, params, suffix=None):
+ """
+ Write out gof values to file as well as some meta info in the first lines
+ """
+ head = "# Parameters"
+ for p in params:
+ head += " %s"%p
+ head += "\n"
+
+ import os
+ if not os.path.exists(outdir):
+ os.makedirs(outdir)
+ fname = "covariances.txt" if suffix is None else "covariances-%s.txt"%suffix
+ outname = os.path.join(outdir, fname)
+ with open(outname, "w") as f:
+ f.write(head)
+
+ for c in covs:
+ temp =""
+ for i in c:
+ temp+=" %e"%i
+ f.write(temp.strip()+"\n")
+GOFS=[]
+MINS=[]
+COVS=[]
+if len(args)==4:
+ with open(args[3]) as f:
+ for l in f:
+ GOFS.append(float(l.strip()))
+
+else:
+ NBINS=0
+ PNAMES=[]
+ for i in xrange(opts.NSAMPLES):
+ thisSEED=opts.NSAMPLES*(1 + opts.SUBRUN) + i + 1
+ IDICT, PARAMS, RUNS = mkSmearedIpolSet(IHISTOS, USEDPARAMS, thisSEED)
+ PNAMES=PARAMS[0]
+ # Write the ipol into a temp file
+ temp = prof.writeIpol("temp", IDICT, PARAMS, RUNS, "", "")
+ del IDICT
+ # Also smear the input data
+ DHISTOS_smeared = smearDataHistos(DHISTOS, thisSEED)
+ gof, minimum, covariance, nbins = mkSingleTune(DHISTOS_smeared, temp.name, None, matchers, opts.FILTER)
+ GOFS.append(gof)
+ MINS.append(minimum)
+ COVS.append(covarianceToList(covariance, PNAMES))
+ NBINS=nbins
+
+ # Delete temp file
+ import os
+ os.remove(temp.name)
+
+ writeToFile(GOFS, NBINS, PNAMES, opts.OUTDIR,"gof", opts.SUBRUN)
+ writeToFile(MINS, NBINS, PNAMES, opts.OUTDIR,"results", opts.SUBRUN)
+ writeToFile(COVS, NBINS, PNAMES, opts.OUTDIR,"covariances", opts.SUBRUN)
+ # writeMINS(MINS, opts.OUTDIR, PNAMES, opts.SUBRUN)
+ # writeCOVS(COVS, opts.OUTDIR, PNAMES, opts.SUBRUN)
+
+ # from IPython import embed
+ # embed()
+ import sys
+ sys.exit(1)
+ # with open("%s-%s"%(opts.PREFIX,opts.OUTFILE), "w") as f:
+ # for g in GOFS:
+ # f.write("%f\n"% g)
+
+
+
+mkGoFPlot(GOFS)
+
+IHISTOS, METADATA = prof.read_ipoldata(IFILE)
+PNAMES = METADATA["ParamNames"].split()
+DBINS, IBINS, MAXERRS = prof.prepareBins(DHISTOS, IHISTOS, None, matchers, opts.FILTER)
+F = prof.mk_fitfunc("prof.simpleGoF", PNAMES, "profGoF", ["DBINS", "IBINS", "MAXERRS"])
+#F=mkFitFunc(DHISTOS, IFILE,None, matchers, opts.FILTER)
+
+import yoda
+if opts.RESULT is not None:
+ # Read the tuning result
+ P_min, OTH = prof.readResult(opts.RESULT)
+ C_param = prof.getParamCov(OTH) # That's the minimiser covariance
+ import professor2 as prof
+ # S, T_fwd are the return values if linalg.eig(COV)
+ # with S being the eigenvalues and T_fwd being composed of
+ # the eigenvectors
+ T_fwd, S, T_back = prof.eigenDecomposition(C_param)
+
+ # Get the gof target according to the required percentile
+ import numpy as np
+ gof_target = np.percentile(GOFS, opts.PERCENTILE)
+
+ # Calculate the eigen tunes (points in parameter space)
+ E_plus = mkEigentunes(T_fwd, T_back, P_min, F, gof_target)
+ E_minus = mkEigentunes(T_fwd, T_back, P_min, F, gof_target, plus=False)
+ ETs = dict(E_plus+E_minus)
+
+ # Get the corresponding ipol histos
+ EThists={}
+ for k, ET in ETs.iteritems():
+ thisEThists = mkScatters(IHISTOS, ET)
+ sgn = "+" if k > 0 else "-"
+ yoda.writeYODA(thisEThists, "%s-Eigentunes_%.1f_%i%s.yoda"%(opts.PREFIX,opts.PERCENTILE, int(abs(k)), sgn))
+ EThists[k]=thisEThists
+
+ # And for convenience corresponding envelopes
+ H_min = mkScatters(IHISTOS, P_min)
+ envelopes = mkEnvelopes(H_min, EThists)
+ for k, v in envelopes.iteritems():
+ yoda.writeYODA(v, "%s-EigentunesComb_%.1f_%i.yoda"%(opts.PREFIX,opts.PERCENTILE, k))
+
+
+ # This is a (conservative) symmetric real parameter covariance matrix reflecting the eigentunes
+ C_param_new = mkNewParamCorrelation(T_fwd, T_back, P_min, F, gof_target)
+
+ # This is now error propagation
+ yoda.writeYODA(mkErrorPropagationHistos(IHISTOS, P_min,C_param_new, True), "%s-ih_comberr_eigcov_%.1f.yoda"%(opts.PREFIX,opts.PERCENTILE))
+ yoda.writeYODA(mkErrorPropagationHistos(IHISTOS, P_min,C_param , True), "%s-ih_comberr_mincov_%.1f.yoda"%(opts.PREFIX,opts.PERCENTILE))
+ # from IPython import embed
+ # embed()
diff --git a/bin/prof2-errors b/bin/prof2-errors
new file mode 100755
--- /dev/null
+++ b/bin/prof2-errors
@@ -0,0 +1,410 @@
+#! /usr/bin/env python
+
+"""\
+%prog <ipolfile>=ipol.dat [opts]
+
+Interpolate histo bin values as a function of the parameter space by loading
+run data and parameter lists from run directories in $runsdir (often "mc")
+
+TODO:
+ * Use weight file position matches to exclude some bins, as well as path matching
+ * Handle run combination file/string (write a hash of the run list into the ipol filename?)
+ * Support asymm error parameterisation
+"""
+
+import optparse, os, sys
+op = optparse.OptionParser(usage=__doc__)
+op.add_option("--ierr", dest="IERR", default="symm", help="Whether to interpolate MC errors: none, mean, median, symm (default: %default)") #< add rel, asymm
+op.add_option("-n", dest="NSAMPLES", type=int, default=1, help="Number of samples")
+op.add_option("-j", dest="MULTI", type=int, default=1, help="Number of threads to use")
+op.add_option("--minos", dest="MINOS", default=False, action="store_true", help="Run Minos algorithm after minimisation")
+op.add_option("--wfile", dest="WFILE", default=None, help="Path to a weight file, used to restrict ipol building to a subset of bins (default: %default)")
+op.add_option("--limits", dest="LIMITS", default=None, help="Simple text file with parameter limits and fixed parameters")
+op.add_option("-g", "--gradient", dest="GRADIENT", action="store_true", default=False, help="Run minimisation with analytic gradient (EXPERIMENTAL!)")
+op.add_option("-s", "--strategy", dest="STRATEGY", default=1, type=int, help="Set Minuit strategy [0 fast, 1 default, 2 slow]")
+op.add_option("-r", "--result", dest="RESULT", default=None, help="Minimisation result to use for eigentunes calculation")
+op.add_option("--filter", dest="FILTER", action="store_true", default=False, help="Filter out data bins that have 0 error")
+op.add_option("--perc", dest="PERCENTILE", type=float, default=95, help="Percentile for eigentunes")
+op.add_option("-v", "--debug", dest="DEBUG", action="store_true", default=False, help="Turn on some debug messages")
+op.add_option("-q", "--quiet", dest="QUIET", action="store_true", default=False, help="turn off messages")
+op.add_option("-o", "--outfile", dest="OUTFILE", default="gofs.txt", help="Output file name for the gof measures")
+op.add_option("-p", "--prefix", dest="PREFIX", default="BURST", help="Output file name prefix")
+op.add_option("-O", "--outdir", dest="OUTDIR", default=None, help="Output directory")
+op.add_option("-S", "--subrun", dest="SUBRUN", default=0, type=int, help="Subrun in case of distributed computing")
+opts, args = op.parse_args()
+
+
+
+## Get mandatory arguments
+if len(args) < 1:
+ print "Argument missing... exiting\n\n"
+ op.print_usage()
+ sys.exit(1)
+INDIR = args[0]
+IFILE = args[1]
+
+## Load the Professor machinery
+import professor2 as prof
+if not opts.QUIET:
+ print prof.logo
+
+
+IHISTOS, IMETA = prof.read_ipoldata(IFILE)
+# PARAMS = prof.read_params(RUNSDIR)
+# # Throw away those points not used in the original ipol
+# USEDRUNS = IMETA["Runs"].split()
+# from collections import OrderedDict
+# USEDPARAMS = OrderedDict()
+# for k, v in PARAMS.iteritems():
+ # if k in USEDRUNS:
+ # USEDPARAMS[k]=v
+
+## Weight file parsing to select a histos subset
+if opts.WFILE:
+ matchers = prof.read_pointmatchers(opts.WFILE)
+ for hn in IHISTOS.keys():
+ if not any(m.match_path(hn) for m in matchers.keys()):
+ del IHISTOS[hn]
+ elif opts.DEBUG:
+ print "Observable %s passed weight file path filter" % hn
+ print "Filtered observables by path, %d remaining" % len(IHISTOS)
+
+
+REFDIR = args[2]
+DHISTOS_raw = prof.read_all_histos(REFDIR)
+DHISTOS ={}
+
+# Throw away all data histos not needed
+for k in IHISTOS.keys():
+ DHISTOS[k]=DHISTOS_raw[k]
+del DHISTOS_raw
+
+matchers = prof.read_pointmatchers(opts.WFILE) if opts.WFILE else None
+
+
+
+
+
+def mkGoFPlot(X, dof, OUTDIR):
+ import numpy as np
+ g_68 = np.percentile(X, 68.8)
+ g_95 = np.percentile(X, 95)
+ g_99 = np.percentile(X, 99.9)
+ import pylab
+ pylab.axvspan(min(X),g_68, label="68.8 pct", facecolor='b', alpha=0.3)
+ pylab.axvspan(g_68, g_95, label="95 pct", facecolor='r', alpha=0.3)
+ pylab.axvspan(g_95, g_99, label="99.9 pct", facecolor='g', alpha=0.3)
+ pylab.hist(X, 50, histtype="step")
+
+ # from scipy.stats import chi2
+ # x = np.linspace(chi2.ppf(0.01, dof), chi2.ppf(0.99, dof), 100)
+ # pylab.plot(x, chi2.pdf(x, dof), 'r-', lw=5, alpha=0.6, label='chi2 pdf dof=%i'%dof)
+
+ pylab.legend()
+ pylab.savefig("%s/myprofgof.pdf"%OUTDIR)
+
+ pylab.clf()
+ pylab.hist(X, 50, cumulative=True, histtype="step")
+ pylab.savefig("%s/myprofgofcumulative.pdf"%OUTDIR)
+
+def mkResDistPlot(X, pname, OUTDIR):
+ import pylab
+ pylab.clf()
+ pylab.hist(X, 50, histtype="step")
+ pylab.xlabel(pname)
+ pylab.savefig("%s/distr_%s.pdf"%(OUTDIR, pname))
+
+
+def mySolve(center_t, direction_t, TRAFO, GOFdef, target):
+ exec GOFdef in globals() # Note globals!
+
+ def getVal(a):
+ temp_t = center_t + a*direction_t
+ temp = TRAFO * temp_t.transpose()
+ temp_r = temp.transpose().tolist()[0]
+ return profGoF(*temp_r) - target
+
+ def getP(a):
+ temp_t = center_t + a*direction_t
+ temp = TRAFO * temp_t.transpose()
+ return temp.transpose().tolist()[0]
+
+ from scipy.optimize import fsolve
+ x=fsolve(getVal,1)
+ return getP(x)
+
+def mkNewParamCorrelation(T_trans, T, point, GOFdef, target):# center_t, TRAFO, GOFdef, target, oldCOV, TRAFO_trans):
+ exec GOFdef in globals() # Note globals!
+
+ from numpy import matrix
+ rv = matrix(point.values())
+ center_t = (T_trans * rv.transpose()).transpose()
+
+
+ DIM=len(point.values())
+ from scipy.optimize import fsolve
+ from numpy import zeros,diag
+
+ def getVal(a, direction):
+ temp_t = center_t + a*direction
+ temp = T * temp_t.transpose()
+ temp_r = temp.transpose().tolist()[0]
+ return profGoF(*temp_r) - target
+
+ def getX(a, direction):
+ return center_t + a*direction
+
+ newdiag=[]
+ for i in xrange(DIM):
+ ev = zeros(DIM)
+ ev[i] = 1
+ temp = fsolve(lambda x:getVal(x, ev), 1)
+ temp2 = fsolve(lambda x:getVal(x, -ev), 1)
+ if temp > temp2:
+ newdiag.append(temp)
+ else:
+ newdiag.append(temp2)
+
+ N = diag([x[0] for x in newdiag])
+ return T_trans*N*T
+
+
+def mkEigentunes(T_trans, T, point, GOFdef, target, plus=True):
+ """
+ COV ... real symmetric covariance matrix
+ point ... could be any point in the true parameter space but should be
+ the minimisation result i.e. the center of COV
+ """
+ from numpy import sqrt, zeros, matrix
+ # Trsnform the minimisation result into the other coordinate system
+ rv = matrix(point.values())
+ rv_trans = (T_trans * rv.transpose()).transpose()
+
+ ret = []
+
+
+ # Construct all base vectors (in rotated system) with pos and neg directions
+ dim = len(point.values())
+ EVS = []
+ for i in xrange(dim):
+ ev = zeros(dim) # A zero vector in len(S) dimensions
+ # Set one of the coordinates to 1 or -1
+ ev[i] = 1 if plus else -1
+ EVS.append(ev)
+
+ # Get the eigentunes
+ for num, ev in enumerate(EVS):
+ thisEigentune = mySolve(rv_trans, ev, T, GOFdef, target)
+ ret.append([int(ev[num]*(num+1)), thisEigentune])
+
+ return ret
+
+def calcHistoCov(h, COV_P, result):
+ """
+ Propagate the parameter covariance onto the histogram covariance
+ using the ipol gradients.
+ """
+ IBINS=h.bins
+ from numpy import zeros
+ COV_H = zeros((h.nbins, h.nbins))
+ from numpy import array
+ for i in xrange(len(IBINS)):
+ GRD_i = array(IBINS[i].grad(result))
+ for j in xrange(len(IBINS)):
+ GRD_j = array(IBINS[j].grad(result))
+ pc =GRD_i.dot(COV_P).dot(GRD_j)
+ COV_H[i][j] = pc
+ return COV_H
+
+def mkErrorPropagationHistos(IHISTOS, point, COV, combine=False):
+ covs={}
+ properrs = {}
+ ipolerrs = {}
+ ipolvals = {}
+ from numpy import sqrt
+ for k, v in IHISTOS.iteritems():
+ covs[k] = calcHistoCov(v, COV, point)
+ properrs[k] = sqrt(0.5*covs[k].diagonal())
+ ipolerrs[k] = [b.err for b in v.toDataHisto(point).bins]
+ ipolvals[k] = [b.val for b in v.toDataHisto(point).bins]
+
+ scatters=[]
+ for k, v in IHISTOS.iteritems():
+ T=v.toDataHisto(point)
+ for i in xrange(T.nbins):
+ T.bins[i].errs=properrs[k][i] if combine is False else sqrt(properrs[k][i]**2 + ipolerrs[k][i]**2)
+ scatters.append(T.toScatter2D())
+ return scatters
+
+
+def mkScatters(ipolH, ppoint):
+ scatters_e =[]
+ for k in sorted(ipolH.keys()):
+ v = ipolH[k]
+ T=v.toDataHisto(ppoint)
+ scatters_e.append(T.toScatter2D())
+ return scatters_e
+
+
+
+def mkEnvelopes(central, etunes):
+ ret = {}
+ for i in xrange(1, len(etunes.keys())-2): # TODO minus 3???
+ ret[i] = []
+ Hplus = etunes[i]
+ Hminus = etunes[-i]
+ for num_h, h in enumerate(central):
+ temp = h.clone()
+ for num_p, p in enumerate(temp.points):
+ yplus = Hplus[num_h].points[num_p].y
+ yminus = Hminus[num_h].points[num_p].y
+ if yplus > p.y:
+ eplus = yplus - p.y
+ eminus = p.y - yminus
+ else:
+ eplus = yminus - p.y
+ eminus = p.y - yplus
+ p.yErrs = (eminus, eplus)
+ ret[i].append(temp)
+ return ret
+
+def mkTotvelopes(central, etunes):
+ ret = []
+ for num_h, h in enumerate(central):
+ temp = h.clone()
+ allThis = [x[num_h] for x in etunes.values()]
+ for num_p, p in enumerate(temp.points):
+ dybin = [x.points[num_p].y - p.y for x in allThis]
+ pos = [x for x in dybin if x>=0]
+ neg = [x for x in dybin if x<0]
+ eplus = max(pos) if len(pos) > 0 else 0
+ eminus = abs(min(neg)) if len(neg) > 0 else 0
+ p.yErrs = (eminus, eplus)
+ ret.append(temp)
+ return ret
+
+def mkAddvelopes(central, etunes, addLinear=False):
+ ret = []
+ for num_h, h in enumerate(central):
+ temp = h.clone()
+ allThis = [x[num_h] for x in etunes.values()]
+ for num_p, p in enumerate(temp.points):
+ dybin = [x.points[num_p].y-p.y for x in allThis]
+ pos = [x for x in dybin if x>=0]
+ neg = [x for x in dybin if x<0]
+ from math import sqrt
+ if addLinear:
+ eplus = sum(pos) if len(pos) > 0 else 0
+ eminus = sum([abs(x) for x in neg]) if len(neg) > 0 else 0
+ else:
+ eplus = sqrt(sum([x*x for x in pos])) if len(pos) > 0 else 0
+ eminus = sqrt(sum([x*x for x in neg])) if len(neg) > 0 else 0
+ p.yErrs = (eminus, eplus)
+ ret.append(temp)
+ return ret
+
+
+def readFromFile(fname):
+ nbins=0
+ params=[]
+ ret = []
+ with open(fname) as f:
+ for line in f:
+ l=line.strip()
+ if l.startswith("#"):
+ if "Nbins" in l:
+ nbins = int(l.split()[-1])
+ elif "Parameters" in l:
+ params = l.split()[2:]
+ else:
+ ret.append(map(float, l.split()))
+ return ret, nbins, params
+
+
+import glob
+g_files = glob.glob("%s/*gof*.txt"%INDIR)
+c_files = glob.glob("%s/*covariance*.txt"%INDIR)
+m_files = glob.glob("%s/*results*.txt"%INDIR)
+
+
+import numpy as np
+NBINS, PNAMES=readFromFile(g_files[0])[1:]
+GOFS=np.array([readFromFile(x)[0] for x in g_files]).flatten()
+MINS=np.array([readFromFile(x)[0] for x in m_files]).reshape((len(GOFS), len(PNAMES)))
+
+
+
+
+OUTDIR=INDIR if opts.OUTDIR is None else opts.OUTDIR
+if not os.path.exists(OUTDIR):
+ os.makedirs(OUTDIR)
+mkGoFPlot(GOFS, NBINS-len(PNAMES), OUTDIR)
+
+# Distribution histograms for the minimisation result parameters
+for num, p in enumerate(PNAMES):
+ mkResDistPlot(MINS[:,num], p, OUTDIR)
+
+
+IHISTOS, METADATA = prof.read_ipoldata(IFILE)
+DBINS, IBINS, MAXERRS = prof.prepareBins(DHISTOS, IHISTOS, None, matchers, opts.FILTER)
+# F=mkFitFunc(DHISTOS, IFILE,None, matchers, opts.FILTER)
+F = prof.mk_fitfunc("prof.simpleGoF", PNAMES, "profGoF", ["DBINS", "IBINS", "MAXERRS"])
+
+import yoda
+if opts.RESULT is not None:
+ # Read the tuning result
+ P_min, OTH = prof.readResult(opts.RESULT)
+ C_param = prof.getParamCov(OTH) # That's the minimiser covariance
+ import professor2 as prof
+ # S, T_fwd are the return values if linalg.eig(COV)
+ # with S being the eigenvalues and T_fwd being composed of
+ # the eigenvectors
+ T_fwd, S, T_back = prof.eigenDecomposition(C_param)
+
+ # Get the gof target according to the required percentile
+ import numpy as np
+ gof_target = np.percentile(GOFS, opts.PERCENTILE)
+
+ # Calculate the eigen tunes (points in parameter space)
+ E_plus = mkEigentunes(T_fwd, T_back, P_min, F, gof_target)
+ E_minus = mkEigentunes(T_fwd, T_back, P_min, F, gof_target, plus=False)
+ ETs = dict(E_plus+E_minus)
+
+ for k, v in ETs.iteritems():
+ print "Eigentune", k
+ print v
+
+ # Get the corresponding ipol histos
+ EThists={}
+ for k, ET in ETs.iteritems():
+ thisEThists = mkScatters(IHISTOS, ET)
+ sgn = "+" if k > 0 else "-"
+ yoda.writeYODA(thisEThists, "%s/Eigentunes_%.1f_%i%s.yoda"%(OUTDIR,opts.PERCENTILE, int(abs(k)), sgn))
+ EThists[k]=thisEThists
+
+ # And for convenience corresponding envelopes
+ H_min = mkScatters(IHISTOS, P_min)
+ envelopes = mkEnvelopes(H_min, EThists)
+ for k, v in envelopes.iteritems():
+ yoda.writeYODA(v, "%s/EigentunesComb_%.1f_%i.yoda"%(OUTDIR,opts.PERCENTILE, k))
+
+ # This is the envelope of the eigentunes
+ totvelopes = mkTotvelopes(H_min, EThists)
+ # This is the deltas added in quadrature
+ quadvelopes = mkAddvelopes(H_min, EThists)
+ # This is the deltas added in lineary
+ linvelopes = mkAddvelopes(H_min, EThists, addLinear=True)
+
+ yoda.writeYODA(totvelopes, "%s/Totvelopes_%.1f.yoda"%(OUTDIR,opts.PERCENTILE))
+ yoda.writeYODA(quadvelopes, "%s/Quadvelopes_%.1f.yoda"%(OUTDIR,opts.PERCENTILE))
+ yoda.writeYODA(linvelopes, "%s/Linvelopes_%.1f.yoda"%(OUTDIR,opts.PERCENTILE))
+
+
+ # This is a (conservative) symmetric real parameter covariance matrix reflecting the eigentunes
+ C_param_new = mkNewParamCorrelation(T_fwd, T_back, P_min, F, gof_target)
+
+ # This is now error propagation
+ yoda.writeYODA(mkErrorPropagationHistos(IHISTOS, P_min,C_param_new, combine=False), "%s/ih_comberr_eigcov_%.1f.yoda"%(OUTDIR,opts.PERCENTILE))
+ yoda.writeYODA(mkErrorPropagationHistos(IHISTOS, P_min,C_param , combine=False), "%s/ih_comberr_mincov_%.1f.yoda"%(OUTDIR,opts.PERCENTILE))
+ # from IPython import embed
+ # embed()
diff --git a/bin/prof2-jackknife b/bin/prof2-jackknife
new file mode 100755
--- /dev/null
+++ b/bin/prof2-jackknife
@@ -0,0 +1,362 @@
+#! /usr/bin/env python
+
+"""\
+%prog <runsdir> [<ipolfile>=ipol.dat] [opts]
+
+Compute and plot the relative deviations of the given interpolation from the
+corresponding sample points in the given <runsdir>, for both bin values and bin
+errors. Histograms are made for each histogram independently (summed over bins and MC
+points), and for the sum over all histos and bins, and are written to PDFs and
+to a res.yoda data file.
+
+This can be useful for determining whether the interpolation is providing a
+sufficiently faithful approximation to the true observables, and choosing the
+most appropriate order of fit polynomial. It can also be useful to detect
+overfitting, by comparing residuals distributions between the set of runs used
+to make the fit (via prof-ipol) and an equivalent set not included in the fit.
+
+TODO:
+ * Allow control over the output filenames / directory.
+ * Add an option for value-only (no err) residuals
+ * Support runcomb filtering
+"""
+
+
+from __future__ import print_function
+
+import matplotlib, os
+matplotlib.use(os.environ.get("MPL_BACKEND", "Agg"))
+
+import optparse, os, sys
+op = optparse.OptionParser(usage=__doc__)
+#op.add_option("--ifile", dest="IFILE", default="ipol.dat", help="file from which to read the bin interpolations (default: %default)")
+op.add_option("--pname", "--pfile", dest="PNAME", default="params.dat", help="name of the params file to be found in each run directory (default: %default)")
+op.add_option("-j", dest="MULTI", type=int, default=1, help="Number of threads to use")
+op.add_option("--summ", dest="SUMMARY", default=None, help="Summary description to be written to the ipol output file")
+op.add_option("--wfile", dest="WFILE", default=None, help="Path to a weight file, used to restrict ipol building to a subset of bins (default: %default)")
+op.add_option("--limits", dest="LIMITS", default=None, help="Simple text file with parameter limits and fixed parameters")
+op.add_option("--no-plots", dest="NO_PLOTS", action="store_true", default=False, help="don't write histogram PDFs (default: %default)")
+op.add_option("--tot-only", dest="ONLY_TOT", action="store_true", default=False, help="only make total residuals histograms, not per-observable (default: %default)")
+op.add_option("--logx", dest="LOG_BINS", action="store_true", default=False, help="use a symmetric log binning for better resolution around 0 (default: %default)")
+op.add_option("--log", dest="LOG", action="store_true", default=False, help="Compare with log values")
+op.add_option("-o", dest="OUTDIR", default=".", help="Plot output folder")
+op.add_option("-p", dest="PREFIX", default="jackknife", help="Prefix for Ipol file storage")
+op.add_option("-d", dest="DROPOUT", type=int,default=10, help="Number of drop-out points")
+op.add_option("-n", dest="NITERATIONS", type=int, default=1, help="Number of iterations")
+op.add_option("-m", dest="MINORDER", type=int, default=1, help="Lowest order of polynomials to consider")
+op.add_option("-M", dest="MAXORDER", default=100, type=int, help="Maximal order to consider")
+op.add_option("-v", "--debug", dest="DEBUG", action="store_true", default=False, help="turn on some debug messages")
+op.add_option("-q", "--quiet", dest="QUIET", action="store_true", default=False, help="turn off messages")
+opts, args = op.parse_args()
+
+## Get mandatory arguments
+if len(args) < 1:
+ print("Argument missing... exiting\n\n", file=sys.stderr)
+ op.print_usage()
+ sys.exit(1)
+RUNSDIR = args[0]
+
+
+## Load the Professor machinery
+import professor2 as prof
+if not opts.QUIET:
+ print(prof.logo)
+
+## No point in going further if YODA isn't available, too
+try:
+ import yoda
+except ImportError:
+ print("YODA is required by this tool... exiting", file=sys.stderr)
+ sys.exit(1)
+
+
+
+## Load MC run histos and params
+# TODO: add runcomb file parsing to select a runs subset
+try:
+ PARAMS, HISTOS = prof.read_all_rundata(RUNSDIR, opts.PNAME)
+ RUNS, PARAMNAMES, PARAMSLIST = prof.mk_ipolinputs(PARAMS)
+except Exception, e:
+ print(e, file=sys.stderr)
+ sys.exit(1)
+
+
+## Weight file parsing to select a histos subset
+if opts.WFILE:
+ matchers = prof.read_pointmatchers(opts.WFILE)
+ for hn in HISTOS.keys():
+ if not any(m.match_path(hn) for m in matchers.keys()):
+ del HISTOS[hn]
+ elif opts.DEBUG:
+ print ("Observable %s passed weight file path filter" % hn)
+ print ("Filtered observables by path, %d remaining" % len(HISTOS))
+HNAMES = HISTOS.keys()
+
+## If there's nothing left to interpolate, exit!
+if not HNAMES:
+ print ("No observables remaining... exiting")
+ sys.exit(1)
+
+## Robustness tests and cleaning: only retain runs that contain every histo
+bad, badnum = [], []
+for irun, run in enumerate(RUNS):
+ for hn in HNAMES:
+ if not HISTOS[hn].has_key(run):
+ bad.append(run)
+ badnum.append(irun)
+ break
+if opts.LIMITS != None:
+ limits, fixed = prof.read_limitsandfixed(opts.LIMITS)
+ for irun, run in enumerate(RUNS):
+ for inam, nam in enumerate(PARAMNAMES):
+ if PARAMSLIST[irun][inam] <= limits[nam][0] or PARAMSLIST[irun][inam] >= limits[nam][1]:
+ if not run in bad:
+ bad.append(run)
+ badnum.append(irun)
+ break
+if bad:
+ print ("Found %d bad runs in %d total... removing" % (len(bad), len(RUNS)))
+ goodr, goodp = [], []
+ for irun, run in enumerate(RUNS):
+ if not irun in badnum:
+ goodr.append(run)
+ goodp.append(PARAMSLIST[irun])
+ RUNS = goodr
+ PARAMSLIST = goodp
+
+## If there's nothing left to interpolate, exit!
+if not RUNS:
+ print ("No valid runs remaining... exiting")
+ sys.exit(1)
+
+# # Testing
+# RUNS=[1 for i in xrange(500)]
+# PARAMNAMES=[i for i in xrange(2)]
+
+# Dimension of the parameter space
+DIM=len(PARAMNAMES)
+
+# The jackknifed number of runs --- determine all possible orders
+JACKNUM = len(RUNS) - opts.DROPOUT
+ORDERS=[]
+
+o_temp=opts.MINORDER
+while True:
+ n_temp = prof.numCoeffs(DIM, o_temp)
+ if n_temp > JACKNUM or o_temp > opts.MAXORDER:
+ break
+ ORDERS.append(o_temp)
+ o_temp+=1
+
+def mkIpolFilename(iteration, order):
+ iname=opts.PREFIX+"_I_%i_O_%i.dat"%(iteration, order)
+ return iname
+
+
+
+def mkIpols(RUNS, PARAMSLIST, ORDER, IFILE):
+
+ IHISTOS = {}
+
+ import zlib
+
+ def worker(q, rdict):
+ "Function to make bin ipols and store ipol persistency strings for each histo"
+ while True:
+ if q.empty():
+ break
+ hn = q.get()
+ histos = HISTOS[hn]
+ ih = prof.mk_ipolhisto(histos, RUNS, PARAMSLIST, ORDER, "none", None)
+ if ih is None:
+ print ("Ignoring", hn)
+ else:
+ s = ""
+ for i, ib in enumerate(ih.bins):
+ s += "%s#%d %.5e %.5e\n" % (hn, i, ib.xmin, ib.xmax)
+ s += " " + ib.ival.toString("val") + "\n"
+ if ib.ierrs:
+ s += " " + ib.ierrs.toString("err") + "\n"
+ rdict[hn] = zlib.compress(s, 9) #< save some memory
+ del s
+ del ih #< pro-actively clear up memory
+ del histos
+
+
+ print ("\n\nParametrising...\n")
+ import time, multiprocessing
+ time1 = time.time()
+
+ ## A shared memory object is required for coefficient retrieval
+ from multiprocessing import Manager
+ manager = Manager()
+ tempDict = manager.dict()
+
+ ## The job queue
+ q = multiprocessing.Queue()
+ map(lambda x:q.put(x), HNAMES)
+
+ ## Fire away
+ workers = [multiprocessing.Process(target=worker, args=(q, tempDict)) for i in range(opts.MULTI)]
+ map(lambda x:x.start(), workers)
+ map(lambda x:x.join(), workers)
+
+ ## Finally copy the result dictionary into the object itself
+ for k in tempDict.keys():
+ IHISTOS[k] = tempDict[k]
+
+ ## Timing
+ time2 = time.time()
+ print ('Parametrisation took %0.2f s' % ((time2-time1)))
+
+
+ ## Write out meta info
+ # TODO: Move the format writing into the prof2 Python library
+ with open(IFILE, "w") as f:
+ if opts.SUMMARY is not None:
+ f.write("Summary: %s\n" % opts.SUMMARY)
+ #f.write("DataDir: %s\n" % os.path.abspath(RUNSDIR))
+ f.write("ProfVersion: %s\n" % prof.version())
+ f.write("Date: %s\n" % prof.mk_timestamp())
+ f.write("DataFormat: binned 2\n") # This tells the reader how to treat the coefficients that follow
+ # Format and write out parameter names
+ pstring = "ParamNames:"
+ for p in PARAMNAMES:
+ pstring += " %s" % p
+ f.write(pstring + "\n")
+ # Dimension (consistency check)
+ f.write("Dimension: %i\n" % len(PARAMNAMES))
+ # Interpolation validity (hypercube edges)
+ minstring = "MinParamVals:"
+ for v in prof.mk_minvals(PARAMSLIST):
+ minstring += " %f" % v
+ f.write(minstring + "\n")
+ maxstring = "MaxParamVals:"
+ for v in prof.mk_maxvals(PARAMSLIST):
+ maxstring += " %f" % v
+ f.write(maxstring + "\n")
+ f.write("DoParamScaling: 1\n")
+ # Number of inputs per bin
+ f.write("NumInputs: %i\n" % len(PARAMSLIST))
+ s_runs = "Runs:"
+ for r in RUNS:
+ s_runs +=" %s"%r
+ f.write("%s\n"%s_runs)
+ f.write("---\n")
+
+ ## Write out numerical data for all interpolations
+ s = ""
+ for hname in sorted(IHISTOS.keys()):
+ ih = IHISTOS[hname]
+ s += zlib.decompress(ih)
+
+ # Open file for write/append
+ with open(IFILE, "a") as f:
+ f.write(s)
+ print ("\nOutput written to %s" % IFILE)
+
+
+
+
+
+
+
+def preparePulls(runs, params, it, order):
+ from numpy import zeros, ones, mean, std, square, sqrt
+ IFILE=mkIpolFilename(it, order)
+ if not os.path.exists(IFILE):
+ mkIpols(runs, params, o, IFILE)
+
+ IHISTOS, IMETA = prof.read_ipoldata(IFILE)
+ print("Loaded ipol histos from", IFILE)
+ USEDRUNS=IMETA["Runs"].split()
+ KNIFERUNS=[r for r in RUNS if not r in USEDRUNS]
+
+ PULL = {}
+ for hn in sorted(HNAMES):
+ # for hn in sorted(IHISTOS.keys()):
+ IH = IHISTOS[hn]
+ nbins= IH.nbins
+ pull = 20*ones((len(KNIFERUNS),nbins))
+ for num, run in enumerate(KNIFERUNS):
+ P = PARAMS[run]
+ temp = IH.toDataHisto(P)
+ for i in xrange(nbins):
+ mcval = HISTOS[hn][run].bins[i].val
+ ipval = temp.bins[i].val
+ # if ipval<0:
+ # ipval=1e20
+ # if mcval>0:
+ pull[num][i] = (ipval - mcval)/ipval
+ # else:
+ # pull[num][i] = 0
+ PULL[hn]= pull
+ return PULL
+
+
+xc = [r for r in prof.sampling.xrandomUniqueCombinations(RUNS, len(RUNS)-opts.DROPOUT, opts.NITERATIONS)]
+
+ALL ={}
+import numpy
+for IT in xrange(opts.NITERATIONS):
+ thisRC=xc[IT]
+ # Filtering
+ thisRUNS, thisPARAMSLIST = [], []
+ for num, r in enumerate(RUNS):
+ if r in thisRC:
+ thisRUNS.append(r)
+ thisPARAMSLIST.append(PARAMSLIST[num])
+
+ for o in ORDERS:
+ ttt = preparePulls(thisRUNS, thisPARAMSLIST, IT, o)
+ if ALL.has_key(o):
+ # from IPython import embed
+ # embed()
+ # sys.exit(1)
+ for k, v in ALL[o].iteritems():
+ ALL[o][k] = numpy.append(v, ttt[k], axis=0)
+ else:
+ ALL[o] = ttt
+
+# Got through observables
+BEST={}
+HNAMES=sorted(ALL[ALL.keys()[0]].keys())
+for hn in HNAMES:
+ if ALL[ORDERS[0]][hn].ndim == 1:
+ nbins = 1
+ else:
+ nbins=ALL[ORDERS[0]][hn].shape[1]
+ nentries=ALL[ORDERS[0]][hn].shape[0]
+ temp_best=[]
+ # Iterate over bins
+ for i in xrange(nbins):
+ temp_pull = []
+ for o in ORDERS:
+ if nbins>1:
+ temp_pull.append(sum([x*x for x in ALL[o][hn][:,i]])/nentries)
+ else:
+ temp_pull.append(sum([x*x for x in ALL[o][hn]])/nentries)
+ gof = [numpy.sqrt(temp_pull[x]) for x in xrange(len(ORDERS))]
+ best = min(gof)
+ bestorderforbin = ORDERS[gof.index(best)]
+ temp_best.append(bestorderforbin)
+ import pylab
+ pylab.clf()
+ if any([x<=0 for x in gof]):
+ pylab.plot(ORDERS, gof)
+ else:
+ pylab.semilogy(ORDERS, gof)
+ pylab.title("pull %s bin#%i"%(hn ,i))
+ pylab.savefig("pull_%s_%i_%i.pdf"%(hn.replace("/","_"),opts.NITERATIONS, i))
+ pylab.clf()
+ pylab.plot(ORDERS[0:-1], [gof[j+1]/gof[j] for j in xrange(len(gof)-1)])
+ print ([gof[j+1]/gof[j] for j in xrange(len(gof)-1)])
+ pylab.title("pull ratio %s bin#%i"%(hn ,i))
+ pylab.savefig("ratio_%s_%i_%i.pdf"%(hn.replace("/","_"),opts.NITERATIONS, i))
+ BEST[hn] = temp_best
+
+print (BEST)
+from IPython import embed
+embed()
+sys.exit(1)
+

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:31 PM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806047
Default Alt Text
(47 KB)

Event Timeline