// =======================================================================================
// F2_numerics.h
// Time-stamp: <2004-02-25 15:37:39 aknoblau>
// Definition of classes for numerical computations
// programmed in May 2001 by Andreas Knoblauch
// Department of Neural Information Processing
// University of Ulm, Germany
// Version 1.0 alpha
// last documented change: May 6, 2001
// =======================================================================================
// types/classes/interfaces:
// TFunctionLibrary
// TMLookUpTable
// ---------------------------------------------------------------------------------------
// global objects:
// ---------------------------------------------------------------------------------------
// related modules:
// ---------------------------------------------------------------------------------------
// notes:
// ---------------------------------------------------------------------------------------
// history:
// Version 1.0 alpha programmed on May 6, 2001 by Andreas Knoblauch
// ---------------------------------------------------------------------------------------
// bugs:
// ---------------------------------------------------------------------------------------
// 2 do:
// =======================================================================================
#ifndef F2_numerics_H
#define F2_numerics_H
#include
#include
#include "F2_types.h"
#include "F2_layout.h"
#include "F2_parser.h"
/*
// **************************************************************************************************************************
// **************************************************************************************************************************
// TKernelLibrary
// Library class for important kernel functions
// Status: version 1.0 alpha
// last change: May 30, 2001
// **************************************************************************************************************************
class TKernelLibrary {
public:
static inline TFloat* setGaussian(TCuboidLT layout, TFloat* kernel, int* center, TFloat* sig, TFloat* phi, // general case
TFloat scale, TFloat offset, TFloat normalize) {
int N,nDim,*mid;
TFloat *a,*sin_phi,*cos_phi,*x_;
TFloat sum, sumEuklid;
N = layout->nDim;
nDim = layout->nDim;
mid = new int[nDim];
a = new TFloat[nDim];
sin_phi = new TFloat[nDim-1];
cos_phi = new TFloat[nDim-1];
x_ = new TFloat[nDim-1];
for(i=0;isize[i]&1)||(layout->size[i]<0)) {
clog << "TKernelLibrary::setGaussian() (F2_numerics) : layout contains non-odd or negative dimension!\n layout= " << layout << "\n";
throw TFelixException("TKernelLibrary::setGaussian() (F2_numerics) : layout contains non-odd or negative dimension!");
}
mid[i]=layout->size[i]/2;
a [i]=4.0/(((TFloat)layout->size[i]-1.0)*((TFloat)layout->size[i]-1.0));
}
delete[] mid;
delete[] a;
delete[] sin_phi;
delete[] cos_phi;
delete[] x_;
clog << "F2_numerics: setGaussian() general case not yet implemented!!!!!!\n";
throw TFelixException("F2_numerics: setGaussian() general case not yet implemented!!!!!!");
return kernel; // not yet implemented!!!!!!!!!!
}
static inline TFloat* setGaussian(int KY, int KX, TFloat* kernel, int CY, int CX, int sigY, int sigX, TFloat phi, // 2D case
TFloat scale, TFloat offset, TFloat normalize) {
int i,j, mid_x,mid_y;
float a,b,integral,sin_phi,cos_phi,x,y,x_,y_,mx_,my_;
float euklid;
mid_x = KX/2;
mid_y = KY/2;
a = 4.0/(((float)KX-1.0)*((float)KX-1.0)); // Ellipsenparameter der Gleichung ax^2+by^2=1
b = 4.0/(((float)KY-1.0)*((float)KY-1.0)); // siehe dazu AB1/69 Seite 88
integral = 0.0;
euklid = 0.0;
sin_phi = (float)sin(phi);
cos_phi = (float)cos(phi);
for(i=0;i
static inline TNum min(TNum x1, TNum x2) {
return (x1
static inline TNum max(TNum x1, TNum x2) {
return (x1
static inline TNum sign(TNum x) {
return (x==0) ? 0 : (x>0) ? 1 : -1;
}
template
static inline TInt round(TNum x) {
return (TInt)(x+0.5);
}
template
static inline TNum mod(TNum a, TNum b) { // zur Modulo-Rechnung in C siehe auch AB1/92 S.146
return ( (a)>=0 ? (a)%(b) : ( ((-(a))%(b)) ? ((a) + ((-(a))/(b)+1)*(b)) : 0 ) );
}
template
static inline TNum fmod(TNum a, TNum b) { // zur Modulo-Rechnung in C siehe auch AB1/92 S.146
return ( (a)>=0 ? fmod((a),(b)) : ( (fmod((-(a)),(b))>EPS) ? ((a) + floor((-(a))/(b)+1.0)*(b)) : (0) ) );
}
static inline TFloat info(TFloat p1) {
if((p1<=0)||(p1>=1)) return 0;
return (-p1*log(p1)-(1-p1)*log(1-p1))/TConst::LN2;
}
static inline TFloat transinfo(TFloat p1, TFloat p01, TFloat p10, int norm) { // p1 = pr[x=1], p01 = pr[y=1|x=0], p10 = pr[y=0|x=1]
TFloat PB[2][2]; // conditional probability PB[i][j] = pr[y=j|x=i] (i->j)
TFloat PV[2][2]; // compound probability PV[i][j] = pr[x=i,y=j] (i&&j)
TFloat py; // pr[y=1]
TFloat I_Y,I_YgivenX; // T(X,Y) = I(Y)-I(Y|X)
int i,j;
py = (1.0-p1)*p01 + p1*(1.0-p10);
PB[0][0] = 1.0-p01; PB[0][1] = p01; PB[1][0] = p10; PB[1][1] = 1.0-p10;
PV[0][0] = (1.0-p1)*PB[0][0]; PV[0][1] = (1.0-p1)*PB[0][1]; PV[1][0] = p1*PB[1][0]; PV[1][1] = p1*PB[1][1];
if((py<=0)||(py>=1)) I_Y=0.0;
else I_Y = (-py*log(py)-(1.0-py)*log(1.0-py))/TConst::LN2;
for(i=0,I_YgivenX=0.0;i<2;i++)for(j=0;j<2;j++)
if(PB[i][j]>0) I_YgivenX += (-PV[i][j]*log(PB[i][j]));
I_YgivenX /= TConst::LN2;
return norm?((I_Y>0)?(I_Y-I_YgivenX)/I_Y:0):(I_Y-I_YgivenX);
}
static inline TFloat transinfo(TFloat p1, TFloat p01, TFloat p10) {
return transinfo(p1,p01,p10,0);
}
static inline TFloat transinfo01(TFloat p1, TFloat p01, TFloat p10) {
return transinfo(p1,p01,p10,1);
}
static inline TFloat binomialCoeff(TFloat n, TFloat k) { // binomial coefficient
double v=1;
long long int n_=(long long int)(n+0.5), k_=(long long int)(k+0.5);
while(k_>0) v*=((double)(n_--))/((double)(k_--));
return (TFloat)v;
}
static inline TFloat pdfBinomial(TFloat x, TFloat n, TFloat p1) { // Binomial probability density function
if(((TFloat)(int)x)!=x) return 0;
if((x<0)||(x>n)) return 0;
return binomialCoeff(n,x)*pow(p1,x)*pow(1-p1,n-x);
}
static inline TFloat cpdfBinomial(TFloat x, TFloat n, TFloat p1) { // Binomial cumulative probability distribution function
double v;
if(x<0) return 0;
if(x>=n) return 1;
long long int k_=(long long int)(x+0.5),n_=(long long int)(n+0.5),i;
if(x<=(n/2)) {
i=0;v=0;
while(i<=k_) v+=pdfBinomial((TFloat)(i++),(TFloat)n_,p1);
}else{
i=k_+1;v=1;
while(i>k_) v-=pdfBinomial((TFloat)(i--),(TFloat)n_,p1);
}
return (TFloat)v;
}
static inline TFloat gauss1D_(TFloat x, TFloat ex, TFloat sig) {
return SQRT2PI_INV/sig*exp(-(x-ex)*(x-ex)/2.0/sig/sig);
}
static inline TFloat gauss1D(TFloat* par) {
return gauss1D_(par[0],par[1],par[2]);
}
static inline TFloat gauss1D_interval(TFloat* par) {
TFloat x=par[0];
return ((x>=par[3])&&(x<=par[4]))?gauss1D(par):0;
}
static inline TFloat sqrt1D(TFloat* par) {
return (TFloat)sqrt(par[0]);
}
static inline TFloat fermi1D_(TFloat x, TFloat beta) {
return 1.0/(1.0+exp(-(beta*x)));
}
static inline TFloat fermi1D(TFloat* par) {
return fermi1D_(par[0],par[1]);
}
static inline TFloat semiPower_(TFloat x, TFloat exp) {
return (x<=0)?0:pow(x,exp);
}
static inline TFloat semiPower(TFloat* par) {
return semiPower_(par[0],par[1]);
}
static TFloat* getParBuffer(TFloat i1); // delivers TFloat array {i1}
static TFloat* getParBuffer(TFloat i1, TFloat i2); // delivers TFloat array {i1,i2}
static TFloat* getParBuffer(TFloat i1, TFloat i2, TFloat i3); // delivers TFloat array {i1,i2,i3}
static TFloat* getParBuffer(TFloat i1, TFloat i2, TFloat i3, TFloat i4); // delivers TFloat array {i1,i2,i3,i4}
static TFloat* getParBuffer(TFloat i1, TFloat i2, TFloat i3, TFloat i4, TFloat i5); // delivers TFloat array {i1,i2,i3,i4}
static TFloat* getParBuffer(TFloat i1, TFloat i2, TFloat i3, TFloat i4, TFloat i5, TFloat i6); // delivers TFloat array {i1,i2,i3,i4}
protected:
static const int lenParBuffer=5; // length of parBuffer
static TFloat parBuffer[lenParBuffer]; // buffer for parameters
};
typedef TFunctionLibrary TFunct;
// **************************************************************************************************************************
// **************************************************************************************************************************
/**
* TMLookUpTable
* look-up tables for looking-up or linear interpolation of mathematical functions
*
* @short look-up tables for looking-up or linear interpolation of mathematical functions
* @author Andreas Knoblauch
* @version 1.0 alpha,
* last change: May 30, 2001
*/
template
class TMLookUpTable : public TMNumMarray, public TFunctionLibrary {
// **************************************************************************************************************************
public:
string id; // identifier, usually a name referring to the tabled function
string& scopeID_TMLookUpTable;
int& nTable; // size of the look-up table; a reference to TMMarray::N
TNum*& table; // look-up table; a reference to TMMarray::data
TNum* derivs; // of size N-1; derivs[i]=(table[i+1]-table[i])/dx
TNum* x; // of size N; x[i]=x1+x*dx
TNum x1; // independent variable corresponding to table[0]
TNum dx; // resolution of the axis of the independent variable
TNum dx_inv; // inverse of dx
int nArgs; // number of arguments of function (including independent variable)
TNum (*fun) (TNum*); // function
TNum* args; // Arguments for function fun
TMLookUpTable(); // constructors
TMLookUpTable(const char* id_arg, int nTab_arg, TNum x1_arg, TNum dx_arg, int nArgs_arg, TNum(*fun_arg)(TNum*), TNum* pars);
TMLookUpTable(const char* id_arg, int nTab_arg, TNum x1_arg, TNum dx_arg, TNum* table_arg);
TMLookUpTable(char* fname);
TMLookUpTable(TParser& parser);
TMLookUpTable(const TMLookUpTable & copy); // copy constructor
~TMLookUpTable(); // destructor
void setTable(TNum *pars); // set table according to fun and pars
void set_derivs_x(); // (re-)allocate and set derivs and x
void operator=(const TMLookUpTable & copy); // copy function
void integrate(int leadingZero); // if leadingZero a zero is inserted (if necessary) at the beginning of table
void cutForInterval(TNum minVal, TNum maxVal); // restrict tab values to interval [minVal;maxVal]
void invertTable(int nTab_arg, TNum x1_arg, TNum dx_arg); // inverts the look-up-table
inline TNum invert(TNum y); // estimates x in equation y=fun(x) by linear interpolation; tab invertible!
void normalize(TNum sum_target); // normalizes that sum over table is sum_target
void inverseDistribution(int nTab_arg); // interpretes the look-up-table as density and computes inverse distribution
inline TNum interpolate(TNum x_arg); // computes interpolated function value
inline TNum interpolate_safe(TNum x_arg); // with index checking
inline TNum nearest(TNum x_arg); // returns the look-up-table entry nearest to x
inline TNum nearest_safe(TNum x_arg); // with index checking
inline TNum exact(TNum x_arg); // returns exact function value of x (if fun is allocated)
void writeToFile(const char* fname); // write table to file fname in format FELIX
void writeToFile(const char* fname, TFileFormat format); // write table to file fname with a specified format,e.g. to plot it in Matlab
friend ostream & operator<< (ostream & os, const TMLookUpTable & par); // output operator
friend TParser & operator>> (TParser & parser, TMLookUpTable & par); // input operator for TParser
private:
static const char* scopeString;
};
template
const char* TMLookUpTable::scopeString="TMLookUpTable";
template
ostream & operator<<(ostream & os, const TMLookUpTable & par); // output operator
template
TParser & operator>>(TParser & parser, TMLookUpTable & par); // input operator
// **************************************************************************************************************************
// **************************************************************************************************************************
// template operators
// for some classes it is necessary to allow generic operators, e.g. binary operators like >> for float types
// Status: version 1.0 alpha
// last change: May 30, 2001
// **************************************************************************************************************************
// shiftRight functions instead of operator>>
inline TUnsigned8 shiftRight(TUnsigned8 u) {return u>>1;}
inline TUnsigned16 shiftRight(TUnsigned16 u) {return u>>1;}
inline TUnsigned32 shiftRight(TUnsigned32 u) {return u>>1;}
inline TUnsigned8 shiftRight(TUnsigned8 u1, TUnsigned8 u2) {return u1>>u2;}
inline TUnsigned16 shiftRight(TUnsigned16 u1, TUnsigned8 u2) {return u1>>u2;}
inline TUnsigned32 shiftRight(TUnsigned32 u1, TUnsigned8 u2) {return u1>>u2;}
template
inline T shiftRight(T f) {
clog << "F2_numerics: shiftRight() for non-unsigned types not allowed!\n"
<< "Probably you used a float or signed int type for container of packed bits!\n";
throw TFelixException("F2_numerics: shiftRight() for non-unsigned types not allowed!");
}
template
inline T shiftRight(T f1, TUnsigned8 u2) {return shiftRight(f1);}
// include inline methods
#define INLINE 1
#include "F2_numerics.cpp"
#undef INLINE
#endif /* F2_numerics_H */
| Generated by: aknoblau on synfire on Sat May 1 14:32:16 2004, using kdoc 2.0a54. |