_math.h

Go to the documentation of this file.

/*  _________________________________________________________________________
 *
 *  UTILIB: A utility library for developing portable C++ codes.
 *  Copyright (c) 2001, Sandia National Laboratories.
 *  This software is distributed under the GNU Lesser General Public License.
 *  For more information, see the README file in the top UTILIB directory.
 *  _________________________________________________________________________
 */

#ifndef ___math_h
#define ___math_h

#if defined(__cplusplus) && defined(ANSI_HDRS)
#include <cmath>
#else
#include <math.h>
#endif
#include "utilib_dll.h"
#include "_generic.h"

#ifndef _MSC_VER

#include <values.h>

#else

#if defined(__cplusplus) && defined(ANSI_HDRS)
#include <climits>
#include <cfloat>
#else
#include <limits.h>
#include <float.h>
#endif
#define MAXINT INT_MAX
#define MAXDOUBLE DBL_MAX

#endif

#ifdef SOLARIS
#include <sunmath.h>
#endif


#ifdef _MSC_VER                                 

#define M_PI 3.14159265358979323846 

#define M_E exp(1.0)
#endif

#ifndef PI
#define PI M_PI
#endif

#ifndef MAXINT
#define MAXINT (int)(~((int)0) ^ (1 << (sizeof(int)*8-1)))
#endif


#ifndef USING_STL
#ifdef min
#undef min
#endif

template <class T, class V>
inline UTILIB_API T min(T x, V y)
{return ((x)<(y)? (x) : (y));}


#ifdef max
#undef max
#endif

template <class T, class V>
inline UTILIB_API T max(T x, V y)
{return ((x)>(y)? (x) : (y));}


template <class T>
inline UTILIB_API void swap(T& i, T& j)
{ T tmp = i; i = j; j = tmp;}
#else
#include <algorithm>
#endif



#ifdef __cplusplus
extern "C" {
#endif

UTILIB_API int calc_filesize(char* filename);

UTILIB_API int round _((double x));

UTILIB_API void setup_bufexp _((int tabsz, double xmin, double xmax));
UTILIB_API double bufexp _((double x));

#if defined(PARAGON) | defined(nCUBE2) | defined(RS6K) | defined(SGI) | defined(COUGAR) 

int isinf(double x);
int isnan(double x);
double exp2(double x);
#endif
#if defined(PARAGON) | defined(nCUBE2) | defined(COUGAR)

void bzero(char* b, int length);
#endif

#ifdef __cplusplus
};
#endif

#ifdef BUFFERED_EXP
#define Exp(x)  bufexp(x)
#else
#define Exp(x)  exp(x)
#endif


#if defined(__BasicArray_h) && !defined(BasicArray_DEFN)
#define __math_h_NO_ARRAYS
#else
#ifndef __math_h_NO_ARRAYS
#include "DUniform.h"
#include "Uniform.h"
#include "Basic2DArray.h"
#include "IntVector.h"
#include "DoubleVector.h"
#include "BasicArray.h"
#include "default_rng.h"

class BitArray;
#endif
#endif



#ifndef __math_h_NO_ARRAYS


template <class T>
UTILIB_API T sum(const BasicArray<T>& vec)
{
T ans=0;

for (size_type i=0; i<vec.size(); i++)
  ans += vec[i];
return(ans);
}


template <class T>
UTILIB_API void shuffle(BasicArray<T>& vec, RNG* rng, size_type num)
{
Uniform urnd(rng);
T temp;

for (size_type i = 0; i<num; i++) {
  size_type j = Discretize(urnd(), i, vec.size()-1);
  temp = vec[j];
  vec[j] = vec[i];
  vec[i] = temp;
  }
}


UTILIB_API void shuffle(BitArray& vec, RNG* rng=&default_rng);


template <class T>
inline UTILIB_API void shuffle(BasicArray<T>& vec, RNG* rng=&default_rng)
{shuffle(vec,rng,vec.size());}


template <class T>
UTILIB_API void subshuffle(BasicArray<T>& vec, RNG* rng, size_type start, 
                size_type end)
{
DUniform urnd(rng);
T temp;

for (size_type i = start; i<end; i++) {
  size_type j = urnd(i,end-1);
  temp = vec[j];
  vec[j] = vec[i];
  vec[i] = temp;
  }
}


template <class T>
UTILIB_API double length(const BasicArray<T>& a)
{
double ans=0.0;

for (size_type i=0; i<a.size(); i++)
  ans += a[i]*a[i];

return sqrt(ans);
}


template <class T>
UTILIB_API double dist(const BasicArray<T>& a, const BasicArray<T>& b)
{
double ans=0.0;

for (size_type i=0; i<a.size(); i++)
  ans += (a[i]-b[i])*(a[i]-b[i]);

return sqrt(ans);
}


template <class T>
UTILIB_API size_type argmin(const BasicArray<T>& vector)
{
size_type i, temp=0;

if (vector.size() == 1)
   return(0);

for (i=1; i<vector.size(); i++)
  if (vector[i] < vector[temp])  temp=i;

return(temp);
}


template <class T>
UTILIB_API size_type argmax(const BasicArray<T>& vector)
{
size_type i, temp=0;

if (vector.size() == 1)
   return(0);

for (i=1; i<vector.size(); i++)
  if (vector[i] > vector[temp])  temp=i;

return(temp);
}


size_type argmedian(double* x, size_type n, int* ws, RNG* rng);

               
template <class T>
UTILIB_API size_type argmedian(const BasicArray<T>& vector, RNG* rng=&default_rng)
{
T* tmp_array = (T*)&(vector[0]);
return argmedian( tmp_array, vector.size(), 0, rng );
}


template <class T>
UTILIB_API size_type argmedian(const BasicArray<T>& vec, BasicArray<int>& ws, 
                        RNG* rng=NULL)
{ return argmedian(vec.data(), vec.size(), ws.data(), rng); }


template <class T>
inline UTILIB_API T min(const BasicArray<T>& vec) 
{ return vec[argmin(vec)]; }


template <class T>
inline UTILIB_API T max(const BasicArray<T>& vec) 
{ return vec[argmax(vec)]; }


template <class T>
inline UTILIB_API T median(const BasicArray<T>& vec) 
{ return vec[argmedian(vec)]; }


template <class T>
UTILIB_API BasicArray<T> concat(const BasicArray<T>& a, const BasicArray<T>& b)
{
BasicArray<T> ans(a.size() + b.size());

for (size_type i=0; i<a.size(); i++)
  ans[i] = a[i];
for (size_type j=a.size(); j<(a.size()+b.size()); j++)
  ans[j] = b[j-a.size()];

return(ans);
}


template <class T>
double mean(const BasicArray<T>& array)
{
double temp=0.0;

for (size_type i=0; i<array.size(); i++)
  temp += array[i];

return( temp/array.size() );
}


template <class T>
double trimmed_mean(BasicArray<T>& array, const double percent)
{
double temp=0.0;

sort(array);

size_type boundary = (size_type) floor(array.size() * percent);

for (size_type i=boundary; i<(array.size()-boundary); i++)
  temp += array[i];

return( temp/(array.size() - 2.0*boundary) );
}


template <class T>
double var(const BasicArray<T>& array, const int sampleflag=TRUE)
{
double array_mean = mean(array);

double temp=0.0;
for (size_type i=0; i<array.size(); i++)
  temp += (array[i] - array_mean) * (array[i] - array_mean);

if ((sampleflag == FALSE) || (array.size() == 1))
   return( temp / array.size() );
else
   return( temp / (array.size() - 1) );
}


template <class T>
double mad(BasicArray<T>& array, BasicArray<double>& work)
{
double meanval = mean(array);

for (size_type i=0; i<work.size(); i++)
  work[i] = fabs(array[i] - meanval);

return median(work);
}


template <class T>
double mad(BasicArray<T>& array)
{
BasicArray<double> tmparray(array.size());
return mad(array,tmparray);
}


template <class T>
UTILIB_API void lapply(BasicArray<T>& vec, double (*func)(T val))
{
for (size_type i=0; i<vec.size(); i++)
  vec[i] = (*func)(vec[i]);
}


template <class T>
UTILIB_API double inner_product(const BasicArray<T>& v1,
                                const BasicArray<T>& v2)
{
double ans = 0.0;

for (size_type i=0; i<v1.size(); i++)
  ans += v1[i]*v2[i];

return ans;
}


#include "DoubleMatrix.h"

template <class T>
NumArray<double> mean(const Basic2DArray<T>& matrix, const int stats_flag)
{
NumArray<double> temp(matrix.ncols());

temp = 0.0;

for (size_type j=0; j<matrix.nrows(); j++)
  for (size_type i=0; i<matrix.ncols(); i++)
    temp[i] += matrix[j][i];

if (stats_flag == TRUE)
     temp /= ((double) (matrix.nrows() - 1));
else
     temp /= ((double) matrix.nrows());

return( temp );
}


template <class T>
DoubleVector var(const Basic2DArray<T>& mat, const int sampleflag=TRUE)
{
DoubleVector array_mean;
return( var(mat, array_mean, sampleflag) );
}


template <class T>
DoubleVector var(const Basic2DArray<T>& mat, BasicArray<double>& array_mean, 
                                                        const int sampleflag=TRUE)
{
DoubleVector ans(mat.ncols());
array_mean &= mean(mat);

ans = 0.0;
for (size_type i=0; i<mat.nrows(); i++)
  for (size_type j=0; j<ans.size(); j++)
    ans[j] += ((mat[i][j] - array_mean[j]) * (mat[i][j] - array_mean[j]));

if ((sampleflag == FALSE) || (mat.nrows() == 1))
   ans /= (double) (mat.nrows());
else
   ans /= (double)(mat.nrows()-1);
return ans;
}


template <class T>
UTILIB_API int min(const Basic2DArray<T>& mat)
{
T ans,tmp;

ans = mat[0][0];
for (size_type i=0; i<mat.nrows(); i++)
  for (size_type j=0; j<mat.ncols(); j++) {
    tmp = mat[i][j];
    if (tmp < ans) ans=tmp;
    }

return ans;
}


template <class T>
UTILIB_API T max(const Basic2DArray<T>& mat)
{
T ans,tmp;

ans = mat[0][0];
for (size_type i=0; i<mat.nrows(); i++)
  for (size_type j=0; j<mat.ncols(); j++) {
    tmp = mat[i][j];
    if (tmp > ans) ans=tmp;
    }

return ans;
}


template <class T>
UTILIB_API T sum(const Basic2DArray<T>& mat)
{
T ans;

ans = 0;
for (size_type i=0; i<mat.nrows(); i++)
  for (size_type j=0; j<mat.ncols(); j++)
    ans += mat[i][j];

return ans;
}


template <class T>
void rowscale(Basic2DArray<T>& a)
{
double tmp;
size_type j,i;

for (i=0; i<a.nrows(); i++) {
  tmp = 0.0;
  for (j=0; j<a.ncols(); j++)
    tmp += a[i][j];
  if (tmp != 0.0)
     for (j=0; j<a.ncols(); j++)
       a[i][j] /= tmp;
  }
}

UTILIB_API int cholesky (DoubleMatrix& A, DoubleMatrix& G, int n);

#endif // __math_h_NO_ARRAYS


#endif