function [A parents beta sigma]=BNGaussian(data,order,maxpar,BF,nu0,sigma0,stepwise)
% BN and DBN learning for linear Gaussian networks
%
% [A parents beta sigma]=BNGaussian(data,order,maxpar,BF,nu0,sigma0,stepwise)
% data: data matrix (genes,time points)
% order: variables order (null if dynamic)
% BF: Bayes factor
% nu0,sigma0: prior parameters (for other prior parameters (R0 and beta0)
% see function scoregene)
% stepwise=1 if I want to do stepwise search, otherwise stepwise=0
% maxpar: max number parents per node
% namedir: string containg name of results directory
%
% A: A(i,j)=1 if i is a parent of j
% parents: cell array with parents of each variable
% beta: cell array with regression parameters for each var (same order as
% parents)
% sigma: vector with learned conditionl variances
%
% By Laboratory for Biomedical Informatics
% Dipartimento di Informatica e Sistemistica - Universita' degli Studi di
% Pavia, Italy

if (strcmp(order,''))
    dynamic=true;
else
    dynamic=false;
end

nvar=size(data,1);
ntot=size(data,2);



if (dynamic)
    % dynamic net
    par_data=data(:,1:ntot-1)';
    ydata=data(:,2:ntot)';
else
    % static net
    par_data=data';
    ydata=data';
end

%K2 search

%%%%%%INPUT variables%%%%%%%%%%%
%data
%number of nodes: nvar

if(dynamic)
    %% (dyn) number of samples
    n=ntot-1;
else
    %% (stat) number of samples
    n=ntot;
end

%% %%%%%%%%%RESULTS%%%%%%%%%%%%%%%%%%%%%

%1) cell array with parents for each gene
parents={};

%2) matrix summarizing relationships
%A(i,j)=1 if i is a parent of j
A=zeros(nvar,nvar);

X={};
y={};

%scores
P=[];

%beta vectors
beta={};
%tau values
tau=[];
%sigma^2 values
sigma=[];

%partial scores for each gene
Ppartial={};


%Bayes factors
BFlist={};

%barra di scorrimento
global bar;
set(bar,'maximum',nvar);
set(bar,'StringPainted',1);
set(bar,'String','Please wait...');

for i=1:nvar

    set(bar,'value',i);
    calin=java.util.Date();
   
    %could take it away but it is shorter
    y{i}=ydata(:,i);

    w=1;

    %parents for gene i: at the beginning I have no parents
    par=[];

    X{i}=ones(n,1);

    %calculate score of node i
    score=scoregene(y{i},X{i},nu0,sigma0);
    P(i)=score.likelihood;
    beta{i}=score.beta;
    tau(i)=score.tau;
    sigma(i)=score.sigma;

    Ppartial{i}(1)=P(i);

    %set first BF to 1, by convention
    BFlist{i}(1)=1;

    %control variable to exit loop of gene i
    done=0;

    while (~done && length(par)<maxpar )
        %potential parents for node i
        %NB: setdiff orders the resulting vector
        if(dynamic)
            pps=setdiff(1:nvar,par);
        else
            pos=find(order==i);
            pps=setdiff(order(1:(pos-1)),par);
        end
        if (length(pps)>=1)

            %could take away this initialization
            pscore=[];
            be={};
            ta=[];
            si=[];
            Xp={};
            Xadd={};

            for k=1:length(pps)
                %matrix X containing parent data
                Xadd{k}=par_data(:,pps(k));
                Xp{k}=[X{i} Xadd{k}];

                %calculate score of model in which I add pps(k) to the
                %parents
                sg=scoregene(y{i},Xp{k},nu0,sigma0);
                pscore(k)=sg.likelihood;
                be{k}=sg.beta;
                ta(k)=sg.tau;
                si(k)=sg.sigma;

            end

            %choose parent with max score
            [maxscore maxind]=max(pscore);
            %use maxind(1) in case the max is not unique
            mind=maxind(1);
            new_parent=pps(mind);


            %stepwise:
            if (stepwise==1 && (length(par)>=2) )
                %add pot.parents to gene i parents, consider previous parents
                %of gene i and remove them one by one
                %calculate score and see if, with this removal, score increases
                Xrem={};
                pscore_rem=[];
                parents_rem={};
                brem={};
                trem=[];
                srem=[];

                for s=1:length(par)
                    %take away one parent, add potential parent and
                    %calculate score (I do not use setdiff because it
                    %reorders the elements)
                    parents_rem{s}=par;
                    parents_rem{s}(s)=[];
                    parents_rem{s}=[parents_rem{s} new_parent];
                    parents_rem{s};

                    %matrix parent data
                    Xrem{s}=Xp{mind};
                    %take away colum (s+1) (which corresponds to parent
                    %s)
                    Xrem{s}(:,s+1)=[];

                    sr=scoregene(y{i},Xrem{s},nu0,sigma0);
                    pscore_rem(s)=sr.likelihood;
                    brem{s}=sr.beta;
                    trem(s)=sr.tau;
                    srem(s)=sr.sigma;

                end
                %max score with parent removal
                [maxscore_rem maxind_rem]=max(pscore_rem);
                mind_rem=maxind_rem(1);

                if (maxscore_rem>=maxscore)
                    par_new=parents_rem{mind_rem};
                    maxscore=maxscore_rem;
                    Xnew=Xrem{mind_rem};
                    taunew=trem(mind_rem);
                    betanew=brem{mind_rem};
                    sigmanew=srem(mind_rem);
                else
                    par_new=[par new_parent];
                    Xnew=Xp{mind};
                    taunew=ta(mind);
                    betanew=be{mind};
                    sigmanew=si(mind);

                end

            else  %if stepwise=0 or if length(par)<2 I do not try parent removal
                par_new=[par new_parent];
                Xnew=Xp{mind};
                taunew=ta(mind);
                betanew=be{mind};
                sigmanew=si(mind);
            end

            if (maxscore-P(i)>log(BF))
                par=par_new;
                P(i)=maxscore;
                X{i}=Xnew;
                beta{i}=betanew;
                tau(i)=taunew;
                sigma(i)=sigmanew;

                w=w+1;
                Ppartial{i}(w)=maxscore;

                BFlist{i}(w)=exp(Ppartial{i}(w)-Ppartial{i}(w-1));


            else
                done=1;
            end

        else
            done=1;  %if there are no more potential parents, exit cycle
        end

    end %end WHILE (cycle on variable parents)

    parents{i}=par;
    A(parents{i},i)=1;

    %calculate approximate time to the end of algorithm
    calfine=java.util.Date();
    time(i)=calfine.getTime()-calin.getTime(); %time for variable i in milliseconds
    wait=round((mean(time))*(nvar-i)/1000); %approximate wait time in seconds
    if (wait<=120)
        set(bar,'String',strcat('Please wait...',num2str(wait),' seconds left'));
    else
        set(bar,'String',strcat('Please wait...',num2str(ceil(wait/60)),' minutes left'));
    end

end %end FOR (loop on variables)


%%
function score=scoregene(y,X,nu0,sigma0)
% By Fulvia Ferrazzi
% Laboratory for Biomedical Informatics
% Dipartimento di Informatica e Sistemistica - Universita' degli Studi di Pavia
% Copyright 2006

p=size(X,2)-1;
n=length(y);

%prior parameters
beta0=zeros(p+1,1);
R0=eye(p+1);
%%%%%%%%%%%%%%%%%%%

alfa1=nu0/2;
alfa2=2/(nu0*sigma0);

%%%%%%%%%
alfa1n=nu0/2+n/2;
Rn=R0+X'*X;
betan=inv(Rn)*(R0*beta0+X'*y);
inv_alfa2n=(-betan'*Rn*betan+y'*y+beta0'*R0*beta0)/2+1/alfa2;
alfa2n=1/inv_alfa2n;
nun=nu0+n;
sigman=2/(nun*alfa2n);

%%%%%%%

%pD=1/((2*pi)^(n/2))*(sqrt(det(R0))/sqrt(det(Rn)))*(gamma(nun/2)/gamma(nu0/2))*((nu0*sigma0/2)^(nu0/2))/((nun*sigman/2)^(nun/2));

taun=alfa1n*alfa2n;
beta_new=betan;
sigma_new=nun*sigman/(nun-2);

lpD=-n/2*log(2*pi)+0.5*(log(det(R0))-log(det(Rn)))+gammaln(nun/2)-gammaln(nu0/2)+nu0/2*log(nu0*sigma0/2)-nun/2*log(nun*sigman/2);

%%%%%%%

score=struct('likelihood',lpD,'tau',taun,'beta',beta_new,'sigma',sigma_new);