WSINDy_CellCluster/classify_models.m at main · MathBioCU/WSINDy_CellCluster · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
classify_inputs_txt=fileread('classify_inputs.m');

J_fx = length(fx_fcn_cell);
J_fv = length(fv_fcn_cell);
J_hx = length(hx_fcn_cell);
J_hv = length(hv_fcn_cell);
J_dx = length(dx_fcn_cell);
J_dv = length(dv_fcn_cell);

test_tinds = 1:floor(length(tobs)*test_tinds_frac);

nux = nu_learned*nufac_x;
nuv = nu_learned*nufac_v;
Xtest=Xscell_obs{1}(:,:,test_tinds);
Vtest=Vscell_obs{1}(:,:,test_tinds);
tobs_test=tobs(test_tinds);
[N,d,M]=size(Xtest);

species_inds={};
species_models={};
valpairs_cell={};
species_indsInModels={};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% learn models

subinds = 1:length(ninds);

while ~isempty(subinds)

    %%% cluster models
    ge = compute_errs(ninds(subinds),simdat(subinds),Xscell_obs,Vscell_obs,neighbs(subinds),errfun);
    [Wsmat,Mod,inds_keep,ge_score,models,inds_pat,ninds_rep]=...
        clustermodels(alpha,beta,gamma,normGE,normALL,tol1,tol2,ninds(subinds),simdat(subinds),neighbs(subinds),algout(subinds),J_fv,J_fx,J_hv,J_hx,J_dv,J_dx,ge);

    %%% compute top model
    if ~isempty(inds_pat)

        species_indsInModels = [species_indsInModels,{{Wsmat(:,inds_pat{1}{1}),inds_pat{1}{1},Mod,ninds_rep}}];

        W=mean(Wsmat(:,inds_pat{1}{1}),2);
        %%% could use median here
        W(log10(abs(W))<max(log10(abs(W)))-logcutoff)=0;

        [f_learned,h_learned,d_learned] = gen_force_fcn_xv(W,fx_fcn_cell,fv_fcn_cell,hx_fcn_cell,hv_fcn_cell,dx_fcn_cell,dv_fcn_cell);

        %%% get cells to validate
        valid_cells=ninds(subinds);
        num_gen=length(valid_cells);

        %%% Validate model on remaining cells
        [valpairs,~] = test_neighbs(Xtest,Vtest,tobs,nu_learned,nufac_x,nufac_v,f_learned,h_learned,d_learned,valid_cells,opts,subdt,avg_v0,test_tinds_frac,verbose,1);

        valpairs_cell=[valpairs_cell,{valpairs}];

        %%% compute log val. errs
        relerr = compute_errs(valid_cells,valpairs,{Xtest},{Vtest},[],errfun);
        frelerr = log10(relerr);

        %%% if >100(halt_prob)% of cells have error > 100(accept_err)%, check for 2-species
        if all([sum(relerr>accept_err)/length(relerr)>halt_prob length(species_inds)<max_species length(relerr)>halt_num])
            %%% fit to gaussian
            gm1=fitgmdist(frelerr(:),1);
            idk1 = cluster(gm1,frelerr(:));

            %%% if 2-guassian mixture has lower mean BIC, 2-species, if not,
            %%% forms single species.
            idks=repmat(idk1*0,1,num_gm_tries);
            bic = zeros(1,num_gm_tries);
            for j=1:num_gm_tries
                gm=fitgmdist(frelerr(:),2,'RegularizationValue',10^-6);
                idk = cluster(gm,frelerr(:));
                bic(j)=gm.BIC;
                [~,ii]=min(gm.mu);
                idks(:,j)= idk==ii;
            end
            if mean(bic) < gm1.BIC
                disp('multi-species')
                subinds=subinds(mean(idks,2)>0.5);
            else
                disp('mono-species')
            end
            species_inds=[species_inds,{subinds}];
            species_models=[species_models,{{W,f_learned,h_learned,d_learned}}];
            subinds = oppinds(cell2mat(species_inds),length(ninds));
        else
            low_err_inds = find(relerr<accept_err)';
            if ~isempty(low_err_inds)
                species_inds=[species_inds,{low_err_inds}];
                species_models=[species_models,{{W,f_learned,h_learned,d_learned}}];
            end
            species_inds=[species_inds,{oppinds(cell2mat(species_inds),length(ninds))}];
            species_models=[species_models,{{}}];
            subinds = [];
        end
    else
        species_inds=[species_inds,{oppinds(cell2mat(species_inds),length(ninds))}];
        species_models=[species_models,{{}}];
        subinds = [];
    end
end