/*
**  MaxOLS
**
**  Purpose:
**    Explain maximization is a fully different manner
**
**  Date:
**    18/8/2011
**
**  Author:
**    Charles Bos
*/
#include <oxstd.h>
#include <oxfloat.h>
#include <oxdraw.h>
// #include <packages/gnudraw/gnudraw.h>
// #include <packages/oxutils/oxutils.h>
#import <maximize>

static decl s_vY, s_mPar;

/*
**  Initialise(const avY, const vP, const iN)
**
**  Purpose:
**    ...
**
**  Inputs:
**    ...
**
**  Return value:
**    ..
*/
Initialise(const avY, const vP, const iN)
{
  decl dSigma, dMu;
  
  [dMu, dSigma]= {vP[0], vP[1]};
  
  avY[0]= dSigma*rann(1, iN) + dMu;
  
  print ("Y: ", "%c", {"Mean", "Variance"}, meanr(avY[0])~varr(avY[0]));
  Draw(0, avY[0]);
  DrawDensity(1, avY[0], "y", 1, 0, 1);
  DrawQQ(2, avY[0], "y", QQ_N, 0, 1);
  ShowDrawWindow();
}


/*
**  AvgLnLiklNorm(const vP, const adLnPdf, const avScore, const amHess)
**
**  Purpose:
**    Calculate LL of normal model
**
**  Inputs:
**    vP    2 x 1 vector of mu and sigma
**    s_vY  1 x iN vector of observations
**
**  Outputs:
**    adLnPdf   double, average LL
**
**  Return value:
**    br        boolean, TRUE if all well
*/
AvgLnLiklNorm(const vP, const adLnPdf, const avScore, const amHess)
{
  decl dMu, dSigma, vE, vLL;
  
  adLnPdf[0]= M_NAN;
  [dMu, dSigma]= {vP[0], vP[1]};

  vE= s_vY - dMu;
  vLL= -log(dSigma) - 0.5*sqr(vE/dSigma);
  adLnPdf[0]= meanr(vLL);
  
  s_mPar~= (vP|adLnPdf[0]);
  
  return !ismissing(adLnPdf[0]);
}

/*
**  AvgLnLiklNormS2(const vP, const adLnPdf, const avScore, const amHess)
**
**  Purpose:
**    Calculate LL of normal model, in terms of sigma^2
**
**  Inputs:
**    vP    2 x 1 vector of mu and sigma^2
**    s_vY  1 x iN vector of observations
**
**  Outputs:
**    adLnPdf   double, average LL
**
**  Return value:
**    br        boolean, TRUE if all well
*/
AvgLnLiklNormS2(const vP, const adLnPdf, const avScore, const amHess)
{
  decl dMu, dSigma2;
  
  [dMu, dSigma2]= {vP[0], vP[1]};

  return AvgLnLiklNorm(dMu|sqrt(dSigma2), adLnPdf, avScore, amHess);
}

/*
**  AvgLnLiklNormLogS(const vP, const adLnPdf, const avScore, const amHess)
**
**  Purpose:
**    Calculate LL of normal model, in terms of log sigma
**
**  Inputs:
**    vP    2 x 1 vector of mu and log sigma
**    s_vY  1 x iN vector of observations
**
**  Outputs:
**    adLnPdf   double, average LL
**
**  Return value:
**    br        boolean, TRUE if all well
*/
AvgLnLiklNormLogS(const vP, const adLnPdf, const avScore, const amHess)
{
  decl dMu, dLnSigma;
  
  [dMu, dLnSigma]= {vP[0], vP[1]};

  return AvgLnLiklNorm(dMu|exp(dLnSigma), adLnPdf, avScore, amHess);
}

/*
**  PlotGrid(const fnFunc, const mLU, const sTitle)
**
**  Purpose:
**    
**
**  Inputs:
**    
**
**  Return value:
**    
*/
PlotGrid(const fnFunc, const mLU, const sTitle)
{
  decl v1, v2, vP, i1, i2, i, j, ir, dLnPdf, mLL, vMax; 
  
  v1= range(mLU[0][0], mLU[0][1], mLU[0][2]);
  v2= range(mLU[1][0], mLU[1][1], mLU[1][2]);
  print (mLU, v1, v2);
  i1= sizerc(v1);
  i2= sizerc(v2);
  
  mLL= zeros(i1, i2);

  for (i= 0; i < i1; ++i)
    for (j= 0; j < i2; ++j)
      {
        vP= v1[i]|v2[j];
        ir= fnFunc(vP, &dLnPdf, 0, 0);
        mLL[i][j]= dLnPdf;
      }
//   print (sTitle, mLL);      
  
  DrawXYZ(0, v1, v2, mLL', 1);
  DrawTitle(0, sTitle);
  ShowDrawWindow();
  
  vMax= maxc(mLL);
  j= maxcindex(vMax');
  vMax= maxc(mLL');
  i= maxcindex(vMax');
  
  println ("Found a maximum ", mLL[i][j], " at v1= ", v1[i], " and v2= ", v2[j]);

}

main()
{
  decl iN, dMu, dSigma, vP, vY, ir, dLnPdf;
  
  // Magics
  iN= 100;
  dMu= 42;
  dSigma= 3.5;
  
  // Initialisation
  vP= dMu|dSigma;
  Initialise(&vY, vP, iN);
  
  s_vY= vY;
  s_mPar= <>;
  ir= AvgLnLiklNorm(vP, &dLnPdf, 0, 0);
  println ("AvgL returns: ", ir, " with LL ", dLnPdf);
  
//   PlotGrid(AvgLnLiklNorm, <35, 50, .5; 1, 7, .1>, "mu and sigma");
//   PlotGrid(AvgLnLiklNormS2, <35, 50, .5; 1, 49, 1>, "mu and sigma^2");
//   PlotGrid(AvgLnLiklNormLogS, <35, 50, .5; 0, 2, .01>, "mu and log sigma");

  vP= <31; 7>;
  ir= MaxBFGS(AvgLnLiklNorm, &vP, &dLnPdf, 0, 1);
  println ("Max returns ", MaxConvergenceMsg(ir), " at ", vP, 
  " with LL= ", double(iN*dLnPdf), " using ", columns(s_mPar), " function evals");
  
  DrawXMatrix(0, s_mPar[0][], "Mu", s_mPar[1][], "sigma", 2, 2);
  DrawZ(range(1, columns(s_mPar)), "", ZMODE_VALUE);

  vP= <31; 7>;
  s_mPar= <>;
  vP[1]= log(vP[1]);
  ir= MaxBFGS(AvgLnLiklNormLogS, &vP, &dLnPdf, 0, 1);
  println ("Max returns ", MaxConvergenceMsg(ir), " at ", vP, 
  " with LL= ", double(iN*dLnPdf), " using ", columns(s_mPar), " function evals");

  DrawXMatrix(1, s_mPar[0][], "Mu", s_mPar[1][], "sigma (using logs)", 2, 2);
  DrawZ(range(1, columns(s_mPar)), "", ZMODE_VALUE);
  
  ShowDrawWindow();
    
}