#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include "mkl_vml.h"
#include "mex.h"
#include "matrix.h"
#include "mkl_vsl.h"
#include <time.h>
#define SEED time(NULL)


/* sub-functions */
double sum(double *vector, int dim);
double updatephi(int n, double *x, double phi, double accept, double sigma2, double out[3], VSLStreamStatePtr stream) ;


/* main fucntion */
void mexFunction(int nlhs,  mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    
    int n ;
    double  *x, phi, sigma2, *out, accept;
    
    /* for distributions */
    VSLStreamStatePtr stream;
    vslNewStream( &stream, VSL_BRNG_MT19937, SEED );
    
    
    /* make pointers to input data */
    n = (int)mxGetScalar(prhs[0]);
    x = mxGetPr(prhs[1]);
    phi = (double)mxGetScalar(prhs[2]);
    sigma2 = (double)mxGetScalar(prhs[3]);
    accept = (double)mxGetScalar(prhs[4]);
    
    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix( 1, 3, mxREAL);
    out = mxGetPr(plhs[0]);
    
    
    /* phi */
    updatephi(n, x, phi, accept, sigma2, out, stream);
    

    
    
    /* Deleting the stream */
    vslDeleteStream( &stream );
    
    
    return;
}

double updatephi(int n, double *x, double phi, double accept, double sigma2, double out[3], VSLStreamStatePtr stream)
{

    double thetacan[1], phican, u[1], logold, logcan, theta, logf, sumx;
    
    
    /* step 1 */
    theta = log(1.0 + phi) - log(1.0 -phi);
    
    vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, stream, 1, thetacan, theta, sqrt(sigma2) );

    
    /* step 2 */
    phican = (exp(thetacan[0])-1.0)/(exp(thetacan[0])+1.0);
    
    /* step 3 */
    sumx = sum(x, n);
    logcan = sumx*log(phican) + ((double)n - sumx)*log(1.0 - phican) + thetacan[0] -2*log(1.0 +exp(thetacan[0]));
    logold = sumx*log(phi) + ((double)n - sumx)*log(1.0 - phi) + theta -2*log(1.0 +exp(theta));
    logf = logcan - logold;
    
    /* step 4 */
    vdRngUniform( VSL_RNG_METHOD_UNIFORM_STD, stream, 1, u, 0.0, 1.0 );
    
    /* step 5 */
    if(log(u[0])<logf){
        
        out[1] = phican;
        out[2] = accept + 1.0;
    }
    else
    {
        out[1] = phi;
        out[2] = accept;
    }
    
    out[0] = logf;
    
    return(0);
}



/* ============================================================== */
double sum(double *vector, int dim)
{
    int i;
    double sum = 0.0;
    for(i=0; i< dim ;i++){ sum += vector[i]; }
    
    return(sum);
}