python-include/ACF.py at main · k50112113/python-include · GitHub

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
import numpy as np
import os
import random
import copy

def ACF_DIRECT(Data,time,dt_index_range,dt_index_max):
    Data = np.array(Data)
    ACF = []
    Intg_ACF = []

    dt_list = time[:dt_index_max]
    dim = len(Data[0])

    quantity = 0.0
    cntt = 0
    for start_index in range(0,len(Data),dt_index_range):
        print("t0 = %f"%(time[start_index]))
        now_lim_index = start_index+dt_index_range
        if now_lim_index < len(Data):
            Cor = np.zeros((len(dt_list)))
            cntt += 1
            for dt_index in range(len(dt_list)):
                tmp = np.zeros(dim)
                cnt = 0
                t_index = start_index
                while t_index+dt_index < now_lim_index:
                    tmp += Data[t_index]*Data[t_index+dt_index]
                    cnt += 1
                    t_index += 1
                Cor[dt_index] = np.mean(tmp/float(cnt))
            sum = 0
            ACF.append(Cor)
            Intg_ACF.append([])
            for dt_index in range(len(dt_list)):
                sum += Cor[dt_index]
                quantity = (sum-Cor[0]/2.0-Cor[dt_index]/2.0)
                Intg_ACF[len(Intg_ACF)-1].append(quantity)
        else:
            break
    print("Number of profiles = %d"%(cntt))

    return dt_list,ACF,Intg_ACF


def ACF_FFT(Data,time,dt_index_range,dt_index_max): #See Tildesley Page 280
    #Data: N x m
    ACF = []
    Intg_ACF = []
    dt_list = time[:dt_index_max]
    dim = Data.shape[1]
    quantity = 0.0
    cntt = 0
    for start_index in range(0,Data.shape[0],dt_index_range):
        #print("t0 = %f"%(time[start_index]))
        now_lim_index = start_index+dt_index_range
        if now_lim_index < Data.shape[0]:
            Cor = np.zeros(dt_list.shape[0])
            cntt += 1
            for i in range(dim):
                tmpData = Data[start_index:now_lim_index][:,i]
                l = tmpData.shape[0]
                tmpfft = np.fft.fft(np.concatenate((tmpData,np.zeros(l)),axis=0))
                c = np.fft.ifft(tmpfft*np.conjugate(tmpfft)).real
                corr = c[:len(c)//2]
                corr /= np.array([l-x for x in range(l)])
                Cor += corr[:dt_index_max]
            Cor /= dim
            intg = 0
            ACF.append(Cor)
            Intg_ACF.append([])
            for dt_index in range(len(dt_list)):
                intg += Cor[dt_index]
                quantity = (intg-Cor[0]/2.0-Cor[dt_index]/2.0)
                Intg_ACF[len(Intg_ACF)-1].append(quantity)

    print("Number of profiles = %d"%(cntt))
    return np.array(dt_list),np.array(ACF),np.array(Intg_ACF)