% script file to plot IACT as a function of dimension
% for AM anf twalk 
% Target is normal distribution with identity covariance 
% (note, both algorithms invariant under affine transformation in low dim)

addpath /home/colin/Computing/MTCpaper/twalk
addpath /home/colin/Computing/MTCpaper/am
addpath /home/colin/Computing/uwerr

ds = [1 2 4 8 16]; %dimensions

nsamp = 1e5;
nd = length(ds); 
it1 = zeros(1,nd); %to store iact for 1st coord in twalk
itl = zeros(1,nd); %to store iact for log likelihood in twalk
ia1 = zeros(1,nd); %to store iact for 1st state in am
ial = zeros(1,nd); %to store iact for log likelihood in am

% setup for twalk
ft = @(x,p) -sum(x.^2/2); % twalk takes log target (p is parameter passed to f)

% setup for am
clear model data params options
%
model.ssfun = @(x,data)  sum(x.^2);     % am takes minus 2 times log target (sum of squares)
% options for |mcmcrun| with am
options.nsimu   = nsamp;
options.method  = 'am';  % use am, default is DRAM
options.verbosity = 0; % shut the f**k up
data = [];


for cnt = 1:length(ds)
  d = ds(cnt);
  sig = eye(d); % don't actually need this -- use as intial covariance in am

% twalk
  x0 = zeros(d,1);     % initial state for chains
  x1 = 0.1* ones(d,1); % second starting value for twalk 
 
  [xx,lt] = twalkparam(ft,nsamp,x0,x1,sig);
  [value,dvalue,ddvalue,tauint,dtauint,Qval] = UWerr(xx(1,:)',1.5,length(xx),0);
  it1(cnt) = tauint*2; %to get stats definition of IACT
  [value,dvalue,ddvalue,tauint,dtauint,Qval] = UWerr(lt(1,:)',1.5,length(lt),0);
  itl(cnt) = tauint*2; %to get stats definition of IACT

% am needs some parameters and other stuff
  model.N = d; % AM does not appear to use this to set dimension -- perhaps for recording only
  options.qcov = sig;
  data.lam = inv(sig); % intial covariance, perhaps needed in AM
  data.x0 = zeros(1,d); % mean vector
  % Create the compulsory prameter structure in a loop -- it appears that dimension, etc, are set by this
  clear params
  for i=1:d
    %            'name'                   initial
    params{i} = {sprintf('\\theta_{%d}',i), x0(i)};
  end
  [results,chain,sigchain,sschain] = mcmcrun(model,data,params,options);

  [value,dvalue,ddvalue,tauint,dtauint,Qval] = UWerr(chain(:,1),1.5,length(chain),0);
  ia1(cnt) = tauint*2; %to get stats definition of IACT (of first component)
  [value,dvalue,ddvalue,tauint,dtauint,Qval] = UWerr(sschain,1.5,length(sschain),0);
  ial(cnt) = tauint*2; %to get stats definition of IACT (of log likelihood)
end