# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "dream_user.h"
# include "pdflib.h"
# include "rnglib.h"
# include "problem1_covariance.h"

int main ( );
void test01 ( int par_num, int sample_num );
double *r8mat_covariance ( int m, int n, double x[] );
void r8mat_print ( int m, int n, double a[], char *title );
void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, char *title );

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

int main ( )

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for PROBLEM1_MAIN.

  Discussion:

    The coding of PROBLEM1 is tricky enough that I want to be able to
    try it out independently of the DREAM code.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt
*/
{
  int chain_num;
  int cr_num;
  int gen_num;
  int pair_num;
  int par_num;
  int sample_num;

  timestamp ( );
  printf ( "\n" );
  printf ( "PROBLEM1_MAIN\n" );
  printf ( "  C version\n" );
/*
  Initialize the random number generator library.
*/
  initialize ( );
/*
  By calling PROBLEM_SIZE, we implicitly set up the covariance as well.
*/
  problem_size ( &chain_num, &cr_num, &gen_num, &pair_num, &par_num );

  sample_num = 10000;
  test01 ( par_num, sample_num );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "PROBLEM1_MAIN\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

  return 0;
}
/******************************************************************************/

void test01 ( int par_num, int sample_num )

/******************************************************************************/
/*
  Purpose:

    TEST01 calls the sampling function.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt
*/
{
  double *cov_sample;
  int i;
  int j;
  double z;
  double *zp;
  double *zp_array;
  double *zp_ave;
  double *zp_max;
  double *zp_min;

  printf ( "\n" );
  printf ( "TEST01\n" );
  printf ( "  Call PRIOR_SAMPLE many times.\n" );
  printf ( "  Compare statistics to PDF parameters.\n" );
  printf ( "  Note that the covariance estimate can be very bad\n" );
  printf ( "  unless the matrix is strongly diagonal.\n" );
  printf ( "\n" );
  printf ( "  Parameter dimension is %d\n", par_num );
  printf ( "  Number of samples is %d\n", sample_num );
/*
  Compute N multinormal samples.
*/
  zp_array = ( double * ) malloc ( par_num * sample_num * sizeof ( double ) );
  for ( j = 0; j < sample_num; j++ )
  {
    zp = prior_sample ( par_num );
    for ( i = 0; i < par_num; i++ )
    {
      zp_array[i+j*par_num] = zp[i];
    }
    free ( zp );
  }

  zp_ave = ( double * ) malloc ( par_num * sizeof ( double ) );
  zp_max = ( double * ) malloc ( par_num * sizeof ( double ) );
  zp_min = ( double * ) malloc ( par_num * sizeof ( double ) );

  for ( i = 0; i < par_num; i++ )
  {
    zp_min[i] =   HUGE_VAL;
    zp_max[i] = - HUGE_VAL;
    zp_ave[i] = 0.0;
    for ( j = 0; j < sample_num; j++ )
    {
      z = zp_array[i+j*par_num];
      if ( z < zp_min[i] )
      {
        zp_min[i] = z;
      }
      if ( zp_max[i] < z )
      {
        zp_max[i] = z;
      }
      zp_ave[i] = zp_ave[i] + z;
    }
    zp_ave[i] = zp_ave[i] / ( double ) ( sample_num );
  }
  printf ( "\n" );
  printf ( " Index       Min            Ave              Max             MU\n" );
  printf ( "\n" );
  for ( i = 0; i < par_num; i++ )
  {
    printf ( "  %4d  %14g  %14g  %14g  %14g\n",
      i, zp_min[i], zp_ave[i], zp_max[i], zp_mean[i] );
  }

  cov_sample = r8mat_covariance ( par_num, sample_num, zp_array );

  r8mat_print ( par_num, par_num, cov_sample, "  Sample covariance:" );

  r8mat_print ( par_num, par_num, c, "  PDF covariance:" );

  free ( cov_sample );
  free ( zp_array );
  free ( zp_ave );
  free ( zp_max );
  free ( zp_min );

  return;
}
/******************************************************************************/

double *r8mat_covariance ( int m, int n, double x[] )

/******************************************************************************/
/*
  Purpose:

    R8MAT_COVARIANCE computes the sample covariance of a set of vector data.

  Discussion:

    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt.

  Parameters:

    Input, int M, the size of a single data vectors.

    Input, int N, the number of data vectors.
    N should be greater than 1.

    Input, double X[M*N], an array of N data vectors, each
    of length M.

    Output, double C[M*M], the covariance matrix for the data.
*/
{
  double *c;
  int i;
  int j;
  int k;
  double *x_mean;

  c = ( double * ) malloc ( m * m * sizeof ( double ) );
  for ( j = 0; j < m; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      c[i+j*m] = 0.0;
    }
  }
/*
  Special case of N = 1.
*/
  if ( n == 1 )
  {
    for ( i = 0; i < m; i++ )
    {
      c[i+i*m] = 1.0;
    }
    return c;
  }
/*
  Determine the sample means.
*/
  x_mean = ( double * ) malloc ( m * sizeof ( double ) );
  for ( i = 0; i < m; i++ )
  {
    x_mean[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      x_mean[i] = x_mean[i] + x[i+j*m];
    }
    x_mean[i] = x_mean[i] / ( double ) ( n );
  }
/*
  Determine the sample covariance.
*/
  for ( j = 0; j < m; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      for ( k = 0; k < n; k++ )
      {
        c[i+j*m] = c[i+j*m] 
          + ( x[i+k*m] - x_mean[i] ) * ( x[j+k*m] - x_mean[j] );
      }
    }
  }

  for ( j = 0; j < m; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      c[i+j*m] = c[i+j*m] / ( double ) ( n - 1 );
    }
  }

  free ( x_mean );

  return c;
}
/******************************************************************************/

void r8mat_print ( int m, int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8MAT_PRINT prints an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

    Entry A(I,J) is stored as A[I+J*M]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 May 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows in A.

    Input, int N, the number of columns in A.

    Input, double A[M*N], the M by N matrix.

    Input, char *TITLE, a title.
*/
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    R8MAT_PRINT_SOME prints some of an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Input, double A[M*N], the matrix.

    Input, int ILO, JLO, IHI, JHI, designate the first row and
    column, and the last row and column to be printed.

    Input, char *TITLE, a title.
*/
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }
/*
  Print the columns of the matrix, in strips of 5.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:  ");
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %7d     ", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( m < ihi )
    {
      i2hi = m;
    }
    else
    {
      i2hi = ihi;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
/*
  Print out (up to) 5 entries in row I, that lie in the current strip.
*/
      fprintf ( stdout, "%5d:", i - 1 );
      for ( j = j2lo; j <= j2hi; j++ )
      {
        fprintf ( stdout, "  %14f", a[i-1+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}