-
Notifications
You must be signed in to change notification settings - Fork 3
/
calc_PESC_galaxy_reclengthloop.py
121 lines (102 loc) · 4.49 KB
/
calc_PESC_galaxy_reclengthloop.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
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 17 21:00:45 2017
@author: dschaffner
"""
import numpy as np
import matplotlib.pylab as plt
from loadnpyfile import loadnpyfile
from calc_PE_SC import PE, CH, PE_dist, PE_calc_only
#import Cmaxmin as cpl
from collections import Counter
from math import factorial
import time
#calc_PESC_fluid.py
datadir = 'C:\\Users\\dschaffner\\OneDrive - brynmawr.edu\\Galatic Dynamics Data\\GalpyData_July2018\\'
datadir = 'C:\\Users\\dschaffner\\Dropbox\\From OneDrive\\Galatic Dynamics Data\\GalpyData_July2018\\'
#fileheader = 'IDdatabase_Type_1_4co_data_3000'#1621
#fileheader = 'IDdatabase_Type_2_4co_data_3000'#11370
#fileheader = 'IDdatabase_Type_31_4co_data_3000'#2212
#fileheader = 'IDdatabase_Type_32_4co_data_3000'#4157
#fileheader = 'IDdatabase_Type_4_4co_data_3000'#5637
#fileheader = 'IDdatabase_Type_1_data_3000' #3227 orbits
#fileheader = 'IDdatabase_Type_2_data_4000' #25387 orbits
#fileheader = 'IDdatabase_Type_31_data_3000' #5770 orbits
#fileheader = 'IDdatabase_Type_32_data_3000'#9798 orbits
#fileheader = 'IDdatabase_Type_4_data_4000' #5818 orbits
#fileheader = 'IDdatabase_Type_1_10co_data_3000' #3226 orbits
#fileheader = 'IDdatabase_Type_2_10co_data_1000' #8972 orbits
#fileheader = 'IDdatabase_Type_31_10co_data_3000' #370 orbits
#fileheader = 'IDdatabase_Type_32_10co_data_8000'#920 orbits
#fileheader = 'IDdatabase_Type_4_10co_data_1000' #11511 orbits
#fileheader = 'IDdatabase_Type_1_7co_data_3000'#1621
#fileheader = 'IDdatabase_Type_2_7co_data_3000'#11370
#fileheader = 'IDdatabase_Type_31_7co_data_3000'#2212
#fileheader = 'IDdatabase_Type_32_7co_data_3000'#4157
#fileheader = 'IDdatabase_Type_4_7co_data_3000'#5637
#fileheader = 'IDdatabase_Type_1_8co_data_3000'#1621
#fileheader = 'IDdatabase_Type_2_8co_data_1000'#11370
#fileheader = 'IDdatabase_Type_31_8co_data_1000'#2212
#fileheader = 'IDdatabase_Type_32_8co_data_8000'#4157
#fileheader = 'IDdatabase_Type_4_8co_data_1000'#5637
fileheader = 'CR6_3t_Rg_Full'
npy='.npy'
datafile = loadnpyfile(datadir+fileheader+npy)
"""
import glob
data=glob.glob(datadir+'*data.npy')
prop=glob.glob(datadir+'*prop.npy')
datafile = loadnpyfile(data[1000])
print datafile.shape
propfile = loadnpyfile(prop[1000])
print propfile.shape
"""
num_orbits = int(datafile.shape[0])
#record length
#length=[500,550,600,650,700,750,800,850,900,950,1000,1050,1100,1150,1200,1250]#,1300,1350,1400,1450,1500,1550,1600,1650,1700,1750,1800,1850,1900,1950,2000]
length=[500,550,600,650,700,750,800,850,900,950,1000,1050,1100,1150,1200,1250,1300,1350,1400,1450,1500,1550,1600,1650,1700,1750,1800,1850,1900,1950,2000,2200,2500,3000]
#length=[500,550,600,650,700,750,800,850,900,950,1000]#,1050,1100,1150,1200,1250,1300,1350,1400]
length=[6000]
length=[3500]
#num_orbits = 3000
###Storage Arrays###
#delta_t = 1.0/(40000.0)
#delays = np.arange(2,250) #248 elements
#taus = delays*delta_t
#freq = 1.0/taus
num_delays = 499#249
PEs = np.zeros([num_delays+1])
SCs = np.zeros([num_delays+1])
delay = 1
embed_delay = 5
nfac = factorial(embed_delay)
start_time = time.time()
for length_loop in length:
print 'On Length ',length_loop
for loop_delay in np.arange(1,num_delays+1):
print 'On Delay ',loop_delay
print("--- %s minutes ---" % np.round((time.time() - start_time)/60.0,4))
permstore_counter = []
permstore_counter = Counter(permstore_counter)
tot_perms = 0
num_orbits_computed = 0
num_orbits_skipped = 0
for shot in np.arange(num_orbits):#(1,120):
if (shot%1000)==0: print 'On Orbit: ',shot
if np.min(datafile[shot,1:length_loop])<0.1:
num_orbits_skipped+=1
continue
arr, nperms = PE_dist(datafile[shot,1:length_loop],5,delay=loop_delay)
permstore_counter = permstore_counter+arr
tot_perms = tot_perms+nperms
num_orbits_computed+=1
if num_orbits_computed == 5000: break
PE_tot,PE_tot_Se = PE_calc_only(permstore_counter,tot_perms)
C = -2.*((PE_tot_Se - 0.5*PE_tot - 0.5*np.log2(nfac))
/((1 + 1./nfac)*np.log2(nfac+1) - 2*np.log2(2*nfac)
+ np.log2(nfac))*(PE_tot/np.log2(nfac)))
PEs[loop_delay]=PE_tot/np.log2(nfac)
SCs[loop_delay]=C
filename='PE_SC_'+fileheader+'_'+str(num_delays)+'_delays_'+str(num_orbits_computed)+'orbits_'+str(length_loop)+'_timesteps.npz'
np.savez(datadir+filename,PEs=PEs,SCs=SCs,num_orbits=num_orbits_computed,num_skipped=num_orbits_skipped)#,delta_t=delta_t,taus=taus,delays=delays,freq=freq)
"""