#include <oxstd.h>
#include <oxfloat.h>
#include <oxstd.h>
#include <oxfloat.h>
#include <maximize.h>
#import <solvenle>
#include <oxdraw.h>
#include <oxprob.h>

enum{week,lakes,mcal,lsal,fwal,phal,boal,gtal,npal,caal,frl,gr,pr,pa,triprob
  ,mcg,lsg,fwg,phg,bog,gtg,npg,cag,mcfl,lsfl,fwfl,phfl,bofl,gtfl,npfl
  ,cafl,mcp,lsp,fwp,php,bop,gtp,npp,cap,tqg,dm1,dm2,dm3,dm4
  ,seas1,seas2,seas3,seas4,seas5,seas6,seas7,seas8,seas9,seas10,seas11,seas12}
enum{alpha0,beta0,alpha1,beta1,alpha2,beta21,beta22,beta23,beta24
    ,beta3,sig1,sig2,sig12,gamma}
decl rw;
decl wk;
decl par;
decl par1;
decl like;
decl iLLF;


likelihood(const par, const adFunc, const avScore, const amHessian)
 {
decl y;
decl x;
decl B;
decl G;
decl Sig;
decl prob;
decl likem01, likem1;
decl likem02, likem2;
decl likem;

B = new matrix[2][2];
G = new matrix[2][7];
Sig = new matrix[2][2];

wk = rw[][week];
println("beta3 = ", par[beta3]);
println("par = ", par);
y = log(rw[][tqg])~log(rw[][gr]);
x = ones(wk)~rw[][lakes]~rw[][gr]~rw[][dm1]~rw[][dm2]~rw[][dm3]~rw[][dm4];
// println("tqg", meanc(rw[][tqg]));
// println("lakes", meanc(rw[][lakes]));
// println("gr", meanc(rw[][gr]));
// println("dm1", meanc(rw[][dm1]));
// println("dm2", meanc(rw[][dm2]));
// println("dm3", meanc(rw[][dm3]));
// println("dm4", meanc(rw[][dm4]));
B[0][0] = 1.0;
B[1][1] = 1.0;
B[0][1] = -par[alpha1];
B[1][0] = -par[beta1];
G[0][0] = par[alpha0];
G[0][1] = par[alpha2];
G[0][2] = 0.0;
G[0][3] = 0.0;
G[0][4] = 0.0;
G[0][5] = 0.0;
G[0][6] = 0.0;
G[1][0] = par[beta0];
G[1][1] = 0.0;
G[1][2] = 0.0;
G[1][3] = par[beta21];
G[1][4] = par[beta22];
G[1][5] = par[beta23];
G[1][6] = par[beta24];
Sig[0][0] = par[sig1]*par[sig1];
Sig[0][1] = par[sig12];
Sig[1][0] = Sig[0][1];
Sig[1][1] = par[sig2]*par[sig2];
prob      = exp(par[gamma])/(1.0+exp(par[gamma]));
likem01 = (B*y'-G*x'-par[beta3]*ones(wk)')';
likem1 = diagonal(likem01*invert(Sig)*likem01');

likem02 = (B*y'-G*x')';
likem2 = diagonal(likem02*invert(Sig)*likem02');
likem = prob*exp(-0.5*likem1)+(1.0-prob)*exp(-0.5*likem2); 
likem = log(determinant(Sig)^(-0.5)*determinant(B)*likem);
like = sumc(likem');
println("like =", like );
adFunc[0] = like;
return 1;
}

main()
 {
rw = loadmat("railway.mat");
par1 = loadmat("input.mat");
par = par1[0:gamma];
println("rw = ", rw[0:10][0:5]);
MaxControl(100,1);
MaxBFGS(likelihood,&par,&iLLF,0,TRUE);
println("parameter");
println(par);
// like = likelihood(par);
 }