#include <stdio.h>
#include <math.h>
#include "randomS.c"

void genfield(x,n,m,beta,lambda,iter,seed)  
/************************************************************************* 
Generates a Gaussian MRF via MCMC.  Assumes 3rd order nbhd
where beta[1..6] correspond to h,v,d+-,d++,hh,vv neighbors
respectively.  Because this matrix is transposed relative to
the Splus matrix, the horizontal adjacencies are reversed here.
***********************************************************************/
double x[1000][1000], *beta, lambda;
int n, m, iter;
{
 int ns, i, j, k, by, i1, i2, j1, j2, i0, j0;
 int samp;
 double mean, meanprec, prec, sd, nbrs, sumbeta;
 double rnor();


 /* initialize x */

 sumbeta = 0.0;
 for(i=0;i<6;i++) sumbeta += beta[i];
 sd = 1.0/sqrt(lambda);
 for(i=0;i<n;i++) for(j=0;j<m;j++) x[i][j] = rnor(0.0,sd); 

 for(samp=1;samp<=iter;samp++){

 for(i0=0;i0<n;i0++) for(j0=0;j0<m;j0++){

   mean = 0.0;
   nbrs = 0.0;
    /* h */
   mean += beta[1]*x[i0][(j0-1+m)%m];
   mean += beta[1]*x[i0][(j0+1)%m];
    /* v */
   mean += beta[0]*x[(i0-1+n)%n][j0];
   mean += beta[0]*x[(i0+1+n)%n][j0];
    /* d+- */
   mean += beta[2]*x[(i0-1+n)%n][(j0+1+m)%m];
   mean += beta[2]*x[(i0+1+n)%n][(j0-1+m)%m];
    /* d++ */
   mean += beta[3]*x[(i0-1+n)%n][(j0-1+m)%m];
   mean += beta[3]*x[(i0+1+n)%n][(j0+1+m)%m];
    /* hh */
   mean += beta[5]*x[i0][(j0-2+m)%m];
   mean += beta[5]*x[i0][(j0+2+m)%m];
    /* vv */
   mean += beta[4]*x[(i0-2+n)%n][j0];
   mean += beta[4]*x[(i0+2+n)%n][j0];

   sd = 1.0/sqrt(lambda);
   
   x[i0][j0] = rnor(mean,sd);
 }
 } /* samp = 1,..,iter */
   
}

/*****************************************************************

 This is a driver program for sgenfield().  You can uncomment this
 and skip all the Splus interface stuff.  You have to specify parameters
 here in this part of the program before calling sgenfield().

*****************************************************************/

/*
void main()
{
  double x[10000], beta[6], lambda[1];
  int i, iter, seed, m, n;
  void sgenfield();

  printf("in main...\n");
  n=100; m=100;
  iter = 100; seed = 2;
  for(i=0;i<6;i++) beta[i] = 1*1.0/12.0;
  
  beta[0] = .595;
  beta[1] = -.209;
  beta[2] = .198;
  beta[3] = .105;
  beta[4] = -.179;
  beta[5] = -.011;
  
  lambda[0] = 1.0/1.0*1;

  printf("in main: calling sgenfield...\n");
  sgenfield(x,&n,&m,beta,lambda,&iter,&seed);

  for(i=0;i<n*m;i++) printf("%6.2f\n",x[i]);
  printf("done!\n");
}
*/

  

void sgenfield(x,n,m,beta,lambda,iter,seed)
/*************************************************************
Splus driver for the genfield function
*************************************************************/
double *x, *beta, *lambda;
int *n, *m, *iter, *seed;
{
 int i, j, k;
 double xmat[1000][1000];
 /* double **Sarray2(); */

 printf("in sgenfield...\n");
 /* set the seed */
 setseedS(*seed);

 printf("in sgenfield: seed set...\n");
 /* note that S uses col major, C uses row major */
 /* xmat = Sarray2(*m,*n); */
 for(k=0;k<m[0] * n[0];k++){
   i = k%n[0];
   j = k/n[0];
   xmat[i][j] = x[k];
 }

 printf("in sgenfield: array allocated...\n");

 /* run the simulator */
 printf("xmat[0][0] = %f\n",xmat[0][0]);
 genfield(xmat,*m,*n,beta,*lambda,*iter);
 /* printf("xmat[0][0] = %f\n",xmat[0][0]); */
 /* for(i=0;i<m[0]*n[0];i++) x[i] = xmat[0][i]; */
 for(k=0;k<m[0] * n[0];k++){
   i = k/n[0];
   j = k%n[0];
   x[k] = xmat[i][j];
 }

}

/*
double **Sarray2(n1,n2)
***********************************************************************
  allocates space for a 2 index array (matrix) 0..n1-1, 0...n2-1
***********************************************************************
int n1,n2;
{
  double **a;
  int i, j;
 
  
  a = (double **) S_alloc(n1 , sizeof(double *));
  
  a = (double **) calloc(n1 , sizeof(double *));
  if(a == NULL) {
    printf("allocation error in Sarray\n");
  }
   
  a[0] = (double *) S_alloc(n1 * n2 , sizeof(double));
   
  a[0] = (double *) calloc(n1 * n2 , sizeof(double));
  if(a == NULL) {
    printf("allocation error in Sarray\n");
  }
  for(i=1;i<n1;i++) a[i] = a[i-1] + n2;

  return a;
}
*/