EPIFIL_LFData/DiffEnt_PVals.m at master · EdwinMichaelLab/EPIFIL_LFData · 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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150


%% Description: to calculate differential entropy and save p-values

%% Input: Data files

%% Output: Differential entropy and p-values

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model={'EPIFIL','TRANSFIL','LYMFASIM'}; % models
Scenario= 1:4;
Site = {'Kirare','Alagramam','Peneng'};
Code = {'A','B','C'}; % A = SSA, B = IND, C = PNG

for itype=1:length(model)

    % new file to save p values
    filename0 = sprintf('%s_DiffEntPVals.xlsx',model{itype});

    % initialize matrix
    entrp0 = zeros(length(Scenario)+1,length(Scenario)+1);

    for n =1:length(Site)

        NumRounds=[];

        % model only
        filename = sprintf('%s_modelonly_%s.xlsx',model{itype},Site{n});
        if exist(filename,'file')
            [num,txt,~]=xlsread(filename);
            fclose('all');
            colnum=2; % mf prev outputs start

            % number of parameter sets for each scenario
            id = 1:length(txt(:,1))-1;

            % calculate timelines to WHO 1% mf threshold
            k = 1;
            for i = id(1):id(end)
                T = find((num(i,colnum:50+colnum-1)) < 1);
                if ~isempty(T)
                    NumRounds(k,1) = T(1);
                    k = k+1;
                end
            end

            % original entropy
            entrp0(n,1) = Entropy(NumRounds(1:length(id),1));
        end

        % model+data
        filename = sprintf('%s_%s.xlsx',model{itype},Site{n});

        if exist(filename,'file')
            [num,txt,~]=xlsread(filename);
            fclose('all');
            colnum=2; % mf prev outputs start

            for iscen=1:length(Scenario)

                % number of parameter sets for each scenario, looks for
                % scenario code such as 'A1'
                if itype == 2 % these files used incorrect codes
                    id = find(strcmp(txt(:,2),sprintf('%s%s',Code{1},int2str(Scenario(iscen)))))-1;
                elseif itype == 3 % these files used incorrect codes and have an extra row in num
                    id = find(strcmp(txt(:,2),sprintf('%s%s',Code{1},int2str(Scenario(iscen)))));
                else
                    id = find(strcmp(txt(:,2),sprintf('%s%s',Code{n},int2str(Scenario(iscen)))))-1;
                end

                % calculate timelines to WHO 1% mf threshold
                k = 1;
                if itype == 2 % TRANSFIL
                    for i = id(1):id(end)
                        T = find((num(i,colnum:50+colnum-1)) < 0.01); % prev out of 1
                        if ~isempty(T)
                            NumRounds(k,iscen+1) = T(1);
                            k = k+1;
                        end
                    end
                else % EPIFIL and LYMFASIM
                    for i = id(1):id(end)
                        T = find((num(i,colnum:50+colnum-1)) < 1); % prev out of 100
                        if ~isempty(T)
                            NumRounds(k,iscen+1) = T(1);
                            k = k+1;
                        end
                    end
                end

                % record original entropy
                entrp0(n,iscen+1) = Entropy(NumRounds(:,iscen+1));

            end

            % original differential entropy
            dse0 = zeros(length(entrp0(1,:)),length(entrp0(1,:)));
            for i = 1:length(entrp0)
                for j = 1:length(entrp0)
                    dse0(i,j) = entrp0(n,i)-entrp0(n,j);
                end
            end

            % permutations
            nPerm = 20000;
            dse1 = zeros(nPerm,1);
            pVal = zeros(length(entrp0(1,:)),length(entrp0(1,:)));
            for i = 1:length(entrp0)
                for j = 1:length(entrp0)
                    x = nonzeros([NumRounds(:,i);NumRounds(:,j)]);

                    a = length(nonzeros(NumRounds(:,i))); % how many to allocate to list 1
                    b = length(nonzeros(NumRounds(:,j))); % how many to allocate to list 2

                    parfor m = 1:nPerm
                        % permute the list
                        id = randperm(length(x));
                        NumRounds1 = x(id(1:a),:);
                        NumRounds2 = x(id(a+1:a+b),:);

                        % new entropies
                        entrp1 = Entropy(NumRounds1);
                        entrp2 = Entropy(NumRounds2);
                        dse1(m) = entrp1-entrp2;

                    end

                    % calculate p value
                    pVal(i,j) = length(find(abs(dse1)>abs(dse0(i,j))))/nPerm;
                end
            end

            % save to model-specific excel files
            if n == 1
                xlswrite(filename0,pVal,1,'A1');
            elseif n == 2
                xlswrite(filename0,pVal,1,'A7');
            elseif n == 3
                xlswrite(filename0,pVal,1,'A13');
            end

        end

        % trying to save memory...
        fclose('all');
        system('taskkill /F /IM EXCEL.EXE');
        clearvars num txt x dse1
    end
end