-
Notifications
You must be signed in to change notification settings - Fork 0
/
statistical_analysis.py
106 lines (92 loc) · 4.11 KB
/
statistical_analysis.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
'''
Classes which help define an Time series or quench as well as a class which defines multiple quenches
Carter Francis
'''
from compute_vornoi import *
from load import load_traj, reduce
import numpy as np
import matplotlib.pyplot as plt
# a Time series object is a combination of the information from a traj and lammps output file. Has information about
# the vornoi polyhedron, temperature, volume, density...
class Timeseries:
def __init__(self, ts, temp, aa, astd, av, vstd, vi, vf,aas, astds, vas, vstds,den):
self.timestep = ts
self.sampling = 100
self.temp = temp
self.surface_area = aa
self.area_std = astd
self.average_volume = av
self.volume_std = vstd
self.density = den
# determining the unique vornoi indices and giving them an index...
self.indexes, self.positions = np.unique(np.reshape(vi, (10*self.sampling, 15)),
axis=0, return_inverse=True)
self.all_average_area, self.all_area_std = self.get_areas(aas, astds)
self.all_average_volume, self.all_volume_std = self.get_volumes(vas, vstds)
self.all_top_ten_freq = self.get_top(vf)
# note this is kind of weird implementation because it is unwrapped to help with the calculations.
# it just helps by giving all the possible Vornoi indices an index... (I should maybe do that myself.
def get_index_stats(self, property):
unwrapped = np.reshape(property, (10 * self.sampling)) # unwrap
resampled = []
for index in range(1, len(self.indexes)):
where = np.where(index == self.positions)[0] # once again the pythons weird arrays strike again...
is_in = np.floor_divide(where, 10)
index_prop = []
count = 0
for i in range(0, self.sampling):
if i in is_in:
index_prop.append(float(unwrapped[where[count]]))
count += 1
else:
index_prop.append(0)
resampled.append(index_prop)
return resampled
def get_top(self,freq,plot=False):
resampled_freq =self.get_index_stats(freq)
print(np.sum(resampled_freq, axis =1))
return resampled_freq
def get_areas(self, aa, astd, plot=False):
area = self.get_index_stats(aa)
astd = self.get_index_stats(astd)
if plot:
for a in area:
plt.scatter(self.temp, a)
plt.show()
return area, astd
def get_volumes(self,vas, vstds, plot=False):
volumes = self.get_index_stats(vas)
vstd = self.get_index_stats(vstds)
if plot:
for v in volumes:
plt.scatter(self.temp, v)
plt.show()
return volumes, vstd
def plot_top(self):
for r in self.all_top_ten_freq:
plt.scatter(self.temp,r)
plt.show()
def get_info_from_index(self, index):
# python is dumb as hell sometimes
print(self.indexes)
pos = [num for num,ind in enumerate(self.indexes) if all(ind == index )][0]
return self.all_average_area[pos], self.all_area_std[pos], self.all_average_volume[pos],\
self.all_volume_std[pos],self.all_top_ten_freq[pos],self.temp
# a Sample is just a list of time series with added information about quench rates.. It might be overkill to have this
# as a separate class but it works for now.
class Sample:
def __init__(self,timeseries, quench_rates):
self.timeseries = timeseries # list of timeseries objects
self.qenchrates = quench_rates
def compare_index(self, index):
avg_area,area_std, avg_vol, vol_std, freq, temp = [],[],[],[],[],[]
for quench in self.timeseries:
holder = quench.get_info_from_index(index)
avg_area.append(holder[0])
area_std.append(holder[1])
avg_vol.append(holder[2])
vol_std.append(holder[3])
freq.append(holder[4])
temp.append(holder[5])
return avg_area, area_std, avg_vol, vol_std, freq, temp