-
Notifications
You must be signed in to change notification settings - Fork 0
/
OffsetDistribution.py
executable file
·142 lines (98 loc) · 4 KB
/
OffsetDistribution.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
Calculate the distribution of fractions in the centre pixel for a psf of
given size as the psf is moved away from the centre pixel
'''
import sys
from jg.subs import progressbarClass
import argparse
from scipy.integrate import dblquad
import numpy as np
from Config import *
import matplotlib.pyplot as plt
import tables
def Gaussian2D(y, x, fwhm, offset):
sigma = fwhm / 2.35
arg1 = (x-offset[0])**2
arg2 = (y - offset[1])**2
return np.exp(-(arg1 + arg2) / (2. * sigma**2))
class App(object):
"""docstring for App"""
def __init__(self, args):
super(App, self).__init__()
self.args = args
self.fwhm = FWHM
self.N = self.args.niter
self.xRange = [-0.5, 0.5]
self.yRange = [-0.5, 0.5]
self.run()
def Integrate(self, lims, offset):
return dblquad(Gaussian2D, lims[0], lims[1], lambda x: lims[2], lambda x: lims[3], args=(self.fwhm, offset))
def run(self):
# The total integral of a 2d Gaussian to infinity in each direction
total = self.Integrate((-np.Inf, np.Inf, -np.Inf, np.Inf), (0., 0.))[0]
fractions = []
pb = progressbarClass(self.N)
counter = 0
# Main loop
while counter < self.N:
# Pick two random coordinates
x = (self.xRange[1] - self.xRange[0]) * np.random.ranf() + self.xRange[0]
y = (self.yRange[1] - self.yRange[0]) * np.random.ranf() + self.yRange[0]
# Integrate the new offset Gaussian in the centre pixel and take the fraction
Fraction = self.Integrate((-0.5, 0.5, -0.5, 0.5), (x, y))[0] / total
# Add the result to the list
fractions.append(Fraction)
pb.progress(counter+1)
counter += 1
med_val = np.median(fractions)
med_val_err = np.std(fractions)
# Log the fractions
# Makes the plot easier to read
fractions = np.log10(fractions)
# Create the histogram
vals, edges = np.histogram(fractions, bins=50, range=(-1, 0))
# Counting errors
errs = np.sqrt(vals)
# Width of a bin
binWidth = np.diff(edges)[0]
# Normalise the histogram
# I realise that the numpy.histogram function inclues a density parameter which
# achieves the same thing but I need the errors from the raw values
Integral = float(np.sum(vals))
# Divide the values by the integral of the system
normalisedVals = vals / Integral
normalisedErrs = errs / Integral
# Get the bin centres
centres = edges[:-1] + binWidth / 2.
# Print a line at the most probable value
MostProbable = 10 ** centres[normalisedVals==normalisedVals.max()][0]
fig = plt.figure()
ax = fig.add_subplot(111)
# Plot the histogram
ax.plot(10 ** centres, normalisedVals, 'k-', drawstyle='steps-mid')
ax.set_xlabel(r'Fraction')
ax.set_ylabel(r'Probability')
ax.set_xscale('log')
ticks = [0.1, 0.2, 0.5, 1]
ax.set_xticks(ticks)
ax.set_xticklabels(ticks)
ax.axvline(MostProbable, color='k', ls=':', zorder=-10)
ax.axvline(med_val, color='g')
ax.axvline(med_val - med_val_err / 2., color='b')
ax.axvline(med_val + med_val_err / 2., color='b')
with tables.open_file('out.h5', 'w') as outfile:
outfile.create_array('/', 'x', 10 ** centres)
outfile.create_array('/', 'y', normalisedVals)
outfile.root._v_attrs.most_probable = med_val
outfile.root._v_attrs.sd = med_val_err
plt.show()
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument("-d", "--device", help="PGPLOT device",
required=False, type=str, default="1/xs")
parser.add_argument("-N", "--niter", help="Number of iterations (use > 5000)",
required=False, type=int, default=5000)
args = parser.parse_args()
app = App(args)