-
Notifications
You must be signed in to change notification settings - Fork 3
/
calc_PESC_noise.py
77 lines (66 loc) · 2.86 KB
/
calc_PESC_noise.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
# -*- 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\\'
fileheader = 'noise_data_8000'#4157
npy='.npy'
num_orbits = int(1000)
#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=[8000]
#length=[500,700,900,1000,1200,1500,1800,2000,2500,3000]
#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
datafile=np.random.uniform(1,10,size=8001)
arr, nperms = PE_dist(datafile[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)