%%% Monte Carlo Dynamic Classifier tool - MCDC tool
% http://www.ismvideo.org/NDE0UP/MCDCtool/en/index.html

% Copyright (C) 2014 Gareth W. Peters, Tomoko Matsui, Francois Septier, and Akira Tamamori

% This program is free software; you can redistribute it and/or modify it under the terms
% of the GNU General Public License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.

% This program is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE.  See the GNU General Public License for more details.

% You should have received a copy of the GNU General Public License along with this
% program; if not, see <http://www.gnu.org/licenses/>.

%
% LorenzModelPMCMC2Estimation.m
%

rng('shuffle');

%
% Generate state and observation time series from Lorenz model
%

T      = 1000;
Sigma  = 10;
Rho    = 28;
Beta   = 8/3;
Alpha0 = 1;
q      = 1e-1;
r      = 1e-1;
dt     = 1e-2;

[x, y] = LorenzModel(T, Sigma, Rho, Beta, Alpha0, q, r, dt);

save('logs/LorenzModelPMCMC2_Data.mat');

xaux = y(:, 1);

grids = { ...
  [-30:2:30], ...
  [-30:2:30], ...
  linspace(min(xaux), max(xaux), 31) ...
};

stateKernelGens = { ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10)), ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10))  ...
};

obsKernelGens = {
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10)), ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10)), ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10))  ...
};

stateMeanFuncs = { ...
  @(x) (Sigma .* x(2,:) .* dt + x(1,:) .* (1 - Sigma .* dt)),  ...
  @(x) (x(1,:) .* (Rho - x(3,:)) .* dt + x(2,:) .* (1 - dt))   ...
};

obsMeanFuncs = { ...
  @(x) ((x(1,:) .* x(2,:) .* dt) + x(3,:) .* (1 - Beta * dt)), ...
  @(x) (Alpha0 * x(1,:) .^ 2),                                 ...
  @(x) (Alpha0 * x(2,:) .^ 2)                                  ...
};

x0 = zeros(1, 2);
N  = 500;
J  = 2000;

gridDimForGramMatrix = 1; % Use x1 for building gram matrix.

aspect = Aspect();
aspect.LogFileName = 'logs/LorenzModelPMCMC2.log';
aspect.MatFileNamePrefix = 'logs/LorenzModelPMCMC2';
aspect.SavesIntermediateMat = true;
aspect.IntermediateMatInterval = 100;

disp('PMCMC2 Estimation using SPAPI');

%%% Choice 1
% kernelGeneratorStrategy = @RBFKernelGeneratorStrategyChoice1;

%%% Choice 2
kernelGeneratorStrategy = @RBFKernelGeneratorStrategyChoice2;

%%% Choice 1
% mcdcSrategy = @MCDCStrategyChoice1;

%%% Choice 2
mcdcStrategy = @MCDCStrategyChoice2;

%%% Choice 3: build gram matrices here.
% for i = 1:size(stateMeanFuncs, 2)
%   stateSigma{i} = eye(length(grids{gridDimForGramMatrix}));
% end
% for i = 1:size(obsMeanFuncs, 2)
%   obsSigma{i} = eye(length(grids{gridDimForGramMatrix}));
% end
% mcdcStrategy = MCDCStrategyChoice3.Prepare(stateSigma, obsSigma);

splineHandle = @GenericSpline;
[ IDX, SKP, OKP, FV, GV, XE, YE, loglik ] = MCDCTrain( ...
  PMCMC2(kernelGeneratorStrategy, mcdcStrategy), ...
  grids, ...
  stateKernelGens, ...
  obsKernelGens, ...
  stateMeanFuncs, ...
  obsMeanFuncs, ...
  gridDimForGramMatrix, ...
  splineHandle, ...
  x0, ...
  xaux', ...
  [], ... % No controls
  y', ...
  N, ...
  J, ...
  1, ...
  aspect ...
);

save('logs/LorenzModelPMCMC2.mat');

disp('Done');
