probability.h

#ifndef PROBABILITY_H
#define PROBABILITY_H

//includes
#include <cmath>
#include <string>
#include <sstream>
#include "mystringtools.h"

//Using linspace interpolation

//matrix indexing:   0 <= m <= M
//                   0 <= n <= N

//density or domain indexing:   -mCenter <= m <= M-nCenter
//                              -nCenter <= n <= N-nCenter

class Density{
    public:
        Density(int n,int m):
                                matrixPtr(NULL),
                                rowPtr(NULL),   //Sum over all Columns
                                colPtr(NULL),
                                expectedNgivenM_ptr(NULL),
                                expectedNNgivenM_ptr(NULL),
                                expectedMgivenN_ptr(NULL),
                                expectedMMgivenN_ptr(NULL),
                                nCenter(n),
                                mCenter(m),
                                totalWeight(1.f), /*inital weight of uniform distribution (equivalent)
                                                to one sampling point */
                                M(2*m+1),
                                N(2*n+1),
                                expectedN(0),
                                expectedM(0),
                                expectedMN(0),
                                expectedNN(0),
                                expectedMM(0)
        {

            matrixPtr = new float*[M];
            for (m = 0; m < M; m++ ){
                dynamicallyExtendVector(&matrixPtr[m],N,0,1.f/float(M*N));
            }

            dynamicallyExtendVector(&rowPtr,N,0,1.f/float(N));
            dynamicallyExtendVector(&colPtr,M,0,1.f/float(M));

            initMoments();
        }

        /*
            Domain Indexing
        */
        float prob_x1_x2(float m, float n) const {
            // std::cout << "m: " << m << std::endl;
            // std::cout << "n: " << n << std::endl;

            int nCeil = n+1;
            int nFloor = n;
            int mCeil = m+1;
            int mFloor = m;

            // std::cout << "nCeil: " << nCeil << std::endl;
            // std::cout << "mCeil: " << mCeil << std::endl;
            // std::cout << "nFloor: " << nFloor << std::endl;
            // std::cout << "mFloor: " << mFloor << std::endl;

            // std::cout << "(nCeil-n): " << nCeil-n << std::endl;
            // std::cout << "(n-nFloor): " << n-nFloor << std::endl;
            // std::cout << "(mCeil-m): " << mCeil-m << std::endl;
            // std::cout << "(m-mFloor): " << m-mFloor << std::endl;


            float prob = matrixPtr[mCeil][nCeil]*(m-mFloor)*(n-nFloor);
            prob +=  matrixPtr[mCeil][nFloor]*(m-mFloor)*(nCeil-n);
            prob +=  matrixPtr[mFloor][nCeil]*(mCeil-m)*(n-nFloor);
            prob +=  matrixPtr[mFloor][nFloor]*(mCeil-m)*(nCeil-n);

            return prob/totalWeight;
        }
        float prob_x1(float m) const {
            if (m <= -mCenter-1 || m >= M - mCenter + 1)
                return 0;
            else {
                int mCeil = m+1;
                int mFloor = m;
                if (m <= -mCenter)
                    return rowPtr[mFloor+mCenter]*(mCeil-m)/totalWeight;
                else if (-mCenter < m && m <= M - mCenter){
                    return ( rowPtr[mCeil+mCenter]*(m-mFloor) + rowPtr[mFloor+mCenter]*(mCeil-m) )/totalWeight;
                } else if (m > M-mCenter )
                    return rowPtr[mCeil+mCenter]*(m-mFloor)/totalWeight;
            }
        }

        float prob_x2(float n) const {
            if (n <= -nCenter-1 || n >= N - nCenter + 1)
                return 0;
            else {
                int nCeil = n+1;
                int nFloor = n;
                if (n <= -nCenter)
                    return rowPtr[nFloor+nCenter]*(nCeil-n)/totalWeight;
                else if (-nCenter < n && n <= N - nCenter){
                    return ( rowPtr[nCeil+nCenter]*(n-nFloor) + rowPtr[nFloor+nCenter]*(nCeil-n) )/totalWeight;
                } else if (n > N-nCenter )
                    return rowPtr[nCeil+nCenter]*(n-nFloor)/totalWeight;
            }
        }

        float prob_x1_given_x2(float m, float nGiven) const {
            float x2_prob = prob_x2(nGiven);

            if (x2_prob > 0.f)
                return prob_x1_x2(m,nGiven)/x2_prob;

            return 0.f;
        }

        float prob_x2_given_x1(float mGiven, float n) const {
            float x1_prob = prob_x1(n);

            if (x1_prob > 0.f)
                return prob_x1_x2(mGiven,n)/x1_prob;
            else

            return 0.f;
        }

        void addPointToDensity(float m, float n, int weight = 1.f){
            // std::cout << "m: " << m << std::endl;
            // std::cout << "n: " << n << std::endl;

            int nCeil = n+1;
            int nFloor = n;
            int mCeil = m+1;
            int mFloor = m;

            // std::cout << "nCeil: " << nCeil << std::endl;
            // std::cout << "mCeil: " << mCeil << std::endl;
            // std::cout << "nFloor: " << nFloor << std::endl;
            // std::cout << "mFloor: " << mFloor << std::endl;

            // std::cout << "(nCeil-n): " << nCeil-n << std::endl;
            // std::cout << "(n-nFloor): " << n-nFloor << std::endl;
            // std::cout << "(mCeil-m): " << mCeil-m << std::endl;
            // std::cout << "(m-mFloor): " << m-mFloor << std::endl;

            addCellToMatrix(mCeil+mCenter,  nCeil+nCenter,  weight*(m-mFloor)*(n-nFloor));
            addCellToMatrix(mCeil+mCenter,  nFloor+nCenter, weight*(m-mFloor)*(nCeil-n));
            addCellToMatrix(mFloor+mCenter, nCeil+nCenter,  weight*(mCeil-m)*(n-nFloor));
            addCellToMatrix(mFloor+mCenter, nFloor+nCenter, weight*(mCeil-m)*(nCeil-n));
        }

        /*
            DOMAIN INDEXING
        */
        void extendDomainToCell(int m, int n) {
            if (n < -nCenter)
                extendDomainHorizontally(n+nCenter);//negative
            else if (n > N-nCenter)
                extendDomainHorizontally(n+nCenter-N);//positive

            if (m < -mCenter)
                extendDomainVertically(m+mCenter);//negative
            else if (m > M-mCenter)
                extendDomainVertically(m+mCenter-M);//positive
        }

        /*
            MATRIX INDEXING
        */
        void addCellToMatrix(int m, int n, float weight = 1.f){
            //fix

            float nCentralized = n-nCenter;
            float mCentralized = m-mCenter;
            float nCentralizedMoment = nCentralized*weight;
            float mCentralizedMoment = mCentralized*weight;


            extendDomainToCell(mCentralized, nCentralized);
            //std::cout << "m: " << m << ", n: " << n << ", weight: " << weight << std::endl;

            rowPtr[n] += weight;
            colPtr[m] += weight;
            totalWeight += weight;
            matrixPtr[m][n] += weight;

            expectedM += mCentralizedMoment;
            expectedN += nCentralizedMoment;

            expectedMN += mCentralized*nCentralizedMoment;
            expectedNN += nCentralized*nCentralizedMoment;
            expectedMM += mCentralized*mCentralizedMoment;

            //new
            expectedNgivenM_ptr[m] += nCentralizedMoment;
            expectedNNgivenM_ptr[m] += nCentralized*nCentralizedMoment;

            expectedMgivenN_ptr[n] += mCentralizedMoment;
            expectedMMgivenN_ptr[n] += mCentralized*mCentralizedMoment;
        }

        void extendDomainHorizontally(int horizontalExtension){
            if (horizontalExtension == 0) return;

            dynamicallyExtendVector(&rowPtr,horizontalExtension,N,0.f);
            dynamicallyExtendVector(&expectedMgivenN_ptr,horizontalExtension,N,0.f);
            dynamicallyExtendVector(&expectedMMgivenN_ptr,horizontalExtension,N,0.f);

            for (int m = 0; m < M; m++){
                dynamicallyExtendVector(&matrixPtr[m],horizontalExtension,N,0.f);
            }
            nCenter += (horizontalExtension < 0) ? std::abs(horizontalExtension) : 0;
            N += std::abs(horizontalExtension);
        }
        void extendDomainVertically(int verticalExtension){
            if (verticalExtension == 0) return;

            dynamicallyExtendVector(&colPtr,verticalExtension,M,0.f);
            dynamicallyExtendVector(&expectedNgivenM_ptr,verticalExtension,M,0.f);
            dynamicallyExtendVector(&expectedNNgivenM_ptr,verticalExtension,M,0.f);

            int newLength = M + std::abs(verticalExtension);
            float ** tempPtr = new float*[newLength];

            if (verticalExtension > 0){
                for (int m = 0; m < M; m++){
                    tempPtr[m] = matrixPtr[m];
                    matrixPtr[m] = NULL;
                }
                for (int m = M; m < newLength; m++){
                    tempPtr[m] = new float[N];
                    dynamicallyExtendVector(&tempPtr[m],newLength,0,0.f);
                }
            } else {

                mCenter += std::abs(verticalExtension);
                for ( int m = 0; m < std::abs(verticalExtension); m++ ){
                    dynamicallyExtendVector(&tempPtr[m],newLength,0,0.f);
                }
                for ( int m = std::abs(verticalExtension); m < newLength; m++ ){
                    tempPtr[m] = matrixPtr[m - (int)std::abs(verticalExtension)];
                    matrixPtr[m] = NULL;
                }
            }

            delete [] matrixPtr;
            matrixPtr = tempPtr;

            M = newLength;
        }

        const char* toString(const char* filler = "\t") const {
            std::stringstream ss("DynamicMatrix instance address: ");

            ss << filler << "DynamicMatrix instance address: " << this << std::endl;
            ss << filler << "\trows M:       " << M << std::endl;
            ss << filler << "\trows N:       " << N << std::endl << std::endl;
            ss << filler << "\trows mCenter: " << mCenter << std::endl;
            ss << filler << "\trows nCenter: " << nCenter << std::endl << std::endl;
            ss << filler << "\ttotal weight: " << totalWeight << std::endl;

            if (M > 0  && M > 0){
                ss << "\n\t" << filler << "matrixPtr[...][...]" << std::endl;

                ss << StringTools<std::string>::matrix2string(matrixPtr,M,mCenter,N,nCenter,"\t\t\t");

                ss << StringTools<std::string>::rowVector2string("rowPtr[...]",rowPtr,N,filler);
                ss << StringTools<std::string>::rowVector2string("expectedMgivenN_ptr[...]",expectedMgivenN_ptr,N,filler);
                ss << StringTools<std::string>::rowVector2string("expectedMMgivenN_ptr[...]",expectedMMgivenN_ptr,N,filler);

                ss << StringTools<std::string>::columnVector2string("colPtr[...]",colPtr,M,filler);
                ss << StringTools<std::string>::columnVector2string("expectedNgivenM_ptr[...]",expectedNgivenM_ptr,M,filler);
                ss << StringTools<std::string>::columnVector2string("expectedNNgivenM_ptr[...]",expectedNNgivenM_ptr,M,filler);

            }

            ss << "\n";

            return ss.str().c_str();
        }

        float getExpectedM() const {
            return expectedM;
        }

        float getExpectedN() const {
            return expectedN;
        }

        float getExpectedMN() const {
            return expectedMN;
        }
        float getExpectedMM() const {
            return expectedMM;
        }
        float getExpectedNN() const {
            return expectedNN;
        }

        float* getColPtr() const {
            return colPtr;
        }

        float* getRowPtr() const {
            return rowPtr;
        }

        float** getMatrixPtr () const {
            return matrixPtr;
        }

        float* getExpectedNgivenM_ptr() const {
            return expectedNgivenM_ptr;
        }

        float* getExpectedMgivenN_ptr() const {
            return expectedMgivenN_ptr;
        }

        float* getExpectedNNgivenM_ptr() const {
            return expectedNNgivenM_ptr;
        }

        float* getExpectedMMgivenN_ptr() const {
            return expectedMMgivenN_ptr;
        }
        int getN() const {
            return N;
        }

        int getM() const {
            return M;
        }

        int getnCenter() const {
            return nCenter;
        }

        int getmCenter() const {
            return mCenter;
        }
        float getTotalWeight() const {
            return totalWeight;
        }


    private:

        //  -Dynamically allocates a new vector and replaces values
        //      appropriately from the old matrix of shorter length,
        //      and appropriately scales the new subdomain as
        //      an appropriate percent of the probability.
        void dynamicallyExtendVector(   float ** ptr,
                                        int horizontalExtension,
                                        std::size_t oldLength,
                                        const float initialValue){

            if (horizontalExtension == 0) return;
            if (*ptr == NULL)
                oldLength = 0;

            int newLength = std::abs(horizontalExtension) + oldLength;

            float * tempPtr = new float[newLength];

            if (horizontalExtension > 0){
                for ( int k = 0; k < oldLength; k++ ){
                    tempPtr[k] = (*ptr)[k];
                }
                for ( int k = oldLength; k < newLength; k++ ){
                    tempPtr[k] = initialValue;
                }
            } else if (horizontalExtension < 0){
                for ( int k = 0; k < std::abs(horizontalExtension); k++ ){
                    tempPtr[k] = initialValue;
                }
                for ( int k = std::abs(horizontalExtension); k < newLength; k++ ){
                    tempPtr[k] = (*ptr)[k-(int)std::abs(horizontalExtension)];
                }
            }

            delete [] (*ptr);
            (*ptr) = tempPtr;
        }

        void initMoments(){
            dynamicallyExtendVector(&expectedNgivenM_ptr,M,0,0.f);
            dynamicallyExtendVector(&expectedMgivenN_ptr,N,0,0.f);

            //Not 0 - must integrate over
            dynamicallyExtendVector(&expectedNNgivenM_ptr,M,0,1.f/float(M));
            dynamicallyExtendVector(&expectedMMgivenN_ptr,N,0,1.f/float(N));

            for (int m = 0; m < M; m++){
                for (int n = 0; n < N; n++){
                    expectedNNgivenM_ptr[m] *= (n-nCenter)*(n-nCenter);
                    expectedMMgivenN_ptr[n] *= (m-mCenter)*(m-mCenter);
                }
            }
        }


        float expectedN;
        float expectedM;
        float expectedMN;
        float expectedNN;
        float expectedMM;

        float * colPtr;
        float * rowPtr;
        float ** matrixPtr;
        float * expectedNgivenM_ptr;
        float * expectedNNgivenM_ptr;
        float * expectedMgivenN_ptr;
        float * expectedMMgivenN_ptr;

        //domain control
        int N;
        int M;
        int nCenter;
        int mCenter;
        float totalWeight;

};

    //Extension
    #ifdef MATLAB_WRAPPER_H
        #include "probability_and_matlab_extensions.h"
    #endif

#endif