No OneTemporary
Actions

Size

29 KB

Subscribers

None

View Options

	diff --git a/debinningMath/dataGen.py b/debinningMath/dataGen.py
	index dfce2dc..79857db 100644
	--- a/debinningMath/dataGen.py
	+++ b/debinningMath/dataGen.py
	@@ -1,255 +1,267 @@
	'''
	dataGen.py
	Functions to generate fake test data.
	Copyright (C) 2015 Abram Krislock

	This program 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.

	This program 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, 110., 1.5, 50., 140.
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program. If not, see http://www.gnu.org/licenses/gpl-3.0.txt.
	'''
	from __future__ import division
	import numpy as np
	import bisect
	import debinningMath.orderStatistics as oStat

	class MyRandom(np.random.RandomState):
	def __init__(self, seedIn):
	np.random.RandomState.__init__(self)
	self.setSeed(seedIn)
	def setSeed(self, seedIn):
	if type(seedIn) in (int, long):
	self.seed(abs(seedIn) % (2**32 - 1))
	else:
	self.seed(abs(hash(str(seedIn))) % (2**32 - 1))

	randomGenerator = MyRandom('default')


	class DataGenContainer(object):
	''' Generates a data sample given an input probability distribution function.
	Arguments for initialization:
	pdf: The probability distribution function used to generate data.
	WARNING: There is no default. pdf must be provided as a function
	with call signature 'pdf(x, pars)' (with x as a numpy array).
	Keyword arguments for initialization:
	pars [=None]: A list of parameters to pass to the pdf.
	nPoints [=100]: The number of points to generate in the sample.
	domain [=np.arange(0., 200., 0.5)]: The discretized x values which
	may be generated.
	rng [=dataGen.randomGenerator]: A flat random number generator with
	call signature 'rng.rand(n)'. It must return a numpy array of
	length n, containing flat random numbers in the range (0, 1).
	'''
	def __init__(self, pdf, **kwargs):
	# Set generating probability distribution function
	self.pdf = pdf
	self.pars = kwargs.get('pars', None)

	# Set data generation parameters
	self.n = kwargs.get('nPoints', 100)
	self.x = np.array(kwargs.get('domain', np.arange(0., 200., 0.5)))
	self.rng = kwargs.get('rng', randomGenerator)

	# initial calculations...
	self.f_unnormed = self.pdf(self.x, self.pars)
	self.F_unnormed = oStat.A(self.f_unnormed, oStat.forwardDiff(self.x))
	self.f = self.f_unnormed / self.F_unnormed[-1]
	self.F = self.F_unnormed / self.F_unnormed[-1]

	def changeFunction(self, f_unnormed_in):
	''' Changes self.f_unnormed to a new numerical function, based upon the same
	x domain as the original. '''
	self.f_unnormed = f_unnormed_in
	# reinitialize
	self.F_unnormed = oStat.A(self.f_unnormed, oStat.forwardDiff(self.x))
	self.f = self.f_unnormed / self.F_unnormed[-1]
	self.F = self.F_unnormed / self.F_unnormed[-1]

	def changeDomain(self, x=0):
	''' Changes the x domain to a given numpy array. If x is an integer,
	chooses a new domain for which the pdf is non-zero, with x points
	discarded at both the beginning and end of the domain. '''
	if type(x) == np.ndarray:
	self.x = x
	else:
	nonzero = self.x[self.f.nonzero()]
	self.x = np.linspace(nonzero[x], nonzero[-1-x], len(self.x))
	# reinitialize
	self.f_unnormed = self.pdf(self.x, self.pars)
	self.F_unnormed = oStat.A(self.f_unnormed, oStat.forwardDiff(self.x))
	self.f = self.f_unnormed / self.F_unnormed[-1]
	self.F = self.F_unnormed / self.F_unnormed[-1]

	def generate(self, sort=False):
	# bisect_left(F, p) locates the left-most entry in F which is >= p
	# Use the unnormed cdf, F_unnormed, to avoid tiny number roundoff error.
	self.data = np.array([self.x[bisect.bisect_left(self.F_unnormed, p)]
	for p in (self.rng.rand(self.n) * self.F_unnormed[-1])])
	if sort:
	self.data.sort()
	return self.data


	class DataGen2dContainer(object):
	''' Generates a data sample given an input probability distribution function.
	Arguments for initialization:
	pdf: The probability distribution function used to generate data.
	WARNING: There is no default. pdf must be provided as a function
	with call signature 'pdf(x, y, pars)' (with x,y as numpy arrays).
	Keyword arguments for initialization:
	pars [=None]: A list of parameters to pass to the pdf.
	nPoints [=100]: The number of points to generate in the sample.
	domain [=( np.arange(0., 200., 0.5), np.arange(0., 200., 0.5) )]:
	The domain of discretized (x,y) values which may be generated.
	rng [=dataGen.randomGenerator]: A flat random number generator with
	call signature 'rng.rand(n)'. It must return a numpy array of
	length n, containing flat random numbers in the range (0, 1).
	'''
	def __init__(self, pdf, **kwargs):
	# Set generating probability distribution function
	self.pdf = pdf
	self.pars = kwargs.get('pars', None)

	# Set data generation parameters
	self.n = kwargs.get('nPoints', 100)
	self.x, self.y = kwargs.get('domain', (np.arange(0., 200., 0.5),
	np.arange(0., 200., 0.5)))
	self.x = np.array(self.x)
	self.y = np.array(self.y)
	self.dx = self.x[1] - self.x[0]
	self.dy = self.y[1] - self.y[0]
	self.xx, self.yy = np.meshgrid(self.x, self.y, indexing='ij')
	self.rng = kwargs.get('rng', randomGenerator)

	# initial calculations...
	self.f_unnormed = self.pdf(self.x, self.y, self.pars)
	self.F_Hunnormed = oStat.A_H(self.f_unnormed, self.dx, self.dy)
	self.F_Vunnormed = oStat.A_V(self.f_unnormed, self.dx, self.dy)
	self.f_H = self.f_unnormed / self.F_Hunnormed[-1][-1]
	self.f_V = self.f_unnormed / self.F_Vunnormed[-1][-1]
	self.F_H = self.F_Hunnormed / self.F_Hunnormed[-1][-1]
	self.F_V = self.F_Vunnormed / self.F_Vunnormed[-1][-1]

	def changeFunction(self, f_unnormed_in):
	''' Changes self.f_unnormed to a new numerical function, based upon the same
	x domain as the original. '''
	self.f_unnormed = f_unnormed_in
	# reinitialize
	self.F_Hunnormed = oStat.A_H(self.f_unnormed, self.dx, self.dy)
	self.F_Vunnormed = oStat.A_V(self.f_unnormed, self.dx, self.dy)
	self.f_H = self.f_unnormed / self.F_Hunnormed[-1][-1]
	self.f_V = self.f_unnormed / self.F_Vunnormed[-1][-1]
	self.F_H = self.F_Hunnormed / self.F_Hunnormed[-1][-1]
	self.F_V = self.F_Vunnormed / self.F_Vunnormed[-1][-1]

	- def generate(self):
	- ''' Generate a data sample from the stored PDF. '''
	+ def generate(self, sort=False):
	+ ''' Generate a data sample from the stored PDF. If sort=True, argsortH and
	+ argsortV are both automatically called. '''
	# bisect_left(F, p) locates the left-most entry in F which is >= p
	# Use the unnormed cdf, F_Vunnormed, to avoid tiny number roundoff error.
	self.data = []
	for p in (self.rng.rand(self.n) * self.F_Vunnormed[-1][-1]):
	i = bisect.bisect_left(self.F_Vunnormed.T[-1], p)
	j = bisect.bisect_left(self.F_Vunnormed[i], p)
	self.data.append((self.x[i], self.y[j]))
	- self.data = np.array(self.data)
	+ self.data = np.array(self.data, dtype=[('x', float), ('y', float)])
	+ if 'r_Harray' in dir(self):
	+ del self.r_Harray
	+ if 'r_Varray' in dir(self):
	+ del self.r_Varray
	+ if sort:
	+ self.argsortH()
	+ self.argsortV()
	return self.data

	- def sortH(self):
	- ''' Return a "horizontally" sorted copy of the generated data. '''
	+ def argsortH(self):
	+ ''' Return an array of indices which would "horizontally" sort the data. '''
	if 'data' not in dir(self):
	self.generate()
	+ self.r_Harray = self.data.argsort(order=['y','x']).argsort()
	+ return self.r_Harray

	- def sortV(self):
	- ''' Return a "vertically" sorted copy of the generated data. '''
	+ def argsortV(self):
	+ ''' Return an array of indices which would "vertically" sort the data. '''
	if 'data' not in dir(self):
	self.generate()
	+ self.r_Varray = self.data.argsort(order=['x','y']).argsort()
	+ return self.r_Varray


	class DataPlusBGContainer(DataGenContainer):
	''' Generates a data sample given signal and background probability
	distribution functions, using the formula:
	pdf(x) = (1 - alpha) * b_pdf + alpha * s_pdf
	Arguments for initialization:
	s_pdf: The "signal" probability distribution function.
	b_pdf: The "background" probability distribution function.
	alpha: The signal-background mixing.
	Keyword arguments for initialization:
	pars [=None]: A list of parameters to pass to the pdf.
	nPoints [=100]: The number of points to generate in the sample.
	domain [=np.arange(0., 200., 0.5)]: The discretized x values which
	may be generated.
	rng [=dataGen.randomGenerator]: A flat random number generator with
	call signature 'rng.rand(n)'. It must return a numpy array of
	length n, containing flat random numbers in the range (0, 1).
	'''
	def __init__(self, s_pdf, b_pdf, alpha, **kwargs):
	# Set generating probability distribution function
	self.s_pdf = s_pdf
	self.b_pdf = b_pdf
	if alpha < 0 or alpha > 1:
	raise ValueError('alpha should be between 0 and 1')
	self.alpha = alpha
	self.pars = kwargs.get('pars', None)

	# Set data generation parameters
	self.n = kwargs.get('nPoints', 100)
	self.x = np.array(kwargs.get('domain', np.arange(0., 200., 0.5)))
	if not type(self.x) == np.ndarray:
	raise TypeError('Use a list or array to specify the x domain')
	self.rng = kwargs.get('rng', randomGenerator)

	# initial calculations...
	# signal
	self.fs_unnormed = self.s_pdf(self.x, self.pars)
	self.Fs_unnormed = oStat.A(self.fs_unnormed, oStat.forwardDiff(self.x))
	self.fs = self.fs_unnormed / self.Fs_unnormed[-1]

	# background
	self.fb_unnormed = self.b_pdf(self.x, self.pars)
	self.Fb_unnormed = oStat.A(self.fb_unnormed, oStat.forwardDiff(self.x))
	self.fb = self.fb_unnormed / self.Fb_unnormed[-1]

	# combined
	self.f_unnormed = (1. - self.alpha) * self.fb + self.alpha * self.fs
	self.F_unnormed = oStat.A(self.f_unnormed, oStat.forwardDiff(self.x))

	# normed
	self.f = self.f_unnormed / self.F_unnormed[-1]
	self.F = self.F_unnormed / self.F_unnormed[-1]

	def changeDomain(self, x=0):
	''' Changes the x domain to a given numpy array. If x is an integer,
	chooses a new domain for which the pdf is non-zero, with x points
	discarded at both the beginning and end of the domain. '''
	if type(x) == np.ndarray:
	self.x = x
	else:
	nonzero = self.x[self.f.nonzero()]
	self.x = np.linspace(nonzero[x], nonzero[-1-x], len(self.x))
	# reinitialize
	# signal
	self.fs_unnormed = self.s_pdf(self.x, self.pars)
	self.Fs_unnormed = oStat.A(self.fs_unnormed, oStat.forwardDiff(self.x))
	self.fs = self.fs_unnormed / self.Fs_unnormed[-1]

	# background
	self.fb_unnormed = self.b_pdf(self.x, self.pars)
	self.Fb_unnormed = oStat.A(self.fb_unnormed, oStat.forwardDiff(self.x))
	self.fb = self.fb_unnormed / self.Fb_unnormed[-1]

	# combined
	self.f_unnormed = (1. - self.alpha) * self.fb + self.alpha * self.fs
	self.F_unnormed = oStat.A(self.f_unnormed, oStat.forwardDiff(self.x))

	# normed
	self.f = self.f_unnormed / self.F_unnormed[-1]
	self.F = self.F_unnormed / self.F_unnormed[-1]
	diff --git a/debinningMath/plotting.py b/debinningMath/plotting.py
	index 4ce9ebb..aaf41a2 100644
	--- a/debinningMath/plotting.py
	+++ b/debinningMath/plotting.py
	@@ -1,247 +1,256 @@
	'''
	plotting.py
	Handy utilities for debinning plots with matplotlib.
	Copyright (C) 2015 Abram Krislock

	This program 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.

	This program 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 this program. If not, see http://www.gnu.org/licenses/gpl-3.0.txt.
	'''
	import numpy as np
	import bisect
	from matplotlib import pyplot as plt
	import debinningMath.orderStatistics as oStat
	from debinningCore import settings


	class HistoContainer(object):
	''' A histogram as a container.
	Arguments:
	binEdges: A numpy array containing the locations of the bin edges.
	data: Initial data to fill the histogram. Data values may also be
	added one at a time with addValue(value).
	normed: Whether or not to normalize the histogram when plotting.
	'''
	def __init__(self, binEdges, data=None, normed=False):
	self.binEdges = binEdges
	self.binContents = np.zeros(len(binEdges)-1)
	self.area = None
	self.maxi = None
	if type(data) in [list, tuple, np.ndarray]:
	for value in data:
	self.addValue(value)
	self.normed = normed

	def __len__(self):
	return len(self.binContents)

	def __call__(self, value):
	self.addValue(value)

	def integral(self):
	''' Returns the total integration over the entire histogram. '''
	if self.area == None:
	self.area = np.sum(self.binContents * oStat.forwardDiff(self.binEdges))
	return self.area

	def maximum(self):
	''' Returns the content of the histogram bin with maximum content. '''
	if self.maxi == None:
	self.maxi = np.max(self.binContents)
	if self.maxi > 0:
	return self.maxi
	else:
	return np.infty

	def normalize(self, norm=np.infty):
	''' Normalize the histogram contents. '''
	if norm == 'fwhm':
	self.binContents /= (self.maximum() / 2)
	elif norm < np.infty:
	self.binContents /= norm
	else:
	self.binContents /= self.integral()
	# demand that integral() and maximum() must recalculate
	self.maxi = None
	self.area = None

	def fwhm(self):
	''' Obtain the Full Width at Half Maximum for this histogram. '''
	self.normalize('fwhm')
	if self.maximum() == np.infty:
	return (-np.infty, np.infty)
	lower = np.infty
	upper = -np.infty
	for bindex in xrange(len(self)):
	if self.binContents[bindex] >= 1.:
	if self.binEdges[bindex+1] > upper:
	upper = self.binEdges[bindex+1]
	if self.binEdges[bindex] < lower:
	lower = self.binEdges[bindex]
	return (lower, upper)

	def plot(self, **kwargs):
	''' Plots the result, filling under the 'histogram'.
	Passes **kwargs on to the matplotlib.pyplot.fill_between function.
	Returns a proxy matplotlib.patches.Rectangle for legend control.
	'''
	if self.normed:
	self.normalize()
	x = np.sort(np.concatenate((self.binEdges, self.binEdges)))
	y = np.array([[value, value] for value in self.binContents / self.integral()])
	y = np.concatenate(([1e-6], y.flatten(), [1e-6]))
	plt.fill_between(x, y, 1e-6, **kwargs)
	return Rectangle((0, 0), 1, 1, fill=True, **kwargs)

	def addValue(self, value):
	''' Insert a value, or set of values, into the histogram. '''
	if type(value) in [list, tuple, np.ndarray]:
	for v in value:
	self.addValue(v)
	else:
	# demand that integral() and maximum() must recalculate
	self.area = None
	self.maxi = None
	# bisect_left(x, datum) locates the left-most entry in x which is >= datum
	index = bisect.bisect_left(self.binEdges, value) - 1
	if index == -1:
	if settings.debug:
	raise ValueError('HistoContainer warning: Underflow')
	elif index >= len(self.binContents):
	if settings.debug:
	raise ValueError('HistoContainer warning: Overflow')
	else:
	self.binContents[index] += 1


	class Histo2dContainer(object):
	''' A 2d histogram as a container.
	Arguments:
	binEdgesX: A numpy array containing the x-values of the bin edges.
	binEdgesY: A numpy array containing the y-values of the bin edges.
	normed: Whether or not to normalize the histogram when plotting.
	Data may be added one ordered pair at a time with "addValue((x,y))".
	'''
	def __init__(self, binEdgesX, binEdgesY, normed=False):
	self.binEdgesX = binEdgesX
	self.binEdgesY = binEdgesY
	self.binContents = np.zeros((len(binEdgesX)-1, len(binEdgesY)-1))
	self.volume = None
	self.maxi = None
	self.normed = normed

	def __len__(self):
	return len(self.binContents)

	def __call__(self, orderedPair):
	self.addValue(orderedPair)

	def integral(self):
	''' Returns the total integration over the entire histogram. '''
	if self.volume == None:
	dx, dy = np.meshgrid(np.diff(self.binEdgesX), np.diff(self.binEdgesY),
	indexing='ij')
	self.volume = np.sum(self.binContents * dx * dy)
	return self.volume

	def maximum(self):
	''' Returns the content of the histogram bin with maximum content. '''
	if self.maxi == None:
	self.maxi = np.max(self.binContents)
	if self.maxi > 0:
	return self.maxi
	else:
	return np.infty

	def normalize(self, norm=np.infty):
	''' Normalize the histogram contents. '''
	if norm == 'fwhm':
	self.binContents /= (self.maximum() / 2)
	elif norm < np.infty:
	self.binContents /= norm
	else:
	self.binContents /= self.integral()
	# demand that integral() and maximum() must recalculate
	self.maxi = None
	self.volume = None

	def plot(self, axis, **kwargs):
	''' Plots the result, filling under the 'histogram'.
	Arguments:
	axis: An existing matplotlib 3d axis object on which to plot.
	Passes **kwargs on to the matplotlib.pyplot.plot_surface function.
	Returns a proxy matplotlib.patches.Rectangle for legend control.
	'''
	if self.normed:
	self.normalize()
	+ cmap = kwargs.pop('cmap', None)
	+ if cmap != None:
	+ plotmax = np.max(self.binContents)
	+ plotmin = np.min(self.binContents)
	for i in np.arange(len(self.binEdgesX)-1):
	for j in np.arange(len(self.binEdgesY)-1):
	x = np.linspace(self.binEdgesX[i], self.binEdgesX[i+1], 10)
	y = np.linspace(self.binEdgesY[j], self.binEdgesY[j+1], 10)
	xx, yy = np.meshgrid(x, y, indexing='ij')
	zz = np.ones_like(xx) * self.binContents[i,j]
	- axis.plot_surface(xx, yy, zz, **kwargs)
	+ if cmap != None:
	+ kwargs.pop('color', None)
	+ colorfloat = (self.binContents[i,j] - plotmin) / (plotmax - plotmin)
	+ axis.plot_surface(xx, yy, zz, color=cmap(colorfloat), **kwargs)
	+ else:
	+ axis.plot_surface(xx, yy, zz, **kwargs)

	def addValue(self, orderedPair):
	''' Insert an ordered pair, (x,y), into the histogram. '''
	# demand that integral() and maximum() must recalculate
	self.volume = None
	self.maxi = None
	# bisect_left(x, datum) locates the left-most entry in x which is >= datum
	i = bisect.bisect_left(self.binEdgesX, orderedPair[0]) - 1
	j = bisect.bisect_left(self.binEdgesY, orderedPair[1]) - 1
	if i == -1:
	if settings.debug:
	raise ValueError('HistoContainer warning: X Underflow')
	if j == -1:
	if settings.debug:
	raise ValueError('HistoContainer warning: Y Underflow')
	if i >= len(self.binEdgesX)-1:
	if settings.debug:
	raise ValueError('HistoContainer warning: X Overflow')
	if j >= len(self.binEdgesY)-1:
	if settings.debug:
	raise ValueError('HistoContainer warning: Y Overflow')
	if all([i >= 0, j >= 0, i < len(self.binEdgesX)-1, i < len(self.binEdgesY)-1]):
	self.binContents[i,j] += 1



	def verticalColorBand(histo, x, dx, colors=None):
	''' Plot the result of a debinning Monte Carlo as colored uncertainty bands:
	This function adds one histogram, painted as a vertical band of color.
	The band is a rectangle of width [x-0.5dx, x+0.5dx].
	'''
	if not type(histo) == HistoContainer:
	raise TypeError('Please provide a HistoContainer instance for "histo"')
	if colors == None:
	colors = np.array([(1.,1.,0.), (1.,0.,0.), (1.,1.,1.), (0.,0.,1.)])

	# Special normalization to show FWHM transition.
	histo.normalize('fwhm')

	# Horizontal width of the band.
	base = np.array([x-dx/2, x+dx/2])

	for bindex in xrange(len(histo)):
	upper = np.array([histo.binEdges[bindex+1], histo.binEdges[bindex+1]])
	lower = np.array([histo.binEdges[bindex], histo.binEdges[bindex]])
	value = histo.binContents[bindex]
	if value >= 1.:
	color = tuple((value - 1.) * colors[0] + (2. - value) * colors[1])
	alp = 1.
	else:
	color = tuple((value) * colors[2] + (1. - value) * colors[3])
	alp = 0.35
	if value > 0.1:
	plt.fill_between(base, upper, lower,
	edgecolor='none', facecolor=color, alpha=alp, zorder=2)
	diff --git a/examples/nfrDemo.py b/examples/nfrDemo.py
	index d4abc1a..ba66997 100644
	--- a/examples/nfrDemo.py
	+++ b/examples/nfrDemo.py
	@@ -1,89 +1,89 @@
	'''
	nfrDemo.py
	Demonstration of the "rth of n" probability distribution.

	Copyright (C) 2015 Abram Krislock

	This program 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.

	This program 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 this program. If not, see http://www.gnu.org/licenses/gpl-3.0.txt.
	'''

	### Imports
	from __future__ import division
	import numpy as np

	from debinningCore import settings
	import debinningMath.dataGen as dataGen
	import debinningMath.sampleFunctions as funcs
	import debinningMath.orderStatistics as oStat

	from matplotlib import pyplot as plt
	###

	nPulls = settings.nfrDemoPulls
	gen = dataGen.DataGenContainer(funcs.endpointPlusBG, nPoints=nPulls)
	fit = dataGen.DataGenContainer(funcs.endpointBG, nPoints=nPulls)
	bgf_plot = fit.f_unnormed / gen.F_unnormed[-1]

	pullData = {index:[] for index in xrange(nPulls)}
	for iteration in xrange(2000):
	for i,v in enumerate(gen.generate(True)):
	pullData[i].append(v)

	ithPDF = {}
	for i in xrange(1, nPulls+1):
	- ithPDF[i-1] = np.array([gen.f[j] * gen.F[j]*(i-1) (1 - gen.F[j])**(nPulls-i)
	- * oStat.nKr(nPulls, i) for j in xrange(len(gen.x))])
	+ ithPDF[i-1] = (gen.f * gen.F*(i-1) (1 - gen.F)**(nPulls-i)
	+ * oStat.nKr(nPulls, i))

	def nfrPlot(i=0):
	fig, axes = plt.subplots(figsize=(9,5), facecolor='#ffffff')
	plot_f, = plt.plot(gen.x, gen.f, '--', color='#0088ff', alpha=0.6,
	label=r's+bg PDF', linewidth=3)
	plot_bgf, = plt.plot(gen.x, bgf_plot, '-.', color='#0088ff', alpha=0.6,
	label=r'bg PDF', linewidth=3)
	plot_hist = ( plt.hist(pullData[i], bins=np.arange(0, 201, 5), color='#66a61e',
	alpha=0.4, normed=True,
	label=r'$%i$ of $%i$ data' % (i+1, nPulls)) )[2][0]
	plot_nfr, = plt.plot(gen.x, ithPDF[i], '#964a01', linewidth=3,
	label=r'${}_{%i} f_{%i}$' % (nPulls, i+1))

	axes.set_xticks(np.arange(0, 201, 50))
	axes.set_xticks(np.arange(0, 201, 10), minor=True)
	axes.set_xticklabels([r'$%i$' % int(num) for num in np.arange(0, 201, 50)],
	fontsize=20)

	axes.set_ylim(bottom=-0.0001, top=0.0401)
	axes.set_yticks(np.arange(0, 0.04, 0.01))
	axes.set_yticklabels([r'$%0.3f$' % num for num in np.arange(0, 0.04, 0.01)],
	fontsize=20)
	axes.set_yticks(np.arange(0, 0.04, 0.01/5), minor=True)

	axes.set_xlabel('Physical Observable', labelpad=5, fontsize=22)
	axes.set_ylabel('Probability', labelpad=5, fontsize=22)

	axes.tick_params(length=10, width=1.2, labelsize=22, zorder=10, pad=10)
	axes.tick_params(which='minor', length=5, width=1.2, zorder=11)
	axes.minorticks_on()

	plt.legend((plot_hist, plot_nfr, plot_f, plot_bgf),
	(plot_hist.get_label(), plot_nfr.get_label(), plot_f.get_label(),
	plot_bgf.get_label()),
	shadow=True, loc='best', prop={'size':26},
	labelspacing=0.25, handlelength=1.2, handletextpad=0.35)

	if __name__ == '__main__':
	print 'Now running the "Being rth of n" demonstration...'
	plt.ion()
	for i in xrange(nPulls):
	nfrPlot(i)
	print 'Done.\n\n'
	diff --git a/examples/nfrVrHDemo.py b/examples/nfrVrHDemo.py
	new file mode 100644
	index 0000000..2d55baf
	--- /dev/null
	+++ b/examples/nfrVrHDemo.py
	@@ -0,0 +1,70 @@
	+#!/usr/bin/python
	+'''
	+nfrVrHDemo.py
	+Demonstration of the 2d "r_V and r_H of n" probability distribution.
	+
	+Copyright (C) 2015 Abram Krislock
	+
	+This program 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.
	+
	+This program 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 this program. If not, see http://www.gnu.org/licenses/gpl-3.0.txt.
	+'''
	+
	+### Imports
	+from __future__ import division
	+import numpy as np
	+from mpl_toolkits.mplot3d import axes3d
	+from matplotlib import pyplot as plt
	+from matplotlib import cm
	+from time import time
	+import argparse
	+
	+from debinningCore import settings
	+import debinningMath.dataGen as dataGen
	+import debinningMath.sampleFunctions as funcs
	+import debinningMath.orderStatistics as oStat
	+from debinningMath.plotting import Histo2dContainer as Hist2d
	+###
	+
	+def nfrVrHPlot(r_V, r_H):
	+ t = time()
	+ nPulls = settings.nfrDemoPulls
	+ gen = dataGen.DataGen2dContainer(funcs.footballBumpBG, nPoints=nPulls,
	+ domain=(np.arange(0, 200.1, 1), np.arange(0, 200.1, 1)))
	+
	+ hist = Hist2d(np.arange(0, 200.1, 10.), np.arange(0, 200.1, 10.), True)
	+ for iteration in xrange(2000000):
	+ for i,v in enumerate(gen.generate(True)):
	+ if gen.r_Varray[i] == r_V-1 and gen.r_Harray[i] == r_H-1:
	+ hist.addValue(v)
	+
	+ nfrVrH = (oStat.nKr(nPulls, r_V) * oStat.nKr(nPulls, r_H) * gen.F_V**(r_V-1)
	+ * gen.F_H*(r_H-1) (1 - gen.F_V)**(nPulls-r_V)
	+ * (1 - gen.F_H)*(nPulls-r_H) gen.f_V)
	+ print 'Calculations finished after', (time() - t), 'seconds.\n'
	+
	+ fig = plt.figure(figsize=(22,14))
	+ ax = fig.add_subplot(111, projection='3d')
	+ ax.plot_wireframe(gen.xx, gen.yy, nfrVrH, rstride=5, cstride=5, cmap=cm.jet)
	+ ax.set_xlabel('X')
	+ ax.set_ylabel('Y')
	+ ax.set_zlabel('Z')
	+ hist.plot(ax, cmap=cm.jet, alpha=0.4)
	+ plt.show()
	+
	+if __name__ == '__main__':
	+ parser = argparse.ArgumentParser(description='Being r_V and r_H of n')
	+ parser.add_argument('r_V', help='r_V is an integer between 1 and n', type=int)
	+ parser.add_argument('r_H', help='r_H is an integer between 1 and n', type=int)
	+ args = parser.parse_args()
	+ nfrVrHPlot(args.r_V, args.r_H)
	+ print 'Done.\n\n'
	diff --git a/examples/test2d.py b/examples/test2d.py
	deleted file mode 100644
	index b2103ae..0000000
	--- a/examples/test2d.py
	+++ /dev/null
	@@ -1,33 +0,0 @@
	-from debinningMath import dataGen
	-from debinningMath import sampleFunctions
	-from debinningMath import plotting
	-
	-from mpl_toolkits.mplot3d import axes3d
	-import matplotlib.pyplot as plt
	-import numpy as np
	-
	-gen = dataGen.DataGen2dContainer( sampleFunctions.footballBumpBG, nPoints=10000,
	- domain=(np.arange(0, 200.1, 1),
	- np.arange(0, 200.1, 1)) )
	-
	-fig = plt.figure(figsize=(22,14))
	-ax = fig.add_subplot(111, projection='3d')
	-ax.plot_wireframe(gen.xx, gen.yy, gen.f_V, rstride=5, cstride=5)
	-ax.set_xlabel('X')
	-ax.set_ylabel('Y')
	-ax.set_zlabel('Z')
	-
	-gen.generate()
	-hist = plotting.Histo2dContainer(np.arange(0, 200.1, 10.),
	- np.arange(0, 200.1, 10.), True)
	-for orderedPair in gen.data:
	- hist.addValue(orderedPair)
	-hist.plot(ax, color='r', alpha=0.4, edgecolors='r')
	-
	-plt.show()
	-
	-# YAY! This works. I can verify that data generation is now working as expected
	-# using my current A_H and A_V operators. However, they are slow, so I should
	-# rethink what they actually do, and try to optimize them.
	-#
	-# Once they are optimized for a decent amount of speed, I will test my nfrVrH idea.

File Metadata

Mime Type: text/x-diff
Expires: Thu, Apr 24, 6:31 AM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4881375
Default Alt Text: (29 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions