

function [Run] = PROC_cpBGe_regularized(data, steps, p, k)

global DATA_ALL;

[DATA_ALL] = SHIFT_DATA(data);

global T_0;
global my_0;
global v;
global alpha;
global Prior;
global g_vec;
global g0_vec;
global G_vec;
global G0_vec;

[T_0, my_0, v, alpha, Prior, g_vec, g0_vec, G_vec, G0_vec] = SET_HYPERPARAMETERS(DATA_ALL, p,k);

[Run] = INITIALISE_cpBGe_regularized(DATA_ALL, steps);
[Run] = START_cpBGe_regularized(DATA_ALL, steps, Run);

return;

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

function [DATA_SHIFTED] = SHIFT_DATA(data)

[n, t] = size(data);

n_plus = n+1; % Covariance matrices will be of size (n+1)x(n+1)

for node_i=1:n
            % For each variable X_i consider
            % data_1[X_1(t-1),...,X_(n)(t-1),...,X_1(t-slice),...,X_(n)(t-slice), X_(i)(t)]      
            obs                  = 1:(t-1);  
            data_new             = zeros(n_plus,t-1);    
            data_new(1:n,obs)    = data(1:n,1:(t-1));       
            data_new(n+1,obs)    = data(node_i,2:t);   
            DATA_SHIFTED{node_i} = data_new;      
end

return;

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


function [T_0, my_0, v, alpha, Prior, g_vec, g0_vec, G_vec, G0_vec] = SET_HYPERPARAMETERS(DATA_ALL, p_transition, k_success)

[n_plus, n_obs] = size(DATA_ALL{1});

n_nodes = n_plus-1;

Prior = zeros(n_nodes+1,1); 

% THESE VALUES MAY BE MODIFIED BY A USER: v, alpha, my_0, and B 
v     = 1;
alpha = n_plus+2;
my_0  = 0 * ones(n_plus,1);
v_vec = 1 * ones(n_plus,1); 
B     = zeros(n_plus,n_plus);

% COMPUTE PRECISION-MATRIX W
W =  1/v_vec(1);

for i = 1:(n_plus-1)
    b_vec = B(1:i,(i+1));
    W = [W + (b_vec * b_vec')/v_vec(i+1) , (-1)*b_vec/v_vec(i+1) ; (-1)*b_vec'/v_vec(i+1), 1/v_vec(i+1)];
end

% COMPUTE T_0: 
T_0   = (v*(alpha-n_plus-1))/(v+1) * inv(W);   

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

g_vec   = zeros(1,100);
g0_vec  = zeros(1,100);

for t=k_success:100
    g_vec(t) = prod(1:(t-1))/(prod(1:(k_success-1))*prod(1:(t-k_success))) * p_transition^(k_success) * (1-p_transition)^(t-k_success);    
end

for t=1:100
    g0_vec_t = 0;
    for i=1:min([t,k_success])
       g0_vec_t = g0_vec_t + (1/k_success) * prod(1:(t-1))/(prod(1:(i-1))*prod(1:(t-i))) * p_transition^i * (1-p_transition)^(t-i);
    end
    g0_vec(1,t) = g0_vec_t;
end

G_vec = cumsum(g_vec);

G0_vec = cumsum(g0_vec);
 
return;

function [Run] = INITIALISE_cpBGe_regularized(DATA_ALL, steps)

[n_plus, m]= size(DATA_ALL{1});

n_nodes = length(DATA_ALL);

DAG = zeros(n_nodes,n_nodes);

MATRIX = ones(n_nodes,m);

NODE_VEC = 1:n_nodes;

for node_i=1:n_nodes   
    data = DATA_ALL{node_i};
    node_i_comp = NODE_VEC(node_i);
    for component=1:max(MATRIX(node_i_comp,:))      
        DATA{node_i}{component} = data(:,find(MATRIX(node_i_comp,:)==component));
    end  
end

[log_score] = COMPUTE_LOG_SCORE_MIX(DATA, DAG, MATRIX, NODE_VEC);

for i=1:steps+1
    Run.dag{i}         = 0;
    Run.matrix{i}      = 0;
    Run.log_score(i)   = 0;
    Run.Node_Vec{i}    = 0;
end    
    
Run.dag{1}         = DAG;
Run.matrix{1}      = MATRIX;
Run.Node_Vec{1}    = NODE_VEC;
Run.steps(1)       = 1;
Run.log_score      = log_score;

return;


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

function [log_score] = COMPUTE_LOG_SCORE_MIX(DATA, DAG, MATRIX, NODE_VEC)

% The hyperparameters:
global T_0;
global v;
global alpha;
global my_0;
global g_vec;
global g0_vec;
global G_vec;
global G0_vec;

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

n_nodes = length(DAG);

counts         = hist(NODE_VEC,1:n_nodes);
ind_non_empty  = find(counts);

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

G_vec_zero = [0,G_vec];

[n_nodes, m]= size(MATRIX);


log_prob_breaks = 0;

for ind=ind_non_empty
        
    break_vec         = find(MATRIX(ind,2:end) - MATRIX(ind,1:(end-1)));

    if (length(break_vec)==0)
       log_prob_breaks_i = log(1-G0_vec(m-1)); 
       
    elseif (length(break_vec)==1)
        log_prob_breaks_i = log(g0_vec(break_vec(1))) + log(1-G_vec_zero(m-break_vec(1)));       
    else 
        log_prob_breaks_i = log(g0_vec(break_vec(1))) + log(1-G_vec_zero(m-break_vec(end)));
        for i=2:(length(break_vec))
            log_prob_breaks_i = log_prob_breaks_i + log(g_vec(break_vec(i)-break_vec(i-1)));
        end
    end
   
    log_prob_breaks = log_prob_breaks + log_prob_breaks_i ;
end


log_prob_data = 0;

% Compute the local score for each node:

for i_node=1:n_nodes
    
    k_i = length(DATA{i_node}); % k_i is the number of mixture components for the i-th node
    
    log_score_node_i = 0;
    
    for component=1:k_i
    
        [n_plus, n_obs_in_component] =  size(DATA{i_node}{component});

        if(n_obs_in_component==0)
                % do nothing 
        else
                data = DATA{i_node}{component};
    
                % Compute T_m for the current combination of node and
                % component:
                
                parents_of_node  = find(DAG(:,i_node));
                parents_and_node = [parents_of_node; n_plus];
                
                S_m = (n_obs_in_component-1) * cov(data(parents_and_node,:)');
  
                T_m = T_0(parents_and_node,parents_and_node) + S_m + (v*n_obs_in_component)/(v+n_obs_in_component) * (my_0(parents_and_node,1) - mean(data(parents_and_node,:)',1)') * (my_0(parents_and_node,1) - mean(data(parents_and_node,:)',1)')';  

                [log_score_local_node_nominator]   = Gauss_Score_complete(length(parents_and_node), n_obs_in_component, v, alpha, T_0(parents_and_node,parents_and_node), T_m);
                [log_score_local_node_denominator] = Gauss_Score_complete(length(parents_of_node),  n_obs_in_component, v, alpha, T_0(parents_of_node,parents_of_node),   T_m(1:(end-1),1:(end-1)));
                 
                log_score_node_i = log_score_node_i + (log_score_local_node_nominator - log_score_local_node_denominator);
        end
        
    end

    log_prob_data = log_prob_data + log_score_node_i;
    
end

log_score = log_prob_breaks + log_prob_data;

return

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

function [log_Score_i] = Gauss_Score_complete(n_nodes_in_sub, n_obs, v, alpha, T_0_sub, T_m_sub)

if n_nodes_in_sub==0
    log_Score_i = 0;
else
      
    sum_1 = (-1)*(n_nodes_in_sub)*(n_obs)/2 * log(2*pi) + (n_nodes_in_sub/2) * log(v/(v+n_obs));

    sum_3 = (alpha/2) * log(det(T_0_sub)) + (-1)*(alpha+n_obs)/2 * log(det(T_m_sub));

    log_value_a = (alpha         *n_nodes_in_sub/2)*log(2) + (n_nodes_in_sub*(n_nodes_in_sub-1)/4)*log(pi);
    log_value_b = ((alpha+n_obs) *n_nodes_in_sub/2)*log(2) + (n_nodes_in_sub*(n_nodes_in_sub-1)/4)*log(pi);
 
    for i=1:n_nodes_in_sub
        log_value_a = log_value_a + gammaln(( alpha   +1-i)/2);
        log_value_b = log_value_b + gammaln(((alpha+n_obs)+1-i)/2);
    end
   
    log_Score_i = sum_1 + (-1)*(log_value_a - log_value_b) + sum_3; 

end

return 
    
function [Run] = START_cpBGe_regularized(DATA_ALL, steps, Run)
                     
last           = Run.steps(1);          % index of last sampled dags
DAG            = Run.dag{last};         % the current DAG
MATRIX         = Run.matrix{last};      % the current allocation vector
NODE_VEC       = Run.Node_Vec{last}; 

% Output on screen: 
fprintf('\n###########################################################\n')
fprintf('The Gibbs simulation for the regularized cpBGe model has been started \n')
fprintf('###########################################################\n')

CONSTANT       = exp(300);
n_nodes        = length(DAG);
counts         = hist(NODE_VEC,1:n_nodes);
ind_non_empty  = find(counts);

DAG_full = ones(n_nodes,n_nodes);

for ind=ind_non_empty
    nodes_in_comp            = find(NODE_VEC==ind); 
    [Q_VEC]                  = Dynamic_Q_MIX(DATA_ALL, nodes_in_comp, DAG_full, CONSTANT);
    [Vector]                 = Simulate_changepoints_MIX(Q_VEC, DATA_ALL, nodes_in_comp, DAG_full, CONSTANT);
    MATRIX(ind,:)            = Vector;
end 





%%% Start of the MCMC simulation
for i = last:(steps)
    i
    [DAG, MATRIX, NODE_VEC, log_score] = MCMC_INNER_MIX(DAG, MATRIX, NODE_VEC, DATA_ALL, CONSTANT);
    Run.dag{i+1}             = DAG;
    Run.matrix{i+1}          = MATRIX;
    Run.Node_Vec{i+1}        = NODE_VEC;
    Run.steps(1)             = i+1;
    Run.log_score(i+1)       = log_score;
end

return;

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

function [DAG, MATRIX, NODE_VEC, log_score] = MCMC_INNER_MIX(DAG, MATRIX, NODE_VEC, DATA_ALL, CONSTANT)

[thrash, n_nodes] = size(DAG);

for child_node=1:n_nodes
        
        data = DATA_ALL{child_node};
       
        clear DATA;
        
        child_node_comp = NODE_VEC(child_node);
        
        for component=1:max(MATRIX(child_node_comp,:))      
            DATA{child_node}{component} = data(:,find(MATRIX(child_node_comp,:)==component));
        end
               
        [SCORES] = SCORE_ALL_PARENTS(DATA, child_node, n_nodes);
                
        cumulative = cumsum(SCORES.exp_scores)/sum(SCORES.exp_scores);
               
        x = rand;
            
        indicis = find(cumulative>=x);
        index   = indicis(1); 
       
        new_parents_of_node = nonzeros(SCORES.parents(index,:));
       
        DAG(:,child_node)=0;  
        DAG(new_parents_of_node,child_node)=1;
                  
        clear DATA;
end 
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
         
counts         = hist(NODE_VEC,1:n_nodes);
ind_non_empty  = find(counts);

for ind=ind_non_empty
    nodes_in_comp            = find(NODE_VEC==ind); 
    [Q_VEC]                  = Dynamic_Q_MIX(DATA_ALL, nodes_in_comp, DAG, CONSTANT);
    [Vector]                 = Simulate_changepoints_MIX(Q_VEC, DATA_ALL, nodes_in_comp, DAG, CONSTANT);
    MATRIX(ind,:)            = Vector;
end 

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

% BIRTH/DEATH MOVES

P_Death=0.5;
    
for steps=1:10
    
counts   = hist(NODE_VEC,1:n_nodes);
ind_eq1  = find(counts==1);
ind_gr1  = find(counts>1);
ind_geq1 = find(counts>=1);

n_eq1       = length(ind_eq1);
n_gr1       = length(ind_gr1);
n_geq1      = length(ind_geq1);
    
    x = rand;
    
    if (x<=P_Death) % -> Death move
      
        if (n_eq1==0)
            % do nothing
        else
        
        indicis_1 = randperm(n_eq1); % the single node with label_1 WILL GET label_2 -> i.e. the component with label_1 "dies"
        label_1   = ind_eq1(indicis_1(1));
        
        indicis_2 = randperm(n_geq1);
        label_2   = ind_geq1(indicis_2(1));
        
        if (label_1==label_2)
            label_2 = ind_geq1(indicis_2(2)); % this ensures that label_1 is different from label_2
        end
        
        single_node_with_label_1 = find(NODE_VEC==label_1); % the node with label_1
        nodes_with_label_2       = find(NODE_VEC==label_2); % the node(s) with label_2
       
        NODE_VEC_NEW                           = NODE_VEC; 
        NODE_VEC_NEW(single_node_with_label_1) = label_2; % the single node with label_1 has just obtained label_2 
       
        nodes_with_label_2_new = find(NODE_VEC_NEW==label_2); % the nodes with label_2 after death move
        
        [Q_new]    = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
        [Vector_2] = Simulate_changepoints_MIX(Q_new, DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
        
        MATRIX_NEW            = MATRIX;
        MATRIX_NEW(label_2,:) = Vector_2;
              
      
        
        counts_new   = hist(NODE_VEC_NEW,1:n_nodes);
        ind_gr1_new  = find(counts_new>1);
        n_gr1_new    = length(ind_gr1_new);
       
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [Q_old1]   = Dynamic_Q_MIX(DATA_ALL, single_node_with_label_1, DAG, CONSTANT);
        [Q_old2]   = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2, DAG, CONSTANT);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        log_forward  = log(1/n_eq1)     + log(1/(n_geq1-1))                      - log(Q_new(1)); 
        log_backward = log(1/n_gr1_new) + log(1/length(nodes_with_label_2_new))  - log(Q_old1(1)) - log(Q_old2(1)) + log(CONSTANT);
        log_hastings = log_backward - log_forward; 
        
        A = exp(log_hastings);
        
        y = rand;
        

        if (y<A)
           MATRIX    = MATRIX_NEW;
           NODE_VEC  = NODE_VEC_NEW;
        end
        

        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    else  % -> Birth move
        
        if(n_gr1==0)
            % do nothing
        else
            
        indicis_1 = randperm(n_gr1); % one node with label_1, namely node_out, WILL GET the label_2 
        label_1   = ind_gr1(indicis_1(1));
        
        ind_empty  = find(counts==0);
        label_2    = ind_empty(1);

        nodes_with_label_1_old   = find(NODE_VEC==label_1);
        n_nodes_with_label_1_old = length(nodes_with_label_1_old);
        indicis                  = randperm(n_nodes_with_label_1_old);
        node_out                 = nodes_with_label_1_old(indicis(1));
        
        NODE_VEC_NEW           = NODE_VEC;
        NODE_VEC_NEW(node_out) = label_2;
         
        nodes_with_label_1_new = find(NODE_VEC_NEW==label_1);
        
        [Q_new1]    = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
        [Vector_1]  = Simulate_changepoints_MIX(Q_new1, DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
    
        
        [Q_new2]    = Dynamic_Q_MIX(DATA_ALL, node_out, DAG, CONSTANT);
        [Vector_2]  = Simulate_changepoints_MIX(Q_new2, DATA_ALL, node_out, DAG, CONSTANT);
        
        
        MATRIX_NEW             = MATRIX;
        MATRIX_NEW(label_1,:)  = Vector_1;
        MATRIX_NEW(label_2,:)  = Vector_2;
        
       
        counts_new   = hist(NODE_VEC_NEW,1:n_nodes);
        ind_eq1_new  = find(counts_new==1);
        n_eq1_new    = length(ind_eq1_new);
        
        ind_geq1_new  = find(counts_new>=1);
        n_geq1_new    = length(ind_geq1_new);
        
     
        [Q_old1]   = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_old, DAG, CONSTANT);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        log_forward  = log(1/n_gr1)     + log(1/n_nodes_with_label_1_old)  - log(Q_new1(1)) - log(Q_new2(1)) + log(CONSTANT); 
        log_backward = log(1/n_eq1_new) + log(1/n_geq1_new)                - log(Q_old1(1));
        log_hastings = log_backward - log_forward; 
        
        A = exp(log_hastings);

        y = rand;
        
        if (y<A)
           MATRIX    = MATRIX_NEW;
           NODE_VEC  = NODE_VEC_NEW;
        end
        
        end
    
    end
    
   
    % RE-ALLOCATION MOVE
    
 
    counts   = hist(NODE_VEC,1:n_nodes);

    ind_gr1  = find(counts>1);
    ind_geq1 = find(counts>=1);

    n_gr1       = length(ind_gr1);
    n_geq1      = length(ind_geq1);
  
    if (n_gr1==0 | n_geq1==1)
        % do nothing
    else
     indicis_1 = randperm(n_gr1); % one node with label_1, namely node_out, WILL GET the label_2 
     label_1   = ind_gr1(indicis_1(1));
        
     indicis_2  = randperm(n_geq1);
     label_2    = ind_geq1(indicis_2(1));
     
     if (label_1==label_2)
         label_2 = ind_geq1(indicis_2(2));
     end

     nodes_with_label_1_old   = find(NODE_VEC==label_1);
     nodes_with_label_2_old   = find(NODE_VEC==label_2);
     
     n_nodes_with_label_1_old = length(nodes_with_label_1_old);
     indicis                  = randperm(n_nodes_with_label_1_old);
     node_out                 = nodes_with_label_1_old(indicis(1));
        
     NODE_VEC_NEW           = NODE_VEC;
     NODE_VEC_NEW(node_out) = label_2;
         
     nodes_with_label_1_new = find(NODE_VEC_NEW==label_1);
     nodes_with_label_2_new = find(NODE_VEC_NEW==label_2);
     
     [Q_new1]    = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
     [Vector_1]  = Simulate_changepoints_MIX(Q_new1, DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
       
     [Q_new2]    = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
     [Vector_2]  = Simulate_changepoints_MIX(Q_new2, DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
        
        
     MATRIX_NEW             = MATRIX;
     MATRIX_NEW(label_1,:)  = Vector_1;
     MATRIX_NEW(label_2,:)  = Vector_2;
        
   
        
     counts_new   = hist(NODE_VEC_NEW,1:n_nodes);
     ind_gr1_new  = find(counts_new>1);
     n_gr1_new    = length(ind_gr1_new);
        
     ind_geq1_new  = find(counts_new>=1);
     n_geq1_new    = length(ind_geq1_new);
        
    
     
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     [Q_old1]   = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_old, DAG, CONSTANT);
     [Q_old2]   = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2_old, DAG, CONSTANT);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     log_forward  = log(1/n_gr1)     + log(1/(n_geq1-1))    - log(Q_new1(1))  - log(Q_new2(1)); 
     log_backward = log(1/n_gr1_new) + log(1/(n_geq1_new-1))- log(Q_old1(1))  - log(Q_old2(1));
     log_hastings = log_backward - log_forward; 
     
        
     A = exp(log_hastings);
      

     y = rand;
        
     if (y<A)
        MATRIX    = MATRIX_NEW;
        NODE_VEC  = NODE_VEC_NEW;
     end
        
     end
       
end
    
   
DATA = [];   

for node=1:n_nodes  
    node_label = NODE_VEC(node);
    data = DATA_ALL{node};
    for component=1:max(MATRIX(node_label,:))      
        DATA{node}{component} = data(:,find(MATRIX(node_label,:)==component));
    end  
end

[log_score] = COMPUTE_LOG_SCORE_MIX(DATA, DAG, MATRIX, NODE_VEC);
clear DATA

        
return

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

function [SCORES] = SCORE_ALL_PARENTS(DATA, i_node, n_nodes)

Counter = 1;    
Structure.parents(Counter,1:3) = [0 0 0];
log_likelihood  = COMPUTE_LIKELIHOOD_LOCAL_2(DATA, [], i_node); 
Structure.scores(Counter) = log_likelihood;
Counter = Counter +1;

% Local Prior scores of parent-sets of size 1:
    for parent_1 = 1:n_nodes    
        
            Structure.parents(Counter,1:3) = [0, 0, parent_1];
            log_likelihood            = COMPUTE_LIKELIHOOD_LOCAL_2(DATA, [parent_1], i_node);
            Structure.scores(Counter) = log_likelihood;
            Counter = Counter +1;
    end
       
% Local Prior scores of parent-sets of size 2:
    for parent_1 = 1:(n_nodes-1)
        for parent_2 = (parent_1+1):n_nodes
           
                Structure.parents(Counter,1:3) = [0, parent_1, parent_2];
                log_likelihood            = COMPUTE_LIKELIHOOD_LOCAL_2(DATA, [parent_1;parent_2], i_node);
                Structure.scores(Counter) = log_likelihood; 
                Counter = Counter +1;
        
        end
    end

% LocalPrior scores of parent-sets of size 3: 

for parent_1 = 1:(n_nodes-2)
     for parent_2 = (parent_1+1):(n_nodes-1)
       for parent_3 = (parent_2+1):n_nodes    
           
                Structure.parents(Counter,1:3) = [parent_1, parent_2, parent_3];
                log_likelihood            = COMPUTE_LIKELIHOOD_LOCAL_2(DATA, [parent_1;parent_2;parent_3], i_node);
                Structure.scores(Counter) = log_likelihood;
                Counter = Counter +1;
          
       end
     end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sort with respect to the local Bayesian scores
[thrash, indicis]    = sort(Structure.scores,'descend');

SCORES.parents       = Structure.parents(indicis,:);
SCORES.exp_scores    = exp(Structure.scores(indicis));

clear Structure;


return

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

function [log_score] = COMPUTE_LIKELIHOOD_LOCAL_2(DATA, parents_of_node, i_node)

% Some hyperparameters:
global T_0;
global v;
global alpha;
global my_0;

% Compute the local score for each node:
k_i = length(DATA{i_node}); % k_i is the number of mixture components for the i-th node

log_score_node_i = 0;
    
for component=1:k_i
    
    [n_plus, n_obs_in_component] =  size(DATA{i_node}{component});
           
    data = DATA{i_node}{component};
    
    % Compute T_m for the current combination of node and component:
    parents_and_node = [parents_of_node; n_plus];
                
    S_m = (n_obs_in_component-1) * cov(data(parents_and_node,:)');
  
    T_m = T_0(parents_and_node,parents_and_node) + S_m + (v*n_obs_in_component)/(v+n_obs_in_component) * (my_0(parents_and_node,1) - mean(data(parents_and_node,:)',1)') * (my_0(parents_and_node,1) - mean(data(parents_and_node,:)',1)')';  
           
    [log_score_local_node_nominator]   = Gauss_Score_complete(length(parents_and_node), n_obs_in_component, v, alpha, T_0(parents_and_node,parents_and_node), T_m);
    [log_score_local_node_denominator] = Gauss_Score_complete(length(parents_of_node),  n_obs_in_component, v, alpha, T_0(parents_of_node,parents_of_node),   T_m(1:(end-1),1:(end-1)));
           
    log_score_node_i = log_score_node_i + (log_score_local_node_nominator - log_score_local_node_denominator); 
                
end
    
log_score = log_score_node_i;    

return

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

function [Vector] = Simulate_changepoints_MIX(Q_vec, DATA_ALL, nodes_in_comp, DAG, CONSTANT)

data_1    = DATA_ALL{1};
[n,m_new] = size(data_1); % NOTE: m_new = m-1

Vector = zeros(1,m_new); %1st column: cp at t=0 -> [1,2,3,...,m_new]
                         %2nd column: cp at t=1 -> [2,3,4,...,m_new]
                                   
                         %(m_new)-th column: cp at t=m_new-1 -> [m_new]
                                  
Vector(:,1) = 1; % changpoint at t=0 

for t=1:(m_new-1) % corresponds to possible cps at 0:(m_new-2)
                  % Note that m_new-1 is the LAST possible cp anyway
   
       n_t = 0;    
        
       cps_at_t = find(Vector(1,:)==1);
       last_cp = cps_at_t(end); % last_cp out of {1,...,m_new}
       if (last_cp==t)
           n_t         = n_t+1;
       end 
 
   if (n_t==1)
         
       [Prob] = compute_prob_MIX(Q_vec, DATA_ALL, nodes_in_comp, DAG, t, CONSTANT);
       
       Prob = Prob/sum(Prob); % I guess this normalisation is required!!!! 
    
       index_no_cp = length(Prob); % check will refer to the last probability in Prob -> no further changepoint
        
       [OUTPUT] = CARPENTER(1,Prob); % A sample of size 1 from the discrete distribution Prob with realisations in {1,...,m_new)
                                       % 1 means cp at 1 -> Vector(1,2)
                                       % 2 means cp at 2 -> Vector(1,3)
                                       % m_new-1 means cp at m_new-1 ->
                                       % Vector(1,m_new)
                                       % m_new means no further changepoint
                                      
       indicis_for_cp = find(OUTPUT~=index_no_cp);
       OUTPUT         = OUTPUT(indicis_for_cp);
    
       for j=1:length(OUTPUT)
            Vector(1,OUTPUT(j)+1)=1; % OUTPUT(j)+1 because the first column of Matrix corresponds to cp at 0
       end  
   end
    
end

Vector = cumsum(Vector,2);

return


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

function [Prob] = compute_prob_MIX(Q_vec, DATA_ALL, nodes_in_comp, DAG, t, CONSTANT)

% t = 1,2,3,....,m_new-1
% corresponds to the old changepoint
% tau_old = 0,1,2,...,m_new-2

global g_vec;
global g0_vec;
global G_vec;
global G0_vec;

tau_old = t-1;

data_1    = DATA_ALL{1};
[n,m_new] = size(data_1);

Prob = zeros(1,m_new); % success time 1,2,3,...,m_new-1 AND m_new means no further changepoints

if (tau_old>0)

    for tau=(tau_old+1):(m_new-1)
        Prob(1,tau) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, (tau_old+1):tau, DAG)+log(Q_vec(tau+1))+log(g_vec(tau-tau_old))-log(Q_vec(tau_old+1)));    
    end
    
    Prob(1,m_new) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, (tau_old+1):m_new, DAG)+log(1-G0_vec(m_new-tau_old-1))+log(CONSTANT)-log(Q_vec(tau_old+1)));
    
    
else % tau_old==0
    
    for tau=1:(m_new-1)
        Prob(1,tau) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:tau, DAG)+log(Q_vec(tau+1))+log(g0_vec(tau))-log(Q_vec(1)));    
    end

    Prob(1,m_new) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:m_new, DAG)+log(1-G0_vec(m_new-1))+log(CONSTANT)-log(Q_vec(1)));

end
     
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




function [log_score] = COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, interval, DAG)

global T_0;
global v;
global alpha;
global my_0;

log_score = 0;

for node = nodes_in_comp
    data_node = DATA_ALL{node};
    data_node = data_node(:,interval);
    [n_plus,m_new]  = size(data_node);

    parents = find(DAG(:,node));
    parents_and_node = [parents; n_plus];

    data_node = data_node(parents_and_node,:);

    S_m = (m_new-1) * cov(data_node');

    T_m = T_0(parents_and_node,parents_and_node) + S_m + (v*m_new)/(v+m_new) * (my_0(parents_and_node,1) - mean(data_node',1)') * (my_0(parents_and_node,1) - mean(data_node',1 )')';  

    [log_score_local_nominator]   = Gauss_Score_complete(length(parents_and_node), m_new, v, alpha, T_0(parents_and_node,parents_and_node), T_m);
    [log_score_local_denominator] = Gauss_Score_complete(length(parents),          m_new, v, alpha, T_0(parents,parents)                  , T_m(1:(end-1),1:(end-1)));
                 
    log_score = log_score + (log_score_local_nominator - log_score_local_denominator);
    
end
        
return

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


function [Q_VEC] = Dynamic_Q_MIX(DATA_ALL, nodes_in_comp, DAG, CONSTANT)

data_1    = DATA_ALL{1};
[n,m_new] = size(data_1); % NOTE: m_new = m-1

Q_VEC = zeros(1,m_new);   % Q_VEC(1,m_new) corresponds to Q(m_new)=P(m_new:m_new|cp at m_new-1)
                          % Q_VEC(1,1)     corresponds to Q(1)=P(1:m_new| cp at 0)
                  
% Determine Q(m_new)
[BN_score] = COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, m_new:m_new, DAG);

Q_VEC(1,m_new) = exp(BN_score +log(CONSTANT));

for i=1:(m_new-1)
    t = m_new-i;  
    % i=1       gives t=m_new-1 
    % i=m_new-1 gives t=1
    Q_t        = RECURSION_MIX(Q_VEC, t, DATA_ALL, nodes_in_comp, DAG, CONSTANT);
    Q_VEC(1,t) = Q_t;
end

return

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

function [Q_t] = RECURSION_MIX(Q_VEC, t, DATA_ALL, nodes_in_comp, DAG, CONSTANT)

global g_vec;
global g0_vec;
global G_vec;
global G0_vec;

% t ranges from: m_new-1,...,1

if (t>1)
  
    data_1     = DATA_ALL{1};
    [n, m_new] = size(data_1);

    Q_t = 0;

    for s=t:(m_new-1)    
        Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, t:s, DAG)+log(Q_VEC(1,s+1))+log(g_vec(s+1-t))); 
    end
    
    Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, t:m_new, DAG)+log(1-G_vec(m_new-t))+log(CONSTANT));
else % t==1
   
    data_1     = DATA_ALL{1};
    [n, m_new] = size(data_1);

    Q_t = 0;

    for s=1:(m_new-1)    
        Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:s, DAG)+log(Q_VEC(1,s+1))+log(g0_vec(s)));  
    end

    Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:m_new, DAG)+log(1-G0_vec(m_new-1))+log(CONSTANT));
      
end
    
      
return

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

function [OUTPUT] = CARPENTER(n,Prob)

x_vec = exprnd(1,1,n+1);

S = sum(x_vec);

u_vec      = cumsum(x_vec/S);

Q = 0;
U = u_vec(1,1);
j = 1;
i = 1;

OUTPUT = [];

while (i<n+1)
    if (U<(Q+Prob(1,j))) 
        OUTPUT = [OUTPUT,j];
        U = u_vec(1,i+1);
        i = i+1;   
    else
        Q = Q + Prob(1,j);
        j = j+1;
    end         
end

return






