-
Notifications
You must be signed in to change notification settings - Fork 0
/
distrubition.py
114 lines (90 loc) · 3.82 KB
/
distrubition.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
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 8 10:49:17 2016
@author: zhiyiwu
"""
import matplotlib.pyplot as plt
import numpy as np
from dcpyps import dataset
from scipy.optimize import minimize
def funcexp(params, x):
y = np.zeros(len(x))
for i in range(int(len(params)/2)):
y += x*params[i]**(-1)*np.exp(x/-params[i])*params[i+1]
return y
def pdf(x, tau):
return tau**(-1)*np.exp(x/-tau)
def sumpdf(x, tau_list, area_list):
y = 0
for index in range(len(tau_list)):
y += area_list[index]*pdf(x, tau_list[index])
return y
def minfunc(params, x):
tau_list = params[:len(params)/2]
area_list = params[len(params)/2:]
return -sum(np.log(sumpdf(x, tau_list, area_list)))
def con(params):
area_list = params[int(len(params)/2):]
return sum(area_list < 1) + sum(area_list > 0) + sum(area_list)!=1
def plot(ax, filename, tres, exp, os):
tres /= 1e6
screcord = dataset.SCRecord([filename,], tres=tres)
if os == 'open':
x = np.array(screcord.opint)[np.array(screcord.oppro) != 8]
xlim = 100
else:
x = np.array(screcord.shint)[np.array(screcord.shpro) != 8]
xlim = 1000
x*=1000
n=len(x)
if n <= 300:
nbdec = 5.0
elif n <= 1000:
nbdec = 8.0
elif n <= 3000:
nbdec = 10.0
else:
nbdec = 12.0
bins = np.arange(np.log10(tres*1000), np.log10(xlim) + np.log10(10)/nbdec, np.log10(10)/nbdec)
bins = 10**bins
hist, bin_edges = np.histogram(x, bins = bins)
plothist = np.hstack((0, np.repeat(hist, 2), 0)).astype(float)
plotbin_edges = np.repeat(bin_edges, 2)
plothist = np.sqrt(plothist)
ax.plot(plotbin_edges, plothist, 'k')
res = minimize(minfunc, np.hstack((exp['value'], np.array(exp['area'])/100)), args=(x), constraints = {'type':'eq', 'fun': con})
exp['value'] = res.x[:len(res.x)/2]
exp['area'] = res.x[len(res.x)/2:]
x = 10**np.linspace(np.log10(0.01), np.log10(xlim), 1000)
sum_y=np.zeros(1000)
for i in range(len(exp['value'])):
y = x*exp['value'][i]**(-1)*np.exp(x/-exp['value'][i])*exp['area'][i]*2.30259*np.log10(10)/nbdec*n
ax.plot(x,np.sqrt(y), 'k--')
sum_y+=y
ax.plot(x, np.sqrt(sum_y), 'k')
ax.set_yticklabels([str(i**2) for i in ax.get_yticks()])
ax.set_xscale('log', basex=10)
ax.set_xlim([0.01, xlim])
if os == 'open':
ax.set_xticks([0.01,0.1,1,10,100])
ax.set_xticklabels(['0.01','0.1','1','10','100'], fontsize=11, ha = 'left')
else:
ax.set_xticks([0.01,0.1,1,10,100, 1000])
ax.set_xticklabels(['0.01','0.1','1','10','100', '1000'], fontsize=11, ha = 'left')
fig = plt.figure(figsize=(16,8))
ax1 = fig.add_subplot(2,4,1)
ax2 = fig.add_subplot(2,4,5)
ax3 = fig.add_subplot(2,4,2)
ax4 = fig.add_subplot(2,4,6)
ax5 = fig.add_subplot(2,4,3)
ax6 = fig.add_subplot(2,4,7)
ax7 = fig.add_subplot(2,4,4)
ax8 = fig.add_subplot(2,4,8)
plot(ax1, '/Users/zhiyiwu/D.SCN', 30, {'value': [0.046, 0.235, 2.27], 'area': [49.3, 45.8, 4.9]}, 'open')
plot(ax2, '/Users/zhiyiwu/D.SCN', 30, {'value': [0.0128, 0.706, 11.6, 97.4], 'area': [24.5, 7.6, 35.9,32]}, 'close')
plot(ax3, '/Users/zhiyiwu/F.SCN', 30, {'value': [0.106, 0.410, 3.49], 'area': [42, 51.9, 6.1]}, 'open')
plot(ax4, '/Users/zhiyiwu/F.SCN', 30, {'value': [0.0206, 1.96, 9.11,159.8], 'area': [24.5, 36.5, 37.4, 1.5]}, 'close')
plot(ax5, '/Users/zhiyiwu/I.SCN', 20, {'value': [0.0889, 0.3, 3.01], 'area': [43.4, 53.5, 3.1]}, 'open')
plot(ax6, '/Users/zhiyiwu/I.SCN', 20, {'value': [0.0227, 2.21, 10.3, 111.9], 'area': [11.3, 49.2, 38, 1.5]}, 'close')
plot(ax7, '/Users/zhiyiwu/M.SCN', 20, {'value': [0.0302, 0.208, 2.52], 'area': [17.8, 77, 5.1]}, 'open')
plot(ax8, '/Users/zhiyiwu/M.SCN', 20, {'value': [0.0141, 1.87, 10.1], 'area': [7.6, 78, 14.1]}, 'close')