text:
/*! \file GMatrix.h
 *  \brief this file defines the \ref Math::GMatrix class that is part of the Math namespace
 */

#ifndef GMATRIX_H
#define GMATRIX_H

#include <iostream>
#include <complex>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <sstream>
#include <stdexcept>
#include <Precision.h>
#include <Coordinate.h>
#include <Matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

using namespace std;
namespace Math
{

/*! \class Math::GMatrix
    \brief General n x m matrix
    General matrix class that eventually should replace \ref Math::Matrix class.
    will have to adjust error checks for this class and adjust inverse
    and make most calls inline to increase speed.
*/
class GMatrix
{
    private:

    int row,col;
    Double_t *matrix;

    protected:

    ///useful comparison functions which can be used with qsort
    int compare_Double(const void* a,const void* b) {
        Double_t arg1 = *((Double_t*) a);
        Double_t arg2 = *((Double_t*) b);
        if( arg1 < arg2 ) return -1;
        else if( arg1==arg2 ) return 0;
        else return 1;
    }

    public:

    /// \name Constructors & Destructors
    //@{
    /// Constructor (where data can be is initilized) otherwise set to zero
    GMatrix(int r, int c, Double_t *data=NULL)
    {
        if (r*c>0){
            row=r;col=c;
            matrix=new Double_t[row*col];
            if (data!=NULL){
                for (int i = 0; i < row; i++)
                    for (int j = 0; j < col; j++)
                        matrix[i*col+j] = data[i*col+j];
            }
            else{
                for (int i = 0; i < row; i++)
                    for (int j = 0; j < col; j++)
                        matrix[i*col+j] = 0.;
            }
        }
        else {
            row=0;col=0;matrix=NULL;
            printf("Error, net size 0 with row,col=%d,%d\n",row,col);
        }
    }

    /// Constructor based on Coordinate.
    GMatrix(const Coordinate &c)
    {
        row=3;col=1;
        matrix=new Double_t[row*col];
        for (int i = 0; i < row; i++) matrix[i] = c[i];
    }

    /// Constructor Matrix
    GMatrix(const Matrix &m)
    {
        row=3;col=3;
        matrix=new Double_t[row*col];
        for (int i=0; i<row; i++)
            for (int j=0;j<col;j++) matrix[i*col+j] = m(i,j);
    }

    /// Copy Constructor
    GMatrix(const GMatrix &m)
    {
        row=m.row;col=m.col;
        matrix=new Double_t[row*col];
        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                matrix[i*col+j] = m.matrix[i*col+j];
    }

    /// Identity Matrix of dimension N
    GMatrix(int N)
    {
        row=N;col=N;
        matrix=new Double_t[row*col];
        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                matrix[i*col+j] = (i==j);
    }

    ~GMatrix(){
        if (row>0&&col>0) delete[] matrix;
    }
    //@}

    /// \name Get & Set functions
    //@{
    ///get row dimensions
    int Row()const {return row;}
    ///get col dimensions
    int Col()const {return col;}
    ///get rank
    int Rank()const {return min(row,col);}
    //@}

    /// \name Overloaded Operators
    //@{
    /// How to access a matrix element: M(i, j)
    Double_t& operator () (int i, int j) {
        //if (j<col&&i<row)
        return matrix[j + col * i];
        //else {printf("Error, (i,j)=(%d, %d) exceed matrix dimesions of %d %d \n", i,j,row,col); return 0;}
    }
    const Double_t& operator() (int i, int j) const
    {
        //if (j<col&&i<row)
        return matrix[j + col * i];
        //else {printf("Error, (i,j)=(%d, %d) exceed matrix dimesions of %d %d \n", i,j,row,col); return 0;}
    }

    std::runtime_error incompatible_dimensions_error(const GMatrix &m)
    {
        std::ostringstream os;
        os << "Error, incompatible dimensions: "
           << "(" << row << "," << col << ") + "
           << "(" << m.row << "," << m.col << ")";
        return std::runtime_error(os.str());
    }

    /// Matrix addition
    GMatrix operator + (const GMatrix& m)
    {
        if (row==m.row&&col==m.col){
            Double_t data[row*col];
            for (int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                    data[i*col+j] = matrix[i*col+j]+m.matrix[i*col+j];
            return GMatrix(row,col,data);
        }
        throw incompatible_dimensions_error(m);
    }

    /// Matrix subtraction
    GMatrix operator - (const GMatrix& m)
    {
        if (row==m.row&&col==m.col){
            Double_t data[row*col];
            for (int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                    data[i*col+j] = matrix[i*col+j]-m.matrix[i*col+j];
            return GMatrix(row,col,data);
        }
        throw incompatible_dimensions_error(m);
    }

    /// Multiply by a scalar: a * M
    friend GMatrix operator * (Double_t a, const GMatrix& m);

    /// Multiply by a scalar: M * a
    GMatrix operator * (Double_t a)
    {
        Double_t data[row*col];
        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                data[i*col+j] = matrix[i*col+j]*a;
        return GMatrix(row,col,data);
    }

    /// Multiply by another Gmatrix
    GMatrix operator * (const GMatrix& m)
    {
        if (col==m.row){
            Double_t data[row*m.col];
            for (int i = 0; i < row; i++)
                for (int j = 0; j < m.col; j++){
                    data[i*m.col+j]=0.;
                    for (int k=0;k<col;k++)
                        data[i*m.col+j]+=matrix[i*col+k]*m.matrix[k*m.col+j];
                }
            return GMatrix(row,m.col,data);
        }
        else {
            printf("Error, incompatible dimensions (%d,%d) + (%d,%d)\n",row,col,m.row,m.col);
            return GMatrix(0,0);
        }

    }
    /// Asignment constructor
    GMatrix& operator = (const GMatrix m)
    {
        delete[] matrix;
        row=m.row;
        col=m.col;
        matrix=new Double_t[row*col];
        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                matrix[i*col+j] = m.matrix[i*col+j];
        return *this;
    }
    //@}

    /// \name Mathematical Operations
    //@{

    /// Trace the matrix.
    Double_t Trace() const
    {
        Double_t trace=1.0;
        for (int i = 0; i < row; i++)trace*=matrix[i*col+i];
        return trace;
    }

    /// return the Transpose of the matrix.
    GMatrix Transpose() const
    {
        Double_t data[col*row];
        for (int i = 0; i < col; i++)
            for (int j = 0; j < row; j++)
                data[i*row+j]=(*this)(i,j);
        return GMatrix(col,row,data);
    }

    /// Transpose matrix in place (changes the matrix permanently!)
    void TransposeInPlace()
    {
        Double_t data[col*row];
        int temp;
        for (int i = 0; i < col; i++)
            for (int j = 0; j < row; j++)
                data[i*row+j]=(*this)(i,j);
        temp=row;row=col;col=temp;
        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                matrix[i*col+j]=data[i*col+j];
    }
    ///Get the diagonals
    GMatrix Diag(){
        GMatrix D(row,col);
        for (int i = 0; i < col; i++) D(i,i)=matrix[i*col+i];
        return D;
    }
    /// Return Pivot matrix P of A, where P*A=B, such that B(i,i)!=0
    GMatrix Pivot() const
    {
        Double_t a[row*col];
        GMatrix P(row,col);
        for (int j=0;j<row;j++) P(j,j)=1;
        if (Trace()==0) {
            for (int j=0;j<row;j++) for (int k=0;k<row;k++) a[j*col+k]=matrix[j*col+k];
            for (int j=0;j<row;j++) if (a[j*col+j]*a[j*col+j]<1e-32) {
                int k=j+1,ipivot;
                if (k==row)k=0;
                ipivot=(P(k,k)==1);
                while (!(a[j*col+k]*a[j*col+k]>1e-32&&ipivot)){
                    k++;
                    if (k==row)k=0;
                    ipivot=(P(k,k)==1);
                }
                P(j,j)=0;P(k,k)=0;
                P(j,k)=1;P(k,j)=1;
                for (int l=0;l<row;l++) {
                    Double_t temp=a[l*col+j];
                    a[l*col+j]=a[l*col+k];
                    a[l*col+k]=temp;
                }
            }
        }
        return P;
    }
    /// Return Pivot matrix P of A, where P*A=B, such that B(i,i)!=0
    GMatrix PivotInPlace()
    {
        GMatrix P(row,col);
        for (int j=0;j<row;j++) P(j,j)=1;
        if (Trace()==0) {
            for (int j=0;j<row;j++) if (matrix[j*col+j]*matrix[j*col+j]<1e-32) {
                int k=j+1,ipivot;
                if (k==row)k=0;
                ipivot=(P(k,k)==1);
                while (!(matrix[j*col+k]*matrix[j*col+k]>1e-32&&ipivot)){
                    k++;
                    if (k==row)k=0;
                    ipivot=(P(k,k)==1);
                }
                P(j,j)=0;P(k,k)=0;
                P(j,k)=1;P(k,j)=1;
                for (int l=0;l<row;l++) {
                    Double_t temp=matrix[l*col+j];
                    matrix[l*col+j]=matrix[l*col+k];
                    matrix[l*col+k]=temp;
                }
            }
        }
        return P;
    }

    /// return a SubMatrix of the matrix from row ra to  row rb and col ca to cb
    GMatrix SubMatrix(int ra, int rb, int ca, int cb) const
    {
        int r=(rb-ra+1), c=(cb-ca+1);
        Double_t data[r*c];
        for (int i = ra; i <= rb; i++)
            for (int j = ca; j <=cb; j++)
                data[(i-ra)*c+(j-ca)]=(*this)(i,j);
        return GMatrix((rb-ra+1),(cb-ca+1),data);
    }

    /// return a submatrix ignore column r, c
    GMatrix ExtractMatrix(int r, int c) const
    {
        Double_t data[(row-1)*(col-1)];
        int tempi=0;
        for (int i = 0; i < row; i++)
        for (int j = 0; j <col; j++)
            if (i!=r&&j!=c) data[tempi++]=(*this)(i,j);
        return GMatrix((row-1),(col-1),data);
    }

    /// Check if matrix square symmetric
    int isSymmetric() const
    {
        if (row!=col) return 0;
        for (int i=0;i<row;i++)
        for (int j=i+1;j<col-1;j++)
            if (((*this)(i,j)-(*this)(j,i))*((*this)(i,j)-(*this)(j,i))/(*this)(i,j)/(*this)(i,j)>1e-8) return 0;
        return 1;
    }
    /// Check if matrix is zero matrix
    int isZero() const
    {
        Double_t sm=0;
        for (int i=0;i<row;i++)
            for (int j=0;j<col;j++)
                sm+=matrix[i*row+j]*matrix[i*row+j];
        if (sm==0.0) return 1;
        else return 0;
    }

    /// Determinant of matrix, for now recursively implemented.
    Double_t Det() const
    {
        Double_t total =0.;
        int tempi;
        if (row==col) {
            if (row==2) return matrix[0]*matrix[3]-matrix[1]*matrix[2];
            else {
                for (int i = 0; i < col; i++) {
                    Double_t data[(row-1)*(row-1)];
                    tempi=0;
                    for (int j=1;j<row;j++)
                        for (int k=0;k<col;k++)
                            if (k!=i) data[tempi++]=(*this)(j,k);
                    total+=matrix[i]*(GMatrix(row-1,col-1,data)).Det();
                }
                return total;
            }
        }
        else {
            printf("not a square matrix %d, %d returning 0\n",row,col);
            return 0.;
        }
    }

    /*! \brief Lower-Upper decompistion

        LU decomposition of matrix stored in such a fashion that
        matrix return is (for 3x3):
        \f[
        \left(\begin{array}{ccc}
            L11 & U12 &U13\\
            L21 & L22 & U23\\
            L31 & L32 & L33
        \end{array}\right)
        \f]
        L is lower triangular matrix, U is upper unit triangular matrix. \n
        Crout's algorithm is used to determine the L and U matrices. System of equations to solve are : \n
        \f[
        \begin{array}{lcl}
            i<j : & l_{i1}*u_{1j}+l_{i2}*u_{2j}+...+l_{ii}*u_{ij}&=a_{ij} \\
            i=j : & l_{i1}*u_{1j}+l_{i2}*u_{2j}+...+l_{ii}*u_{jj}&=a_{ij} \\
            i>j : & l_{i1}*u_{1j}+l_{i2}*u_{2j}+...+l_{ij}*u_{ij}&=a_{ij}
            \end{array}
        \f]
        with condition that \f$ l_{ii}=1 \f$. \n
        If matrix singular and no decomposition, print error and return zero matrix.
        Here matrix is singular if diagonal element is zero. \n
        Must alter this to pivot method so that subroutine reorder matrix if necessary and diagonals
        can be zero.
    */
    GMatrix LUDecomp() const
    {
        if (row!=col) {
            printf("only implemented for square, returning zeros\n");
            return GMatrix(row,col);
        }
        else {
            Double_t A[row*col];
            Int_t i, j, k, p, p_k, p_row, p_col;
            for (j=0;j<row;j++) for (k=0;k<col;k++) A[j*col+k]=matrix[j*col+k];
            for (k = 0, p_k = 0; k < row; p_k += row, k++) {
                for (i = k, p_row = p_k; i < row; p_row += row, i++) {
                    for (p = 0, p_col = 0; p < k; p_col += row, p++)
                        A[p_row + k] -= A[p_row + p] * A[p_col + k];
                }
                if ( A[p_k + k] == 0.0 ) {
                    printf("Singular matrix, returning zero matrix\n");
                    return GMatrix(row,col);
                }
                for (j = k+1; j < row; j++) {
                    for (p = 0, p_col = 0; p < k; p_col += row,  p++)
                        A[p_k + j] -= A[p_k + p] * A[p_col + j];
                    A[p_k + j] /= A[p_k + k];
                }
            }
            return GMatrix(row,col,A);
        }
    }
    /*! \brief Calculate the Inverse of the lower triangular matrix

        This routine calculates the inverse of the lower triangular matrix L.
        The superdiagonal part of the matrix is not addressed.
        The algorithm follows: \n
        Let M be the inverse of L, then \f$ L M = I \f$, \n
        \f$ M[i][i] = 1.0 / L[i][i] \f$ for i = 0, ..., row-1, and \n
        \f$ M[i][j] = -[(L[i][j] M[j][j] + ... + L[i][i-1] M[i-1][j])] / L[i][i], \f$
        for i = 1, ..., row-1, j = 0, ..., i - 1.
    */
    GMatrix Lower_Triangular_Inverse()
    {
        Double_t L[row*col];
        int i, j, k, p_i, p_j, p_k;
        Double_t sum;
        for (j=0;j<row;j++) for (k=0;k<=j;k++) L[j*col+k]=matrix[j*col+k];
        for (j=0;j<row;j++) for (k=j+1;k<col;k++) L[j*col+k]=0;
        //Invert the diagonal elements of the lower triangular matrix L.
        for (k = 0, p_k = 0; k < row; p_k += (row + 1), k++) {
            if (L[p_k] == 0.0) return GMatrix(row,row);//error, returning zero matrix
            else L[p_k] = 1.0 / L[p_k];
        }
        //Invert the remaining lower triangular matrix L row by row.
        for (i = 1, p_i = 0 + row; i < row; i++, p_i += row) {
            for (j = 0, p_j = 0; j < i; p_j += row, j++) {
                sum = 0.0;
                for (k = j, p_k = p_j; k < i; k++, p_k += row)
                    sum += L[p_i + k] * L[p_k + j];
                L[p_i + j] = - L[p_i + i] * sum;
            }
        }
        return GMatrix(row,row,L);
    }
    /*! \brief Calculate the Inverse of the lower triangular matrix

        This routine calculates the inverse of the unit upper triangular matrix U.
        The subdiagonal part of the matrix is not addressed.
        The diagonal is assumed to consist of 1's and is not addressed.
        The algorithm follows: \n
        Let M be the inverse of U, then \f$ U M = I \f$, \n
        \f$ M[i][j] = -( U[i][i+1] M[i+1][j] + ... + U[i][j] M[j][j] ), \f$
        for i = row-2, ... , 0,  j = row-1, ..., i+1.
    */
    GMatrix Unit_Upper_Triangular_Inverse()
    {
        Double_t U[row*col];
        int i, j, k, p_i, p_k;
        for (j=0;j<row;j++) for (k=j+1;k<col;k++) U[j*col+k]=matrix[j*col+k];
        for (j=0;j<row;j++) for (k=0;k<=j;k++) U[j*col+k]=(j==k);
        //         Invert the superdiagonal part of the matrix U row by row where
        //         the diagonal elements are assumed to be 1.0.
        for (i = row - 2, p_i = 0 + row * (row - 2); i >=0; p_i -= row, i-- ) {
            for (j = row - 1; j > i; j--) {
                U[p_i + j] = -U[p_i + j];
                for (k = i + 1, p_k = p_i + row; k < j; p_k += row, k++ )
                    U[p_i + j] -= U[p_i + k] * U[p_k + j];
            }
        }
        return GMatrix(row,row,U);
    }

    /// \brief Adjugate (classical adjoint) of matrix
    /// This is simply the transpose of the cofactor of the matrix
    GMatrix Adjugate() const
    {
        GMatrix cofactorT(row,col);
        for (int i=0;i<row;i++)
            for (int j=0;j<col;j++)
                cofactorT(j,i)=(ExtractMatrix(i,j)).Det()*pow(-1.,i+j);
        return cofactorT;
    }

    /// Inverse of matrix
    GMatrix Inverse() const
    {
        if (row!=col) {
            printf("only implemented for square, returning zeros\n");
            return GMatrix(row,col);
        }
        GMatrix LU(LUDecomp());
        Double_t d=1.0;
        for (int i=0;i<row;i++) d*=LU(i,i);
        if (d  == 0.0) return GMatrix(row, col);
        else {
            GMatrix Uinv=LU.Unit_Upper_Triangular_Inverse();
            GMatrix Linv=LU.Lower_Triangular_Inverse();
            return Uinv*Linv;
        }
    }

    /// Inverse of matrix using a pivot (needed in case diagonals are zero)
    GMatrix InversewithPivot() const
    {
        if (row!=col) {
            printf("only implemented for square, returning zeros\n");
            return GMatrix(row,col);
        }
        GMatrix A(*this);
        GMatrix Pivot(A.PivotInPlace());
        GMatrix LU(A.LUDecomp());
        Double_t d=1.0;
        for (int i=0;i<row;i++) d*=LU(i,i);
        if (d  == 0.0) return GMatrix(row, col);
        else {
            GMatrix Uinv=LU.Unit_Upper_Triangular_Inverse();
            GMatrix Linv=LU.Lower_Triangular_Inverse();
            return Uinv*Linv*Pivot;
        }
    }

    /// helpful rotation definition used in \ref Math::GMatrix::Jacobi (see Numerical recipes)
    inline void ROTATE(GMatrix &m, int i, int j, int k, int l, Double_t sinrot, Double_t tanrot) const {
        Double_t tempa=m(i,j),tempb=m(k,l);
        m(i,j)=tempa-sinrot*(tempb+tempa*tanrot);m(k,l)=tempb+sinrot*(tempa-tempb*tanrot);
    }
    /// Applies jacobi transformations to a copy of a matrix and returns effectively the eigenvalues of the matrix
    GMatrix Jacobi(Double_t tol=1e-2) const
    {
        GMatrix m(*this);
        Double_t eigen[row], b[row], z[row];
        GMatrix Prot(row);//identity
        int nrot=0;
        Double_t crot,srot,trot2,thres,t,theta,sm,temp1,temp2;

        //initialize vectors to diagonal elements of matrix
        for (int i=0;i<row;i++){ eigen[i]=b[i]=m(i,i);z[i]=0.;}

        do {
            //first check to see if off diagonal elements are zero
            // could check to within some precision relative to trace
            sm=0.;
            for (int ip=0;ip<row-1;ip++)
            for (int iq=ip+1;iq<col;iq++)
                sm+=sqrt(m(ip,iq)*m(ip,iq));
            if (sm==0.) return GMatrix(row,1,eigen);
            //set the threshold to be 0 until several rotations have been applied
            nrot<4? thres=tol*sm/(Double_t)(row*row): thres=0.;
            //now start applying Jacobi rotations
            for (int ip=0;ip<row-1;ip++){
                for (int iq=ip+1;iq<col;iq++){
                    //get absolute value of off diagonal element and check if it is small in comparison to
                    //threshold factor of 1/5 is from numerical recipes
                    temp1=sqrt(m(ip,iq)*m(ip,iq))*0.2;
                    //if smaller than tolerance set to 0. after for sweeps and skip the rotation.
                    if (nrot > 4 && (sqrt(eigen[ip]*eigen[ip])+temp1)==(sqrt(eigen[ip]*eigen[ip])+tol)
                        && (sqrt(eigen[iq]*eigen[iq])+temp1)==(sqrt(eigen[iq]*eigen[iq])+tol))
                        m(ip,iq)=0.0;
                    else if (thres<temp1*5.0) {
                        temp2=eigen[iq]-eigen[ip];
                        //find the solution to t^2+2*t*theta-1=0
                        //where theta=(cos^2-sin^2)/(2*s*c)=(m(iq,iq)-m(ip,ip))/(2*m(ip,iq))
                        //note that nr has ad hoc check if theta^2 would overflow the computer
                        /*if ((sqrt(temp2*temp2)+temp1) == sqrt(temp2*temp2)) t=(m(ip,iq))/temp2;
                        else {
                            theta=0.5*temp2/m(ip,iq);
                            t=1.0/(sqrt(theta*theta)+sqrt(1.0+theta*theta));
                            if (theta < 0.0) t = -t;
                        }*/
                        //instead I just determine the value of t
                        theta=0.5*temp2/m(ip,iq);
                        t=1.0/(sqrt(theta*theta)+sqrt(1.0+theta*theta));
                        if (theta < 0.0) t = -t;
                        //calculate cos and sin of rotation angle and tan(rot/2)
                        crot=1.0/sqrt(1.+t*t); srot=t*crot; trot2=srot/(1.0+crot);
                        //not alter values in matrix
                        temp2=t*m(ip,iq);
                        z[ip]-= temp2; z[iq] += temp2;
                        eigen[ip]-= temp2; eigen[iq] += temp2;
                        m(ip,iq)=0.0;
                        //now rotate the ipth values for 0<=j<ip
                        for (int j=0;j<ip;j++) {ROTATE(m,j,ip,j,iq,srot,trot2);}
                        //rotate the values for ip<j<iq
                        for (int j=ip+1;j<iq;j++) {ROTATE(m,ip,j,j,iq,srot,trot2);}
                        //rotate the values for iq<j<=row
                        for (int j=iq+1;j<row;j++) {ROTATE(m,ip,j,iq,j,srot,trot2);}
                        //store the rotations in Prot
                        for (int j=0;j<row;j++) {ROTATE(Prot,j,ip,j,iq,srot,trot2);}
                    }
                }
            }
            //now update vectors used in rotation and reinitialize b
            for (int ip=0;ip<row;ip++) {
                b[ip] += z[ip];
                eigen[ip]=b[ip];
                z[ip]=0.0;
            }
        }while (nrot++<10);
        printf("Too many iterations in routine jacobi, giving zeros\n");
        for (int i=0;i<row;i++) eigen[i]=0.;
        return GMatrix(row,1,eigen);
    }
    /// same as \ref Math::GMatrix::Jacobi above with different arguments and also returns Prot
    void Jacobi(GMatrix &eigenval, GMatrix &Prot, Double_t tol=1e-2) const
    {
        GMatrix m((*this));
        Double_t eigen[row], b[row], z[row];
        for (int i=0;i<row;i++) Prot(i,i)=1.;
        int nrot=0;
        Double_t crot,srot,trot2,thres,t,theta,sm,temp1,temp2;

        //initialize vectors to diagonal elements of matrix
        for (int i=0;i<row;i++){eigen[i]=b[i]=m(i,i);z[i]=0.;}

        do {
            //first check to see if off diagonal elements are zero
            // could check to within some precision relative to trace
            sm=0.;
            for (int ip=0;ip<row-1;ip++)
            for (int iq=ip+1;iq<col;iq++)
                sm+=sqrt(m(ip,iq)*m(ip,iq));
            if (sm==0.) break;
            //set the threshold to be 0 until several rotations have been applied
            nrot<4? thres=tol*sm/(Double_t)(row*row): thres=0.;
            //now start applying Jacobi rotations
            for (int ip=0;ip<row-1;ip++){
                for (int iq=ip+1;iq<col;iq++){
                    //get absolute value of off diagonal element and check if it is small in comparison to
                    //threshold factor of 1/5 is from numerical recipes
                    temp1=sqrt(m(ip,iq)*m(ip,iq))*0.2;
                    //if smaller than tolerance set to 0. after for sweeps and skip the rotation.
                    if (nrot > 4 && (sqrt(eigen[ip]*eigen[ip])+temp1)==(sqrt(eigen[ip]*eigen[ip])+tol)
                        && (sqrt(eigen[iq]*eigen[iq])+temp1)==(sqrt(eigen[iq]*eigen[iq])+tol))
                        m(ip,iq)=0.0;
                    else if (thres<temp1*5.0) {
                        temp2=eigen[iq]-eigen[ip];
                        //find the solution to t^2+2*t*theta-1=0
                        //where theta=(cos^2-sin^2)/(2*s*c)=(m(iq,iq)-m(ip,ip))/(2*m(ip,iq))
                        //note that nr has ad hoc check if theta^2 would overflow the computer
                        /*if ((sqrt(temp2*temp2)+temp1) == sqrt(temp2*temp2)) t=(m(ip,iq))/temp2;
                        else {
                            theta=0.5*temp2/m(ip,iq);
                            t=1.0/(sqrt(theta*theta)+sqrt(1.0+theta*theta));
                            if (theta < 0.0) t = -t;
                        }*/
                        //instead I just determine the value of t
                        theta=0.5*temp2/m(ip,iq);
                        t=1.0/(sqrt(theta*theta)+sqrt(1.0+theta*theta));
                        if (theta < 0.0) t = -t;
                        //calculate cos and sin of rotation angle and tan(rot/2)
                        crot=1.0/sqrt(1.+t*t); srot=t*crot; trot2=srot/(1.0+crot);
                        //not alter values in matrix
                        temp2=t*m(ip,iq);
                        z[ip]-= temp2; z[iq] += temp2;
                        eigen[ip]-= temp2; eigen[iq] += temp2;
                        m(ip,iq)=0.0;
                        //now rotate the ipth values for 0<=j<ip
                        for (int j=0;j<ip;j++) {ROTATE(m,j,ip,j,iq,srot,trot2);}
                        //rotate the values for ip<j<iq
                        for (int j=ip+1;j<iq;j++) {ROTATE(m,ip,j,j,iq,srot,trot2);}
                        //rotate the values for iq<j<=row
                        for (int j=iq+1;j<row;j++) {ROTATE(m,ip,j,iq,j,srot,trot2);}
                        //store the rotations in Prot
                        for (int j=0;j<row;j++) {ROTATE(Prot,j,ip,j,iq,srot,trot2);}
                    }
                }
            }
            //now update vectors used in rotation and reinitialize b
            for (int ip=0;ip<row;ip++) {
                b[ip] += z[ip];
                eigen[ip]=b[ip];
                z[ip]=0.0;
            }
        }while (nrot++<10);
        if (nrot>=10) {
            printf("Too many iterations in routine jacobi, giving zeros\n");
            for (int i=0;i<row;i++) {
                eigenval(i,0)=0.;
                for (int j=0;j<row;j++)
                Prot(i,j)=(i==j);
            }
        }
        for (int i=0;i<row;i++) eigenval(i,0)=eigen[i];
    }

    /*! \brief Calculate eigenvalues of matrix of rank N. Stored as a GMatrix, e1, e2 ...eN.

        The user must know what to expect from this function, whether real eigenvalues
        or complex ones.
        Have to implement the Francis QR algorithm or use successive Jacobi transformations
        to get eigenvalues and eigenvectors (note Jacobi slower than QR)
        but to ensure real eigenvalues, check that the matrix is symmetric
    */
    GMatrix Eigenvalues(Double_t tol=1e-2) const
    {
        if (isSymmetric()==0||isZero()==1) {
            printf("only implemented for nonzero symmetric matrices, returning null matrix\n");
            return GMatrix(0,0);
        }
        else {
            //using Jacobi transformations
            //issue is then must order eigen values
            GMatrix eigen(Jacobi(tol));
            //now order the eigen vaules e0>e1>e2..>eN
            //qsort(eigen.matrix,row,sizeof(Double_t),compare_Double);
            return eigen;
        }
    }

    // Computes the eigenvectors of the matrix, given the eigenvalues.  Sadly, this function will only work
    // properly on symmetric matrices -- i.e., all real eigenvalues.
    //GMatrix Eigenvectors(const Coordinate& e) const;

    /// Calculate eigenvalues \e and eigenvector of matrix of rank N. Stored as a GMatrix, e1, e2 ...eN.
    void Eigenvalvec(GMatrix &eigenval, GMatrix &eigenvector, Double_t tol=1e-2) const
    {
        if (isSymmetric()==0||isZero()==1) {
            printf("only implemented for nonzero symmetric matrices, returning null matrix\n");
            for (int i=0;i<row;i++) {
                eigenval(i,0)=0.;
                for (int j=0;j<row;j++)
                    eigenvector(i,j)=0.;
            }
        }
        else {
            //using Jacobi transformations
            Jacobi(eigenval,eigenvector,tol);
        }
    }

    //@}
};

}
#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
/*
 * @file modulated.h
 * @brief Cosine modulated analysis and synthesis filter banks.
 * @author John McDonough and Kenichi Kumatani
 */
#ifndef MODULATED_H
#define MODULATED_H

#include <stdio.h>
#include <assert.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include "common/jexception.h"

#include "stream/stream.h"
//#include "btk.h"

#ifdef HAVE_LIBFFTW3
#include <fftw3.h>
#endif

inline int powi(int x, int p)
{
  if(p == 0) return 1;
  if(x == 0 && p > 0) return 0;
  if(p < 0) {assert(x == 1 || x == -1); return (-p % 2) ? x : 1;}

  int r = 1;
  for(;;) {
    if(p & 1) r *= x;
    if((p >>= 1) == 0)	return r;
    x *= x;
  }
}

/**
* \defgroup FilterBanks Filter Banks
* This hierarchy of classes provides the capability to divide
* signal into 'M' subbands and then resynthesize the original time-domain signal.
*/
/*@{*/

// ----- definition for class `BaseFilterBank' -----
//
class BaseFilterBank {
 public:
  virtual ~BaseFilterBank();

  virtual void reset() = 0;

 protected:
  class RealBuffer_ {
  public:
    /*
      @brief Construct a circular buffer to keep samples periodically.
             It keeps nsamp arrays which is completely updated with the period 'nsamp'.
             Each array holds actual values of the samples.
      @param unsigned len [in] The size of each array
      @param unsigned nsamp [in] The period of the circular buffer
    */
    RealBuffer_(unsigned len, unsigned nsamp)
      : len_(len), nsamp_(nsamp), zero_(nsamp_ - 1), samples_(new gsl_vector*[nsamp_])
    {
      for (unsigned i = 0; i < nsamp_; i++)
	samples_[i] = gsl_vector_calloc(len_);
    }
    ~RealBuffer_()
    {
      for (unsigned i = 0; i < nsamp_; i++)
	gsl_vector_free(samples_[i]);
      delete[] samples_;
    }

    const double sample(unsigned timeX, unsigned binX) const {
      unsigned idx = index_(timeX);
      const gsl_vector* vec = samples_[idx];
      return gsl_vector_get(vec, binX);
    }

    void nextSample(const gsl_vector* s = NULL, bool reverse = false) {
      zero_ = (zero_ + 1) % nsamp_;

      gsl_vector* nextBlock = samples_[zero_];

      if (s == NULL) {
	gsl_vector_set_zero(nextBlock);
      } else {
	if (s->size != len_)
	  throw jdimension_error("'RealBuffer_': Sizes do not match (%d vs. %d)", s->size, len_);
	assert( s->size == len_ );
	if (reverse)
	  for (unsigned i = 0; i < len_; i++)
	    gsl_vector_set(nextBlock, i, gsl_vector_get(s, len_ - i - 1));
	else
	  gsl_vector_memcpy(nextBlock, s);
      }
    }

    void nextSample(const gsl_vector_float* s) {
      zero_ = (zero_ + 1) % nsamp_;

      gsl_vector* nextBlock = samples_[zero_];

      assert( s->size == len_ );
      for (unsigned i = 0; i < len_; i++)
	gsl_vector_set(nextBlock, i, gsl_vector_float_get(s, i));
    }

    void nextSample(const gsl_vector_short* s) {
      zero_ = (zero_ + 1) % nsamp_;

      gsl_vector* nextBlock = samples_[zero_];

      assert( s->size == len_ );
      for (unsigned i = 0; i < len_; i++)
	gsl_vector_set(nextBlock, i, gsl_vector_short_get(s, i));
    }

    void zero() {
      for (unsigned i = 0; i < nsamp_; i++)
	gsl_vector_set_zero(samples_[i]);
      zero_ = nsamp_ - 1;
    }

  private:
    unsigned index_(unsigned idx) const {
      assert(idx < nsamp_);
      unsigned ret = (zero_ + nsamp_ - idx) % nsamp_;
      return ret;
    }

    const unsigned				len_;
    const unsigned				nsamp_;
    unsigned					zero_; // index of most recent sample
    gsl_vector**				samples_;
  };

  BaseFilterBank(gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0, bool synthesis = false);

  const unsigned				M_;
  const unsigned				Mx2_;
  const unsigned				m_;
  const unsigned				mx2_;
  const unsigned				r_;
  const unsigned				R_;
  const unsigned				Rx2_;
  const unsigned				D_;
};

// ----- definition for class `NormalFFTAnalysisBank' -----
//
/**
   @class do FFT on time discrete samples multiplied with a window.
*/
class NormalFFTAnalysisBank
  : protected BaseFilterBank, public VectorComplexFeatureStream {
 public:
  NormalFFTAnalysisBank(VectorFloatFeatureStreamPtr& samp,
                        unsigned fftLen,  unsigned r = 1, unsigned windowType = 1,
                        const String& nm = "NormalFFTAnalysisBank");
  ~NormalFFTAnalysisBank();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  unsigned fftlen() const { return N_; }

#ifdef ENABLE_LEGACY_BTK_API
  unsigned fftLen() const { return fftlen(); }
#endif

protected:
  void update_buf_();
  virtual void update_buffer_(int frame_no);

#ifdef HAVE_LIBFFTW3
  fftw_plan                          fftwPlan_;
#endif
  const VectorFloatFeatureStreamPtr  samp_;
  int                                winType_; // 1 = hamming, 2 = hann window
  unsigned                           N_;       // FFT length
  const unsigned	             processing_delay_;
  unsigned			     framesPadded_;
  RealBuffer_			     buffer_;
  gsl_vector*			     convert_;
  RealBuffer_		             gsi_;
  double*			     output_;
  const gsl_vector*		     prototype_;
};

typedef Inherit<NormalFFTAnalysisBank, VectorComplexFeatureStreamPtr> NormalFFTAnalysisBankPtr;

/**
* \defgroup OversampledFilterBank Oversampled Filter Bank
*/
/*@{*/

// ----- definition for class `OverSampledDFTFilterBank' -----
//
class OverSampledDFTFilterBank : public BaseFilterBank {
public:
  ~OverSampledDFTFilterBank();

  virtual void reset();

  protected:
  OverSampledDFTFilterBank(gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0, bool synthesis = false, unsigned delayCompensationType=0, int gainFactor=1 );

  double polyphase(unsigned m, unsigned n) const {
    return gsl_vector_get(prototype_, m + M_ * n);
  }

  unsigned				laN_; /*>! the number of look-ahead */
  const unsigned			N_;
  unsigned				processing_delay_;
  const gsl_vector*			prototype_;
  RealBuffer_				buffer_;
  gsl_vector*				convert_;
  RealBuffer_				gsi_;
  const int                             gain_factor_;
};

/*@}*/

/**
* \defgroup PerfectReconstructionFilterBank Perfect Reconstruction Filter Bank
*/
/*@{*/

// ----- definition for class `PerfectReconstructionFilterBank' -----
//
class PerfectReconstructionFilterBank : public BaseFilterBank {
public:
  ~PerfectReconstructionFilterBank();

  virtual void reset();

protected:
  PerfectReconstructionFilterBank(gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0, bool synthesis = false);

  double polyphase(unsigned m, unsigned n) const {
    return gsl_vector_get(prototype_, m + Mx2_ * n);
  }

  const unsigned				N_;
  const unsigned				processing_delay_;

  const gsl_vector*				prototype_;
  RealBuffer_					buffer_;
  gsl_vector*					convert_;
  gsl_vector_complex*				w_;
  RealBuffer_					gsi_;
};

/*@}*/

/**
* \addtogroup OversampledFilterBank
*/
/*@{*/

// ----- definition for class `OverSampledDFTAnalysisBank' -----
//
class OverSampledDFTAnalysisBank
: protected OverSampledDFTFilterBank, public VectorComplexFeatureStream {
 public:
  OverSampledDFTAnalysisBank(VectorFloatFeatureStreamPtr& samp,
			     gsl_vector* prototype, unsigned M, unsigned m, unsigned r, unsigned delayCompensationType =0,
			     const String& nm = "OverSampledDFTAnalysisBank");
  ~OverSampledDFTAnalysisBank();

  virtual const gsl_vector_complex* next(int frame_no = -5);

  virtual void reset();

  unsigned fftlen() const { return M_; }
  unsigned shiftlen() const { return D_; }

#ifdef ENABLE_LEGACY_BTK_API
  unsigned fftLen()      const { return fftlen(); }
  unsigned nBlocks()     const { return m_; }
  unsigned subSampRate() const { return r_; }
#endif
  bool is_end(){return is_end_;}
  int frame_no() const { return frame_no_; }

  using OverSampledDFTFilterBank::polyphase;

 private:
  void update_buf_();
  bool update_buffer_(int frame_no);

#ifdef HAVE_LIBFFTW3
  fftw_plan					fftwPlan_;
#endif
  const VectorFloatFeatureStreamPtr		samp_;
  double*					polyphase_output_;
  unsigned					framesPadded_;
};

typedef Inherit<OverSampledDFTAnalysisBank, VectorComplexFeatureStreamPtr> OverSampledDFTAnalysisBankPtr;

// ----- definition for class `OverSampledDFTSynthesisBank' -----
//
class OverSampledDFTSynthesisBank
: private OverSampledDFTFilterBank, public VectorFloatFeatureStream {
 public:
  OverSampledDFTSynthesisBank(VectorComplexFeatureStreamPtr& samp,
			      gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0,
			      unsigned delayCompensationType = 0, int gainFactor=1,
			      const String& nm = "OverSampledDFTSynthesisBank");

  OverSampledDFTSynthesisBank(gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0,
			      unsigned delayCompensationType = 0, int gainFactor=1,
			      const String& nm = "OverSampledDFTSynthesisBank");

  ~OverSampledDFTSynthesisBank();

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  using OverSampledDFTFilterBank::polyphase;

  void input_source_vector(const gsl_vector_complex* block){ update_buf_(block); }
  void no_stream_feature(bool flag=true){ no_stream_feature_ = flag; }
#ifdef ENABLE_LEGACY_BTK_API
  void inputSourceVector(const gsl_vector_complex* block){ input_source_vector(block); }
  void doNotUseStreamFeature(bool flag=true){ no_stream_feature(flag); }
#endif

 private:
  bool update_buffer_(int frame_no);
  void update_buf_(const gsl_vector_complex* block);

  const VectorComplexFeatureStreamPtr		samp_;
  bool                                          no_stream_feature_;
#ifdef HAVE_LIBFFTW3
  fftw_plan					fftwPlan_;
#endif
  double*					polyphase_input_;
};

typedef Inherit<OverSampledDFTSynthesisBank, VectorFloatFeatureStreamPtr> OverSampledDFTSynthesisBankPtr;

/*@}*/

/**
* \addtogroup PerfectReconstructionFilterBank
*/
/*@{*/

// ----- definition for class `PerfectReconstructionFFTAnalysisBank' -----
//
class PerfectReconstructionFFTAnalysisBank
: protected PerfectReconstructionFilterBank, public VectorComplexFeatureStream {
 public:
  PerfectReconstructionFFTAnalysisBank(VectorFloatFeatureStreamPtr& samp,
				       gsl_vector* prototype, unsigned M, unsigned m, unsigned r,
				       const String& nm = "PerfectReconstructionFFTAnalysisBank");
  ~PerfectReconstructionFFTAnalysisBank();

  virtual const gsl_vector_complex* next(int frame_no = -5);

  virtual void reset();

#ifdef ENABLE_LEGACY_BTK_API
  unsigned fftLen()	 const { return Mx2_; }
  unsigned nBlocks()	 const { return 4; }
  unsigned subSampRate() const { return 2; }
#endif

  using PerfectReconstructionFilterBank::polyphase;

 protected:
  void update_buf_();
  virtual void update_buffer_(int frame_no);

#ifdef HAVE_LIBFFTW3
  fftw_plan					fftwPlan_;
#endif
  double*					polyphase_output_;
  unsigned					framesPadded_;

  const VectorFloatFeatureStreamPtr		samp_;
};

typedef Inherit<PerfectReconstructionFFTAnalysisBank, VectorComplexFeatureStreamPtr> PerfectReconstructionFFTAnalysisBankPtr;

// ----- definition for class `PerfectReconstructionFFTSynthesisBank' -----
//
class PerfectReconstructionFFTSynthesisBank
: private PerfectReconstructionFilterBank, public VectorFloatFeatureStream {
 public:
  PerfectReconstructionFFTSynthesisBank(VectorComplexFeatureStreamPtr& samp,
					gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0,
					const String& nm = "PerfectReconstructionFFTSynthesisBank");

  PerfectReconstructionFFTSynthesisBank(gsl_vector* prototype, unsigned M, unsigned m, unsigned r = 0,
					const String& nm = "PerfectReconstructionFFTSynthesisBank");

  ~PerfectReconstructionFFTSynthesisBank();

  virtual const gsl_vector_float* next(int frame_no = -5);

  virtual void reset();

  using PerfectReconstructionFilterBank::polyphase;

 private:
  void update_buffer_(int frame_no);
  void update_buffer_(const gsl_vector_complex* block);

  const VectorComplexFeatureStreamPtr		samp_;
#ifdef HAVE_LIBFFTW3
  fftw_plan					fftwPlan_;
#endif
  double*					polyphase_input_;
};

typedef Inherit<PerfectReconstructionFFTSynthesisBank, VectorFloatFeatureStreamPtr> PerfectReconstructionFFTSynthesisBankPtr;

// ----- definition for class `DelayFeature' -----
//
class DelayFeature : public VectorComplexFeatureStream {
 public:
  DelayFeature( const VectorComplexFeatureStreamPtr& samp, float time_delay=0.0, const String& nm = "DelayFeature");
  ~DelayFeature();

  void set_time_delay(float time_delay){ time_delay_ = time_delay; }
  virtual const gsl_vector_complex* next(int frame_no = -5);

  virtual void reset();

private:
  const VectorComplexFeatureStreamPtr samp_;
  float                               time_delay_;
};

typedef Inherit<DelayFeature, VectorComplexFeatureStreamPtr> DelayFeaturePtr;

gsl_vector* get_window(unsigned winType, unsigned winLen);

void write_gsl_format(const String& fileName, const gsl_vector* prototype);

/*@}*/

/*@}*/

#endif

ext: h


text:
/** @file
 * Based on the paper:
 * Quintana, F.A. (1998) ``Nonparametric Bayesian analysis for assessing homogeneity in $ k\times l$
 * contingency tables with fixed right margin totals''.
 * Also as Technical Report PUC/FM-96/7. Journal of the American Statistical Association, 93(443), 1140-1149.
 * Available as compressed postscript.
 * http://www.mat.puc.cl/~quintana/hcttr.ps.gz
 */

#ifndef BIO_CONTINGENCY_HOMOGENEITY_H_
#define BIO_CONTINGENCY_HOMOGENEITY_H_

#include "bio/defs.h"
#include "bio/partition.h"

#include <boost/assign/list_of.hpp>

#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_machine.h>
#include <gsl/gsl_sf_gamma.h>

#include <numeric>
#include <stdexcept>

BIO_NS_START

typedef std::vector<double> gamma_vec_t;
extern gamma_vec_t contingency_gamma_exponential_prior;

template <class LambdaIt>
double
calculate_ln_D(
	LambdaIt lambda_begin,
	LambdaIt lambda_end)
{
	double result = gsl_sf_lngamma(std::accumulate(lambda_begin, lambda_end, double(0)));
	for (LambdaIt i = lambda_begin; lambda_end != i; ++i)
	{
		result -= gsl_sf_lngamma(*i);
	}

	return result;
}

/** Calculate natural logarithm of a term for one partition in Lm(n). */
template <class PartitionIt, class CountIt, class LambdaIt>
double
contingency_calculate_ln_L_term(
	PartitionIt part_it,
	unsigned m,
	CountIt n_begin,
	CountIt n_end,
	LambdaIt lambda_begin,
	LambdaIt lambda_end)
{
	double ln_result = 1.0;

#ifdef BIO_CONTINGENCY_PRINT_PARTITIONS
	//can only do this for smallish numbers
	if (26 >= part_it.M.size())
	{
		std::vector<char> v;
		for (unsigned i = 0; i != part_it.M.size(); ++i)
		{
			v.push_back('a' + i);
		}
		std::cout << *part_it[v] << '\n';
	}
#endif //BIO_CONTINGENCY_PRINT_PARTITIONS

	//for each set in the partition
	for (unsigned i = 0; i < m; ++i)
	{
		//calculate lambda + n
		typedef std::vector<double> lambda_vec_t;
		lambda_vec_t lambda_plus_n;
		{
			std::copy(lambda_begin, lambda_end, std::inserter(lambda_plus_n, lambda_plus_n.begin()));
			unsigned e_i = 0;
			unsigned j = 0;
			//for each index
			for (CountIt n = n_begin; n_end != n; ++j, ++n)
			{
				//is this j in this partition?
				if (i == part_it.M[j])
				{
					for (unsigned k = 0; k != (*n).size(); ++k)
					{
						lambda_plus_n[k] += (*n)[k];
					}
					if (1 < e_i)
					{
						ln_result += gsl_sf_log(double(e_i)); //the factorial term
					}
					++e_i;
				}
			}
		}

		ln_result += calculate_ln_D(lambda_begin, lambda_end);
		ln_result -= calculate_ln_D(lambda_plus_n.begin(), lambda_plus_n.end());
	}

	return ln_result;
}

/** Calculate Lm(n). */
template <class CountIt, class LambdaIt>
double
contingency_calculate_L(
	unsigned m,
	CountIt n_begin,
	CountIt n_end,
	LambdaIt lambda_begin,
	LambdaIt lambda_end)
{
	unsigned set_size = unsigned(n_end - n_begin);
	BOOST_ASSERT(unsigned(lambda_end - lambda_begin) == set_size);

	double result = 0;

	//iterate over all the partitions
	partition::iterator_k part_it(set_size, m);
	try
	{
		while (true)
		{
			const double ln_L_term =
				contingency_calculate_ln_L_term(
					part_it,
					m,
					n_begin,
					n_end,
					lambda_begin,
					lambda_end);

			//check for underflow - in which case it will be 0.
			if (ln_L_term > GSL_LOG_DBL_MIN)
			{
				result += gsl_sf_exp(ln_L_term);
			}

			++part_it;
		}
	}
	catch (std::overflow_error&) //how the partition code determines there are no more partitions!
	{
	}

	return result;
}

template <class Value>
Value
factorial(Value value)
{
	Value result = Value(1);
	while (value != 0)
	{
		result *= value;
		--value;
	}
	return result;
}

typedef double(* contingency_cluster_prior)(double c, void * void_params);

struct ContingencyClusterPrior
{
	virtual double operator()(double c) = 0;
	virtual double get_ln_prior(double c);
};

struct ExponentialPrior : ContingencyClusterPrior
{
	double mean;
	ExponentialPrior(double mean = 1.0);
	double operator()(double c);
	double get_ln_prior(double c);
};

struct GammaParams
{
	unsigned m;
	unsigned k;
	ContingencyClusterPrior * prior;
};

double
gamma_exponential_prior_integrand(double c, void * void_params);

double
contingency_calculate_gamma(
	unsigned m,
	unsigned k,
	ContingencyClusterPrior * prior);

template <class GammaIt, class CountIt, class LambdaIt>
double
contingency_calculate_bayes_factor(
	GammaIt gamma_begin,
	GammaIt gamma_end,
	CountIt n_begin,
	CountIt n_end,
	LambdaIt lambda_begin,
	LambdaIt lambda_end)
{
	const unsigned k = unsigned(n_end == n_begin ? 0 : n_begin->end() - n_begin->begin());

	double result = contingency_calculate_L(1, n_begin, n_end, lambda_begin, lambda_end);

	result *= (1 - factorial(k-1) * (*gamma_begin));

	result /= factorial(k-1);

	double sum = 0;
	for (unsigned i = 2; k + 1 != i; ++i)
	{
		sum += *(gamma_begin + i - 1) * contingency_calculate_L(i, n_begin, n_end, lambda_begin, lambda_end);
	}
	result /= sum;

	return result;
}

/** Contains the data for a contingency table and defines types. */
template <unsigned K, unsigned L>
struct ContingencyTable
{
	typedef boost::array<unsigned, K> row_t;
	typedef boost::array<row_t, L> table_t;
	typedef boost::array<double, L> lambda_vec_t;

	table_t data;

	static unsigned get_k() { return K; }
	static unsigned get_l() { return L; }
};

template <unsigned K, unsigned L>
std::ostream &
operator<<(std::ostream & os, const ContingencyTable<K, L> & table)
{
	for (unsigned l = 0; L != l; ++l)
	{
		for (unsigned k = 0; K != k; ++k)
		{
			os << table.data[l][k] << '\t';
		}
		os << '\n';
	}
	return os;
}

BIO_NS_END

#endif //BIO_CONTINGENCY_HOMOGENEITY_H_

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
/*
 * Copyright (c) 2011-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include <lapacke.h>
#include "dplasma.h"
#include "dplasmatypes.h"
#include "parsec/data_dist/matrix/sym_two_dim_rectangle_cyclic.h"
#include "parsec/data_dist/matrix/two_dim_rectangle_cyclic.h"

/**
 *******************************************************************************
 *
 *  Generic case
 *
 *******************************************************************************
 */
static int
dplasma_zlatms_operator( parsec_execution_stream_t *es,
                         const parsec_tiled_matrix_dc_t *descA,
                         void *_A,
                         PLASMA_enum uplo, int m, int n,
                         void *args )
{
    int tempmm, tempnn, ldam, i;
    double            *cond = (double*)args;
    parsec_complex64_t *A    = (parsec_complex64_t*)_A;
    (void)es;

    tempmm = ((m)==((descA->mt)-1)) ? ((descA->m)-(m*(descA->mb))) : (descA->mb);
    tempnn = ((n)==((descA->nt)-1)) ? ((descA->n)-(n*(descA->nb))) : (descA->nb);
    ldam = BLKLDD( descA, m );

    if (m == n) {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR, lapack_const( uplo ), tempmm, tempnn,
            0., 0., A, ldam);

        /* D(i) = 1 - (i-1)/(N-1)*(1 - 1/COND) */
        if ( m == 0 ) {
            A[0] = 1.;
            i = 1;
        }
        else {
            i = 0;
        }
        if ( descA->n > 1 ) {
            double tmp = 1. / (*cond);
            double alp = ( 1. - tmp ) / ((double)( descA->n - 1 ));
            int minmn = dplasma_imin( tempmm, tempnn );
            for(; i < minmn; i++){
                A[i+i*ldam] = (parsec_complex64_t)( (double)(descA->n-(descA->nb*n+i+1)) * alp + tmp );
            }
        }
    } else {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR, 'A', tempmm, tempnn,
            0., 0., A, ldam);
    }
    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg - Generates a special test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *           - PlasmaGeneral:   Generate a general matrix
 *           - PlasmaSymmetric: Generate a symmetric matrix
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *          On exit, the matrix A generated.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg
 * @sa dplasma_dpltmg
 * @sa dplasma_spltmg
 *
 ******************************************************************************/
int
dplasma_zlatms( parsec_context_t *parsec,
                PLASMA_enum mtxtype, double cond,
                parsec_tiled_matrix_dc_t *A,
                unsigned long long int seed)
{
    two_dim_block_cyclic_t Q, T;
    int nodes, rank, mb, nb, m, n, mt, nt, P;
    int rc;

    /* Init the diagonal of A */
    {
        parsec_taskpool_t *tp;
        double *condptr = malloc(sizeof( double ));
        *condptr = cond;
        tp = parsec_apply_New( PlasmaUpperLower, A, dplasma_zlatms_operator, condptr );
        if ( tp != NULL ) {
            rc = parsec_context_add_taskpool(parsec, tp);
            rc = parsec_context_start( parsec );
            rc = parsec_context_wait( parsec );
            parsec_apply_Destruct(tp);
        }
        else {
            return -1;
        }
    }

    nodes = A->super.nodes;
    rank  = A->super.myrank;
    mb    = A->mb;
    nb    = A->nb;
    m     = A->m;
    n     = A->n;
    mt    = A->mt;
    nt    = A->nt;

    if ( A->dtype & two_dim_block_cyclic_type ) {
        P = ((two_dim_block_cyclic_t*)A)->grid.rows;
    }
    else if ( A->dtype & sym_two_dim_block_cyclic_type ) {
        P = ((sym_two_dim_block_cyclic_t*)A)->grid.rows;
    }
    else {
        P = 1;
    }

    /* Init the random matrix R */
    two_dim_block_cyclic_init( &Q, matrix_ComplexDouble, matrix_Tile,
                               nodes, rank, mb, nb, m, n, 0, 0, m, n, 1, 1, P );
    Q.mat = parsec_data_allocate((size_t)Q.super.nb_local_tiles *
                                (size_t)Q.super.bsiz *
                                (size_t)parsec_datadist_getsizeoftype(Q.super.mtype));
    parsec_data_collection_set_key((parsec_data_collection_t*)&Q, "Q");

    /* Init the T matrix */
    two_dim_block_cyclic_init( &T, matrix_ComplexDouble, matrix_Tile,
                               nodes, rank, 32, nb, mt*32, n, 0, 0, mt*32, n, 1, 1, P );
    T.mat = parsec_data_allocate((size_t)T.super.nb_local_tiles *
                                (size_t)T.super.bsiz *
                                (size_t)parsec_datadist_getsizeoftype(T.super.mtype));
    parsec_data_collection_set_key((parsec_data_collection_t*)&T, "T");

    if ( mtxtype == PlasmaGeneral ) {
        if ( m >= n ) {
            parsec_tiled_matrix_dc_t *subA = tiled_matrix_submatrix( A, 0, 0, n, n );
            parsec_tiled_matrix_dc_t *subQ = tiled_matrix_submatrix( (parsec_tiled_matrix_dc_t *)&Q,
                                                                0, 0, n, n );
            parsec_tiled_matrix_dc_t *subT = tiled_matrix_submatrix( (parsec_tiled_matrix_dc_t *)&T,
                                                                0, 0, nt*32, n );

            /* Multiply on the right by an unitary matrix */
            dplasma_zplrnt( parsec, 0, subQ, seed + 1 );
            dplasma_zgeqrf( parsec, subQ, subT );
            dplasma_zunmqr( parsec, PlasmaRight, PlasmaNoTrans,
                            subQ, subT, subA );

            /* Multiply on the left by an unitary matrix */
            dplasma_zplrnt( parsec, 0,
                            (parsec_tiled_matrix_dc_t *)&Q, seed );
            dplasma_zgeqrf( parsec,
                            (parsec_tiled_matrix_dc_t*)&Q,
                            (parsec_tiled_matrix_dc_t*)&T );
            dplasma_zunmqr( parsec, PlasmaLeft, PlasmaNoTrans,
                            (parsec_tiled_matrix_dc_t*)&Q,
                            (parsec_tiled_matrix_dc_t*)&T, A );

            free(subA); free(subQ); free(subT);
        }
        else {
            parsec_tiled_matrix_dc_t *subA = tiled_matrix_submatrix( A, 0, 0, m, m );
            parsec_tiled_matrix_dc_t *subQ = tiled_matrix_submatrix( (parsec_tiled_matrix_dc_t *)&Q,
                                                                0, 0, m, m );
            parsec_tiled_matrix_dc_t *subT = tiled_matrix_submatrix( (parsec_tiled_matrix_dc_t *)&T,
                                                                0, 0, mt*32, m );

            /* Multiply on the left by an unitary matrix */
            dplasma_zplrnt( parsec, 0, subQ, seed );
            dplasma_zgeqrf( parsec, subQ, subT );
            dplasma_zunmqr( parsec, PlasmaLeft, PlasmaNoTrans,
                            subQ, subT, subA );

            /* Multiply on the right by an unitary matrix */
            dplasma_zplrnt( parsec, 0,
                            (parsec_tiled_matrix_dc_t *)&Q, seed );
            dplasma_zgeqrf( parsec,
                            (parsec_tiled_matrix_dc_t*)&Q,
                            (parsec_tiled_matrix_dc_t*)&T );
            dplasma_zunmqr( parsec, PlasmaRight, PlasmaNoTrans,
                            (parsec_tiled_matrix_dc_t*)&Q,
                            (parsec_tiled_matrix_dc_t*)&T, A );
            free(subA); free(subQ); free(subT);
        }
    }
    else {
        assert( mtxtype == PlasmaHermitian );

        /* Init the unitary matrix */
        dplasma_zplrnt( parsec, 0,
                        (parsec_tiled_matrix_dc_t *)&Q, seed );
        dplasma_zgeqrf( parsec,
                        (parsec_tiled_matrix_dc_t*)&Q,
                        (parsec_tiled_matrix_dc_t*)&T );

        dplasma_zunmqr( parsec, PlasmaLeft, PlasmaNoTrans,
                        (parsec_tiled_matrix_dc_t*)&Q,
                        (parsec_tiled_matrix_dc_t*)&T, A );
        dplasma_zunmqr( parsec, PlasmaRight, PlasmaConjTrans,
                        (parsec_tiled_matrix_dc_t*)&Q,
                        (parsec_tiled_matrix_dc_t*)&T, A );
    }

    free(Q.mat);
    parsec_tiled_matrix_dc_destroy( (parsec_tiled_matrix_dc_t*)&Q );
    free(T.mat);
    parsec_tiled_matrix_dc_destroy( (parsec_tiled_matrix_dc_t*)&T );

    return 0;
}

ext: c

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
#ifndef ComplexMatrix_h
#define ComplexMatrix_h
#include <cmath>
#include <complex>
#include <iostream>
#include <fstream>
#include <iomanip>

#ifndef NOBLAS
    extern "C" {
        #include <cblas.h>
    }

//    #ifdef SINGLEPRECISION
//        #define phi_hemm  cblas_chemm
//        #define phi_symm  cblas_csymm
//        #define phi_her   cblas_cher
//        #define phi_gerc  cblas_cgerc
//        #define phi_copy  cblas_ccopy
//        #define phi_axpy  cblas_caxpy
//        #define phi_dotc_sub cblas_cdotc_sub
//        #define phi_dotu_sub cblas_cdotu_sub
//        #define phi_scal  cblas_zscal
//        #define phi_hemv  cblas_chemv
//        #define phi_gemv  cblas_cgemv
//    #else
//        #define phi_hemm  cblas_zhemm
//        #define phi_symm  cblas_zsymm
//        #define phi_her   cblas_zher
//        #define phi_gerc  cblas_zgerc
//        #define phi_copy  cblas_zcopy
//        #define phi_axpy  cblas_zaxpy
//        #define phi_dotc_sub cblas_zdotc_sub
//        #define phi_dotu_sub cblas_zdotu_sub
//        #define phi_scal  cblas_zscal
//        #define phi_hemv  cblas_zhemv
//        #define phi_gemv  cblas_zgemv
//    #endif

#else

    #ifndef CBLAS_ENUM_DEFINED_H
        #define CBLAS_ENUM_DEFINED_H
        enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102 };
        enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113,
                         AtlasConj=114};
        enum CBLAS_UPLO  {CblasUpper=121, CblasLower=122};
        enum CBLAS_DIAG  {CblasNonUnit=131, CblasUnit=132};
        enum CBLAS_SIDE  {CblasLeft=141, CblasRight=142};
    #endif

#endif

using namespace std;

#ifdef SINGLEPRECISION
typedef float Float;
typedef complex<float> Complex;
#else
typedef double Float;
typedef complex<double> Complex;
#endif

struct Index2D {
    int i;
    int j;
};

//#ifdef NOBLAS
void phi_gemm(const int N, const Complex *alpha, const Complex *A,
              const Complex *B, const Complex *beta, Complex *C);

void phi_hemm(const enum CBLAS_SIDE Side, const int N, const Complex *alpha, const Complex *A,
                 const Complex *B, const Complex *beta, Complex *C);

void phi_symm(const enum CBLAS_SIDE Side, const int N, const Complex *alpha, const Complex *A,
                 const Complex *B, const Complex *beta, Complex *C);

void phi_her(const int N, const Float alpha, const Complex *X,
                Complex *A);

void phi_gerc(const int N, const Complex *alpha, const Complex *X,
                 const Complex *Y, Complex *A);

void phi_copy(const int N, const Complex *X, const int incX,
                 Complex *Y, const int incY);

void phi_axpy(const int N, const Complex *alpha, const Complex *X,
                 const int incX, Complex *Y, const int incY);

void phi_dotc_sub(const int N, const Complex *X, const int incX,
                       const Complex *Y, const int incY, Complex *dotc);

void phi_dotu_sub(const int N, const Complex *X, const int incX,
                       const Complex *Y, const int incY, Complex *dotu);

void phi_scal(const int N, const Complex *alpha, Complex *X, const int incX);

void phi_hemv(const int N, const Complex *alpha, const Complex *A,
                 const Complex *X, const Complex *beta, Complex *Y);

void phi_gemv(const int N, const Complex *alpha, const Complex *A,
                 const Complex *X, const Complex *beta, Complex *Y);

//#endif

//Equal size square matrices only
Complex * outerProduct(const Complex *,const Complex *, const int);

Complex * toLiouville(const Complex *,const int);

void printMatrix(const Complex *, const int);

void printMatrix(const Complex *, const Index2D *, const int);

void writeMatrixToFile(const Complex *, const int, ofstream &);

#endif

max_stars_repo_path: PHI code/phi-source/src/ComplexMatrix.h

max_issues_repo_path: PHI code/phi-source/src/ComplexMatrix.h

max_forks_repo_path: PHI code/phi-source/src/ComplexMatrix.h


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_hpr (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int N = 2;
   float alpha = 0.1f;
   float Ap[] = { -0.273f, -0.499f, -0.305f, -0.277f, 0.238f, -0.369f };
   float X[] = { 0.638f, -0.905f, 0.224f, 0.182f };
   int incX = -1;
   float Ap_expected[] = { -0.26467f, 0.0f, -0.30718f, -0.245116f, 0.360607f, 0.0f };
   cblas_chpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr(case 1418) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr(case 1418) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 2;
   float alpha = 0.1f;
   float Ap[] = { -0.273f, -0.499f, -0.305f, -0.277f, 0.238f, -0.369f };
   float X[] = { 0.638f, -0.905f, 0.224f, 0.182f };
   int incX = -1;
   float Ap_expected[] = { -0.26467f, 0.0f, -0.30718f, -0.308884f, 0.360607f, 0.0f };
   cblas_chpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr(case 1419) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr(case 1419) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 2;
   float alpha = 0.1f;
   float Ap[] = { -0.273f, -0.499f, -0.305f, -0.277f, 0.238f, -0.369f };
   float X[] = { 0.638f, -0.905f, 0.224f, 0.182f };
   int incX = -1;
   float Ap_expected[] = { -0.26467f, 0.0f, -0.30718f, -0.245116f, 0.360607f, 0.0f };
   cblas_chpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr(case 1420) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr(case 1420) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 2;
   float alpha = 0.1f;
   float Ap[] = { -0.273f, -0.499f, -0.305f, -0.277f, 0.238f, -0.369f };
   float X[] = { 0.638f, -0.905f, 0.224f, 0.182f };
   int incX = -1;
   float Ap_expected[] = { -0.26467f, 0.0f, -0.30718f, -0.308884f, 0.360607f, 0.0f };
   cblas_chpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr(case 1421) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr(case 1421) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int N = 2;
   double alpha = 1;
   double Ap[] = { 0.265, 0.362, -0.855, 0.035, 0.136, 0.133 };
   double X[] = { -0.278, -0.686, -0.736, -0.918 };
   int incX = -1;
   double Ap_expected[] = { 1.64942, 0.0, -0.020644, -0.214692, 0.68388, 0.0 };
   cblas_zhpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr(case 1422) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr(case 1422) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 2;
   double alpha = 1;
   double Ap[] = { 0.265, 0.362, -0.855, 0.035, 0.136, 0.133 };
   double X[] = { -0.278, -0.686, -0.736, -0.918 };
   int incX = -1;
   double Ap_expected[] = { 1.64942, 0.0, -0.020644, 0.284692, 0.68388, 0.0 };
   cblas_zhpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr(case 1423) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr(case 1423) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 2;
   double alpha = 1;
   double Ap[] = { 0.265, 0.362, -0.855, 0.035, 0.136, 0.133 };
   double X[] = { -0.278, -0.686, -0.736, -0.918 };
   int incX = -1;
   double Ap_expected[] = { 1.64942, 0.0, -0.020644, -0.214692, 0.68388, 0.0 };
   cblas_zhpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr(case 1424) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr(case 1424) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 2;
   double alpha = 1;
   double Ap[] = { 0.265, 0.362, -0.855, 0.035, 0.136, 0.133 };
   double X[] = { -0.278, -0.686, -0.736, -0.918 };
   int incX = -1;
   double Ap_expected[] = { 1.64942, 0.0, -0.020644, 0.284692, 0.68388, 0.0 };
   cblas_zhpr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr(case 1425) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr(case 1425) imag");
     };
   };
  };

}


text:
/**
 * @file videofeature.h
 * @brief Audio-visual speech recognition front end.
 * @author Munir Georges, John McDonough, Friedrich Faubel
 */
#include "btk.h"

#ifdef AVFORMAT

#ifdef OPENCV

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sort_vector.h>

#ifdef HAVE_LIBFFTW3
#include <fftw3.h>
#endif

#undef HAVE_CONFIG_H

#include <opencv/cv.h>
#include <opencv/highgui.h>

#define HAVE_CONFIG_H

#include "stream/stream.h"
#include "common/mlist.h"

#ifndef VIDEOFEATURE_H
#define VIDEOFEATURE_H

// ----- definition for class VideoFeature -----
//
class VideoFeature : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  VideoFeature(int mode, unsigned width, unsigned height, const String& nm = "VideoFeature");

  ~VideoFeature();

  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

  void read(const String& fileName, int from, int _to);

  int width() const { return _width; }

  int height() const { return _height; }

  int frameNumber() const { return _frameNumber;}
  int Frames();

 private:
  CvCapture*					_capture;
  IplImage*					_frame;
  IplImage*					_gray;
  IplImage*					_hsv;

  unsigned					_width;
  unsigned					_height;

  gsl_vector_float*				_R;
  gsl_vector_float*				_G;
  gsl_vector_float*				_B;

  String					_filename;
  unsigned					_from;
  unsigned					_to;
  int						_mode;

  bool						_endOfSamples;
  string					_fileName;
  int						_frameNumber;
};

typedef Inherit<VideoFeature, VectorFloatFeatureStreamPtr> VideoFeaturePtr;

// ----- definition for class 'ImageROI' -----
//
class ImageROI : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageROI(const VectorFloatFeatureStreamPtr& src, unsigned width, unsigned height, unsigned x=0, unsigned y=0, unsigned w=2, unsigned h=3, const String& nm = "ImageROI");

  ~ImageROI();
  void setROI(int x,int y);

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

  unsigned width() const { return _width; }

  unsigned height() const { return _height; }

  unsigned x_pos() const { return _x_pos; }

  unsigned y_pos() const { return _y_pos; }

  unsigned w() const { return _w; }

  unsigned h() const { return _h; }

 private:
  VectorFloatFeatureStreamPtr			_src;
  bool						_endOfSamples;
  const unsigned				_width;
  const unsigned				_height;
  gsl_vector_float*				_srcVec;
  const unsigned				_w;
  const unsigned				_h;
  unsigned					_x_pos;
  unsigned					_y_pos;
};

typedef Inherit<ImageROI, VectorFloatFeatureStreamPtr> ImageROIPtr;

// ----- definition for class ImageSmooth -----
//
class ImageSmooth : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageSmooth(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, int smoothtype, int param1, int param2, const String& nm = "ImageSmooth");

  ~ImageSmooth();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  bool						_endOfSamples;
  unsigned					_width;
  unsigned					_height;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;
  IplImage*					_imgf;
  int						_smoothtype;
  int						_param1;
  int						_param2;
};

typedef Inherit<ImageSmooth, VectorFloatFeatureStreamPtr> ImageSmoothPtr;

// ----- definition for class ImageMorphology -----
//
class ImageMorphology : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageMorphology(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, int type, int param, const String& nm = "ImageMorphology");

  ~ImageMorphology();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  bool						_endOfSamples;
  unsigned					_width;
  unsigned					_height;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;
  IplImage*					_imgf;
  int						_type;
  int						_param;
};

typedef Inherit<ImageMorphology, VectorFloatFeatureStreamPtr> ImageMorphologyPtr;

// ----- definition for class ImageMorphologyEx -----
//
class ImageMorphologyEx : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageMorphologyEx(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, int type, int param, const String& nm = "ImageMorphologyEx");

  ~ImageMorphologyEx();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  bool						_endOfSamples;
  unsigned					_width;
  unsigned					_height;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;
  IplImage*					_imgf;
  IplImage*					_temp;
  int						_type;
  int						_param;
};

typedef Inherit<ImageMorphologyEx, VectorFloatFeatureStreamPtr> ImageMorphologyExPtr;

// ----- definition for class Canny -----
//
class Canny : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  Canny(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, int param0, int param1, int param2, const String& nm = "Canny");

  ~Canny();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  bool						_endOfSamples;
  unsigned					_width;
  unsigned					_height;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;
  IplImage*					_imgf;
  int						_param0;
  int						_param1;
  int						_param2;
};

typedef Inherit<Canny, VectorFloatFeatureStreamPtr> CannyPtr;

// ----- definition for class ImageThreshold -----
//
class ImageThreshold : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageThreshold(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, float param0, int param1, int param2, const String& nm = "Canny");

  ~ImageThreshold();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  bool						_endOfSamples;
  unsigned					_width;
  unsigned					_height;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;
  IplImage*					_imgf;
  int						_param0;
  int						_param1;
  int						_param2;
};

typedef Inherit<ImageThreshold, VectorFloatFeatureStreamPtr> ImageThresholdPtr;

/*Do not use this funktion within python - there is a problem with the window-handle or threading....*/
// ----- definition for class imageshow -----
//
class ImageShow : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageShow(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, const String& nm = "ImageShow");

  ~ImageShow();

  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  unsigned					_width;
  unsigned					_height;
  IplImage*					_img;
  gsl_vector_float*				_srcVec;
  bool						_endOfSamples;
};

typedef Inherit<ImageShow, VectorFloatFeatureStreamPtr> ImageShowPtr;

// ----- definition for class 'SaveImage' -----
//
class SaveImage {
 public:
  // sz := number of elements in _vector
  SaveImage(unsigned width, unsigned height);

  ~SaveImage();

  void save(const gsl_vector_float* V, const String& filename);
  void savedouble(const gsl_vector* V, const String& filename);

 private:
  unsigned					_width;
  unsigned					_height;
  IplImage*					_img;
  const String					_filename;
};

typedef refcount_ptr<SaveImage> SaveImagePtr;

// ----- definition for class 'ImageDetection' -----
//
class ImageDetection : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageDetection(VectorFloatFeatureStreamPtr& src, unsigned width, unsigned height, unsigned w, unsigned h,
		 const String& filename, double scale_factor, int min_neighbors, int flags, int min_sizeX, int min_sizeY,
		 const String& nm = "ImageDetection");

  ~ImageDetection();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  double _linearnext(gsl_vector * v);

  VectorFloatFeatureStreamPtr			_src;
  const unsigned				_width;
  const unsigned				_height;
  unsigned					_w;
  unsigned					_h;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;

  const String					_filename;
  CvHaarClassifierCascade*			_pCascade;
  CvMemStorage*					_pStorage;
  CvSeq*					_pRectSeq;
  double					_scale_factor;
  int						_min_neighbors;
  int						_flags;
  CvSize					_min_size;

  int						_x_pos;
  int						_y_pos;

  const int					_his;
  gsl_vector*					_x_his;
  gsl_vector*					_y_his;
  int						_tmp_x;
  int						_tmp_y;
  bool						_endOfSamples;
};

typedef Inherit<ImageDetection, VectorFloatFeatureStreamPtr> ImageDetectionPtr;

// ----- definition for class FaceDetection -----
//
class FaceDetection : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  FaceDetection(VectorFloatFeatureStreamPtr& src, unsigned width, unsigned height, int region,
		const String& filename_eye, double scale_factor_eye, int min_neighbors_eye, int flags_eye, int min_sizeX_eye, int min_sizeY_eye,
		const String& filename_nose, double scale_factor_nose, int min_neighbors_nose, int flags_nose, int min_sizeX_nose, int min_sizeY_nose,
		const String& filename_mouth, double scale_factor_mouth, int min_neighbors_mouth, int flags_mouth, int min_sizeX_mouth, int min_sizeY_mouth,
		const String& nm = "FaceDetection");

  ~FaceDetection();

  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  bool						_endOfSamples;
  VectorFloatFeatureStreamPtr			_src;
  gsl_vector_float*				_srcVec;
  unsigned					_width;
  unsigned					_height;
};

typedef Inherit<FaceDetection, VectorFloatFeatureStreamPtr> FaceDetectionPtr;

// ----- definition for class ImageCentering -----
//
class ImageCentering : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  ImageCentering(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height,	const String& nm = "ImageCentering");

  ~ImageCentering();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  bool                                          _endOfSamples;
  VectorFloatFeatureStreamPtr	                _src;
  gsl_vector_float*                             _srcVec;
  IplImage*					_img;
  IplImage*					_gmi;
  IplImage*					_out;
  unsigned                                      _width;
  unsigned                                      _height;
};

typedef Inherit<ImageCentering, VectorFloatFeatureStreamPtr> ImageCenteringPtr;

// ----- definition for class LinearInterpolation -----
//
class LinearInterpolation : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  LinearInterpolation(const VectorFloatFeatureStreamPtr& src, double fps_src, double fps_dest,
		      const String& nm = "LinearInterpolation");
  ~LinearInterpolation();

  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  VectorFloatFeatureStreamPtr			_src;
  const double					_DeltaTs;
  const double					_DeltaTd;

  int						_sourceFrameX;
  gsl_vector_float*				_x_n;
  bool						_endOfSamples;
};

typedef Inherit<LinearInterpolation, VectorFloatFeatureStreamPtr> LinearInterpolationPtr;

// ----- definition for class OpticalFlowFeature -----
//
class OpticalFlowFeature : public VectorFloatFeatureStream {
 public:
  OpticalFlowFeature(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, const char* filename, const String& nm = "OpticalFlowFeature");

  ~OpticalFlowFeature();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  float StringToFloat(const std::string& str_input);
  int StringToInt(const std::string& str_input);

  void disalloc_matrix(float **matrix, long nx, long ny);
  void alloc_matrix(float ***matrix, long nx, long ny);
  float ** _gray;
  float ** _ux;
  float ** _uy;
  typedef float fptype;
  typedef int itype;
  fptype **f1_c1;         /* in     : 1st image, channel 1*/
  fptype **f2_c1;         /* in     : 2nd image, channel 1*/
  fptype **f1_c2;         /* in     : 1st image, channel 2*/
  fptype **f2_c2;         /* in     : 2nd image, channel 2*/
  fptype **f1_c3;         /* in     : 1st image, channel 3*/
  fptype **f2_c3;         /* in     : 2nd image, channel 3*/
  fptype **u;             /* in+out : u component of flow field */
  fptype **v;             /* in+out : v component of flow field */
  itype  _nx;             /* in     : size in x-direction on current grid */
  itype  _ny;             /* in     : size in y-direction on current grid */
  itype  _bx;             /* in     : boundary size in x-direction */
  itype  _by;             /* in     : boundary size in y-direction */
  fptype hx;              /* in     : grid size in x-dir. on current grid */
  fptype hy;              /* in     : grid size in y-dir. on current grid */
  itype  m_type_d;        /* in     : type of data term */
  itype  m_type_s;        /* in     : type of smoothness term */
  fptype m_gamma_ofc1;    /* in     : weight of grey value constancy */
  fptype m_gamma_ofc2;    /* in     : weight of gradient constancy */
  fptype m_gamma_ofc3;    /* in     : weight of Hessian constancy */
  fptype m_gamma_gradnorm;/* in     : weight of gradient norm constancy */
  fptype m_gamma_laplace; /* in     : weight of Laplacian constancy */
  fptype m_gamma_hessdet; /* in     : weight of Hessian det. constancy */
  itype  m_function_d;    /* in     : type of robust function in data term */
  itype  m_function_s;    /* in     : type of robust function in smoothness term */
  fptype m_epsilon_d;     /* in     : parameter data term */
  fptype m_epsilon_s;     /* in     : parameter smoothness term */
  fptype m_power_d;       /* in     : exponent data term */
  fptype m_power_s;       /* in     : exponent smoothness term */
  fptype m_rhox;          /* in     : integration scale in x-direction*/
  fptype m_rhoy;          /* in     : integration scale in y-direction*/
  fptype m_rhox_F;        /* in     : FLOW integration scale in x-direction*/
  fptype m_rhoy_F;        /* in     : FLOW integration scale in y-direction*/
  fptype m_sigma_F;       /* in     : FLOW presmoothing scale */
  fptype m_alpha;         /* in     : smoothness weight */
  itype  m_color;         /* in     : color flag */
  itype  n_solver;        /* in     : general solver */
  itype  n_max_rec_depth; /* in     : max rec depth */
  itype  n_mg_res_prob;   /* in     : resample problem */
  itype  n_theta;         /* in     : discretisation parameter */
  itype  n_mg_cycles;     /* in     : number of cycles */
  itype  n_mg_solver;     /* in     : type of multigrid solver */
  itype  n_mg_pre_relax;  /* in     : number of pre relaxation steps */
  itype  n_mg_post_relax; /* in     : number of post relaxation steps */
  itype  n_mg_rec_calls;  /* in     : number of recursive calls */
  itype  n_mg_centered;   /* in     : intergrid transfer: 0-cell/1-vertex */
  itype  n_mg_mat_order_r;/* in     : matrix multigrid restriction order */
  itype  n_mg_cas_order_r;/* in     : cascadic multigrid restriction order */
  itype  n_mg_cas_order_p;/* in     : cascadic multigrid prolongation order */
  itype  n_mg_cor_order_r;/* in     : correcting multigrid restriction order */
  itype  n_mg_cor_order_p;/* in     : correcting multigrid prolongation order */
  itype  n_iter_out;      /* in     : outer iterations */
  itype  n_iter_in;       /* in     : inner iterations */
  fptype n_omega;         /* in     : overrelaxation parameter */
  fptype n_tau;           /* in     : time step size */
  fptype n_warp_eta;      /* in     : reduction factor */
  itype  n_warp_steps;    /* in     : warping steps per level */
  itype  n_warp_max_rec_depth; /*in : max warp rect depth */
  itype  n_warp_scaling;  /* in     : flag for scaling  */
  itype  n_warp_rem;      /* in     : flag for removing out of bounds estimates */
  itype  n_warp_centered; /* in     : warping intergrid transfer: 0-cell/1-vertex */
  itype  n_warp_order_r;  /* in     : correcting warping restriction order */
  itype  n_warp_order_p;  /* in     : correcting warping prolongation order */
  fptype n_tensor_eps;    /* in     : motion tensor normalisation factor */
  fptype **utruth;        /* in     : x-component of numerical ground truth */
  fptype **vtruth;        /* in     : y-component of numerical ground truth */
  fptype **uref;          /* in     : x-component of problem ground truth */
  fptype **vref;          /* in     : y-component of problem ground truth */
  fptype **error_mag;     /* in     : error magnitude map */
  itype  i_rel_res_flag;  /* in     : flag for relative residual computation */
  fptype *i_rel_res;      /* out    : relative residual */
  itype  i_rel_err_flag;  /* in     : flag for relative error computation */
  fptype *i_rel_err;      /* out    : relative error */
  itype  i_time_flag;     /* in     : flag for runtime computation */
  itype  *i_time_sec;     /* out    : consumed seconds */
  itype  *i_time_usec;    /* out    : consumed microseconds */
  itype  info_flag;       /* in     : flag for information */
  itype  info_step;       /* in     : stepsize for information */
  itype  write_flag;      /* in     : flag for writing out */
  itype  write_step;      /* in     : stepsize for writing out */
  itype  frame_nr;        /* in     : current frame number */

  bool						_endOfSamples;
  VectorFloatFeatureStreamPtr			_src;
  gsl_vector_float*				_srcVec;
  unsigned					_width;
  unsigned					_height;
};

typedef Inherit<OpticalFlowFeature, VectorFloatFeatureStreamPtr> OpticalFlowFeaturePtr;

// ----- definition for class SnakeImage -----
//
class SnakeImage : public VectorFloatFeatureStream {
 public:
  SnakeImage(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, float alpha, float beta, float gamma, const String& nm = "SnakeImage");

  ~SnakeImage();

  gsl_vector_float* data();
  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  bool                          		_endOfSamples;
  VectorFloatFeatureStreamPtr			_src;
  gsl_vector_float*				_srcVec;
  IplImage*					_img;
  unsigned                      		_width;
  unsigned                      		_height;
  float                       		        _alpha;
  float                       		        _beta ;
  float                       		        _gamma;

};

typedef Inherit<SnakeImage, VectorFloatFeatureStreamPtr> SnakeImagePtr;

// ----- definition for class PCAFeature -----
//
class PCAFeature : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  PCAFeature(const VectorFloatFeatureStreamPtr& src,unsigned width,unsigned height, const char *filename, const char *filename_mean, int n, int k, const String& nm = "PCAFeature");

  ~PCAFeature();

  virtual const gsl_vector_float* next(int frameX = -5);

  virtual void reset();

 private:
  void _floatToDouble(gsl_vector_float *src, gsl_vector * dest);
  void _doubleToFloat(gsl_vector *src, gsl_vector_float * dest);
  void _normalizeRows();

  VectorFloatFeatureStreamPtr			_src;

  const unsigned				_width;
  const unsigned				_height;
  const unsigned				_N;
  const unsigned				_M;

  gsl_matrix*					_evec;
  gsl_vector*					_mean;
  gsl_vector_float*				_srcVec;

  bool						_endOfSamples;
};

typedef Inherit<PCAFeature, VectorFloatFeatureStreamPtr> PCAFeaturePtr;

// ----- definition for class IPCAFeature -----
//
class IPCAFeature : public VectorFloatFeatureStream {
 public:
  // sz := number of elements in _vector
  IPCAFeature(const VectorFloatFeatureStreamPtr& src, unsigned width, unsigned height, const String& filename,  const String& filename_mean, int n, int k, const String& nm = "IPCAFeature");
  ~IPCAFeature();

  virtual const gsl_vector_float* next(int frameX = -5);
  const gsl_vector* get_mean() const { return _mean; }
  const gsl_matrix* get_vec(int i) const { return _evec; }

  virtual void reset();

 private:
  void _floatToDouble(gsl_vector_float* src, gsl_vector* dest);
  void _doubleToFloat(gsl_vector* src, gsl_vector_float* dest);
  void _normalizeRows();

  VectorFloatFeatureStreamPtr			_src;

  const unsigned				_width;
  const unsigned				_height;
  const unsigned				_N;
  const unsigned				_M;

  gsl_matrix*					_evec;
  gsl_vector*					_mean;
  gsl_vector_float*				_srcVec;

  bool						_endOfSamples;
};

typedef Inherit<IPCAFeature, VectorFloatFeatureStreamPtr> IPCAFeaturePtr;

// ----- definition for class PCAEstimator -----
//
class PCAEstimator{
 public:
  // sz := number of elements in _vector
  PCAEstimator(unsigned n,unsigned m);

  ~PCAEstimator();

  void clearaccu();
  void accumulate(gsl_vector_float *V);
  void estimate();
  void save(const String& filename);

  const gsl_vector* get_mean()   const { return _mean; }
  const gsl_matrix* get_vec()    const { return _evec; }
  const gsl_vector* get_eigval() const { return _eval; }

 private:
  void     _floatToDouble(gsl_vector_float* src, gsl_vector* dest);
  void     _doubleToFloat(gsl_vector* src, gsl_vector_float* dest);
  unsigned					_N;
  unsigned					_M;

  unsigned					_in_data;
  gsl_vector*					_sum;    //sum per row
  gsl_vector*					_mean;

  gsl_matrix*					_data;   // original data, one image in each row
  gsl_matrix*					_mdata;  // data without row mean
  gsl_matrix*					_mdatat; // data'

  gsl_matrix*					_cov;    // cov = 1/M * (mdata' * data)
  gsl_vector*					_eval;   //eigen values
  gsl_matrix*					_evec;   //eigen vector
};

typedef refcount_ptr<PCAEstimator> PCAEstimatorPtr;

// ----- definition for class PCAModEstimator -----
//
class PCAModEstimator {
 public:
  // sz := number of elements in _vector
  PCAModEstimator(unsigned n, unsigned m);

  ~PCAModEstimator();

  const gsl_vector* get_mean()    const { return _mean; }
  const gsl_matrix* get_data()    const { return _data; }
  const gsl_matrix* get_mdata()   const { return _mdata; }
  const gsl_vector* get_eigval()  const { return _eval; }
  const gsl_matrix* get_eigface() const { return _eigface; }
  const gsl_matrix* get_evecc()   const { return _evecc; }
  const gsl_matrix* get_evec()    const { return _evec; }
  const gsl_matrix* get_cov()     const { return _cov; }

  void clearaccu();
  void accumulate(gsl_vector_float* V);
  void estimate();
  void save(const char* filename1, const char* filename2);

 private:
  void _floatToDouble(gsl_vector_float* src, gsl_vector* dest);
  void _doubleToFloat(gsl_vector* src, gsl_vector_float* dest);
  void _normalizeRows();
  void _normalizeColumn();

  const unsigned				_N;
  const unsigned				_M;

  unsigned					_indata;
  gsl_matrix*					_data;   // original data, one image in each row
  gsl_matrix*					_mdata;  // data without row mean
  gsl_vector*					_sum;    // sum per row
  gsl_vector*					_mean;   // sum per row
  gsl_matrix*					_mdatat; // data'

  gsl_matrix*					_cov;    // cov = 1/M * (mdata' * data)
  gsl_vector*					_eval;   // eigen values
  gsl_matrix*					_evec;   // eigen vector
  gsl_matrix*					_evecc;  // eigen vector
  gsl_matrix*					_evecct; // eigen vector
  gsl_matrix*					_eigface;// eigen vector
};

typedef refcount_ptr<PCAModEstimator> PCAModEstimatorPtr;

#endif

#endif

#endif


text:
#ifndef wf_H_INCLUDED
#define wf_H_INCLUDED

//#define DEBUG

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include "alphas.h"

//*
// Global variables for UGD for nucleus 1 and 2
//*

//run at fixed x
//#define fixedx (0.001)

//number of tabulated points in Y and kT for wavefunction
#define NY 151
#define Nkt 101
#define kT_min (0.02)   // phi=0 below this min kt
#define kT_max (18.)    // phi has a power lower extrapolation
                        // beyond this max kt

#define x0 (0.01)       // separates small/large x
#define XMIN (x0*exp(-15.)) //smallest x tabulated

// Global variables for UGD for nucleus 1 and 2
int wfTAG;
double uGD1[NY][Nkt];
double uGD2[NY][Nkt];
double Y[NY];
double kT[NY][Nkt];
gsl_interp_accel *accWF1[NY];
gsl_spline *WF1[NY];
gsl_interp_accel *accWF2[NY];
gsl_spline *WF2[NY];
gsl_interp_accel *accLargeX;
gsl_spline *LargeX;
//global variable for wfTAG=2
#define QS02 (0.1)
double QS21, QS22;

void ReadInWF(int A1, int A2, int TAG);
double wf(int nucleus, double x, double kT) ;
double wfsimple(int nucleus, double x, double kT) ;
void PrintWF(double x, const char *out) ;

#endif

max_stars_repo_name: kdusling/mpc

max_issues_repo_name: kdusling/mpc

max_forks_repo_name: kdusling/mpc