Extened README file, and modified h2mmtest.py to process FRETbursts data structure

Author: harripd

README.md

@@ -8,4 +8,16 @@ This repo is written in Python by Paul David Harris, with the objective of imple
 
 The primary functions of this library are the EM_H2MM, and the equivalent EM_H2MM_par, which differ only in that EM_H2MM_par is accelerated by parallel processing.
 These functions use Maximum Likelihood Estimators to optimize a hidden Markov Model of data produced primarily by diffusing smFRET measurements- which are composed of arrival times of single photons at two or more detectors. This method is built off of well established hidden Markov modeling (HMM) but extends it to efficiently handel variable times between observations.
-Inputs to the functions are a h2mm_model object, which is composed of the prior or initial state vector, the transition probability matrix, and the observational likelihood matrix.
+Inputs to the functions are:
+    h2mm_model object: an object of the h2mm_model class specific to H2MMpythonlib, which contains three numpy arrays:
+        prior: initial state distribution
+        trans: the transition probability matrix
+        obs: the emmision proability matrix, giving the probability of a photon arriving at a given detector in a given state
+    ArrivalColors: a list of numpy int arrays correspondig to the colors (detectors)(indexed from 0) of the photons
+    ArrivalTimes: a list of numpy int arrays, 1 array per burst, corresponding to the arrival times of the photons
+The output of EM_H2MM and EM_H2MM_par is another h2mm_model object, which contains the optimized parameters.
+
+This algorithm is both powerful and risky to use.
+The main difficulty arrises due to the strong potential of overfitting.
+Yet the potential to detect and quantify hidden dynamics is too good to pass up.
+I ask any researcher to apply this method cautiously, but hopefully.

h2mmtest.py

@@ -9,40 +9,44 @@
 from H2MM import *
 from fretbursts import *
 
-def H2MM_eval(h_model,data):
+def H2MM_eval(h_model,data):    
     DAemDex_mask = data.get_ph_mask(ph_sel=Ph_sel(Dex='DAem'))
     AemDex_mask = data.get_ph_mask(ph_sel=Ph_sel(Dex='Aem'))
     AemDex_mask_red = AemDex_mask[DAemDex_mask]
     d_bursts = data.mburst[0]
-    h_mod = EM_H2MM(h_model,data_ph_short,data_time_short)
-    
+    ph_d = data.get_ph_times(ph_sel=Ph_sel(Dex='DAem'))
+    d_bursts_red = d_bursts.recompute_index_reduce(ph_d)
+    data_color = []
+    data_time = []
+    for start, stop in zip(d_bursts.istart,d_bursts.istop+1):
+        temp_ph_mask = AemDex_mask_red[start:stop]
+        temp_color = np.zeros((temp_ph_mask.shpae),dtype=int)
+        temp_color[temp_ph_mask] = 1
+        data_time.append(ph_d[start:stop])
+        data_color.append(temp_color)
+    h_mod = EM_H2MM_par(h_model,data_color,data_time)
+    return h_mod
 
-trans2 = np.array([[0.9999,0.0001],[0.001,0.999]])
-prior2 = np.array([0.3,0.7])
-obs2 = np.array([[0.1,0.9],[0.4,0.6]])
-h_mod2 = h2mm_model(prior2,trans2,obs2)
-trans3 = np.array([[0.98, 0.01, 0.01],[0.01, 0.98, 0.01],[0.01, 0.01, 0.98]])
-prior3 = np.array([0.2, 0.6, 0.2])
-obs3 = np.array([[0.1, 0.9],[0.5, 0.5],[0.7, 0.3]])
-h_mod3 = h2mm_model(prior3, trans3, obs3)
-data_time_short = np.zeros(12,dtype=int)
-
-for i in range(1,12):
-    if i != 4:
-        data_time_short[i] = i + data_time_short[i-1]
+if __name__ == "__main__" :
+    filename = './data/minus8TA_minus6NTD_RPo_T25C_G150uW_R100uW_1.hdf5'
+    if os.path.isfile(filename):
+        print('Perfect, I found the file!')
     else:
-        data_time_short[i] = i + data_time_short[i-1] + 3
-
-data_time_short = np.array([0, 4, 6, 13, 19, 21, 50, 54, 55, 59, 70, 72])
-data_ph_short = np.zeros(12,dtype=int)
-for i in range(1,len(data_time_short)):
-    data_ph_short[i] = data_time_short[i] % 2
-
-deltas_short = np.diff(data_time_short)
-R_short = ph_factors(deltas_short)
-deltas_short = np.diff(data_time_short)
-R_short = ph_factors(deltas_short)
-transmat_t_short = CalculatePowerofTransMatrices(h_mod2,R_short.R)
-transmat_t3_short = CalculatePowerofTransMatrices(h_mod3,R_short.R)
-
-h_mod = EM_H2MM(h_mod2,[data_ph_short],[data_time_short])
+        print('ERROR: file does not exist')
+    d = loader.photon_hdf5(filename)
+    d.add(det_donor_accept=(0, 1), 
+          alex_period=4000, 
+          D_ON=(2100, 3900), 
+          A_ON=(150, 1900),
+          offset=700)
+    loader.usalex_apply_period(d)
+    d.calc_bg(fun=bg.exp_fit,time_s=50.1, tail_min_us='auto', F_bg=1.7)
+    d.burst_search(m=10, F=6, ph_sel=Ph_sel(Dex='DAem'))
+    d.fuse_bursts(ms=0)
+    d_all = Sel(d, select_bursts.naa, th1=50)
+    d_all = Sel(d_all, select_bursts.size, th1=50)
+    prior2 = np.array([0.1, 0.9])
+    trans2 = np.array([[0.998, 0.002],[0.0001, 0.9999]])
+    obs2 = np.array([[0.3, 0.7],[0.8, 0.2]])
+    h_mod2i = h2mm_model(prior2,trans2,obs2)
+    h_model = H2MM_eval(h_mod,d_all)