/**************************************************************************************************************************************
/* mexFunction [best_t,M,tti,t_samps]=box_D_func(f4,[ts ts(end)+1],n_ints,px1,px,in,int,s,length(s_i)/p1,length(ns_i)/p0,m0,samp_t-1);
/* LP 12/2/02
/*
/* C code to quickly find global maximum of M(t) over circle, as described in Paninski, `03
/*
/* inputs: all inherited from phi_func.m
/* out: best_t = argmax M(t)
/*      M = max(M) = M(best_t)
/*      tti = M(t_samps)
/*      t_samps = some randomly chosen sampling times, grabbed for and "max-distance" "info-cov" choice for e_0 
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mex.h"

void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray*prhs[])
{
  long ii,lf4,mm,nm,i,n_samps;
  int iy,samp,k;
  double *f4,*ts,*n_ints,*px1,*px,*in,*ints,*s,*ls_i,*lns_i,*m0,*samp_ti;
  double *best_ti, *M, *tti, *t_samps;
  double try_m,old_p1_x,new_p1_x,eps=.00001; /* eps is a noise-floor term, to avoid some minor problems caused by underflow errors */
  
  /* input vars */
  f4=mxGetPr(prhs[0]);
  ts=mxGetPr(prhs[1]);
  n_ints=mxGetPr(prhs[2]);
  px1=mxGetPr(prhs[3]);
  px=mxGetPr(prhs[4]);
  in=mxGetPr(prhs[5]);
  ints=mxGetPr(prhs[6]);
  s=mxGetPr(prhs[7]);
  ls_i=mxGetPr(prhs[8]);
  lns_i=mxGetPr(prhs[9]);
  m0=mxGetPr(prhs[10]);
  samp_ti=mxGetPr(prhs[11]);
  
  /* output vars */
  mm = mxGetM(prhs[11]);
  nm = mxGetN(prhs[11]);
  n_samps=mm*nm-1;
  plhs[0]=mxCreateDoubleMatrix(1,1,mxREAL);
  best_ti=mxGetPr(plhs[0]);
  plhs[1]=mxCreateDoubleMatrix(1,1,mxREAL);
  M=mxGetPr(plhs[1]);
  plhs[2]=mxCreateDoubleMatrix(1,n_samps,mxREAL);
  tti=mxGetPr(plhs[2]);
  plhs[3]=mxCreateDoubleMatrix(1,n_samps,mxREAL);
  t_samps=mxGetPr(plhs[3]);
  M[0]=m0[0]; best_ti[0]=0;
  
  try_m=m0[0]; k=0; samp=0;
  mm = mxGetM(prhs[0]);
  nm = mxGetN(prhs[0]);
  lf4=mm*nm;
  /* main loop */
  for(i=0;i<lf4;i++) {
    ii=f4[i]-1;
    iy=ii-n_ints[0]*floor(ii/n_ints[0]);
    if(px[iy]>0) {old_p1_x=px1[iy]/px[iy];}
    else {old_p1_x=0;}
    in[ii]=-in[ii];
    px[iy]=px[iy]+in[ii];
    if(ceil((ii+1)/n_ints[0])<=s[0]) {px1[iy]=px1[iy]+in[ii];}
    if(px[iy]>0) {new_p1_x=px1[iy]/px[iy];}
    else {new_p1_x=0;}
    if(ints[iy]<=s[0]) {try_m=try_m+(new_p1_x-old_p1_x)/ls_i[0];}
    else {try_m=try_m+(old_p1_x-new_p1_x)/lns_i[0];}
    if(i==samp_ti[k]) {
      tti[samp]=try_m;
      t_samps[samp]=ts[i];
      samp++;
      k++;
    }
    if(try_m>M[0] && ts[i+1]>ts[i]+eps) {
      M[0]=try_m;
      best_ti[0]=(ts[i+1]+ts[i])/2;
    }
  }
}