-
Notifications
You must be signed in to change notification settings - Fork 0
/
SFA_bias.py
73 lines (49 loc) · 1.18 KB
/
SFA_bias.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
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy as sp
import pandas as pd
import seaborn as sns
from bisect import bisect_left
import math
def closest_idx(array, value):
"""
Find index of the element of 'array' which is closest to 'value'.
"""
if type(array) is np.ndarray:
pos = bisect_left(array, value)
else:
pos = bisect_left(array.numpy(), value)
if pos == 0:
return 0
elif pos == len(array):
return -1
else:
return pos-1
X = np.loadtxt('')
X = (X - np.min(X)) / (np.max(X) - np.min(X))
T = len(X)
bias = X[:, -2]
t = X[:, 0]
temp = 300.
kb = 0.008314
beta = 1./(kb*temp)
logweight = beta*bias
dt = np.round(t[1]-t[0], 3)
d_tprime = np.copy(np.exp(logweight)*dt)
tprime = np.cumsum(d_tprime)
t_lenth = math.floor(tprime[-1])
x_ = []
for i in range(t_lenth):
idx = closest_idx(tprime, i)
x_.append(X[idx])
X = np.array(x_)
X = (X - np.mean(X, axis=0))
X = X.T
C = np.dot(X, X.T)
dX = np.diff(X)/dt
dC = np.dot(dX, dX.T)
eigenvalues, eigenvectors = sp.linalg.eig(dC, C)
order = eigenvalues.argsort()
eigenvalues = eigenvalues[order]
eigenvectors = eigenvectors[:, order]