fixes sampling based PA and general cleanup

ICB-DCM · Mar 6, 2017 · 5d387c8 · 5d387c8
1 parent e917cf3
commit 5d387c8
Show file tree

Hide file tree

Showing 15 changed files with 229 additions and 596 deletions.
diff --git a/examples/example1/estimation_toy.m b/examples/example1/estimation_toy.m
@@ -120,7 +120,7 @@
 
 clc
 % Options for multi-start optimization
-options.fmincon = optimset('algorithm','interior-point',...
+options.localOptimizerOptions = optimset('algorithm','interior-point',...
     'GradObj','on',...
     'MaxIter',1000,...
     'TolFun',1e-8,...
@@ -130,7 +130,9 @@
 options.n_starts = 20;
 options.comp_type = 'sequential';
 options.mode = 'visual';
-options.calc_profiles = 'true';
+options.calc_profiles = true;
+options.parameter_index = [3,4];
+options = PestoOptions(options);
 parameters_MEM.true = xi;
 
 parameters_MEM = getMultiStarts(parameters_MEM,@(theta) logLMEMOIR(theta,Data,Model),options);

diff --git a/examples/example2/estimation_decay.m b/examples/example2/estimation_decay.m
@@ -132,21 +132,24 @@
 
 clc
 % Options for multi-start optimization
-options.fmincon = optimset('algorithm','interior-point',...
+options.localOptimizerOptions = optimset('algorithm','interior-point',...
     'GradObj','on',...
     'MaxIter',1000,...
     'TolFun',1e-6,...
     'TolX',1e-8,...
     'display','iter',... %        'display','final',...
     'MaxFunEvals',1000*parameters_MEM.number);
-options.n_starts = 20;
+options.n_starts = 5;
 options.comp_type = 'sequential';
 options.mode = 'visual';
-options.calc_profiles = 'true';
+options.calc_profiles = true;
+options.parameter_index = [3,4];
+options = PestoOptions(options);
 parameters_MEM.true = xi;
 
-parameters_MEM = getMultiStarts(parameters_MEM,@(theta) logLMEMOIR(theta,Data,Model),options);
-parameters_MEM = getProfiles(parameters_MEM,@(theta) logLMEMOIR(theta,Data,Model),options);
+options_logL.plot = false;
+parameters_MEM = getMultiStarts(parameters_MEM,@(theta) logLMEMOIR(theta,Data,Model,options_logL),options);
+parameters_MEM = getParameterProfiles(parameters_MEM,@(theta) logLMEMOIR(theta,Data,Model,options_logL),options);
 save(['./project/results/' model '/' datafile '/' date_now '/MS_MEM/S=[' strrep(num2str(S),'  ','_') ']_' filename '.mat'],...
     'parameters_MEM',...
     'Data',...

diff --git a/examples/example2/project/models/experiments_decay_diag.m b/examples/example2/project/models/experiments_decay_diag.m
@@ -10,12 +10,12 @@
     Model.exp{s}.N = 10;
     Model.exp{s}.sigma_noise = 1;
     Model.exp{s}.noise_on = 1;
-    Model.exp{s}.t = [0:0.01:1]';
+    Model.exp{s}.t = [0:0.1:1]';
 
     Model.exp{s}.ind_phi = [1,2];
     Model.exp{s}.sym.phi = [beta(1)+b(1);
         beta(2)+b(2)];
-    Model.exp{s}.sym.sigma_noise = sym(0.01);
+    Model.exp{s}.sym.sigma_noise = sym(0.18);
     Model.exp{s}.sym.sigma_time = sym.empty(0,1);
 
     Model.exp{s}.time_model = 'normal';

diff --git a/getSCTLfits.m b/getSCTLfits.m
@@ -15,6 +15,7 @@
     optionsMultistart.comp_type = 'sequential';
     optionsMultistart.mode = 'visual';
     optionsMultistart.obj_type = 'log-posterior';
+    optionsMultistart.MCMC.nsimu_run = 1e6;
 
     if(isfield(Data{iData},'SCTL'))
         if(~isfield(Data{iData}.SCTL,'T'))
@@ -28,6 +29,8 @@
             objective_function = @(phi) objective_phi_wrapper(phi,Model,Data,iData,icl);
 
             parameters_SCTL{iD,icl} = getMultiStarts(parameters, objective_function, optionsMultistart);
+%             parameters_SCTL{iD,icl} = getParameterProfiles(parameters_SCTL{iD,icl}, objective_function, optionsMultistart);
+            parameters_SCTL{iD,icl} = getParameterSamples(parameters_SCTL{iD,icl}, objective_function, optionsMultistart);
         end
 
        iD = iD+1; 

diff --git a/getSimulationPA.m b/getSimulationPA.m
@@ -15,27 +15,37 @@
     else
         op_SP.approx = 'sp';
     end
+
+    logFlag = strcmp(op_SP.approx,'sp');
 %     SP = testSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.PA.time,phi,Data{s}.condition),xi,Model.exp{s},op_SP);
-    SP = getSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.PA.time,phi,Data{s}.condition),xi,Model.exp{s},op_SP);
+    SP = getSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.PA.time,phi,Data{s}.condition,logFlag),xi,Model.exp{s},op_SP);
 
 %     [g,g_fd_f,g_fd_b,g_fd_c] = testGradient(xi,@(x)getSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.PA.time,phi,Data{s}.condition),x,Model.exp{s},op_SP),1e-4,'my','dmydxi')
 
-    % we log transformed in simulateForSP, now we need to compute the mean
-    % from a lognormal density
-    tmp = arrayfun(@(x) diag(squeeze(SP.Cy(x,:,:))),1:size(SP.Cy,1),'UniformOutput',false);
-    my = exp(SP.my+transpose([tmp{:}])/2);
-    my(isnan(my)) = 0;
+    % we log transformed in simulateForSP (only if we used sigma point approximation)
+    % now we need to compute the mean from a lognormal density
+    if(logFlag)
+        tmp = arrayfun(@(x) diag(squeeze(SP.Cy(x,:,:))),1:size(SP.Cy,1),'UniformOutput',false);
+        my = exp(SP.my+transpose([tmp{:}])/2);
+        my(isnan(my)) = 0;
+    else
+        my = SP.my;
+    end
 
     if(nderiv>0)
-        nt = size(SP.dCydxi,1);
-        np = size(SP.dCydxi,4);
-        ny = size(SP.dCydxi,2);
-        dtmpdxi = arrayfun(@(x,y) diag(squeeze(SP.dCydxi(x,:,:,y))),repmat(1:nt,[np,1]),...
-            repmat(transpose(1:np),[1,nt]),'UniformOutput',false);
-        dmydxi = bsxfun(@times,my,SP.dmydxi) ...
-            + bsxfun(@times,my,permute(reshape([dtmpdxi{:}]/2,...
-            [ny,np,nt]),[3,1,2]));
-        dmydxi(isnan(dmydxi)) = 0;
+        if(logFlag)
+            nt = size(SP.dCydxi,1);
+            np = size(SP.dCydxi,4);
+            ny = size(SP.dCydxi,2);
+            dtmpdxi = arrayfun(@(x,y) diag(squeeze(SP.dCydxi(x,:,:,y))),repmat(1:nt,[np,1]),...
+                repmat(transpose(1:np),[1,nt]),'UniformOutput',false);
+            dmydxi = bsxfun(@times,my,SP.dmydxi) ...
+                + bsxfun(@times,my,permute(reshape([dtmpdxi{:}]/2,...
+                [ny,np,nt]),[3,1,2]));
+            dmydxi(isnan(dmydxi)) = 0;
+        else
+            dmydxi = SP.dmydxi;
+        end
     else
         dmydxi = zeros([size(my,1),size(my,2),length(xi)]);
     end

diff --git a/logLMEMOIR.m b/logLMEMOIR.m
@@ -269,67 +269,68 @@
 
         %% Single cell time-lapse data - Statistics
         if isfield(Data{s},'SCTLstat')
-            % Simulation using sigma points
-            op_SP.options.nderiv = options.nderiv;
-            op_SP.req = [0,0,0,1,1,1,0];
-            op_SP.type_D = Model.type_D;
-            if(extract_flag)
-                SP = testSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.SCTLstat.time,phi,Data{s}.condition),xi,Model.exp{s},op_SP);
-            else
-                SP = getSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.SCTLstat.time,phi,Data{s}.condition),xi,Model.exp{s},op_SP);
-            end
-
-            % Evaluation of likelihood, likelihood gradient and hessian
-
-            % Mean
-            logL_mz = - 0.5*sum(nansum(((Data{s}.SCTLstat.mz - SP.mz)./Data{s}.SCTLstat.Sigma_mz).^2,1),2);
-            if options.nderiv >= 2
-                dlogL_mzdxi = permute(nansum(bsxfun(@times,(Data{s}.SCTLstat.mz - SP.mz)./Data{s}.SCTLstat.Sigma_mz.^2,SP.dmzdxi),1),[2,1]);
-                if options.nderiv >= 3
-                    wdmz_SP = bsxfun(@times,1./Data{s}.SCTLstat.Sigma_mz,SP.dmzdxi);
-                    %                     wdmz_SP = reshape(wdmz_SP,[numel(SP.mz),size(SP.dmdxizdxi,3)]);
-                    ddlogL_mzdxi2 = -wdmz_SP'*wdmz_SP;
-                end
-            end
-
-
-            % Covariance
-            logL_Cz = - 0.5*sum(nansum(nansum(((Data{s}.SCTLstat.Cz - SP.Cz)./Data{s}.SCTLstat.Sigma_Cz).^2,1),2),3);
-            if options.nderiv >= 2
-                dlogL_Czdxi = squeeze(nansum(nansum(bsxfun(@times,(Data{s}.SCTLstat.Cz - SP.Cz)./Data{s}.SCTLstat.Sigma_Cz.^2,SP.dCzdxi),1),2));
-                if options.nderiv >= 3
-                    wdCzdxi = bsxfun(@times,1./Data{s}.SCTLstat.Sigma_Cz,SP.dCzdxi);
-                    wdCzdxi = reshape(wdCzdxi,[numel(SP.Cz),size(SP.dCzdxi,3)]);
-                    ddlogL_Czdxi2 = -wdCzdxi'*wdCzdxi;
-                end
-            end
-
-            % Summation
-            logL = logL + logL_mz + logL_Cz;
-            if options.nderiv >=2
-                dlogLdxi = dlogLdxi + dlogL_mzdxi + dlogL_Czdxi;
-                if options.nderiv >=3
-                    ddlogLdxidxi = ddlogLdxidxi + ddlogL_mzdxi2 + ddlogL_Czdxi2;
-                end
-            end
-
-            % Visulization
-            if options.plot
-                Sim_SCTLstat.mz = SP.mz;
-                Sim_SCTLstat.Cz = SP.Cz;
-                Model.exp{s}.plot(Data{s},Sim_SCTLstat,s);
-            end
-
-            P{s}.SCTLstat.SP = SP;
-
+            error('currently not properly implemented');
+%             % Simulation using sigma points
+%             op_SP.options.nderiv = options.nderiv;
+%             op_SP.req = [0,0,0,1,1,1,0];
+%             op_SP.type_D = Model.type_D;
+%             if(extract_flag)
+%                 SP = testSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.SCTLstat.time,phi,Data{s}.condition),xi,Model.exp{s},op_SP);
+%             else
+%                 SP = getSigmaPointApp(@(phi) simulateForSP(Model.exp{s}.model,Data{s}.SCTLstat.time,phi,Data{s}.condition),xi,Model.exp{s},op_SP);
+%             end
+%             
+%             % Evaluation of likelihood, likelihood gradient and hessian
+%             
+%             % Mean
+%             logL_mz = - 0.5*sum(nansum(((Data{s}.SCTLstat.mz - SP.mz)./Data{s}.SCTLstat.Sigma_mz).^2,1),2);
+%             if options.nderiv >= 2
+%                 dlogL_mzdxi = permute(nansum(bsxfun(@times,(Data{s}.SCTLstat.mz - SP.mz)./Data{s}.SCTLstat.Sigma_mz.^2,SP.dmzdxi),1),[2,1]);
+%                 if options.nderiv >= 3
+%                     wdmz_SP = bsxfun(@times,1./Data{s}.SCTLstat.Sigma_mz,SP.dmzdxi);
+%                     %                     wdmz_SP = reshape(wdmz_SP,[numel(SP.mz),size(SP.dmdxizdxi,3)]);
+%                     ddlogL_mzdxi2 = -wdmz_SP'*wdmz_SP;
+%                 end
+%             end
+%             
+%             
+%             % Covariance
+%             logL_Cz = - 0.5*sum(nansum(nansum(((Data{s}.SCTLstat.Cz - SP.Cz)./Data{s}.SCTLstat.Sigma_Cz).^2,1),2),3);
+%             if options.nderiv >= 2
+%                 dlogL_Czdxi = squeeze(nansum(nansum(bsxfun(@times,(Data{s}.SCTLstat.Cz - SP.Cz)./Data{s}.SCTLstat.Sigma_Cz.^2,SP.dCzdxi),1),2));
+%                 if options.nderiv >= 3
+%                     wdCzdxi = bsxfun(@times,1./Data{s}.SCTLstat.Sigma_Cz,SP.dCzdxi);
+%                     wdCzdxi = reshape(wdCzdxi,[numel(SP.Cz),size(SP.dCzdxi,3)]);
+%                     ddlogL_Czdxi2 = -wdCzdxi'*wdCzdxi;
+%                 end
+%             end
+%             
+%             % Summation
+%             logL = logL + logL_mz + logL_Cz;
+%             if options.nderiv >=2
+%                 dlogLdxi = dlogLdxi + dlogL_mzdxi + dlogL_Czdxi;
+%                 if options.nderiv >=3
+%                     ddlogLdxidxi = ddlogLdxidxi + ddlogL_mzdxi2 + ddlogL_Czdxi2;
+%                 end
+%             end
+%             
+%             % Visulization
+%             if options.plot
+%                 Sim_SCTLstat.mz = SP.mz;
+%                 Sim_SCTLstat.Cz = SP.Cz;
+%                 Model.exp{s}.plot(Data{s},Sim_SCTLstat,s);
+%             end
+%             
+%             P{s}.SCTLstat.SP = SP;
+%             
         end
 
         %% Single cell snapshot data
         if isfield(Data{s},'SCSH')
 
             switch(options.nderiv)
                 case 0
-                    [SP,logL_m,logL_C] = logL_PA(xi, Model, Data, s, options);
+                    [SP,logL_m,logL_C] = logL_SCSH(xi, Model, Data, s, options);
                 case 1
                     [SP,logL_m,logL_C,dlogL_mdxi,dlogL_Cdxi] = logL_PA(xi, Model, Data, s, options);
                 case 2
@@ -357,7 +358,7 @@
                     [SP,logL_m] = logL_PA(xi, Model, Data, s, options);
                 case 1
                     [SP,logL_m,dlogL_mdxi] = logL_PA(xi, Model, Data, s, options);
-                    % [g,g_fd_f,g_fd_b,g_fd_c]=testGradient(xi,@(xi) logL_PA(xi, Model, Data, s, options),1e-3,2,3)
+                    % [g,g_fd_f,g_fd_b,g_fd_c]=testGradient(xi,@(xi) logL_PA(xi, Model, Data, s, options),1e-3,2,3,true)
                 case 2
                     [SP,logL_m,dlogL_mdxi,ddlogL_mdxi2] = logL_PA(xi, Model, Data, s, options);
             end