%% Figure 3. Ichise et al. 20191018
clear

global upper_t dt;
dt=0.5; % [min], time interval: 30 sec
upper_t=101; % Upper time frame: 101, corresponding to 50 min.

%% Parameters
VIS=0; % 1: visualize trajectory of cells 
inner=10; % [ƒÊm]
outer=15; % [ƒÊm]

%% 1: Import NK cells w/o cancer cells
nk0_all=zeros(1,500,3);
len_nk0_all=0;
for file_id=1:3
    disp('NK cells w/o cancers, file_id: #/3');disp(file_id)
    nk0=ImportAndConvertTable(file_id, 0); % 0: nk cells without cancers
    
    %% Applicable to only NK cells
    [nk0, maxt_list0]=RearrangeTime(nk0); % Rearrangement of time origin in the data structure
    [nk0, maxt_list0]=ReplaceMissingTimeToNan(nk0, maxt_list0); % This is critical; the sample number will be significantly decreased if done.
    nk0=ReplaceZeroToNan(nk0, maxt_list0);
    
    %% Updating nk_all and hitprob
    list1=find(maxt_list0>0); % cells id having "all-Nan" are ignored.
    list2=find(maxt_list0==0);
    nk0(list2,:,:)=[];
    maxt_list0(list2)=[];
    
    [l, m, n]=size(nk0);
    nk0_all(len_nk0_all+1:len_nk0_all+l,1:m,1:n)=nk0(1:l,1:m,1:n);
    maxt_list0_all(len_nk0_all+1:len_nk0_all+l)=maxt_list0*dt; % [min]
    
    len_nk0_all=size(nk0_all,1);
end


%% 2: Import NK cells w/ cancers and count NK cells' hit to cancers
nk1_all=zeros(1,500,3);
len_nk1_all=0;
ca_all=zeros(1,500,3);
len_ca_all=0;
for file_id=1:7
    disp('NK cells w/ cancers, file_id: #/7');disp(file_id)
    nk1=ImportAndConvertTable(file_id, 1); % 1: nk cells with cancers
    cancer=ImportAndConvertTable(file_id, 2); % 2: cancer cells
    
    cnt_in=HitCancer(nk1, cancer, inner, outer);
    
    %% 2-1-a. Exclusion of cancers having incompete data sets
    [cancer, maxt_list_c]=RearrangeTime(cancer); % Rearrangement of time origin in the data structure, and the duration time (maxt_list_c)
    [cancer, maxt_list_c]=ReplaceMissingTimeToNan(cancer, maxt_list_c); % This is critical; the sample number will be significantly decreased if done.
    cancer=ReplaceZeroToNan(cancer, maxt_list_c);
    %% 2-1-b. Deletion of NK cells having incompete data sets
    [nk1, maxt_list1]=RearrangeTime(nk1); % Rearrangement of time origin in the data structure
    [nk1, maxt_list1]=ReplaceMissingTimeToNan(nk1, maxt_list1); % This is critical; the sample number will be significantly decreased if done.
    nk1=ReplaceZeroToNan(nk1, maxt_list1);
    
    %% 2-2-a. Updating ca_all and hitprob
    list1=find(maxt_list_c>0); % cells id having "all-Nan" are ignored.
    list2=find(maxt_list_c==0);
    cancer(list2,:,:)=[];
    cnt_in(list2)=[];
    maxt_list_c(list2)=[];
    
    [l, m, n]=size(cancer);
    ca_all(len_ca_all+1:len_ca_all+l,1:m,1:n)=cancer(1:l,1:m,1:n);
    cnt_in_all(len_ca_all+1:len_ca_all+l)=cnt_in;
    maxt_list_c_all(len_ca_all+1:len_ca_all+l)=maxt_list_c*dt; % [min]
    len_ca_all=size(ca_all,1);
    
    %% 2-2-b. Updating nk1_all
    list1=find(maxt_list1>0); % cells id having "all-Nan" are ignored.
    list2=find(maxt_list1==0);
    nk1(list2,:,:)=[];
    maxt_list1(list2)=[];
    [l1, m1, n1]=size(nk1);
    nk1_all(len_nk1_all+1:len_nk1_all+l1,1:m1,1:n1)=nk1(1:l1,1:m1,1:n1);
    maxt_list1_all(len_nk1_all+1:len_nk1_all+l1)=maxt_list1*dt; % [min]
    len_nk1_all=size(nk1_all,1);
    
end


%% 3: Normalization of position
n_nk0_all=SetPositionToOrigin(nk0_all);
n_nk1_all=SetPositionToOrigin(nk1_all);
n_ca_all=SetPositionToOrigin(ca_all);

%% 4: Quantify mode of cell motion
[inst_speed0,msd0]=CellMigrationQuantification(n_nk0_all);
[inst_speed1,msd1]=CellMigrationQuantification(n_nk1_all);


%% 5: Graph
GraphMSD(msd0, msd1);
GraphInstSpeed(inst_speed0, inst_speed1);
disp('w/o cancers: mean, median, std');
[nanmean(inst_speed1(:,upper_t-1)) nanmedian(inst_speed1(:,upper_t-1)) nanstd(inst_speed1(:,upper_t-1))]
disp('w/ cancers: mean, median, std');
[nanmean(inst_speed0(:,upper_t-1)) nanmedian(inst_speed0(:,upper_t-1)) nanstd(inst_speed0(:,upper_t-1))]


%% 6: Hit probability
CalculationHitProbability(cnt_in_all, maxt_list_c_all);

%% 7: Trajectory -- NOTE: a few minutes need to output
if VIS==1
    VisualizeTrajectory(n_nk0_all,0);
    VisualizeTrajectory(n_nk1_all,1);
    VisualizeTrajectory(n_ca_all,2);
end

