%%% 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/>.

%
% KitagawaModelPMCMC2Estimation.m
%

rng('shuffle');

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

T  = 1000;
a1 = 0.5;
b1 = 28;
c1 = 8;
a2 = 0.6;
b2 = 30;
c2 = 10;
d1 = 0.05;
d2 = 0.06;
d3 = 0.07;
d4 = 0.08;
q  = 1e-1;
r  = 1e-1;

[x, y] = KitagawaModel(T, a1, b1, c1, a2, b2, c2, d1, d2, d3, d4, q, r);

save('logs/KitagawaModelPMCMC2_Data.mat');

grids = { ...
  [-30:2:30], ...
  [-30:2:30]  ...
};

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) (a1 .* x(1,:) + b1 .* x(1,:) / (1 + x(1,:) .^ 2)), ...
  @(x) (a2 .* x(2,:) + b2 .* x(2,:) / (1 + x(2,:) .^ 2))  ...
};

obsMeanFuncs = { ...
  @(x) (d1 .* x(1,:) .^ 2 + d2 .* x(2,:) .^ 2), ...
  @(x) (d3 .* x(1,:) .^ 2                    ), ...
  @(x) (                    d4 .* x(2,:) .^ 2)  ...
};

x0 = zeros(1, 2);
u  = repmat(cos(1.2 * [1:T]), 2, 1)';
N  = 500;
J  = 10000;

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

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

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, ...
  [], ... % No auxiliary states
  u', ...
  y', ...
  N, ...
  J, ...
  0, ...
  aspect ...
);

save('logs/KitagawaModelPMCMC2.mat');

disp('Done');
