No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	Index: trunk/unfold.py
	===================================================================
	--- trunk/unfold.py (revision 20)
	+++ trunk/unfold.py (revision 21)
	@@ -1,130 +1,115 @@
	#!/usr/bin/env python
	# -- coding: utf-8 --

	import os
	import sys
	import numpy as np
	from math import *
	from utils import *
	from aru import *
	from kernel import MyKernel

	def GetData(inputFileName):
	from ROOT import TFile
	ntuple = TFile.Open(inputFileName).GetListOfKeys()[0].ReadObj()
	data = np.empty(ntuple.GetEntries())
	for i in xrange(ntuple.GetEntries()):
	ntuple.GetEntry(i)
	data[i] = ntuple.GetArgs()[0]
	return data

	# define input/output, minimal checks
	try:
	inputFilename = sys.argv[1]
	outputFilename = sys.argv[2]
	if not os.path.exists(inputFilename):
	raise SystemExit("Input file %s does not exist."%inputFilename)
	if os.path.exists(outputFilename):
	raise SystemExit("Output file %s does exist, I don't dare to overwrite."%outputFilename)
	except IndexError:
	raise SystemExit("Wrong number of arguments. Usage: unfold.py <inputFile> <outputFile>")

	# inputFile contains the raw data as a TNtupleD with a single branch
	print "reading data..."
	data = GetData(inputFilename)

	# replace implementation of MyKernel by a proper kernel
	kernel = MyKernel()

	# x-range to be unfolded, may be set to other values
	xmin = np.min(data)
	xmax = np.max(data)

	# customize number of knots (use a sufficiently large number) and their positions
	nKnots = 20
	xknots = np.linspace(xmin,xmax,nKnots)

	# main program, encapsulated in ARU object
	print "unfolding..."
	aru = ARU(kernel, ToVector(xknots))
	-c_coefs, c_cov = aru.Unfold(ToVector(data))
	+coefs, cov = aru.Unfold(ToVector(data))
	unfolded = aru.GetUnfoldedSpline()
	folded = aru.GetFoldedSpline()
	-c_refCoefs = aru.GetReferenceCoefficients()
	+refCoefs = aru.GetReferenceCoefficients()

	-coefs = Extract(c_coefs)
	-cov = Extract(c_cov)
	+coefs = Extract(coefs)
	+cov = Extract(cov)

	# writing the output file with contents
	# - TNtupleD with variables
	# x : sample value
	# yref : final reference distribution (for checking)
	# yfitted : folded distribution (the actual fit to the data)
	# yunfolded: unfolded distribution (the solution)
	# ysigma : uncertainty of yunfolded
	# - TVectorD, coefficients of the solution
	# - TMatrixDSym, covariance matrix of the solution
	print "Writing result..."

	-class Sigma(object):
	- def __init__(self,unfolded, coefs, cov):
	- self.unfolded = unfolded
	- self.coefs = coefs
	- self.cov = cov
	- def __call__(self,x):
	- n = self.unfolded.GetSize()
	- bs = [ unfolded.Basis(k)(x) for k in xrange(n) ]
	- return np.sqrt( np.dot(np.dot(bs,cov),bs) )
	-
	xs = np.linspace(xmin,xmax,200)

	-yref = np.vectorize(lambda x: unfolded(c_refCoefs,x))(xs)
	-yfitted = np.vectorize(lambda x: folded (c_coefs,x))(xs)
	-yunfolded = np.vectorize(lambda x: unfolded(c_coefs,x))(xs)
	-sigma = Sigma(unfolded, coefs, cov)
	-ysigma = np.vectorize(sigma)(xs)
	+yref = np.vectorize(lambda x: unfolded(refCoefs,x))(xs)
	+yfitted = np.vectorize(lambda x: folded (coefs,x))(xs)
	+yunfolded = np.vectorize(lambda x: unfolded(coefs,x))(xs)
	+ysigma = GetSigma(unfolded,cov,xs)

	from ROOT import TNtupleD, TFile, TMatrixDSym, TVectorD
	nt = TNtupleD("nt","","x:yref:yfitted:yunfolded:ysigma")
	for i in xrange(len(xs)):
	nt.Fill(xs[i],yref[i],yfitted[i],yunfolded[i],ysigma[i])

	rootCoefs = TVectorD(len(coefs))
	rootCov = TMatrixDSym(len(coefs))
	for i in xrange(len(coefs)):
	rootCoefs[i] = coefs[i]
	for j in xrange(len(coefs)):
	rootCov[i][j] = cov[i][j]
	f = TFile.Open(outputFilename,"RECREATE")
	nt.Write()
	rootCoefs.Write()
	rootCov.Write()
	f.Close()

	try:
	# draw everything, if matplotlib is installed
	from matplotlib import pyplot as plt

	- from matplotlib.patches import Polygon
	- poly = Polygon(MakeSigmaPatch(unfolded, c_coefs, cov, xknots), closed=True)
	- poly.set(ec="b",fc="b",alpha=0.1,fill=True,zorder=0)
	- plt.gca().add_patch(poly)
	-
	plt.plot(xs,yref,"y",lw=2,label="$g(x)$",zorder=1)

	plt.plot(xs,yfitted,"r",lw=1,label="$f(y)$",zorder=2)

	plt.plot(xs,yunfolded,"b",lw=2,label="$b(x)$",zorder=2)
	+ plt.fill_between(xs,yunfolded-ysigma,yunfolded+ysigma,edgecolor="b",facecolor="b",alpha=0.1,zorder=0)

	xh, w, werr = GetScaledHistogram(data,bins=20,range=(xmin,xmax))
	plt.errorbar(xh,w,werr,fmt="ko",capsize=0,label="data",zorder=3)

	plt.ylim(ymin=0)
	plt.xlim(xmin,xmax)
	plt.xlabel("x")
	plt.ylabel(r"$N_\mathrm{data} \times p.d.f.$")
	plt.legend(loc="upper left")
	plt.show()

	except ImportError:
	print "Install matplotlib to get plots of the solution."
	Index: trunk/tests/blobel.py
	===================================================================
	--- trunk/tests/blobel.py (revision 20)
	+++ trunk/tests/blobel.py (revision 21)
	@@ -1,87 +1,82 @@
	#!/usr/bin/env python
	# -- coding: utf-8 --

	import sys
	import numpy as np
	from matplotlib import pyplot as plt
	from utils import *
	from aru import *
	from kernel import BlobelKernel

	kernel = BlobelKernel()

	print "generate"
	try:
	np.random.seed(int(sys.argv[1]))
	except IndexError:
	np.random.seed(1)

	nMC = 5000

	def Model(x):
	b = (1.0,10.0,5.0)
	g = (2.0,0.2,0.2)
	c = (0.4,0.8,1.5)
	result = 0.0
	for k in xrange(3):
	result += b[k]g[k]2/((x-c[k])2+g[k]*2)
	return result

	rxs = []
	while True:
	xy = np.random.rand(2)
	x = xy[0]*2
	y = xy[1]*12
	if y < Model(x):
	rxs.append(x)
	if len(rxs) == nMC:
	break

	data = []
	for x in rxs:
	y = kernel.Y(x) + 0.1*np.random.randn()
	z = np.random.rand()
	if z < kernel.Efficiency(y):
	data.append(y)

	xknots = np.linspace(0,2,20)

	aru = ARU(kernel, ToVector(xknots))
	coefs, cov = aru.Unfold(ToVector(data))
	unfolded = aru.GetUnfoldedSpline()
	folded = aru.GetFoldedSpline()
	refCoefs = aru.GetReferenceCoefficients()

	# draw everything

	xs = np.linspace(xknots[0],xknots[-1],200)

	from scipy.integrate import quad
	ymodel = nMC*Model(xs)/quad(Model,0,2)[0]
	yprior = np.vectorize(lambda x: unfolded(refCoefs,x))(xs)
	yfitted = np.vectorize(lambda x: folded (coefs,x))(xs[xs<1.8])
	yunfolded = np.vectorize(lambda x: unfolded(coefs,x))(xs)
	-
	-plt.figure(1)
	-
	-from matplotlib.patches import Polygon
	-poly = Polygon(MakeSigmaPatch(unfolded, coefs, cov, xknots), closed=True)
	-poly.set(ec="b",fc="b",alpha=0.1,fill=True,zorder=0)
	-plt.gca().add_patch(poly)
	+ysigma = GetSigma(unfolded,cov,xs)

	plt.plot(xs,ymodel,"g--",lw=2,label="$t(x)$",zorder=1)

	plt.plot(xs,yprior,"y",lw=3,label="$g(x)$",zorder=1)

	plt.plot(xs[xs<1.8],yfitted,"r",lw=1,label="$f(y)$",zorder=3)

	plt.plot(xs,yunfolded,"b",lw=2,label="$b(x)$",zorder=2)
	+plt.fill_between(xs,yunfolded-ysigma,yunfolded+ysigma,edgecolor="b",facecolor="b",alpha=0.1,zorder=0)

	xh, w, werr = GetScaledHistogram(data,bins=20,range=(folded.GetStart(),folded.GetStop()))
	plt.errorbar(xh,w,werr,fmt="ko",capsize=0,label="data",zorder=4)

	plt.ylim(ymin=0)
	plt.xlabel("x")
	plt.ylabel(r"$\mathrm{d}N / \mathrm{d} x$")
	plt.legend(loc="upper left").get_frame().set_visible(0)
	plt.show()
	Index: trunk/tests/foldspline.py
	===================================================================
	--- trunk/tests/foldspline.py (revision 0)
	+++ trunk/tests/foldspline.py (revision 21)
	@@ -0,0 +1,44 @@
	+#!/usr/bin/env python
	+# -- coding: utf-8 --
	+import numpy as np
	+from matplotlib import pyplot as plt
	+from math import *
	+from pyik.mplext import *
	+from pyik.numpyext import *
	+from aru import *
	+from kernel import *
	+
	+kernel = GaussianKernel(0.1);
	+
	+nKnots = 10
	+xknots = KnotVector(nKnots);
	+for i in xrange(nKnots):
	+ xknots[i] = float(i)/(nKnots-1)
	+
	+unfolded = SplineFunction1D(xknots)
	+folded = FoldedSplineFunction1D(kernel, unfolded)
	+
	+xs = np.linspace(-0.1,1.1,200)
	+
	+selKnot = 1
	+f1 = lambda x : unfolded[selKnot](x)
	+f2 = lambda x : folded[selKnot](x)
	+f1s = np.vectorize(f1)(xs)
	+f2s = np.vectorize(f2)(xs)
	+
	+from scipy.integrate import *
	+
	+f1ia = unfolded[selKnot].Integral()
	+f1ib = quad(f1,unfolded.GetStart(),unfolded.GetStop())[0]
	+f2ia = folded[selKnot].Integral(folded.GetStart(),folded.GetStop())
	+f2ib = quad(f2,folded.GetStart(),folded.GetStop())[0]
	+
	+print f1ia,"[ %g, acc. %i ]"%(f1ib,log10(abs(f1ia-f1ib)))
	+print f2ia,"[ %.12g, acc. %i ]"%(f2ib,log10(abs(f2ia-f2ib)))
	+
	+f1mask = f1s>0
	+f2mask = f2s>0
	+plt.plot(xs[f1mask],f1s[f1mask],"k-")
	+plt.plot(xs[f2mask],f2s[f2mask],"r-")
	+plt.xlim(xknots[0],xknots[xknots.size()-1])
	+plt.show()

	Property changes on: trunk/tests/foldspline.py
	___________________________________________________________________
	Added: svn:executable
	## -0,0 +1 ##
	+*
	\ No newline at end of property
	Index: trunk/tests/regfit.py
	===================================================================
	--- trunk/tests/regfit.py (revision 20)
	+++ trunk/tests/regfit.py (revision 21)
	@@ -1,86 +1,83 @@
	#!/usr/bin/env python
	# -- coding: utf-8 --

	import sys
	import numpy as np
	from matplotlib import pyplot as plt
	from pyik.mplext import auto_adjust
	from utils import *
	from aru import *
	from kernel import *

	print "generate"
	np.random.seed(int(sys.argv[1]))

	nMC = int(sys.argv[2])

	norms = (0.3,0.7)
	means = (0.3,0.7)
	sigmas = (0.1,0.05)
	kernelSigma = 0.1

	# ~norms = (1.,20.)
	# ~means = (0.4,0.6)
	# ~sigmas = (0.05,0.2)
	# ~kernelSigma = 0.1

	model = ToyModel(norms,means,sigmas)
	data = model.rvs(nMC)

	kernel = GaussianKernel(kernelSigma)
	data += kernelSigma*np.random.randn(nMC)

	-xknots = np.linspace(0,1,40)
	+xknots = np.linspace(0,1,20)

	aru = ARU(kernel, ToVector(xknots))
	-c_coefs, c_cov = aru.Unfold(ToVector(data))
	+coefs, cov = aru.Unfold(ToVector(data))
	unfolded = aru.GetUnfoldedSpline()
	folded = aru.GetFoldedSpline()
	-c_refCoefs = aru.GetReferenceCoefficients()
	+refCoefs = aru.GetReferenceCoefficients()

	modelNorm = np.sum((folded.GetStart()<=data)*(data<=folded.GetStop()))

	csum = 0.0
	for k in xrange(unfolded.GetSize()):
	- csum += c_coefs[k]*unfolded[k].Integral()
	+ csum += coefs[k]*unfolded[k].Integral()

	print "integral comparison",modelNorm, csum

	# draw everything

	xs = np.linspace(0,1,200)

	ysmodel = nMC*model(xs)
	-ysprior = np.vectorize(lambda x: unfolded(c_refCoefs,x))(xs)
	-ysfitted = np.vectorize(lambda x: folded (c_coefs,x))(xs)
	-ysunfolded = np.vectorize(lambda x: unfolded(c_coefs,x))(xs)
	-
	-from matplotlib.patches import Polygon
	-poly = Polygon(MakeSigmaPatch(unfolded, c_coefs, c_cov, xknots), closed=True)
	-poly.set(ec="b",fc="b",alpha=0.1,fill=True,zorder=0)
	-plt.gca().add_patch(poly)
	+ysprior = np.vectorize(lambda x: unfolded(refCoefs,x))(xs)
	+ysfitted = np.vectorize(lambda x: folded (coefs,x))(xs)
	+ysunfolded = np.vectorize(lambda x: unfolded(coefs,x))(xs)
	+yssigma = GetSigma(unfolded,cov,xs)

	plt.plot(xs,ysmodel,"g--",lw=2,label="$t(x)$",zorder=1)

	plt.plot(xs,ysprior,"y",lw=3,label="$g(x)$",zorder=1)

	plt.plot(xs,ysfitted,"r",lw=1,label="$f(y)$",zorder=3)

	plt.plot(xs,ysunfolded,"b",lw=2,label="$b(x)$",zorder=2)
	+plt.fill_between(xs, ysunfolded-yssigma, ysunfolded+yssigma, edgecolor="b", facecolor="b", alpha=0.1, zorder=0)

	xh, w, werr = GetScaledHistogram(data,bins=20,range=(0,1))
	plt.errorbar(xh,w,werr,fmt="ko",capsize=0,label="data",zorder=4)

	# ~print "noreg fit"
	# ~c_coefs_noreg = std.vector("double")(nBasis,1)
	# ~objective.SetWeight(0)
	# ~Fit(objective,c_coefs_noreg)
	# ~yunreg = np.vectorize(lambda xx: unfolded(c_coefs_noreg,xx))(x)
	# ~plt.plot(x,yunreg,"b:",lw=1,label="$b_\mathrm{w=0}(x)$",zorder=1)

	plt.ylim(ymin=0)
	plt.xlabel("x")
	plt.ylabel(r"$\mathrm{d}N / \mathrm{d} x$")
	plt.legend(loc="upper left").get_frame().set_visible(0)

	plt.show()
	Index: trunk/src/Integrator.h
	===================================================================
	--- trunk/src/Integrator.h (revision 20)
	+++ trunk/src/Integrator.h (revision 21)
	@@ -1,164 +1,155 @@
	/*
	This file is part of ARU, a software to unfold detector effects
	from data distributions.

	Copyright (C) 2011 Karlsruhe Instiute of Technology,
	contact Hans Dembinski <hans.dembinski@kit.edu>

	ARU is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	ARU is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with ARU. If not, see <http://www.gnu.org/licenses/>.
	*/

	#ifndef _Aru_Integrator_h_
	#define _Aru_Integrator_h_

	#include <cmath>
	#include <vector>
	#include <stdexcept>
	#include <string>
	#include <sstream>
	#include <limits>

	namespace Aru {

	namespace OfflineIntegrator {

	// based on Offline/Utilities/Math/Integrator.h
	// from the Pierre Auger Observatory;
	// LICENSE.Auger applies to everything inside
	// the current namespace

	inline
	double
	PolynomialInterpolation(const size_t n,
	const double px[],
	const double py[],
	const double x,
	double& dy)
	{
	if (!n)
	return 0;

	double minDist = std::fabs(x - px[0]);
	size_t k = 0;
	double c[n];
	double d[n];
	c[0] = py[0];
	d[0] = py[0];
	double oldx = px[0];
	for (size_t i = 1; i < n; ++i) {
	const double newx = px[i];
	oldx = newx;
	const double dist = std::fabs(x - newx);
	if (dist < minDist) {
	minDist = dist;
	k = i;
	}
	c[i] = d[i] = py[i];
	}

	double y = py[k];

	for (size_t m = 1; m < n; ++m) {
	const size_t nm = n - m;
	for (size_t i = 0; i < nm; ++i) {
	const double xa = px[i];
	const double xb = px[i+m];
	const double dx = xb - xa;
	const double cd = (c[i+1] - d[i]) / dx;
	c[i] = (x-xa)*cd;
	d[i] = (x-xb)*cd;
	}
	dy = (2*k < nm) ? c[k] : d[--k];
	y += dy;
	}

	return y;
	}

	/// average of a function represented with equidistant boxes
	template<class Functor>
	inline
	double
	GetBoxAverage(const Functor& functor, const double start, const double step, const size_t n)
	{
	double sum = 0;
	for (size_t i = 0; i < n; ++i) {
	const double x = start + i*step;
	sum += functor(x);
	}
	return sum/n;
	}

	template<class Functor>
	- inline
	- double
	- GetTrapezoidalAverage(const Functor& functor,
	- const double previousApproximation,
	- const double a, const double delta, const size_t level)
	- {
	- const size_t n = 1 << (level - 1); // = 2^(level-1)
	- const double step = delta/n;
	- return 0.5*(previousApproximation +
	- GetBoxAverage(functor, a+0.5*step, step, n));
	- }
	-
	- template<class Functor>
	double
	Integrate(const Functor& functor, const double a, const double b,
	const double accuracy = 1e-8,
	const size_t order = 5, const size_t maxIterations = 15)
	{
	// calculates trapezoidal average of functor as
	// a function of decreasing stepsize h, then
	// extrapolates series to h = 0

	const double delta = b - a;
	const double eps = 0.5*accuracy;

	std::vector<double> tInt(maxIterations+1);
	std::vector<double> tStep(maxIterations+1);

	- double oldInt = tInt[0] = GetBoxAverage(functor, a, delta, 2);
	+ tInt[0] = GetBoxAverage(functor, a, delta, 2);
	tStep[0] = delta;
	+ tInt[1] = 0.5(tInt[0] + functor(a+0.5delta));
	+ tStep[1] = delta/4; // h^2

	double extrapolation = 0;
	- for (size_t i=1; i <= maxIterations; ++i) {
	- oldInt = tInt[i] = GetTrapezoidalAverage(functor, oldInt, a, delta, i);
	- tStep[i] = delta/(1 << (2*i));
	+ for (size_t i=2; i <= maxIterations; ++i) {
	+ const size_t n = 1 << (i - 1); // = 2^(j-1)
	+ const double step = delta/n;
	+ tInt[i] = 0.5(tInt[i-1] + GetBoxAverage(functor, a+0.5step, step, n));
	+ tStep[i] = delta/(1 << (2*i)); // h^2

	const int first = i - order + 1;
	if (first < 0) continue;
	double residual;
	extrapolation =
	PolynomialInterpolation(order, &tStep[first], &tInt[first], 0, residual);
	if (std::isnan(extrapolation))
	throw std::runtime_error("got nan during integration");
	if (std::fabs(residual) < eps*std::fabs(extrapolation))
	return delta*extrapolation;
	}

	if (std::fabs(extrapolation) > accuracy)
	throw std::runtime_error("maximal level reached in integration, result still inaccurate");

	return delta*extrapolation;
	}

	} // NS OfflineIntegrator

	using OfflineIntegrator::Integrate;

	}

	#endif

	Index: trunk/utils.py
	===================================================================
	--- trunk/utils.py (revision 20)
	+++ trunk/utils.py (revision 21)
	@@ -1,430 +1,404 @@
	# -- coding: utf-8 --

	# This file is part of ARU, a software to unfold detector effects
	# from data distributions.
	#
	# Copyright (C) 2011 Karlsruhe Instiute of Technology,
	# contact Hans Dembinski <hans.dembinski@kit.edu>
	#
	# ARU is free software: you can redistribute it and/or modify
	# it under the terms of the GNU General Public License as published by
	# the Free Software Foundation, either version 3 of the License, or
	# (at your option) any later version.
	#
	# Foobar is distributed in the hope that it will be useful,
	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	# GNU General Public License for more details.
	#
	# You should have received a copy of the GNU General Public License
	# along with Foobar. If not, see <http://www.gnu.org/licenses/>.

	import numpy as np
	from math import *
	from scipy.stats import norm

	def Extract(obj):
	"""
	Extracts data from ROOT objects and returns corresponding Numpy arrays.

	Extract can be used with a string argument. In this case, the string is
	interpreted as a filename of a ROOT file. The ROOT file is opened and
	Extract is applied on each object inside the file. The result is returned
	as a dictionary, mapping the key names to the outputs.

	Authors
	-------
	Hans Dembinski <hans.dembinski@kit.edu>
	"""

	import ROOT
	import numpy as np
	from aru import Vector, Matrix

	if isinstance(obj, Vector):
	a = np.empty(obj.size())
	for i in xrange(len(a)):
	a[i] = obj[i]
	return a

	if isinstance(obj, Matrix):
	m = obj.GetShape(0)
	n = obj.GetShape(1)
	a = np.empty((m,n))
	for i in xrange(m):
	for j in xrange(n):
	a[i,j] = obj[i][j]
	return a

	if isinstance(obj, ROOT.std.vector("double")):
	a = np.empty(obj.size())
	for i in xrange(len(a)):
	a[i] = obj[i]
	return a

	if isinstance(obj, str):
	f = ROOT.TFile.Open(obj)
	d = {}
	for key in f.GetListOfKeys():
	d[key.GetName()] = Extract(key.ReadObj())
	return d

	if isinstance(obj, ROOT.TVectorD):
	n = obj.GetNoElements()
	a = np.empty(n)
	for i in xrange(n):
	a[i] = obj[i]
	return a

	if isinstance(obj, ROOT.TMatrixD) or isinstance(obj, ROOT.TMatrixDSym):
	a = np.empty((obj.GetNcols(),obj.GetNrows()))
	for i in xrange(obj.GetNcols()):
	for j in xrange(obj.GetNrows()):
	a[i,j] = obj(i,j)
	return a

	raise ValueError("Extract cannot handle object %s yet, please take a look at pyik.rootext and implement it"%(type(obj)))

	def ToVector(array):
	from aru import Vector
	n = len(array)
	res = Vector(n)
	for i in xrange(n):
	res[i] = array[i]
	return res

	def centers(x):
	"""
	Computes the centers of an array of bin edges.

	Parameters
	----------
	x: array-like
	A 1-d array containing lower bin edges.

	Returns
	-------
	c: array of dtype float
	Returns the centers of the bins.
	hw: array of dtype float
	Returns the half-width of the bins.

	Examples
	--------
	>>> c,hw = centers([0.0,1.0,2.0])
	>>> print c, hw
	[ 0.5 1.5] [ 0.5 0.5]

	Authors
	-------
	Hans Dembinski <hans.dembinski@kit.edu>
	"""

	n = len(x)-1
	c = np.empty(n)
	hw = np.empty(n)
	for i in xrange(len(x)-1):
	hw[i] = 0.5*(x[i+1]-x[i])
	c[i] = x[i] + hw[i]
	return c,hw

	class ToyModel:
	def __init__(self,norms,mus,sigmas):
	self.norms = np.array(norms)
	self.norms /= np.sum(self.norms)
	self.mus = mus
	self.sigmas = sigmas

	def __call__(self,x):
	result = 0.0
	n = len(self.norms)
	for i in xrange(n):
	result += self.norms[i]*norm(self.mus[i],self.sigmas[i]).pdf(x)
	return result

	def rvs(self, nEvents):
	n = len(self.norms)
	ps = np.zeros(n)
	for i in xrange(n):
	ps[i:] += self.norms[i]
	ps /= np.sum(self.norms)
	counts = np.zeros(n,dtype=int)
	for i in xrange(nEvents):
	for j,p in enumerate(ps):
	uniran = np.random.rand()
	if uniran < p:
	counts[j] += 1
	break
	rvs = np.array([])
	for i in xrange(n):
	rvs = np.append(rvs,norm(self.mus[i],self.sigmas[i]).rvs(counts[i]))
	return rvs

	def MinimumMeanSquaredErrorFit(obj, starts):
	from aru import Vector, Matrix
	nBasis = obj.GetBasisSize()

	foldedCache = obj.GetFoldedCache()
	foldedMuCache = obj.GetFoldedMuCache()

	nUnbinned = foldedCache.GetShape(0)
	nBinned = foldedMuCache.GetShape(0)
	nEventsForHistogram = obj.GetNEventsForHistogram()

	fmatrix = obj.GetFMatrix()

	# start values
	c_coefs = Vector(nBasis)
	for i in xrange(nBasis): c_coefs[i] = float(starts[i])
	c_cov = Matrix(nBasis)

	def MeanSquaredError(sqrtw):
	w = sqrtw*sqrtw

	obj.SetWeight(w)
	obj.Minimize(c_coefs)
	obj.CovarianceStat(c_cov, c_coefs)
	obj.UpdateReference(c_coefs)

	result = 0.0
	if w != 0.0:
	for i in xrange(nBasis):
	for j in xrange(i,nBasis):
	result += (1 if i==j else 2)(c_cov[i][j]+c_coefs[i]c_coefs[j])*fmatrix[i][j]

	for i in xrange(nUnbinned):
	folded = 0.0
	for k in xrange(nBasis):
	folded += c_coefs[k]*foldedCache[i][k]
	result -= 2*folded

	for ib in xrange(nBinned):
	hw = obj.GetHistogramWeight(ib)
	if hw<nEventsForHistogram: continue
	foldedMu = 0.0
	for k in xrange(nBasis):
	foldedMu += c_coefs[k]*foldedMuCache[ib][k]
	result -= 2hwfoldedMu/obj.GetHistogramStep(ib)

	print "iterating (%.5e): current weight = %g"%(result,w)

	return result

	from scipy.optimize import brent
	sqrtw = brent(MeanSquaredError, brack=(0,1), tol=1e-2)
	w = sqrtw*sqrtw
	obj.SetWeight(w)
	obj.Minimize(c_coefs)
	return c_coefs

	def GetScaledHistogram(data,bins,range=None):
	ws, xedges = np.histogram(data,bins=bins,range=range)
	if type(bins) is not int:
	bins = len(xedges)-1
	xs = centers(xedges)[0]
	werrs = np.sqrt(ws)
	factor = bins/float(xedges[-1]-xedges[0]) # due to binning
	ws = np.asarray(ws,np.float)
	ws *= factor
	werrs *= factor
	return xs, ws, werrs

	-def MakeSigmaPatch(basis, coef, cov, xknot):
	+def GetSigma(basis, cov, xs):

	- nover = 10
	-
	- xs = []
	- for i in xrange(len(xknot)-1):
	- x0 = xknot[i]
	- x1 = xknot[i+1]
	-
	- dx = (x1-x0)/nover
	-
	- for j in xrange(nover-1):
	- xs.append(x0+j*dx)
	- xs.append(xknot[-1])
	-
	- def CalcSigma(basis,cov,x):
	+ def Sigma(x):
	n = basis.GetSize()
	s2 = 0.0
	for l in xrange(n):
	- for m in xrange(n):
	- s2 += basis[l](x)basis[m](x)cov[l][m]
	+ for m in xrange(l,n):
	+ s2 += (1 if l==m else 2)basis[l](x)basis[m](x)*cov[l][m]
	if s2 < 0:
	+ print s2,"< 0! at",x
	s2 = 0
	return sqrt(s2)

	- y1 = [ basis(coef,x) for x in xs ]
	-
	- y2 = [ CalcSigma(basis,cov,x) for x in xs ]
	-
	- xy = []
	-
	- n = len(xs)
	- # upper
	- for i in xrange(n):
	- xy.append([ xs[i], y1[i]+y2[i] ])
	- # lower
	- for i in xrange(n-1,0,-1):
	- xy.append([ xs[i], y1[i]-y2[i] ])
	-
	- return np.array(xy)
	+ return np.vectorize(Sigma)(xs)

	def JackknifeShifts(obj,c_coefs):
	from aru import TMatrixDSym
	nBasis = obj.GetBasisSize()

	grad0 = Extract(obj.Gradient(c_coefs))
	c_hesse1 = TMatrixDSym(nBasis)
	obj.HesseStat(c_hesse1, c_coefs)
	c_hesse2 = TMatrixDSym(nBasis)
	obj.HesseReg(c_hesse2, c_coefs)
	hesse0 = Extract(c_hesse1) + obj.GetWeight()*Extract(c_hesse2)

	foldedCache = obj.GetFoldedCache
	foldedMuCache = obj.GetFoldedMuCache

	shifts = []
	coefs = Extract(c_coefs)
	for i in xrange(obj.GetUnbinnedSize()):
	folded = 0.0
	for k in xrange(nBasis):
	folded += coefs[k]*foldedCache(k,i)

	grad = np.copy(grad0)
	hesse = np.copy(hesse0)
	for k in xrange(nBasis):
	grad[k] += foldedCache(k,i)/folded
	for l in xrange(k,nBasis):
	delta = foldedCache(k,i)foldedCache(l,i)/folded*2
	hesse[k,l] -= delta
	if k!=l:
	hesse[l,k] -= delta

	u,s,vt = np.linalg.svd(hesse)
	hesseInv = np.zeros((nBasis,nBasis))
	for k in xrange(nBasis):
	for l in xrange(k,nBasis):
	for m in xrange(nBasis):
	hesseInv[k,l] = vt[k,m]*u[l,m]/s[m]
	if k!=l:
	hesseInv[l,k] = hesseInv[k,l]
	shifts.append(-np.dot(hesseInv,grad))

	for ib in xrange(obj.GetBinnedSize()):
	hw = obj.GetHistogramWeight(ib)

	if hw==0:continue

	foldedMu = 0.0
	for k in xrange(nBasis):
	foldedMu += coefs[k]*foldedMuCache(k,ib)

	grad = np.copy(grad0)
	hesse = np.copy(hesse0)
	for k in xrange(nBasis):
	grad[k] += hw*foldedMuCache(k,i)/foldedMu
	for l in xrange(k,nBasis):
	delta = hwfoldedMuCache(k,i)foldedMuCache(l,i)/foldedMu**2
	hesse[k,l] -= delta
	if k!=l:
	hesse[l,k] -= delta

	u,s,vt = np.linalg.svd(hesse)
	hesseInv = np.zeros((nBasis,nBasis))
	for k in xrange(nBasis):
	for l in xrange(k,nBasis):
	for m in xrange(nBasis):
	hesseInv[k,l] = vt[k,m]*u[l,m]/s[m]
	if k!=l:
	hesseInv[l,k] = hesseInv[k,l]
	shifts.append(-np.dot(hesseInv,grad))

	return shifts

	def GetEmpiricalDensity(kernel,unfolded,data):
	xa = unfolded.GetStart()
	xb = unfolded.GetStop()
	result = 0.0
	for y in data:
	x = kernel.X(y)
	if xa <= x and x < xb:
	result += 1.0/kernel.Efficiency(y)
	return result/(xb-xa)

	class ARU(object):
	"""
	Class with a simple interface for unfolding data with ARU.

	Parameters
	----------
	kernel: VKernel
	Instance of the detector kernel, derived from Aru::VKernel.
	xknots: std.vector("double")
	STL vector of doubles containing the knot positions.

	Author
	------
	Hans Dembinski <hans.dembinski@kit.edu>
	"""

	def __init__(self, kernel, xknots):
	from aru import SplineFunction1D,FoldedSplineFunction1D,KnotVector
	self.kernel = kernel
	self.unfolded = SplineFunction1D(KnotVector(ToVector(xknots)))
	self.folded = FoldedSplineFunction1D(kernel, self.unfolded)
	self.referenceCoefs = None
	self.weight = None

	def Unfold(self, data, nEventsForBinning=100):
	"""
	Perform the unfolding.

	Parameters
	----------
	data: std.vector("double")
	STL vector of doubles containing the observations.

	Returns
	-------
	coefs: std.vector("double")
	STL vector of doubles containing the spline coefficients (solution).
	cov: TMatrixDSym
	Covariance matrix of the spline coefficients (solution).
	"""

	from aru import ObjectiveFunction, Vector, Matrix
	nBasis = self.unfolded.GetSize()

	print "creating objective function..."
	objective = ObjectiveFunction(data, self.kernel, self.folded, self.unfolded, 4, nEventsForBinning)

	starts = objective.GetStartValues()

	print "fitting..."
	coefs = MinimumMeanSquaredErrorFit(objective, starts)
	self.refCoefs = Vector(objective.GetReferenceCoefficients())
	self.weight = objective.GetWeight()
	cov = Matrix(nBasis)
	cov2 = Matrix(nBasis)
	objective.CovarianceStat(cov, coefs)
	objective.CovarianceReg(cov2, coefs)
	cov += cov2
	self.objective = objective
	return coefs, cov

	def GetUnfoldedSpline(self):
	return self.unfolded

	def GetFoldedSpline(self):
	return self.folded

	def GetReferenceCoefficients(self):
	if self.refCoefs is None:
	raise StandardError("no reference computed yet")
	return self.refCoefs

	def GetWeight(self):
	if self.weight is None:
	raise StandardError("no weight computed yet")
	return self.weight

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:36 AM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4983094
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions