Completed cross-checks. Cleared out most TODO items.

TODO.txt

@@ -1,13 +1,5 @@
---  double-check line 59 in create_mask: 
-    temp     = ~diff( state_matrix, 1, 2);
-    make sure that this is turning all non-zero values (including negative ones) into zeros, and all zeros into ones
-
 -- double check that R=1/6 works out with the calculations in Vestergaard
 
--- double-check your dx/dy calculation in build_xyl_trackmat.m line 40 
-
 --  add correction for merger/split events with already-dimerized tracks.
 
--- check that changing state_matrix  inside create_mask doesn't affect it outside function.
-
 -- fill out the README

get_diffdat.m

@@ -1,35 +1,37 @@
-function [ Diffdat, dout, dndnp1, D_observations ] =  get_diffdat ( state_matrices_allti, trackdat_xyl, max_state, dt, Nframes ,R )
-% Diffdat is a data structure (list) containing diffusion parameters for each state
-% dout is an x-y array of all the changes in position. Checking for drift
+function [ diffconst_vals, d2out_Plist, dout_all, dndnp1, D_observations ] =  get_diffdat ( state_matrices_allti, trackdat_xyl, max_state, dt, Nframes ,R )
+% diffconst_vals is a data structure (list) containing average diffusion parameters for each state
+% d2out_Plist is an x-y array of all the squared changes in position for each state --likewise for _all 
+% except binned together into a single list. 
 % dndnp1 is an array of correlations between adjacent steps. Checking for drift.
 % D_observations is an array of diffusion constants calculated for each track, regardless of state or duration.
 
 % ----- initialize -----------------
 D_observations  = [];
-dout.x          = [];
-dout.y          = [];
+d2out_Plist     = [];
+dout_all.x      = [];
+dout_all.y      = [];
 dndnp1          = [];
 
 % -------------------------------------------------
-% ---- First collect Diffdat for each polymer state
+% ---- First collect diffconst_vals for each polymer state
 
 Num_comp_tracks = length( trackdat_xyl);
 
 for S= 1:max_state
 
-    dx2               = calculate_dx2(  state_matrices_allti, trackdat_xyl, Nframes, S);
-    dx2_ave           = mean(dx2);
+    [ d2out_Plist{S}.dx2, d2out_Plist{S}.dy2 ] = get_dxy2_list(  state_matrices_allti, trackdat_xyl, Nframes, S);
+    MSD_ave               = mean( ( d2out_Plist{S}.dx2 + d2out_Plist{S}.dy2 ) );
 
-    dxndxnp1          = calculate_dxnnp1( state_matrices_allti, trackdat_xyl, Nframes, S);
-    dxndxnp1_ave      = mean(dxndxnp1);
+    dxndxnp1_list         = get_dxnnp1_list( state_matrices_allti, trackdat_xyl, Nframes, S);
+    dxndxnp1_ave          = mean( dxndxnp1_list );
 
     % this is the diffusion constant, using the formula from Eq. 14 (page 022726-7
     % from Vestergaard et al, Phys. Rev. E. 89, (2014)
-    Diffdat{S}.D       = (dx2_ave/2*(dt)) + (dxndxnp1_ave/dt);
-    Diffdat{S}.sigma2  = R*(dx2_ave) + (2*R-1)*dxndxnp1_ave;
+    diffconst_vals{S}.D       = (MSD_ave/4*(dt)) + (dxndxnp1_ave/dt);
+    diffconst_vals{S}.sigma2  = R*(MSD_ave) + (2*R-1)*dxndxnp1_ave;
 
-    if( Diffdat{S}.sigma2 <1  )
-        disp( strcat("WARNING: sigma^2 =", num2str(Diffdat{S}.sigma2 ), " for state ",num2str( S ) ) )
+    if( diffconst_vals{S}.sigma2 < 1  )
+        disp( strcat("WARNING: sigma^2 =", num2str(diffconst_vals{S}.sigma2 ), " for state ",num2str( S ) ) )
     end
 
 end
@@ -76,8 +78,8 @@
    end   
 
 
-   dout.x  = [dout.x, d.x];
-   dout.y  = [dout.y, d.y];
+   dout_all.x  = [dout_all.x, d.x];
+   dout_all.y  = [dout_all.y, d.y];
    % ----- these will constitute the x-y scatter plot ------------
 
    clear mask;
@@ -104,4 +106,4 @@
 
 end
 
-end
+end

calculate_dxnnp1.m → get_dxnnp1_list.m

@@ -1,4 +1,4 @@
-function [ dnnp1 ] =  calculate_dxnnp1 (state_matrices_allti, trackdat_xyl, Nframes, S)
+function [ dnnp1_list ] =  get_dxnnp1_list (state_matrices_allti, trackdat_xyl, Nframes, S)
 % collect change in position between two adjacent frames
 
 Num_comp_tracks     = length( trackdat_xyl );
@@ -14,7 +14,7 @@
    mask          =  create_mask( state_matrices_allti{ti}, Nframes, S, 2 );
 
    if( size(mask,2) ~= Nframes - 2 || size(mask,1) ~= size( trackdat_xyl(ti).dx, 1) )
-      disp("ERROR: mismatched frame length in calc_dnnp1" )
+      disp("ERROR: mismatched frame length in get_dnnp1_list" )
       return
    end
 
@@ -37,7 +37,7 @@
 
 end
 
-dnnp1 = mean( ( dxnnp1_list + dynnp1_list )/2 );
+dnnp1_list = ( dxnnp1_list + dynnp1_list )/2;
 % if there is any mismatched dimension between x and y,
 % then there is a larger-scale problem that we want flagged.
 

calculate_dx2.m → get_dxy2_list.m

@@ -1,11 +1,11 @@
-function [ d2 ] =  calculate_dx2 (state_matrices_allti, trackdat_xyl, Nframes, S)
+function [ dx2_list, dy2_list ] =  get_dxy2_list (state_matrices_allti, trackdat_xyl, Nframes, S)
 % collect squared change in position between adjacent frames
 
 Num_comp_tracks     = length( trackdat_xyl );
 
 % initialize results
-dx_list = [];
-dy_list = [];
+dx2_list = [];
+dy2_list = [];
 % ---------------------------------
 
 for ti = 1:Num_comp_tracks
@@ -15,7 +15,7 @@
    mask          =  create_mask( state_matrices_allti{ti}, Nframes, S, 1 );
 
    if( ~all( size(mask) == [size( trackdat_xyl(ti).dx, 1), (Nframes-1) ] ) )
-      disp("ERROR: mismatched frame length in calc_d2")
+      disp("ERROR: mismatched frame length in get_d2_list")
       return
    end
 
@@ -32,14 +32,13 @@
        dy2_vec = transpose( dy2_vec );
    end
 
-   dx_list = [ dx_list , dx2_vec ];
-   dy_list = [ dy_list , dy2_vec ];
+   dx2_list = [ dx2_list , dx2_vec ];
+   dy2_list = [ dy2_list , dy2_vec ];
 
 end
 
-d2 = mean( (dx_list+dy_list)/2  );
+% d2_list = (dx_list+dy_list)/2 ;
 % if there is any mismatched dimension between x and y,
 % then there is a larger-scale problem that we want flagged.
 
 end
-

polytrack.m

@@ -1,7 +1,7 @@
 % take TracksFinal type data structure and dt spacing and some label, and
 % output relevant plots
 
-function [ lifetime_list, density, lumen_list,  Diffdat_p, D_observations] =  polytrack ( tracks_input_RAW, Label, dt, px_spacing, R, Nbin)
+function [ lifetime_list, density, lumen_list, d2out_Plist, Diffconst_vals, D_observations] =  polytrack ( tracks_input_RAW, Label, dt, px_spacing, R, Nbin)
 
 % ===================================================
 % dat_in = "E:\NikonTIRF\04-10-18\beta1\141\TrackingPackage\tracks\Channel_1_tracking_result"
@@ -108,18 +108,18 @@
 % same convention as above:
 
 
-[ Diffdat_p, d, dndnp1, D_observations ]  = get_diffdat( state_matrices_allti, trackdat_xyl, max_state, dt, Nframes, R  );
+[ Diffconst_vals, d2out_Plist, dout_all, dndnp1, D_observations ]  = get_diffdat( state_matrices_allti, trackdat_xyl, max_state, dt, Nframes, R  );
 
 figure(7);
 subplot(2,1,1);
-plot(d.x, d.y, '.');
+plot(dout_all.x, dout_all.y, '.');
 xlabel("dx");
 ylabel("dy");
 title( strcat('position change -scatter') );
 
 subplot(2,1,2);
 hist(dndnp1, 2*Nbin);
-xlim([-5, 5]);
+xlim([-0.0005, 0.0005]);
 xlabel("dx_n * dx_{n+1}");
 ylabel("Freq");
 title( strcat('Two-frame drift correlation: ', Label) );
@@ -129,7 +129,7 @@
 xlabel("Observed diffusion constant");
 ylabel("Freq");
 title( strcat('Diffusion constant calculation (like in PNAS 2013): ', Label) );
-xlim([0, 50]);
+xlim([0, 0.0050]);
 
 % for the diffusion constant, consider these functions:
 %  http://tinevez.github.io/msdanalyzer/
@@ -138,4 +138,3 @@
 %
 
 end
-