/*  miscprob.h

    Christian Beckmann & Mark Woolrich, FMRIB Image Analysis Group

    Copyright (C) 1999-2000 University of Oxford  */

/*  CCOPYRIGHT  */

// Miscellaneous maths functions that rely on libprob build ontop of miscmaths


#if !defined(__miscprob_h)
#define __miscprob_h

#include "miscmaths.h"
#include "libprob/libprob.h"
#include "stdlib.h"

using namespace NEWMAT;
// using namespace NEWRAN;

namespace MISCMATHS {
 
  ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1, 
		       const float start = 0, const float end = 1);

  int distribrnd(const ColumnVector& histo);
  
  ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1, 
		       const float mu = 0, const float sigma = 1);

  // returns nsamps*nparams matrix:
  ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp = 1);
  

  ReturnMatrix normpdf(const RowVector& vals, const float mu = 0, const float var = 1);

  ReturnMatrix normpdf(const RowVector& vals, const RowVector& mus, 
		       const RowVector& vars);

  float normpdf(const ColumnVector& val,const ColumnVector& mu,const SymmetricMatrix& sigma);
  
  ReturnMatrix normpdf(const Matrix& val,const ColumnVector& mu,const SymmetricMatrix& sigma);

  ReturnMatrix normcdf(const RowVector& vals, const float mu = 0, const float var = 1);

  ReturnMatrix gammapdf(const RowVector& vals, const float mu = 0, const float var = 1);

  ReturnMatrix gammacdf(const RowVector& vals, const float mu = 0, const float var = 1);

  // returns n! * n matrix of all possible permutations
  ReturnMatrix perms(const int n);

  // ReturnMatrix wishrnd(const SymmetricMatrix&,const int);
  // ReturnMatrix iwishrnd(const SymmetricMatrix&,const int);
  
  class Mvnormrandm
    {
    public:
      Mvnormrandm(){}

      Mvnormrandm(const RowVector& pmu, const SymmetricMatrix& pcovar) :
	mu(pmu),
	covar(pcovar)
	{
	  Matrix eig_vec;
	  DiagonalMatrix eig_val;
	  EigenValues(covar,eig_val,eig_vec);

	  covarw = sqrt(eig_val)*eig_vec.t();
	}

      ReturnMatrix next(int nsamp = 1) const 
	{
	  Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
	  ret.Release();
	  return ret;
	}
    private:      

      RowVector mu;
      SymmetricMatrix covar;

      Matrix covarw;

    };
}
#endif