fix

ssl-ctbn-code/PowerSet/PowerSet.m

-% By Paulo Abelha
-%
-% Returns the powerset of set S
-%
-% S is a cell array
-% P is a cell array of cell arrays
-function [ P ] = PowerSet( S )
-    n = numel(S);
-    x = 1:n;
-    P = cell(1,2^n);
-    p_ix = 2;
-    for nn = 1:n
-        a = combnk(x,nn);
-        for j=1:size(a,1)
-            P{p_ix} = S(a(j,:));
-            p_ix = p_ix + 1;
-        end
-    end
-end
-

ssl-ctbn-code/PowerSet/PowerSet_r.m

-% By Paulo Abelha
-%
-% Returns the powerset of set S
-%
-% S is a cell array
-% P is a cell array of cell arrays
 function [ P ] = PowerSet_r( S ,N)
     N
     x = 1:N;

ssl-ctbn-code/ctbns_base/CTBN_cond_stat_star_sparse_reg_DIMS.m

@@ -17,7 +17,7 @@ if (isempty(node(i).parents)==0)
         M_i=0;       
         pa=ones(t,1);
         ki=1;       
-        for k=node(i).parents           
+        for k = node(i).parents           
               mu_p=MU{k};
               pa=pa.*mu_p(1:t,ps(mn,ki));  
               ki=ki+1;           

ssl-ctbn-code/ctbns_base/CTBN_cond_stat_star_sparse_reg_DIMS_greedy.m

@@ -13,14 +13,13 @@ for k=1:length(node(i).pi)
         %%%%averaging
         sp=size(ps);
         for mn=1:sp(1)
-            
             T_i=0;
             M_i=0;
             pa=ones(t,1);
             ki=1;
             for l=node(i).subsets{k}
                 mu_p=MU{l};
-                pa=pa.*mu_p(1:t,ps(mn,ki));              
+                pa=pa.*mu_p(1:t,ps(mn,ki));
                 ki=ki+1;
             end
             for d=1:D
@@ -32,8 +31,8 @@ for k=1:length(node(i).pi)
                     end
                 end
             end
-            T(mn,:)=T_i(:);
-            M(mn,:,:)=M_i(:,:);
+            T(mn,:)=T_i;
+            M(mn,:,:)=M_i;
         end
     else
         for d=1:D

ssl-ctbn-code/ctbns_base/CTBN_cond_stat_star_sparse_reg_DIMS_greedy_approx.m

+function [T_k,M_k] = CTBN_cond_stat_star_sparse_reg_DIMS_greedy_approx(node,i,t,dt,MU,RHO)
+%CALCULATE COND STATISTICS OF CTBN
+mu=MU{i};
+rho=RHO{i};
+D=node(i).D;
+%%%get effective rate
+for k=1:length(node(i).pi)
+    ps=node_states(node,node(i).subsets{k});
+    sp=size(ps);
+    T=zeros(sp(1),D);
+    M=zeros(sp(1),D,D);
+    if (isempty(node(i).subsets{k})==0)
+        %%%%averaging
+        sp=size(ps);
+        for mn=1:sp(1)
+            T_i=0;
+            M_i=0;
+            pa=ones(t,1);
+            ki=1;
+            for l=node(i).subsets{k}
+                mu_p=MU{l};
+                pa=pa.*mu_p(1:t,ps(mn,ki));
+                ki=ki+1;
+            end
+            for d=1:D
+                T_i(d)=trapz(linspace(0,dt*t,t),mu(1:t,d).*pa);
+                for d_=1:D
+                    if (d~=d_)
+                        q_gm=ctbn_gen_gm_rate_sparse_reg_DIMS_partial(node,MU,i,ps(mn,:),node(i).subsets{k},d,d_);
+                        M_i(d,d_)=trapz(linspace(0,dt*t,t),q_gm(1:t).*pa(1:t).*mu(1:t,d).*(rho(1:t,d_))./(rho(1:t,d)));
+                    end
+                end
+            end
+            T(mn,:)=T_i;
+            M(mn,:,:)=M_i;
+        end
+    else
+        for d=1:D
+            T_i(d)=trapz(linspace(0,dt*t,t),mu(1:t,d));
+            for d_=1:D
+                if (d~=d_)
+                    q_gm=ctbn_gen_gm_rate_sparse_reg_DIMS_partial(node,MU,i,[],[],d,d_);
+                    M_i(d,d_)=trapz(linspace(0,dt*t,t),mu(1:t,d).*q_gm(1:t).*(rho(1:t,d_)./rho(1:t,d)));
+                end
+            end
+        end
+        T(1,:,:)=T_i;
+        M(1,:,:)=M_i;
+    end
+    T_k{k}=T;
+    M_k{k}=M;
+end
+

ssl-ctbn-code/ctbns_base/createMotherCTBN_DIMS.m

@@ -71,6 +71,7 @@ for i=1:L
     
     sps=size(P);
     node(i).pi=pi_0.*ones(1,sps(2));
+    node(i).pi(end)=1;
     node(i).pi=node(i).pi/sum(node(i).pi);
 end
 

ssl-ctbn-code/em/expectation/expectation_greedy/ctbn_expectation_sparse_reg_par_DIMS_greedy_closure.m

+function [MU,RHO,node] = ctbn_expectation_sparse_reg_par_DIMS_greedy_closure(node,dt,Mmax,t0,Z,TZ,thresh)
+%UNTITLED7 Summary of this function goes here
+%   Detailed explanation goes here
+M=cell(1,length(TZ));
+T=cell(1,length(TZ));
+MU=cell(1,length(TZ));
+RHO=cell(1,length(TZ));
+
+for k=1:length(TZ)
+    MU{k}=[];
+    RHO{k}=[];
+    for i=1:length(node)
+        for l=1:length(node(i).pi)
+            T{k}{i}{l}=0;
+            M{k}{i}{l}=0;
+        end
+    end
+end
+parfor k=1:length(TZ)
+    try
+        [mu,rho,~] = P_CVMCTBN_STAR_POST_REG_DIMS_greedy_closure(node,dt,Mmax,t0,Z{k},TZ{k},thresh);
+        MU{k}=mu{Mmax};
+        RHO{k}=rho{Mmax-1};
+        TK=linspace(0,TZ{k}(end),ceil(TZ{k}(end)/dt));
+        for i=1:length(node)
+            [T0,M0] = CTBN_cond_stat_star_sparse_reg_DIMS_greedy_approx(node,i,length(TK)-2,dt,MU{k},RHO{k});
+            for l=1:length(node(i).pi)
+                T{k}{i}{l}=T{k}{i}{l}+T0{l};
+                M{k}{i}{l}=M{k}{i}{l}+M0{l};
+            end
+        end
+    catch
+        MU{k}=[];
+        RHO{k}=[];
+        sprintf('Warning: Could not process trajectory %d',k)
+    end
+end
+for i=1:length(node)
+    for k=1:length(TZ)
+        for l0=1:length(node(i).pi)
+            if T{k}{i}{l0}~=0
+                st=size(T{k}{i}{l0});
+                sm=size(M{k}{i}{l0});
+            end
+        end
+    end
+    for l0=1:length(node(i).pi)
+        T0=zeros(st(1),st(2));
+        M0=zeros(sm(1),sm(2),sm(3));
+        for k_=1:length(TZ)
+            T0=T0+T{k_}{i}{l0};
+            M0=M0+M{k_}{i}{l0};
+        end
+        node(i).trans_k{l0}=M0;
+        node(i).dwell_k{l0}=T0;
+    end
+end
+end
+

ssl-ctbn-code/em/maximization/estimate_pi_sparse_reg_par_DIMS_grad_lam.m

+function [lam] = estimate_pi_sparse_reg_par_DIMS_grad_lam(node,lam,restarts)
+%UNTITLED6 Summary of this function goes here
+%   Detailed explanation goes here
+for i=1:length(node)
+    fun=@(x) pi_llh_DIMS_lam(node,i,x,node(i).pi);
+    x0=rand;
+    A=[];
+    b=[];
+    Aeq=[];
+    beq=[];
+    lb=-10;
+    ub=10;
+    options = optimoptions('fmincon','Display','off','SpecifyObjectiveGradient',false,'Algorithm','sqp');
+    lami = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,[],options);
+    lam(i) = lami;
+end
+
+end
+

ssl-ctbn-code/em/maximization/pi_llh_DIMS_lam.m

+function [F] = pi_llh_DIMS_lam(node,i,lam,pi)
+%UNTITLED5 Summary of this function goes here
+%   Detailed explanation goes here
+node(i).pi=pi;
+%[node] = ctbn_summarize_stats(node);
+%[node] = ctbn_compute_post_rates(node);
+F = -marg_llh_sparse_reg_noexpln_single_DIMS(node,i,lam);
+
+end
+

ssl-ctbn-code/sampler/ctbn_sampler/corrupted_observation_gaussianD.m

@@ -18,7 +18,7 @@ for j=1:N_TRAJ
             Y(d_,:)=normpdf(D_h,d_,SIGMA);
         end
         Z=round(Y,5);
-        Y(Z==0)=10^(-4);
+        Y(Z==0)=10^(-2);
         Y=Y./sum(Y,1);
         D{j}{i}=D_h;
         DATAC{j}{i}=Y;

ssl-ctbn-code/scripts/random_graph_experiment_script.m

@@ -3,7 +3,7 @@ addpath(genpath('./ssl-ctbn-code'))
 L=5; %number of nodes of random graph
 max_par=2;
 num_graphs=30;
-N_TRAJ=10; %number of synthetic trajectories
+N_TRAJ=40; %number of synthetic trajectories
 mworkers=4;
 for graphs=1:num_graphs
     A=zeros(L);

ssl-ctbn-code/scripts/run_greedy_ctbn_ssl.m

@@ -6,9 +6,9 @@
 %        F    : value of objective function at different iterations
 addpath(genpath('./ssl-ctbn-code'))
 %maximal number of parents in greedy search
-K=2;
+K=3;
 %name of experiment
-name=sprintf('my_test_run_%d_%dparents',1,K);
+name=sprintf('my_t0est_run_%d_%dparents',1,K);
 %number of workers running in parallel
 mworkers=4;
 ctbn_gradient_structure_learning_dims_greedy(name,mworkers,K)

ssl-ctbn-code/var_inference/P_CVMCTBN_STAR_POST_REG_DIMS_greedy_closure.m

+function [MU,RHO,F,m] = P_CVMCTBN_STAR_POST_REG_DIMS_greedy_closure(node,dt,M,t0,Z,TZ,thresh)
+%%%%%%%%%%%%%CLUSTER VARIATIONAL STAR-APPROXIMATION FOR CTBNs WITH ERROR MODEL%%%%%
+%%%%%%INPUT:%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%node is NODE STRUCT. DEFINING CTBN
+%D is local dim. (only D=2 possible a.t.m.)
+%T is SIMULATION TIME
+%dt is time-step
+%M is max. number of iterations for iterative ODE solver
+%X0 is initial condition
+%Z is noisy observation of state
+%TZ is time of observation
+%%%%%%OUTPUT:%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%mu(m,t,i,d) is MARGINAL PROB. IN m'th IERATION OF NODE i TO BE IN STATE d
+%AT TIME t
+%rho(m,t,i,d) is LAGRANGE MULTIPLIER " " " " " " " " ".
+%qt(t,i,d) is AVERAGED RATE """""""
+%pst(t,i,d) is AVERAGED CHILD RESPONSE " " " " ".
+
+%%%%%%%%%%%READ OUR SIM. PARAM.%%%%%%%%%%%%%%%%%%
+
+L=length(node);
+T=TZ(end)+dt*t0;
+tau=ceil(TZ(end)/dt)+t0;
+e0=squeeze(Z(:,1,:));
+
+
+xt=linspace(0,T,tau);
+%xt2=linspace(0,T-dt,tau-1);
+
+%%%%%%%%%%SET INITIAL CONDITIONS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%CREATE ANSATZ TRAJECTORY BY PICKING RND CIM%%%
+for i=1:L
+    D=node(i).D;
+    mu_init=zeros(tau,D);
+    for d=1:D
+        mu_init(1,d)=1/D;
+    end
+    
+    TSPAN = linspace(0,T,tau);
+    %%%PICK RANDOM CIM FOR ANSATZ SOLUTION
+    q=zeros(tau,D,D);
+    Q_i=node(i).cellOfCondRM;
+   % c=node(i).allPosStatesOfParents;
+   % sc=size(c);
+    Q=Q_i{1};
+    
+    for d=1:D
+        for d_=1:D
+            q(:,d,d_)=Q(d,d_).*ones(tau,1);
+        end
+    end
+    
+    [~,Y] = ode15s(@(t,y) CVM_CTBN_mu_fastD(t, y,D, xt, q(1:tau,:,:), xt, ones(tau,D)), TSPAN,  mu_init(1,:)');
+    mu_init(1:tau,:)=Y(1:tau,:);
+    for m=1:M
+      MU{m}{i}=mu_init(1:tau,:);
+      RHO{m}{i}=zeros(tau,D);
+    end
+end
+
+
+
+%%%%%%%%%%ACTUAL SIMULATION OF CTBN%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+TSPAN_R = linspace(T,0,tau);
+TSPAN_M = linspace(0,T,tau);
+%options=odeset('RelTol',1e-12,'AbsTol',1e-13);
+%%%%%%%%%%%%%%%%%%%%%%%%ITERATE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+for m=1:M-1
+    %%%%%%%%%%UPDATE LAGRANGE MULT. GIVEN MARG. PROBS%%%%%%%%%%%%%%%%%%%%%%
+    for i=1:L
+        D=node(i).D;
+        psi0=zeros(tau,D);
+        mu=MU{m}{i};
+        rho=RHO{m}{i};
+        
+        [q_am,q_gm] = P_CVMCTBN_EFF_RATES_SPARSE_REG_DIMS_greedy_approx(node,i,m,MU);
+        
+        if m>1
+             %psi0 = P_CVMCTBN_COUP_D(node,i,squeeze(mu(m,:,:,:)),squeeze(rho(m-1,:,:,:)),D);
+            % psi0 = P_CVMCTBN_COUP_DIMS(node,i,m,MU,RHO);
+            [psi0] = P_CVMCTBN_COUP_SPARSE_REG_DIMS_greedy(node,i,m,MU,RHO);
+        end
+      
+        rhoT=ones(1,D);
+        TSPAN_Y=linspace(T,TZ(end),t0);
+        %%%%%%%%%%ACTUAL TIME EVOLUTION OF LAGR. MULT. (BACKWARDS IN TIME)%
+        %options = odeset('Jacobian',@(t,y)J_CVM_CTBN_mu(D,t, xt, q(1:tau,:), xt, squeeze(rho(m,1:tau,i,:))));
+        [~, Y] = ode15s(@(t,y) CVM_CTBN_rho_fast_sparse_reg_DIMS(D,t, y, xt, q_am, q_am, psi0), TSPAN_Y, rhoT');
+        R{length(TZ)+1}=Y;
+        
+        for k=length(TZ):-1:1
+            Y_mem=Y;
+            rhoT=Y(end,:).*(squeeze(Z{i}(:,k))')/sum(Y(end,:).*(squeeze(Z{i}(:,k))'));
+            if k>1
+                TSPAN_Y=linspace(TZ(k),TZ(k-1),ceil((TZ(k)-TZ(k-1))/dt));
+            else
+                TSPAN_Y=linspace(TZ(1),0,ceil(TZ(1)/dt));
+            end
+            
+            %options = odeset('Jacobian',@(t,y)J_CVM_CTBN_mu(D,t, xt, q(1:tau,:), xt, squeeze(rho(m,1:tau,i,:))));
+            [~, Y] = ode15s(@(t,y) CVM_CTBN_rho_fast_sparse_reg_DIMS(D,t, y, xt, q_am,q_gm, psi0), TSPAN_Y, rhoT');
+            R{k}=Y;
+            [msglast, msgidlast] = lastwarn;
+            if isempty(msglast)==0
+                Z{i}(:,k)=1; %remove data-point
+                warning('') %clear last warning
+                msglast=[];
+                sprintf('Warning: Could not process data-point %d of node %d',k,i)
+                
+                rhoT=Y(end,:).*(squeeze(Z{i}(:,k))')/sum(Y(end,:).*(squeeze(Z{i}(:,k))'));
+                if k>1
+                    TSPAN_Y=linspace(TZ(k),TZ(k-1),ceil((TZ(k)-TZ(k-1))/dt));
+                else
+                    TSPAN_Y=linspace(TZ(1),0,ceil(TZ(1)/dt));
+                end           
+                %options = odeset('Jacobian',@(t,y)J_CVM_CTBN_mu(D,t, xt, q(1:tau,:), xt, squeeze(rho(m,1:tau,i,:))));
+                [Tt Y] = ode15s(@(t,y) CVM_CTBN_rho_fast_sparse_reg_DIMS(D,t, y, xt, q_am,q_gm, psi0), TSPAN_Y, rhoT');
+                R{k}=Y;
+            end
+        end
+        for d=1:D
+            A=[];
+            for k=length(TZ)+1:-1:1
+                A=[A ;R{k}(:,d)];
+            end
+           % rho(m,tau:-1:1,i,d)=A(1:tau);
+            rho(tau:-1:1,d)=A(1:tau);
+        end
+        RHO{m+1}{i}=rho;
+        %%%%%%%%%%UPDATE MARG. PROBS GIVEN LAGRANGE MULT. %%%%%%%%%%%%%%%%%
+        %%%%%%%%%%ACTUAL TIME EVOLUTION OF MARG. PROB. (FORWARDS IN TIME)%%
+        %  options = odeset('Jacobian',@(t,y)J_CVM_CTBN_mu(D,t, xt, q(1:tau,:), xt, squeeze(rho(m,1:tau,i,:))));
+       % [Tt Y] = ode15s(@(t,y) CVM_CTBN_mu_fastD(D,t, y, xt, q(1:tau,:,:), xt, RHO{m}{i}, TSPAN_M, mu(m,1,i,:)));
+         [~, Y] = ode15s(@(t,y) CVM_CTBN_mu_fastD(t, y,D, xt, q_gm(1:tau,:,:), xt, rho(1:tau,:)), TSPAN_M, squeeze(mu(1,:)));
+     
+       %mu(m+1,1:tau,i,:)=Y(1:tau,:);
+        MU{m+1}{i}=Y(1:tau,:);
+        
+        
+        if m==M-1
+             node(i).R_am=q_am(1:tau,:,:);
+        end
+        
+    end
+    
+    %%%%%CALCULATE LIKELIHOOD LOWER BOUND TO CHECK FOR CONVERGENCE
+%     if m>1
+%         [F(m),~,~] = CVM_CTBN_likelihood_starD(node,squeeze(mu(m,:,:,:)),squeeze(rho(m-1,:,:,:)),dt);
+%     end
+%     if m>5
+%         stdF=std(F(end:-1:end-4));
+%         if stdF<thresh
+%             break
+%         end
+%         
+%     end
+    
+
+end
+%%%%CONVERGED SOLUTIONS
+%MU=squeeze(mu(m,1:tau-t0,:,:));
+%RHO=squeeze(rho(m-1,1:tau-t0,:,:));
+F=0;
+end
+
+
+