%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% BATCH FOR GLAUBER MODEL %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% General parameters %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

i2mPath;
SizeMax=5;
nbpatterns=1000;
BinSize = 10;
datalen=100000;
N=8;
state=1;
BinIndex = 1;
tauTab = [1.5 2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Initializing storage matrices %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PexpTab = zeros(nbpatterns,SizeMax);
PthTab = zeros(nbpatterns,SizeMax);
PthNoJ1Tab = zeros(nbpatterns,SizeMax);
PthNoJTab = zeros(nbpatterns,SizeMax);
PthNoJ0Tab = zeros(nbpatterns,SizeMax);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Generating the Glauber raster and statistics %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%J1test = random('unif',-0.15,0.15,N,N);
J1test = -0.1 + 0.2*rand(N,N);%or: random('unif',-0.1,0.1,N,N);
htest = -1.05 + 0.05*rand(N,1);%or: random('unif',-1.05,-1.0,N,1);

for k=1:length(tauTab)%Here we generate raster with the Glauber model for two different taus
    %Alternatively, a raster can be loaded by using the LoadRaster and
    %BinRaster functions
    tau = tauTab(k);
    raster{k} = GlauberRaster(J1test,htest,N,datalen,1,tau);%For this parameter tau, we generate a raster
    
    %HERE START THE ANALYSIS WITH THE STATISTICAL MODEL. In this example, the data are in
    %raster{k}, which is a binned raster composed of -1 and 1. See the
    %sub-directory /LoadRaster to load from data instead of using the
    %Glauber model. 
    m{k} = mean(raster{k},2);
    C{k} = raster{k}*transpose(raster{k})/datalen;
    C1{k} = raster{k}(:,1:(datalen-1))*transpose(raster{k}(:,2:datalen))/(datalen-1);
    [mtot,Ctot] = ConvertSpatial(m{k},C{k},C1{k});
    % First guess with the Tanaka approximation
    [htotguess,Jtotguess]=hJinitial2(mtot,Ctot,2*N); 
    [hstatGuess,JstatGuess]=hJinitial2(m{k},C{k},N);       
    
    % Monte Carlo gradient descent
    [htot,Jtot]=GradientMC(mtot,Ctot,htotguess,Jtotguess,10000000,0.1,50,0.007); 
    [hstat,Jstat]=GradientMC(m{k},C{k},hstatGuess,JstatGuess,1000000,0.1,100,0.007);
    [h,J,J1,hr,Jr] = DeConvertSpatial(htot,Jtot,hstat,Jstat); 
    
    %Store the fitted coefficients h and J in some tables for further use. 
    coefsh(:,k,BinIndex) = h;
    coefsJ(:,:,k,BinIndex) = J;
    coefsJ1(:,:,k,BinIndex) = J1;
    coefshr(:,k,BinIndex) = hr;
    coefsJr(:,:,k,BinIndex) = Jr;
    coefshstat(:,k,BinIndex) = hstat;
	coefsJstat(:,:,k,BinIndex) = Jstat;
        
    %Estimating the 3 models performance: predicting probabilities of
    %individual patterns (figure 1)
    for TempSize=1:SizeMax%For patterns of size 1 to SizeMax
    
        %%% Model with spatio-temporal correlations %%%
        rast{1} = raster{k};
        [Pexp,Pth,Features,Plist,RasterSize]=OcPred(rast,hstat,Jstat,hr,Jr,h,J,J1,TempSize,nbpatterns);%See OcPred picks up some patterns, compute their ocurrence empirically, and predict it from the model. 
    	divjs(TempSize,k)=PlotOcFit(Pexp,Pth,Features,BinSize,TempSize,N);% Plot the prediction vs empirical estimations
        
        PexpTab(1:length(Pexp),BinIndex,TempSize,k) = Pexp;% These 3 lines store the results in table for subsequent use in Fig 1.
        PthTab(1:length(Pexp),BinIndex,TempSize,k) = Pth;
        PlistTab{BinIndex,TempSize,k} = Plist;
        FeaturesTab(1:length(Pexp),BinIndex,TempSize,k) = Features;
        ['Divergence, t ' int2str(TempSize),': ', num2str(divjs(1,TempSize,k)) ]
    	close;

        %%% Model with only spatial correlations %%%    
        PthNoJ1 = OcPredFast(Plist,Pexp,hstat,Jstat,zeros(N,1),zeros(N,N),hstat,Jstat,zeros(N,N));
    	divjsNoJ1(TempSize,k)=PlotOcFit(Pexp,PthNoJ1,Features,BinSize,TempSize,N);
    	PthNoJ1Tab(1:length(Pexp),BinIndex,TempSize,k) = PthNoJ1;
        ['Divergence when no Temp, t ' int2str(TempSize),': ', num2str(divjsNoJ1(1,TempSize,k)) ]
        close;

        %%% Model without correlation %%%
        PthNoJ = OcPredFast(Plist,Pexp,hstat,zeros(N,N),zeros(N,1),zeros(N,N),hstat,zeros(N,N),zeros(N,N));
        divjsNoJ(TempSize,k)=PlotOcFit(Pexp,PthNoJ,Features,BinSize,TempSize,N);
        PthNoJTab(1:length(Pexp),BinIndex,TempSize,k) = PthNoJ;
        ['Divergence when no J, t ' int2str(TempSize),': ', num2str(divjsNoJ(1,TempSize,k)) ]
        close;
    
    end

    
end


save([pwd,'/Workspace/GlauberIsingData.mat'])