|
55 | 55 | end |
56 | 56 | nTrials = size(inputs{1},length(size(inputs{1}))); |
57 | 57 |
|
58 | | - for i=1:opts.xtrp |
59 | | - randidx = randperm(nTrials, nTrials); |
60 | | - npartition = [1 2 4]; |
61 | | - PID2 = repmat({zeros(1,nTimepoints)}, 1, length(outputs)); |
62 | | - PID4 = repmat({zeros(1,nTimepoints)}, 1, length(outputs)); |
63 | | - I1_2 = 0; |
64 | | - I2_2 = 0; |
65 | | - I12_2 = 0; |
66 | | - I1_4 = 0; |
67 | | - I2_4 = 0; |
68 | | - I12_4 = 0; |
69 | | - inputs_s = inputs; |
70 | | - idx = repmat({':'}, 1, ndims(inputs{1})); |
71 | | - idx{end} = randidx; |
72 | | - for var = 1:length(inputs) |
73 | | - inputs_s{var} = inputs{var}(idx{:}); |
74 | | - end |
75 | | - for np=2:length(npartition) |
76 | | - for pidx = 1:npartition(np) |
77 | | - inputs_p = partition(inputs_s, npartition(np), pidx,0); |
78 | | - PID_p = feval(corefunc, inputs_p, outputs, plugin_opts); |
79 | | - for outIdx = 1:length(outputs) |
80 | | - if npartition(np)==2 |
81 | | - PID2{outIdx} = PID2{outIdx} + PID_p{outIdx}/2; |
82 | | - elseif npartition(np)==4 |
83 | | - PID4{outIdx} = PID4{outIdx} + PID_p{outIdx}/4; |
| 58 | + if opts.parallel |
| 59 | + noutputs = length(outputs); |
| 60 | + corrected_v_par = zeros(opts.xtrp,nTimepoints); |
| 61 | + I1_corrected_par = zeros(opts.xtrp,nTimepoints); |
| 62 | + I2_corrected_par = zeros(opts.xtrp,nTimepoints); |
| 63 | + I12_corrected_par = zeros(opts.xtrp,nTimepoints); |
| 64 | + parfor i=1:opts.xtrp |
| 65 | + randidx = randperm(nTrials, nTrials); |
| 66 | + npartition = [1 2 4]; |
| 67 | + PID2 = repmat({zeros(1,nTimepoints)}, 1, length(outputs)); |
| 68 | + PID4 = repmat({zeros(1,nTimepoints)}, 1, length(outputs)); |
| 69 | + I1_2 = 0; |
| 70 | + I2_2 = 0; |
| 71 | + I12_2 = 0; |
| 72 | + I1_4 = 0; |
| 73 | + I2_4 = 0; |
| 74 | + I12_4 = 0; |
| 75 | + inputs_s = inputs; |
| 76 | + idx = repmat({':'}, 1, ndims(inputs{1})); |
| 77 | + idx{end} = randidx; |
| 78 | + for var = 1:length(inputs) |
| 79 | + inputs_s{var} = inputs{var}(idx{:}); |
| 80 | + end |
| 81 | + for np=2:length(npartition) |
| 82 | + for pidx = 1:npartition(np) |
| 83 | + inputs_p = partition(inputs_s, npartition(np), pidx,0); |
| 84 | + PID_p = feval(corefunc, inputs_p, outputs, plugin_opts); |
| 85 | + for outIdx = 1:noutputs |
| 86 | + if npartition(np)==2 |
| 87 | + PID2{outIdx} = PID2{outIdx} + PID_p{outIdx}/2; |
| 88 | + elseif npartition(np)==4 |
| 89 | + PID4{outIdx} = PID4{outIdx} + PID_p{outIdx}/4; |
| 90 | + end |
84 | 91 | end |
85 | | - end |
86 | | - if opts.pid_constrained |
87 | | - I1_p = MI({inputs_p{1}, inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
88 | | - I2_p = MI({inputs_p{2}, inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
89 | | - I12_p = MI({cat(1, inputs_p{1}, inputs_p{2}), inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
90 | | - if npartition(np)==2 |
91 | | - I1_2 = I1_2 + I1_p{1}/2; |
92 | | - I2_2 = I2_2 + I2_p{1}/2; |
93 | | - I12_2 = I12_2 + I12_p{1}/2; |
94 | | - elseif npartition(np)==4 |
95 | | - I1_4 = I1_4 + I1_p{1}/4; |
96 | | - I2_4 = I2_4 + I2_p{1}/4; |
97 | | - I12_4 = I12_4 + I12_p{1}/4; |
| 92 | + if opts.pid_constrained |
| 93 | + I1_p = MI({inputs_p{1}, inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
| 94 | + I2_p = MI({inputs_p{2}, inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
| 95 | + I12_p = MI({cat(1, inputs_p{1}, inputs_p{2}), inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
| 96 | + if npartition(np)==2 |
| 97 | + I1_2 = I1_2 + I1_p{1}/2; |
| 98 | + I2_2 = I2_2 + I2_p{1}/2; |
| 99 | + I12_2 = I12_2 + I12_p{1}/2; |
| 100 | + elseif npartition(np)==4 |
| 101 | + I1_4 = I1_4 + I1_p{1}/4; |
| 102 | + I2_4 = I2_4 + I2_p{1}/4; |
| 103 | + I12_4 = I12_4 + I12_p{1}/4; |
| 104 | + end |
98 | 105 | end |
99 | 106 | end |
100 | 107 | end |
101 | | - end |
102 | | - x_extrap = npartition./nTrials; |
103 | | - if strcmp(opts.bias,'qe') |
104 | | - for t = 1:nTimepoints |
105 | | - for outIdx = 1:length(outputs) |
106 | | - y = [plugin_v{outIdx}(t), PID2{outIdx}(t), PID4{outIdx}(t)]; |
107 | | - p = polyfit(x_extrap, y, 2); |
108 | | - corrected_v{outIdx}(1,t) = corrected_v{outIdx}(t) + p(3)/xtrp; |
| 108 | + x_extrap = npartition./nTrials; |
| 109 | + if strcmp(opts.bias,'qe') |
| 110 | + for t = 1:nTimepoints |
| 111 | + for outIdx = 1:noutputs |
| 112 | + y = [plugin_v{outIdx}(t), PID2{outIdx}(t), PID4{outIdx}(t)]; |
| 113 | + p = polyfit(x_extrap, y, 2); |
| 114 | + corrected_v_par(outIdx, i, t) = p(3); |
| 115 | + end |
| 116 | + end |
| 117 | + elseif strcmp(opts.bias,'le') |
| 118 | + for t = 1:nTimepoints |
| 119 | + for outIdx = 1:noutputs |
| 120 | + y = [plugin_v{outIdx}(t), PID2{outIdx}(t)]; |
| 121 | + p = polyfit(x_extrap, y, 1); |
| 122 | + corrected_v_par(outIdx, i, t) = p(2); |
| 123 | + end |
109 | 124 | end |
110 | 125 | end |
111 | | - elseif strcmp(opts.bias,'le') |
112 | | - for t = 1:nTimepoints |
113 | | - for outIdx = 1:length(outputs) |
114 | | - y = [plugin_v{outIdx}(t), PID2{outIdx}(t)]; |
115 | | - p = polyfit(x_extrap, y, 1); |
116 | | - corrected_v{outIdx}(1,t) = corrected_v{outIdx}(t) + p(2)/xtrp; |
| 126 | + if opts.pid_constrained |
| 127 | + if strcmp(opts.bias,'qe') |
| 128 | + for t = 1:nTimepoints |
| 129 | + y = [I1_plugin(t), I1_2(t), I1_4(t)]; |
| 130 | + p = polyfit(x_extrap, y, 2); |
| 131 | + I1_corrected_par(i,t) = p(3); |
| 132 | + y = [I2_plugin(t), I2_2(t), I2_4(t)]; |
| 133 | + p = polyfit(x_extrap, y, 2); |
| 134 | + I2_corrected_par(i,t) = p(3); |
| 135 | + y = [I12_plugin(t), I12_2(t), I12_4(t)]; |
| 136 | + p = polyfit(x_extrap, y, 2); |
| 137 | + I12_corrected_par(i,t) = p(3); |
| 138 | + end |
| 139 | + elseif strcmp(opts.bias,'le') |
| 140 | + for t = 1:nTimepoints |
| 141 | + y = [I1_plugin(t), I1_2(t)]; |
| 142 | + p = polyfit(x_extrap, y, 1); |
| 143 | + I1_corrected_par(i,t) = p(2); |
| 144 | + y = [I2_plugin(t), I2_2(t)]; |
| 145 | + p = polyfit(x_extrap, y, 1); |
| 146 | + I2_corrected_par(i,t) = p(2); |
| 147 | + y = [I12_plugin(t), I12_2(t)]; |
| 148 | + p = polyfit(x_extrap, y, 1); |
| 149 | + I12_corrected_par(i,t) = p(2); |
| 150 | + end |
117 | 151 | end |
118 | 152 | end |
119 | 153 | end |
120 | | - if opts.pid_constrained |
| 154 | + for spar=1:noutputs |
| 155 | + corrected_v{spar}(1,:) = mean(corrected_v_par(spar,:,:),2); |
| 156 | + end |
| 157 | + I1_corrected = mean(I1_corrected_par,1); |
| 158 | + I2_corrected = mean(I2_corrected_par,1); |
| 159 | + I12_corrected = mean(I12_corrected_par,1); |
| 160 | + |
| 161 | + else |
| 162 | + for i=1:opts.xtrp |
| 163 | + randidx = randperm(nTrials, nTrials); |
| 164 | + npartition = [1 2 4]; |
| 165 | + PID2 = repmat({zeros(1,nTimepoints)}, 1, length(outputs)); |
| 166 | + PID4 = repmat({zeros(1,nTimepoints)}, 1, length(outputs)); |
| 167 | + I1_2 = 0; |
| 168 | + I2_2 = 0; |
| 169 | + I12_2 = 0; |
| 170 | + I1_4 = 0; |
| 171 | + I2_4 = 0; |
| 172 | + I12_4 = 0; |
| 173 | + inputs_s = inputs; |
| 174 | + idx = repmat({':'}, 1, ndims(inputs{1})); |
| 175 | + idx{end} = randidx; |
| 176 | + for var = 1:length(inputs) |
| 177 | + inputs_s{var} = inputs{var}(idx{:}); |
| 178 | + end |
| 179 | + for np=2:length(npartition) |
| 180 | + for pidx = 1:npartition(np) |
| 181 | + inputs_p = partition(inputs_s, npartition(np), pidx,0); |
| 182 | + PID_p = feval(corefunc, inputs_p, outputs, plugin_opts); |
| 183 | + for outIdx = 1:length(outputs) |
| 184 | + if npartition(np)==2 |
| 185 | + PID2{outIdx} = PID2{outIdx} + PID_p{outIdx}/2; |
| 186 | + elseif npartition(np)==4 |
| 187 | + PID4{outIdx} = PID4{outIdx} + PID_p{outIdx}/4; |
| 188 | + end |
| 189 | + end |
| 190 | + if opts.pid_constrained |
| 191 | + I1_p = MI({inputs_p{1}, inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
| 192 | + I2_p = MI({inputs_p{2}, inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
| 193 | + I12_p = MI({cat(1, inputs_p{1}, inputs_p{2}), inputs_p{end}}, {'I(A;B)'}, plugin_opts); |
| 194 | + if npartition(np)==2 |
| 195 | + I1_2 = I1_2 + I1_p{1}/2; |
| 196 | + I2_2 = I2_2 + I2_p{1}/2; |
| 197 | + I12_2 = I12_2 + I12_p{1}/2; |
| 198 | + elseif npartition(np)==4 |
| 199 | + I1_4 = I1_4 + I1_p{1}/4; |
| 200 | + I2_4 = I2_4 + I2_p{1}/4; |
| 201 | + I12_4 = I12_4 + I12_p{1}/4; |
| 202 | + end |
| 203 | + end |
| 204 | + end |
| 205 | + end |
| 206 | + x_extrap = npartition./nTrials; |
121 | 207 | if strcmp(opts.bias,'qe') |
122 | 208 | for t = 1:nTimepoints |
123 | | - y = [I1_plugin(t), I1_2(t), I1_4(t)]; |
124 | | - p = polyfit(x_extrap, y, 2); |
125 | | - I1_corrected(t) = I1_corrected(t) + p(3)/xtrp; |
126 | | - y = [I2_plugin(t), I2_2(t), I2_4(t)]; |
127 | | - p = polyfit(x_extrap, y, 2); |
128 | | - I2_corrected(t) = I2_corrected(t) + p(3)/xtrp; |
129 | | - y = [I12_plugin(t), I12_2(t), I12_4(t)]; |
130 | | - p = polyfit(x_extrap, y, 2); |
131 | | - I12_corrected(t) = I12_corrected(t) + p(3)/xtrp; |
| 209 | + for outIdx = 1:length(outputs) |
| 210 | + y = [plugin_v{outIdx}(t), PID2{outIdx}(t), PID4{outIdx}(t)]; |
| 211 | + p = polyfit(x_extrap, y, 2); |
| 212 | + corrected_v{outIdx}(1,t) = corrected_v{outIdx}(t) + p(3)/xtrp; |
| 213 | + end |
132 | 214 | end |
133 | 215 | elseif strcmp(opts.bias,'le') |
134 | 216 | for t = 1:nTimepoints |
135 | | - y = [I1_plugin(t), I1_2(t)]; |
136 | | - p = polyfit(x_extrap, y, 1); |
137 | | - I1_corrected(t) = I1_corrected(t) + p(2)/xtrp; |
138 | | - y = [I2_plugin(t), I2_2(t)]; |
139 | | - p = polyfit(x_extrap, y, 1); |
140 | | - I2_corrected(t) = I2_corrected(t) + p(2)/xtrp; |
141 | | - y = [I12_plugin(t), I12_2(t)]; |
142 | | - p = polyfit(x_extrap, y, 1); |
143 | | - I12_corrected(t) = I12_corrected(t) + p(2)/xtrp; |
| 217 | + for outIdx = 1:length(outputs) |
| 218 | + y = [plugin_v{outIdx}(t), PID2{outIdx}(t)]; |
| 219 | + p = polyfit(x_extrap, y, 1); |
| 220 | + corrected_v{outIdx}(1,t) = corrected_v{outIdx}(t) + p(2)/xtrp; |
| 221 | + end |
| 222 | + end |
| 223 | + end |
| 224 | + if opts.pid_constrained |
| 225 | + if strcmp(opts.bias,'qe') |
| 226 | + for t = 1:nTimepoints |
| 227 | + y = [I1_plugin(t), I1_2(t), I1_4(t)]; |
| 228 | + p = polyfit(x_extrap, y, 2); |
| 229 | + I1_corrected(t) = I1_corrected(t) + p(3)/xtrp; |
| 230 | + y = [I2_plugin(t), I2_2(t), I2_4(t)]; |
| 231 | + p = polyfit(x_extrap, y, 2); |
| 232 | + I2_corrected(t) = I2_corrected(t) + p(3)/xtrp; |
| 233 | + y = [I12_plugin(t), I12_2(t), I12_4(t)]; |
| 234 | + p = polyfit(x_extrap, y, 2); |
| 235 | + I12_corrected(t) = I12_corrected(t) + p(3)/xtrp; |
| 236 | + end |
| 237 | + elseif strcmp(opts.bias,'le') |
| 238 | + for t = 1:nTimepoints |
| 239 | + y = [I1_plugin(t), I1_2(t)]; |
| 240 | + p = polyfit(x_extrap, y, 1); |
| 241 | + I1_corrected(t) = I1_corrected(t) + p(2)/xtrp; |
| 242 | + y = [I2_plugin(t), I2_2(t)]; |
| 243 | + p = polyfit(x_extrap, y, 1); |
| 244 | + I2_corrected(t) = I2_corrected(t) + p(2)/xtrp; |
| 245 | + y = [I12_plugin(t), I12_2(t)]; |
| 246 | + p = polyfit(x_extrap, y, 1); |
| 247 | + I12_corrected(t) = I12_corrected(t) + p(2)/xtrp; |
| 248 | + end |
144 | 249 | end |
145 | 250 | end |
146 | 251 | end |
| 252 | + |
147 | 253 | end |
148 | 254 | if opts.pid_constrained |
149 | 255 | if ~isfield(opts, 'chosen_atom') |
|
232 | 338 | plugin_opts.bias = 'plugin'; |
233 | 339 | plugin_opts.recall = true; |
234 | 340 | nTrials = size(inputs{1},length(size(inputs{1}))); |
| 341 | + |
| 342 | + if opts.parallel |
| 343 | + else |
| 344 | + end |
235 | 345 | for i=1:opts.xtrp |
236 | 346 | randidx = randperm(nTrials, nTrials); |
237 | 347 | npartition = [1 2 4]; |
|
0 commit comments