%% load Experiment
base = "Path/To/Folder/";
result_paths = get_channel_paths(base,false);
[channel_list,cycle_list,fov_list] = load_mothersegger_data(result_paths,false);%load cycles from a specified path (either all (false) or only complete (true))
complete_cycle_ids = get_complete_cycle_ids(cycle_list);
disp('done')
%% load combined file
path = "/Path/to/File/combined_data.mat";
[cycle_list,channel_list,fov_list] = load_combined_data(path);
complete_cycle_ids = get_complete_cycle_ids(cycle_list);
%% get rid of mother cells
mother_cycle_ids = get_mother_cycle_ids(channel_list,cycle_list);
no_mother_ids = complete_cycle_ids(~ismember(complete_cycle_ids,mother_cycle_ids));
%% foci distribution
figure
plot_foci_distribution(cycle_list,complete_cycle_ids,'foci3',100);
figure
plot_foci_violin(cycle_list,complete_cycle_ids,'foci3',0.0668);
%%
all_foci = zeros(0,1);
for cycle_id = complete_cycle_ids
    cycle = cycle_list(cycle_id);
    foci_int = cycle.foci1.intensity;
    foci_int = foci_int(:);
    foci_int = foci_int(~isnan(foci_int));
    all_foci = cat(1,all_foci,foci_int);
end
%%
histogram(all_foci,200)
%% get trajectories
fluor_selector = 1;
foci_selector = 1;
trajectory_min_length = 12;%12;
interpolate_points_in_a_row = 0;
set_mean_to_0 = false;
flip_by_pole = false;
seperate = true;
remove_bb = false; %only for old data sets
cycle_ids = complete_cycle_ids;
trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);
%% show a kymograph of the trajectories
figure
kymograph_foci(trajectory_list)
%% show trajectories
fluor_selector = 1;
show_trajectory_list(cycle_list,channel_list,trajectory_list,fluor_selector);
%% plot position and velocity
figure
relative_position = false;
[pos,dpos] = plot_pos_and_delta_pos(trajectory_list,cycle_list,relative_position);
std_pos = std(pos)*0.0668
std_vel = std(dpos)*0.0668/60
%% plot velocity profile
color = [0, 0.4470, 0.7410];
alpha = 0.1;
edges = -20.5:0.5:20.5;
%edges = -1:0.1:1;
trajectory_list_tmp = trajectory_list;
[x,means,vars,counts] = get_velocity_profile_edges(trajectory_list_tmp,edges);%get_velocity_profile(trajectory_list,n_bins,relative_position);
figure(1)
plot_with_err(x,means,sqrt(vars),color,alpha);
hold on
plot_with_err(x,means,sqrt(vars)./sqrt(counts),color,alpha+0.3);
hold on
plot([x(1),x(end)],[0,0],'--k','Color',[0.5,0.5,0.5,0.5])
hold off
xlim([edges(1),edges(end)])
xlabel('position (pixel)')
xlabel('position (pixel)')%(pixel)')
%ylim([-0.05,0.05])
ylabel('velocity (\mum/sampling time)')%(pixel/sampling time)')
%% os and rp 
figure(3)
plot_os_and_rp(trajectory_list,12)
%% draw four foci kyographs for 1,2,3 and 4 foci case
for foci_selector = 1:4
    fluor_selector = 1;
    trajectory_min_length = 12;
    interpolate_points_in_a_row = 1;
    set_mean_to_0 = false;
    flip_by_pole = true;
    seperate = false;
    remove_bb = false; %only for old data sets
    cycle_ids = complete_cycle_ids;
    trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);
    %trajectory_list_border = get_border_cell_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);
    % foci kymograph
    subplot(2,2,foci_selector)
    kymograph_foci(trajectory_list)
    title(strcat(string(foci_selector)," Foci"))
    colorbar
end
%% Relative Position against relative cell age
plot_relative_kymograph(trajectory_list,100,80)
%% Trajectory MSD
plot_msd_vs_tau(cycle_list,complete_cycle_ids,[1],true,1:60,60,0.0668)
%% Foci Splitting
figure
plot_splitting_event(cycle_list,complete_cycle_ids,1,11,10);
%% Variance Vs Position
plot_variance_vs_position(trajectory_list)
%% fragment length vs time
plot_trajectory_duration_vs_length(trajectory_list,.2)
%% coefs heatmap
plot_coefs_heatmap(cycle_list,complete_cycle_ids)
%% scatter plot of all coefs vs length or rca
figure
[coefs,ls,rcas] = lag_n(trajectory_list,12,1);
scatter(rcas,coefs,8,'MarkerFaceColor','b','MarkerEdgeColor','b',...
        'MarkerFaceAlpha',0.1,'MarkerEdgeAlpha',0)
xlabel('relative cell age')
ylabel('coef')
%% mean speed
tau = 60;
ptmr = 0.0668;
speeds = get_speed(trajectory_list);
speed = mean(speeds)/tau*ptmr;
disp(strcat("mean speed ",string(speed), "(mum/sec)"))
%% mean corr. plus ste plot
edges = 26.5:1:60.5;%ls
sw = 17;
tau = 1;
plot_how_correlation_changes(trajectory_list,edges,sw,tau)
%% sigma and intensity of trajectory list
[intensities,sigmas] = plot_intensity_and_sigma(cycle_list,trajectory_list);
%% --------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%----------------------------------------------------------------Predict D--------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%% get trajectories
fluor_selector = 1;
foci_selector = 1;
trajectory_min_length = 12;
interpolate_points_in_a_row = 1;
set_mean_to_0 = false;
flip_by_pole = false;
seperate = true;
remove_bb = false; %only for old data sets
cycle_ids = complete_cycle_ids;
trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);
%% parameters
ptmr = 0.0668;
%ptmr = 1;
tau = 60;
%dt = 10;
%% predict D from velocity profile
edges = -20.5:1:20.5;
%trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,10,0.4,1,1,false);
[x,means,vars,counts] = get_velocity_profile_edges(rp_list,edges,false);
%
[D,D_prop_uncert,one_over_tau,oov_prop_uncert,f_over_kT,fokT_prop_uncert] = get_D(x*ptmr,tau,means*ptmr/tau,vars*ptmr^2/tau^2,counts,0.7);
disp(strcat("D: ",string(D)," | tau: ",string(1/one_over_tau)," | f/kT: ",string(f_over_kT)))
%% predict tau from vel autocorrelation
range = 0:20;
tau_inital_estimate = 10;
figure
[tau] = predict_tau_from_vel_auto_correlation(os_list,range,tau,tau_inital_estimate);
%% predict tau from position autocorrelation
range = 0:30;
tau_inital_estimate = 100;
[tau_prediction] = predict_tau_from_pos_auto_correlation(trajectory_list,range,tau,tau_inital_estimate);

%% predict D from msd (not acurate)
vel = zeros(0,1);
for tra_it = 1:numel(trajectory_list)
    tra = trajectory_list(tra_it);
    vel = cat(1,vel,diff(tra.positions(:,1,1)));
end
msd = var(vel)/2/1; %pixel^2/min
msd = msd*ptmr^2/tau; %\mu m^2/s
disp(strcat("MSD: ",string(msd)))
%% predict D from corssing events of oscillation
ptmr = 0.0668;
tau = 60;
[rp_list,os_list] = split_trajectory_list_by_coef(trajectory_list,1,12,5);
[D,c_over_kbT] = predict_D_from_crossing_events(os_list,ptmr,tau);
%% --------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------Annotate Data-----------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%% add annotation field to the dataset
anns = cell(numel(trajectory_list,1));
for tra_it = 1:numel(trajectory_list),anns{tra_it}=zeros(1,size(trajectory_list(tra_it).positions,1));end
[trajectory_list.annotation] = anns{:};
%% annotate the data set
trajectory_list = annotate_trajectory_list(cycle_list,channel_list,trajectory_list);
%% save annotated data
save('tl_ann','trajectory_list')
%% --------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%--------------------------------------------------------------Niche Stuff--------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------

%% pole length difference (mother cells have the pole 1)
pole1 = zeros(1,0);
poleMinus1 = zeros(1,0);
for cycle_it = 1:numel(cycle_list)
    cycle = cycle_list(cycle_it);
    if cycle.pole == 1,pole1(numel(pole1)+1)=cycle.lengths(1);end
    if cycle.pole == -1,poleMinus1(numel(poleMinus1)+1)=cycle.lengths(1);end
end
mean(pole1)
mean(poleMinus1) 


%% pole bias 
new_pole = 0;
old_pole = 0;
trajectory_list_tmp = trajectory_list;
for tra_it = 1:numel(trajectory_list_tmp)
    tra = trajectory_list_tmp(tra_it);
    cycle = cycle_list(tra.cycle_id);
    if cycle.pole == 1,new_pole = new_pole+1;end
    if cycle.pole == -1,old_pole = old_pole+1;end
end
new_pole = new_pole
old_pole = old_pole


%% check bias of poles and mother cells
check_poles_and_mother_cells(channel_list,cycle_list,trajectory_list_nsb)

%% get trajectories
fluor_selector = 1;
foci_selector = 1;
trajectory_min_length = 12;
interpolate_points_in_a_row = 1;
set_mean_to_0 = false;
flip_by_pole = false;
seperate = true;
remove_bb = false; %only for old data sets
cycle_ids = complete_cycle_ids;
trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);

%% hist where trajectories come from
%trajectory_list = get_trajectories(complete_cycle_ids,cycle_list,foci_selector,10,0.4,1,1,false);
tra_list_tmp = os_list;
channel_ids = zeros(numel(tra_list_tmp),1);
average_time = zeros(numel(tra_list_tmp),1);
for tra_it = 1:numel(tra_list_tmp)
    tra = tra_list_tmp(tra_it);
    cycle = cycle_list(tra.cycle_id);
    channel_ids(tra_it) = cycle.channel_id;
    average_time(tra_it) = (cycle.times(tra.start)+cycle.times(tra.stop))/2;
end
figure(1)
subplot(1,2,1)

histogram(channel_ids,0.5:1:(numel(channel_list)+0.5))
%histogram(average_time)
sgtitle('where do the oscillating trajectories come from')
xlabel('channel id')
ylabel('number')
subplot(1,2,2)
scatter(channel_ids,average_time,'x')
xlabel('channel id')
ylabel('time (min)')
%% indef view on foci distribution across time and space
plasmid_position = nan(max([channel_list.n_frames]),numel(channel_list));
for channel_it = 1:numel(channel_list)
    channel = channel_list(channel_it);
    for frame_it = 1:channel.n_frames
        cycle_ids = channel.cycle_ids;
        cycle_indices = channel.cycle_indices;
        nof = nan(1,size(cycle_indices,2));%number of foci
        level = nan(1,size(cycle_indices,2));%cell_height
        for cycle_it = 1:size(cycle_indices,2)
            cycle_id = cycle_ids(frame_it,cycle_it);
            cycle_timing = cycle_indices(frame_it,cycle_it);
            if ~isnan(cycle_id)
                cycle = cycle_list(cycle_id);
                if cycle.flags.complete
                    nof(cycle_it) = sum(~isnan(cycle.foci1.intensity(cycle_timing,:)));
                    level(cycle_it) = cycle.bbs(cycle_timing,1,1);
                end
            end
        end
        %[~,pos] = max(level);
        %if(pos)
        plasmid_position(frame_it,channel_it) = min(nof);%mean(nof,'omitnan');%nof(pos);%/size(cycle_indices,2);
        %end
    end
    %plasmid_position(:,channel_it)
end
%
figure(2)
subplot(1,2,1)
imAlpha=ones(size(plasmid_position));
imAlpha(isnan(plasmid_position))=0;
imagesc(plasmid_position,'AlphaData',imAlpha)
set(gca,'color',0*[1 1 1]);
xlabel('channel id')
ylabel('frame')
colorbar()
subplot(1,2,2)
bar(mean(plasmid_position,1,'omitnan'))
xlabel('channel_id')
ylabel('mean min foci')
%% get trajectories
fluor_selector = 1;
foci_selector = 1;
trajectory_min_length = 12;
interpolate_points_in_a_row = 1;
set_mean_to_0 = false;
flip_by_pole = false;
seperate = true;
remove_bb = false; %only for old data sets
cycle_ids = complete_cycle_ids;
trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);

%% collapsing correlation velocity
tau = 60;
figure(1)
for lag_it = 1:8
    y = autocorrelation_vel(trajectory_list,lag_it);
    y = y(1:5*lag_it+1);
    x_all = (0:numel(y)-1)/lag_it;%for plotting all together
    x_indi = (0:numel(y)-1)*tau;%for plotting individual
    tau = tau*lag_it;
    tau = tau_prediction_from_vel(x_indi,y,tau,200);
    x_tau = x_indi(1):0.1:x_indi(end);
    y_tau = velocity_auto_correlation_fun(x_tau,tau,tau);
    subplot(3,3,lag_it+1)
    plot(x_indi/tau,y)
    hold on
    plot(x_tau/tau,y_tau)
    hold off
    title(strcat('autocorr. at lag=',string(lag_it)," tau=",string(tau)))
    xlabel('t (min)')
    ylabel('autocorr.')
    subplot(3,3,1)
    plot(x_all,y)
    if lag_it == 1,hold on;end
end
hold off
xlabel('t/lag')
ylabel('autocorr.')
legend(strcat('lag=',string(1:8)))
title('autocorr. at different lag')
%%
tau = 60;
tau = 100;
figure(2)
offsets_all = [0.001,0.1,0.5,1,2,5,10,100,1000];
offsets_all = [0.001,0.1,1,10,100,1000];
for offset_it = 1:numel(offsets_all)
    dt2 = offsets_all(offset_it)*60;
    %x_indi = (0:numel(y)-1)*dt;%for plotting individual
    x_tau = 0:0.01:dt2*3;
    y_tau = velocity_auto_correlation_fun(x_tau,dt2,tau);
    plot(x_tau/dt2,y_tau)
    if offset_it == 1,hold on;end
end
title(strcat('tau=',string(tau)))
xlabel('t/lag')
ylabel('autocorr.')
legend(strcat('dt=',string(offsets_all)))
hold off
%%
tau = 1;
figure(2)
taus = [0.1,0.5,1,2,5,50,100];
for offset_it = 1:numel(taus)
    tau = taus(offset_it);
    %x_indi = (0:numel(y)-1)*dt;%for plotting individual
    x_tau = 0:0.01:tau*3;
    y_tau = velocity_auto_correlation_fun(x_tau,tau,tau);
    plot(x_tau/tau,y_tau)
    if offset_it == 1,hold on;end
end
title(strcat('lag=',string(tau)))
xlabel('t/lag')
ylabel('autocorr.')
legend(strcat('tau=',string(taus)))
hold off
%% --------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%----------------------------------------------------------------Functions--------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------
%%---------------------------------------------------------------------------------------------------------------------------------------------



function [tau_prediction] = tau_prediction_from_vel(x,y,lag,tau_inital_estimate)
    %velocity_auto_correlation_fun = @(t,dt,tau) (2*exp(-t/tau)-exp(-(abs(t-dt))/tau)-exp(-(t+dt)/tau))./(2-2*exp(-dt/tau));
    %minimize_this_output = @(tau) sum((velocity_auto_correlation_fun(x,lag,tau)-y).^2);
    minimize_this_output = @(tau) sum(abs((velocity_auto_correlation_fun(x,lag,tau)-y)));
    tau_prediction = fminsearch(minimize_this_output,tau_inital_estimate);
    %x = x(1):0.01:x(end);
    %y = velocity_auto_correlation_fun(x,t,tau_prediction);
end

function y = velocity_auto_correlation_fun(t,dt,tau)
    y = (2*exp(-t/tau)-exp(-abs(t-dt)/tau)-exp(-(t+dt)/tau))./(2-2*exp(-dt/tau));
end


function y = autocorrelation_vel(trajectory_list,t)
    res_counter = 0;
    for tra_it = 1:numel(trajectory_list)
        res_counter = res_counter + size(trajectory_list(tra_it).positions,2);
    end
    results = cell(1,res_counter);
    res_pointer = 1;
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        positions = tra.positions;
        for pos_it = size(positions,2)
            pos = positions(:,pos_it,1);
            [corr,~] = get_autocorrelation_vel(pos,t);
            results{res_pointer} = corr;
            res_pointer = res_pointer + 1;
        end
    end
    n = max(cellfun(@numel,results));
    results2 = nan(res_counter,n);
    for res_it = 1:res_counter
        corr = results{res_it};
        results2(res_it,1:numel(corr)) = corr;
    end
    y = mean(results2,1,'omitnan');
    %figure(100)
    %plot(y(1:30))
end

function [corr,lags] = get_autocorrelation_vel(pos,t)
    vel = pos(1+t:end)-pos(1:end-t);
    [corr,lags] = xcorr(vel,'normalized');
    lags = lags(numel(vel):end)/t;
    corr = corr(numel(vel):end);
end

function plot_os_and_rp(trajectory_list,lag)

    trajectory_list_tmp = trajectory_list;%cat(2,os_list,rp_list);
    [coefs,ls,rcas] = lag_n(trajectory_list_tmp,lag,1);
    %edges = 24.5:1:80.5;%ls
    edges = min(ls)-0.5:1:max(ls)+0.5;%ls
    %edges = -0.05:.05:1.05;%rca
    buckets = cell(1,numel(edges));
    for i = 1:numel(edges),buckets{i}=zeros(1,0);end
    for coef_it = 1:numel(coefs)
        t = ls(coef_it);
        %t = rcas(coef_it);
        c = coefs(coef_it);
        index = find((t<edges+1)==1,1,'first');
        if index > 0
            regime_change = buckets{index};
            regime_change(numel(regime_change)+1) = c;
            buckets{index} = regime_change;
        end
    end
    % mean corr. plus ste plot
    means = cellfun(@mean,buckets);
    stds = cellfun(@std,buckets);
    counts = cellfun(@numel,buckets);

    rp_count = cellfun(@(x) sum(x<0),buckets);
    os_count = cellfun(@(x) sum(x>0),buckets);
    edges_microm = edges*0.0668;
    %subplot(4,2,7)
    subplot(1,2,1)
    area(edges_microm,rp_count)
    hold on
    area(edges_microm,os_count,'FaceAlpha',0.85)
    hold off
    %legend('rp','os')
    xlabel('cell length (\mum)')
    ylabel('absolute abundance')
    xlim([edges_microm(1),edges_microm(end)])
    rp_count = cellfun(@(x) sum(x<0)/numel(x),buckets);
    os_count = cellfun(@(x) sum(x>0)/numel(x),buckets);
    %subplot(4,2,8)
    subplot(1,2,2)
    area(edges_microm,os_count+rp_count)
    hold on
    area(edges_microm,os_count)
    hold off
    %legend('rp','os')
    xlabel('cell length (\mum)')
    ylabel('relative abundance')
    xlim([edges_microm(1),edges_microm(end)])
    %sgtitle('Sliding window: 12') 
    sgtitle('')
end

function [trajectory_list_nsb] = get_border_cell_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb)
    not_bad_border_cells = zeros(0,1);
    for cycle_it = 1:numel(cycle_list)
        cycle = cycle_list(cycle_it);
        if(cycle.flags.border_cell && numel(cycle.relation_graph.parent)==1 && ~isnan(cycle.relation_graph.parent(1)) && cycle.duration>5 && ~cycle.flags.complete && ~cycle.flags.inferred_segmentation)
            not_bad_border_cells(numel(not_bad_border_cells)+1,1) = cycle_it;
        end
    end
    cycle_ids = not_bad_border_cells;
    trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);
    %
    trajectory_list_nsb_count = 1;
    trajectory_list_nsb = zeros(0,1);
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        cycle = cycle_list(tra.cycle_id);
        valid_points = cycle.bbs(tra.start:tra.stop,1,1)>20;
        new_start = find(valid_points,1,'first');
        new_stop = find(valid_points,1,'last');
        if(numel(new_start)>0 && new_stop-new_start+1>trajectory_min_length)
            tra.start = tra.start+new_start-1;
            tra.stop = tra.start+new_stop-1;
            tra.positions = tra.positions(new_start:new_stop,:,:);
            tra.lengths = tra.lengths(new_start:new_stop);
            tra.widths = tra.widths(new_start:new_stop);
            tra.relative_cell_age = tra.relative_cell_age(new_start:new_stop);
            if trajectory_list_nsb_count == 1
                trajectory_list_nsb = tra;
            else
                trajectory_list_nsb(trajectory_list_nsb_count) = tra;
            end

            trajectory_list_nsb_count = trajectory_list_nsb_count + 1;
        end
    end
end


function [coefs,ls,rcas,indices] = lag_tra(tra,sw,lag)
    pos = tra.positions(:,1,1);
    n_points = numel(pos);
    if sw == -1
        sw2=n_points-lag-1;
    else
        sw2 = sw;
    end
    n = n_points-sw2-lag;
    coefs = zeros(1,n);
    rcas = zeros(1,n);
    ls = zeros(1,n);
    indices = zeros(1,n);
    if size(pos,1)-lag>2 %enough for p value
        for i = 1:n
            t0 = diff(pos(i:i+sw2));
            t1 = diff(pos(i+lag:i+sw2+lag));
            middle = round(((i)+(i+sw2+1))/2);
            %[corr,p] = corrcoef(t0,t1);
            %coefs(i) = corr(2);
            corr = xcorr(diff(pos(i:i+sw2+lag)),lag,'normalized');
            coefs(i) = corr(end);%/corr(ceil(end/2));
            rcas(i) = tra.relative_cell_age(middle);
            ls(i) = tra.lengths(middle);
            indices(i) = middle;
        end
    end
end

function [coefs,ls,rcas,indices] = lag_tra_pos(tra,sw,lag)
    pos = tra.positions(:,1,1);
    n_points = numel(pos);
    if sw == -1
        sw2=n_points-lag-1;
    else
        sw2 = sw;
    end
    n = n_points-sw2-lag;
    coefs = zeros(1,n);
    rcas = zeros(1,n);
    ls = zeros(1,n);
    indices = zeros(1,n);
    if size(pos,1)-lag>2 %enough for p value
        for i = 1:n
            t0 = pos(i:i+sw2);
            t1 = pos(i+lag:i+sw2+lag);
            middle = round(((i)+(i+sw2+1))/2);
            %[corr,p] = corrcoef(t0,t1);
            %coefs(i) = corr(2);
            corr = xcorr(pos(i:i+sw2+lag),lag,'normalized');
            coefs(i) = corr(end);%/corr(ceil(end/2));
            rcas(i) = tra.relative_cell_age(middle);
            ls(i) = tra.lengths(middle);
            indices(i) = middle;
        end
    end
end

function [D,c_over_kbT] = predict_D_from_crossing_events(trajectory_list,ptmr,dt)

    mid_crossing = zeros(0,1);
    for i = 1:numel(trajectory_list)
        tra = trajectory_list(i);
        pos = tra.positions(:,1,1);
        mc = pos<0;
        mc = diff(mc)~=0; %mid corssing events
        vel = diff(pos);
        mid_crossing = cat(1,mid_crossing,abs(vel(mc)));
    end
    D = var(mid_crossing)/dt/2*ptmr^2;
    c_over_kbT = mean(mid_crossing)/var(mid_crossing)/2/ptmr;
    disp(strcat("predicted D (crossing): ",string(D)))
end

function [tau_prediction] = predict_tau_from_pos_auto_correlation(trajectory_list,range,dt,tau_inital_estimate)

    means = zeros(1,numel(range));
    stds = zeros(1,numel(range));
    counts = zeros(1,numel(range));
    t_l_t = trajectory_list;
    for lag_it = 1:numel(range)
        lag = range(lag_it);
        [coefs,~,~] = lag_n_pos(t_l_t,-1,lag);
        means(lag_it) = mean(coefs);
        stds(lag_it) = std(coefs);
        counts(lag_it) = numel(coefs);
    end
    position_auto_correlation_fun = @(t,dt,tau) exp(-abs(0-t)/tau);
    minimize_this_output = @(x,c) sum((position_auto_correlation_fun(range*dt,dt,x)-means).^2);
    tau_prediction = fminsearch(minimize_this_output,tau_inital_estimate);
    
    xdata = range*dt;
    position_auto_correlation_fun = @(x,xdata) (exp(-abs(0-xdata)/x(1))+x(2))/(1+x(2));
    position_auto_correlation_fun([1,1],2);
    parameters = lsqcurvefit(position_auto_correlation_fun,[100,0],xdata,means);
    
    
    disp(strcat("predicted tau (pos): ",string(parameters(1))," and constant (c): ",string(parameters(2))))

    x = (range(1):0.1:range(end))*dt;
    %y = position_auto_correlation_fun(x,dt,tau_prediction);
    y = position_auto_correlation_fun(parameters,x);
    x = x/dt;
    figure(1)
    plot(range*1,means)
    hold on 
    plot(x,y)
    hold off
    legend('auto corr. exp.','auto corr. prediction')
    title(strcat("predicted tau (pos): ",string(parameters(1))," and constant (c): ",string(parameters(2))))
    xlabel('lag (min)')
    ylabel('corr.')
end


function [tau_prediction] = predict_tau_from_vel_auto_correlation(trajectory_list,range,dt,tau_inital_estimate)
    means = zeros(1,numel(range));
    stds = zeros(1,numel(range));
    counts = zeros(1,numel(range));
    for lag_it = 1:numel(range)
        lag = range(lag_it);
        [coefs,~,~] = lag_n(trajectory_list,-1,lag);
        means(lag_it) = mean(coefs);
        stds(lag_it) = std(coefs);
        counts(lag_it) = numel(coefs);
    end
    position_auto_correlation_fun = @(t,dt,tau) (2*exp(-t/tau)-exp(-(abs(t-dt))/tau)-exp(-(t+dt)/tau))./(2-2*exp(-dt/tau));
    minimize_this_output = @(x) sum((position_auto_correlation_fun(range*dt,dt,x)-means).^2);
    tau_prediction = fminsearch(minimize_this_output,tau_inital_estimate);
    disp(strcat("predicted tau (vel): ",string(tau_prediction)))

    x = (range(1):0.1:range(end))*dt;
    y = position_auto_correlation_fun(x,dt,tau_prediction);
    x = x/dt;
    plot(range*1,means)
    hold on 
    plot(x,y)
    %for i = 140:10:170
    %    y = position_auto_correlation_fun(x,dt,i);
    %    plot(x,y)
    %end
    hold off
    xlim([x(1),x(end)])
    legend('auto corr. exp.','auto corr. prediction')
    title(strcat("predicted tau (vel): ",string(tau_prediction)))
    xlabel('lag (sampling time)')
    ylabel('autocorr.')
end

function plot_trajectory_duration_vs_length(trajectory_list,aplha)
    res = zeros(2,numel(trajectory_list));
    for i = 1:numel(trajectory_list)
        tra = trajectory_list(i);
        mean_length = mean(tra.lengths);
        tra_length = numel(tra.lengths);
        res(:,i) = [mean_length,tra_length];
    end
    x = res(1,:);
    y = res(2,:);
    figure(1)
    scatter(y,x,'MarkerFaceColor','b','MarkerEdgeColor','b',...
        'MarkerFaceAlpha',aplha,'MarkerEdgeAlpha',aplha)
    ylabel('mean fragment length (pixel)')
    xlabel('mean fragment duration (min)')
    title('fragment length vs dration')
end

function [speed_all] = get_speed(trajectory_list)

    speed_all = zeros(0,1);
    for i = 1:numel(trajectory_list)
        tra = trajectory_list(i);
        pos = tra.positions(:,1,1);
        speed = abs(diff(pos));
        speed_all = cat(1,speed_all,speed);
    end
end

function plot_coefs_heatmap(cycle_list,cycle_ids)
    %figure(101)
    results = cell(2,4);
    for foci_selector = 1:4
        fluor_selector = 1;
        %foci_selector = 1;
        trajectory_min_length = 12;
        interpolate_points_in_a_row = 1;
        set_mean_to_0 = false;
        flip_by_pole = true;
        seperate = true;
        remove_bb = false; %only for old data sets
        trajectory_list = get_trajectories(cycle_ids,cycle_list,foci_selector,trajectory_min_length,fluor_selector,interpolate_points_in_a_row,set_mean_to_0,flip_by_pole,seperate,remove_bb);

        tau = 1;
        sw = 12;
        min_length = 6;
        %path = strcat('Tmp/',string(ml_it),'_',string(sw_it),'_os_or_rp.png');
        [rp_list,os_list] = split_trajectory_list_by_coef(trajectory_list,tau,sw,min_length);

        ptmr = 0.0668;
        tmp_list = cat(2,rp_list,os_list);
        lengths = cat(1,tmp_list.lengths)*ptmr;
        coefs = cat(2,tmp_list.coefs)';

        values_x = lengths;
        values_y = coefs;
        edges_x = 1:0.25:8;
        edges_y = -.8:0.05:.6;
        label_x = 'length (\mum)';
        label_y = 'auto corr.';
        plot_title = strcat(string(foci_selector)," Foci (N=",string(numel(values_x)),")");
        subplot(2,2,foci_selector)
        [x,means] = heatmap_with_mean(values_x,values_y,edges_x,edges_y,label_x,label_y,plot_title);
        results{1,foci_selector} = x;
        results{2,foci_selector} = means;
    end
end

function plot_how_correlation_changes(trajectory_list,edges,sw,lag)
    [coefs,ls,rcas] = lag_n(trajectory_list,sw,lag);
    buckets = cell(1,numel(edges));
    for i = 1:numel(edges),buckets{i}=zeros(1,0);end
    for coef_it = 1:numel(coefs)
        t = ls(coef_it);
        %t = rcas(coef_it);
        c = coefs(coef_it);
        index = find((t<edges)==1,1,'first');
        if index > 0
            regime_change = buckets{index};
            regime_change(numel(regime_change)+1) = c;
            buckets{index} = regime_change;
        end
    end
    means = cellfun(@mean,buckets);
    stds = cellfun(@std,buckets);
    counts = cellfun(@numel,buckets);
    figure(2)
    plot([edges(1),edges(end)],[0,0],'--k')
    hold on
    plot_with_err(edges,means,stds./sqrt(counts),[1,0,0],0.5)
    hold off
    xlabel('rca')
    ylabel('corr.')
    title('corr vs length (ste)')
end



function [coefs,ls,rcas] = lag_n(trajectory_list,sw,lag)
    coefs = zeros(1,0);
    ls = zeros(1,0);
    rcas = zeros(1,0);
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        [coefs_t,ls_t,rcas_t] = lag_tra(tra,sw,lag);
        coefs = cat(2,coefs,coefs_t);
        ls = cat(2,ls,ls_t);
        rcas = cat(2,rcas,rcas_t);
    end
end

function [coefs,ls,rcas] = lag_n_pos(trajectory_list,sw,lag)
    coefs = zeros(1,0);
    ls = zeros(1,0);
    rcas = zeros(1,0);
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        [coefs_t,ls_t,rcas_t] = lag_tra_pos(tra,sw,lag);
        coefs = cat(2,coefs,coefs_t);
        ls = cat(2,ls,ls_t);
        rcas = cat(2,rcas,rcas_t);
    end
end



% annotate trajectory list. adds a new field to the tra_list called
% annotations. it contains a vector which can as elements in [-1,0,1] where
% -1 means oscillation, (os), 0 undetermined (un), 1 regular positioning (rp)
function trajectory_list = annotate_trajectory_list(cycle_list,channel_list,trajectory_list)
    current = 1;
    while 1

        tra = trajectory_list(current);
        cycle = cycle_list(tra.cycle_id);
        canvas = display_cycle_with_range2(cycle_list(tra.cycle_id),channel_list,0,'fluor1_subtracted',1,tra.start:1:tra.stop,true);

        subplot(2,1,1)
        imshow(canvas)
        subplot(2,1,2)
        pos = tra.positions;
        p = squeeze(pos(:,1,1))*-1;
        x = 1:numel(p);
        transitions = tra.annotation;
        plot(x,p,'-b')
        hold on
        scatter(x(transitions==-1),p(transitions==-1),'r','fill')
        scatter(x(transitions==0),p(transitions==0),'b','fill')
        scatter(x(transitions==1),p(transitions==1),'g','fill')
        plot(tra.lengths/2,'-k')
        plot(tra.lengths/2*-1,'-k')
        hold off
        xticks(0:5:500)
        xlim([x(1),x(end)])
        %break
        prompt = 'what u want to do: ';
        str = input(prompt,'s');
        if strcmp(str,'exit'),break,end
        if startsWith(str,'os'),tmp=strsplit(str,' ');tra.annotation(str2num(tmp{end-1}):str2num(tmp{end}))=-1;end
        if startsWith(str,'un'),tmp=strsplit(str,' ');tra.annotation(str2num(tmp{end-1}):str2num(tmp{end}))=0;end
        if startsWith(str,'rp'),tmp=strsplit(str,' ');tra.annotation(str2num(tmp{end-1}):str2num(tmp{end}))=1;end
        
        if numel(tra.annotation) > size(tra.positions,1) %if someone puts in a to high ending index this will take care of it
            tra.annotation = tra.annotation(1:size(tra.positions,1));
        end
        
        trajectory_list(current) = tra;


        if strcmp(str,'next'),current=current+1;end
        if strcmp(str,'prev'),current=current-1;end
        if startsWith(str,'goto'),tmp=strsplit(str,' ');current=str2num(tmp{end});end
    end
    disp('it is over')
end

function show_trajectory_list(cycle_list,channel_list,trajectory_list,fluor_selector)
    current = 1;
    %fluor_selector = 1;
    while 1
        
        tra = trajectory_list(current);
        cycle = cycle_list(tra.cycle_id);
        subplot(2,1,1)
        %canvas = display_cycle_with_range(cycle,channel_list,0,{'fluor3_subtracted','fluor2_subtracted',''},fluor_selector,tra.start:1:tra.stop);
        canvas = display_cycle_with_range(cycle,channel_list,0,strcat('fluor',string(fluor_selector),'_subtracted'),fluor_selector,tra.start:1:tra.stop);
        pause(1)
        subplot(2,1,2)
        pos = tra.positions;
        for tra_it = 1:size(pos,2)
            p = squeeze(pos(:,tra_it,1))*-1;
            x = 1:numel(p);
            ptmr = 0.0668;
            plot(x,p*ptmr,'-b')
            if tra_it == 1,hold on;end
            scatter(x,p*ptmr,'b','fill')
        end
        plot(x,tra.lengths/2*ptmr,'-k')
        plot(x,(tra.lengths/2*-1)*ptmr,'-k')
        hold off
        xlabel('time (min)')
        ylabel('position')
        xticks(0:5:500)
        xlim([x(1),x(end)])
        %break
        prompt = 'what u want to do: ';
        str = input(prompt,'s');
        if strcmp(str,'exit'),break,end
        if strcmp(str,'next'),current=current+1;end
        if strcmp(str,'prev'),current=current-1;end
        if startsWith(str,'goto'),tmp=strsplit(str,' ');current=str2num(tmp{end});end
        if startsWith(str,'fluor'),tmp=strsplit(str,' ');fluor_selector=str2num(tmp{end});end
    end
    disp('it is over')
end

%changes the length of all trajectories in the trajectory list by the
%facktor given by scale
function trajectory_list = change_length(trajectory_list,scale)
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        tra.positions = tra.positions*scale;
        tra.lengths = tra.lengths*scale;
        tra.widths = tra.widths*scale;
        trajectory_list(tra_it) = tra;
    end
end



function [trajectory_list] = flip_trajectories(trajectory_list)
    for i = 1:numel(trajectory_list)
        tra = trajectory_list(i);
        tmp = sort(mean(tra.positions(:,:,1),1));
        if sum(tmp <= 0) >= 2
            tra.positions = tra.positions*-1;
            trajectory_list(i) = tra;
        end
    end
end
%{
function [foci_numbers] = plot_foci_distribution(cycle_list,cycle_ids)

    foci_numbers = zeros(0,1);
    for cycle_it = 1:numel(cycle_ids)
        cycle = cycle_list(cycle_ids(cycle_it));
        foci_numbers = cat(1,foci_numbers,sum(~isnan(cycle.foci1.intensity),2));
    end
    histogram(foci_numbers,-.5:1:5.5)
    xlabel('number of foci')
    ylabel('absolute abundance')
end
%}
function [intensities,sigmas] = plot_intensity_and_sigma(cycle_list,trajectory_list)

    intensities = zeros(0,1);
    sigmas = zeros(0,1);
    for i = 1:numel(trajectory_list)
        tra = trajectory_list(i);
        sta = tra.start;
        sto = tra.stop;
        cycle = cycle_list(tra.cycle_id);
        intensity = cycle.foci1.intensity(sta:sto,:);
        intensities = cat(1,intensity(:),intensities);
        sigma = cycle.foci1.sigma(sta:sto,:);
        sigmas = cat(1,sigma(:),sigmas);
    end
    subplot(1,3,1)
    histogram(intensities)
    title('intensities')
    subplot(1,3,2)
    histogram(sigmas)
    title('sigmas')
    subplot(1,3,3)
    hist3(cat(2,intensities,sigmas))
    xlabel('intensities')
    ylabel('sigmas')

end

function check_poles_and_mother_cells(channel_list,cycle_list,trajectory_list)

    [mother_cycle_ids] = get_mother_cycle_ids(channel_list,cycle_list);
    min_l = 1000000;
    max_l = 0;
    figure(1)
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        l = tra.lengths;
        if min_l>min(l),min_l=min(l);end
        if max_l<max(l),max_l=max(l);end
    end
    n = max_l-min_l+1;
    plasmid_position = cell(1,n); %plasmid position based on length
    for i = 1:n
        plasmid_position{i} = zeros(1,0);
    end
    cell_length_based_on_pole = zeros(n,2);
    cell_length_based_on_mother = zeros(n,2);
    for tra_it = 1:numel(trajectory_list)
        tra = trajectory_list(tra_it);
        pole = cycle_list(tra.cycle_id).pole;
        pos = tra.positions(:,1,1);
        l = tra.lengths;
        for i = 1:numel(l)
            tmp = plasmid_position{l(i)-min_l+1};
            tmp(numel(tmp)+1) = pos(i);
            plasmid_position{l(i)-min_l+1} = tmp;
            if pole==1,cell_length_based_on_pole(l(i)-min_l+1,1)=cell_length_based_on_pole(l(i)-min_l+1,1)+1;end
            if pole==-1,cell_length_based_on_pole(l(i)-min_l+1,2)=cell_length_based_on_pole(l(i)-min_l+1,2)+1;end

            if ismember(tra.cycle_id, mother_cycle_ids),cell_length_based_on_mother(l(i)-min_l+1,1)=cell_length_based_on_mother(l(i)-min_l+1,1)+1;end
            if ~ismember(tra.cycle_id, mother_cycle_ids),cell_length_based_on_mother(l(i)-min_l+1,2)=cell_length_based_on_mother(l(i)-min_l+1,2)+1;end
        end
    end
    means = cellfun(@mean,plasmid_position);
    stds = cellfun(@std,plasmid_position);
    counts = cellfun(@numel,plasmid_position);
    subplot(1,3,1)
    plot_with_err(min_l:max_l,means,stds./sqrt(counts),[1,0,0],0.2)
    hold on 
    plot([min_l,max_l],[0,0],'--k')
    xlim([min_l,max_l])
    xlabel('cell length (pixel)')
    ylabel('plasmid position')
    subplot(1,3,2)
    plot(min_l:max_l,counts)
    hold on
    plot(min_l:max_l,cell_length_based_on_pole(:,1))
    plot(min_l:max_l,cell_length_based_on_pole(:,2))
    hold off
    legend('all','pole=1','pole=-1')
    xlim([min_l,max_l])
    xlabel('cell lengths (pixel)')
    ylabel('abundance')
    title('cell length based on pole')
    subplot(1,3,3)
    plot(min_l:max_l,counts)
    hold on
    plot(min_l:max_l,cell_length_based_on_mother(:,1))
    plot(min_l:max_l,cell_length_based_on_mother(:,2))
    hold off
    legend('all','mother','not a mother')
    xlim([min_l,max_l])
    xlabel('cell lengths (pixel)')
    ylabel('abundance')
    title('cell length based on pole')
end

%{

%%
%tl_saved = trajectory_list;

rp_count = 1;
os_count = 1;
un_count = 1;
for tra_it = 1:numel(trajectory_list)
    tra = trajectory_list(tra_it);
    ann = tra.annotation;
    transitions = find(~(diff(tra.annotation)==0));
    transitions = cat(2,0,transitions,size(tra.positions,1));
    
    for seq_it = 1:numel(transitions)-1
        start = transitions(seq_it)+1;
        stop = transitions(seq_it+1);
        identifier = ann(start);
        
        
        trajectory_struct = struct();
        trajectory_struct.cycle_id = tra.cycle_id;
        trajectory_struct.score = tra.score;
        trajectory_struct.start = tra.start+start-1;
        trajectory_struct.stop = tra.start+(stop-start)+start-1;
        trajectory_struct.positions = tra.positions(start:stop,:,:);
        trajectory_struct.relative_cell_age = tra.relative_cell_age(start:stop);
        trajectory_struct.widths = tra.widths(start:stop);
        trajectory_struct.lengths = tra.lengths(start:stop);
        if identifier == -1
            os_list(os_count) = trajectory_struct;
            os_count = os_count + 1;
        elseif identifier == 0
            un_list(un_count) = trajectory_struct;
            un_count = un_count + 1;
        else
            rp_list(rp_count) = trajectory_struct;
            rp_count = rp_count + 1;
        end
        
        
    end
end


%% fix tra list
for i = 1:numel(trajectory_list)
    tra = trajectory_list(i);
    if size(tra.positions,1) ~= numel(tra.annotation)
        i
        [size(tra.positions,1),numel(tra.annotation)]
        if numel(tra.annotation) > size(tra.positions,1) %if someone puts in a to high ending index this will take care of it
            tra.annotation = tra.annotation(1:size(tra.positions,1));
        end
        trajectory_list(i) = tra;
    end
end
%}