1+ import numpy as np ;
2+ import matplotlib .pyplot as plt ;
3+ import scipy .fft as fft ;
4+
5+ class Sequence2param :
6+
7+
8+ def KLDivergence ( P , Q ):
9+ if len (P ) != len (Q ):
10+ return None ;
11+ sum = 0 ;
12+
13+ for i in range (len (P )):
14+ sum = np .dot (P , np .log2 ( np .divide (P , Q )));
15+
16+ return sum ;
17+
18+ def JSDivergence ( P , Q ):
19+ M = (P + Q )/ 2 ;
20+ JS = 0.5 * ( KLDivergence (P ,M ) + KLDivergence (Q , M ));
21+ return JS ;
22+
23+ def __init__ (self , filepath ):
24+ self .AList = np .loadtxt (filepath + "Acceleration.txt" );
25+ self .VList = np .loadtxt (filepath + "LatticeDepth.txt" );
26+
27+ self .AVListIndex = np .loadtxt (filepath + "AVIndex.txt" , dtype = int );
28+ self .MomProb = np .loadtxt (filepath + "MomentumProbability.txt" ); # np array with rows containing momentum probabilities for each [a,V] value pair in AVList
29+ self .a0 = 0.0 ;
30+ self .V0 = 10.0 ;
31+
32+ def BayesianUpdating (self , totalmeasurements , recordstep , totalrecords ):
33+ PossibleMomentumOutcomes = np .array ( [- 10 + 2 * i for i in range (0 ,11 )]); # values of momentum in n\hbar k_L
34+ PossibleOutcomes = range (0 ,11 );
35+ datamom = np .reshape (self .MomProb , (self .AList .size ,self .VList .size ,11 ));
36+
37+ P_actual = np .array (datamom [50 ,25 ,:]); # Fix this hardcoded value
38+
39+ P_actual = P_actual / np .sum (P_actual );
40+ P_simulated = P_actual ; #No errors
41+
42+ Runs = totalmeasurements ; # How many simulated data do we want
43+ outcomes = np .random .default_rng ().choice (PossibleOutcomes ,size = Runs , p = P_simulated );
44+ unique , frequency = np .unique (outcomes , return_counts = True );
45+
46+ PaVprior = np .full ((self .AList .size , self .VList .size ),1 )/ (self .AList .size * self .VList .size );
47+
48+ plotPaV = np .array ([PaVprior ]);
49+
50+ counter = 0 ;
51+ for m in outcomes :
52+ for i in range (self .AList .size * self .VList .size ):
53+ indexpair = self .AVListIndex [i ];
54+ MomentumProbabilities = self .MomProb [i ];
55+
56+ PaVprior [indexpair [0 ], indexpair [1 ]] *= (MomentumProbabilities [m ])
57+ PaVprior /= np .sum (PaVprior )
58+
59+ counter += 1 ;
60+ if counter % recordstep and counter < totalrecords :
61+ plotPaV = np .append (plotPaV ,[PaVprior ], axis = 0 );
62+
63+ return plotPaV ;
64+
65+
66+ def JSDMatrix (self ):
67+ momindices = np .where (self .AVListIndex [:,1 ]== 25 )[0 ];
68+ momproblist = self .MomProb [momindices ];
69+ indices = self .AVListIndex [momindices ];
70+
71+ accindices = indices [:, 0 ];
72+ acc = self .AList [accindices ]
73+
74+ no_of_values = len (accindices );
75+ JSDivergenceMatrix = np .zeros ((no_of_values , no_of_values )) ;
76+
77+ for i in accindices :
78+ for j in accindices :
79+ JSDivergenceMatrix [i ][j ]= JSDivergence (momproblist [i ], momproblist [j ])
80+
81+ return JSDivergenceMatrix ;
0 commit comments