potential.h

gibbs

#ifndef POTENTIAL_H
#define POTENTIAL_H

#include "Distribution.h"
#include "VarSchema.h"
#include "DecisionTree.h"

class Potential
{
    vector<vector<double> > probs;

public:
    Potential() { /* NOP */ }

    Potential(const DecisionTree& dtree, const VarSchema& schema)
        : probs(dtree.maxVal)
    {
        // Create array of values for calling getProb
        double allValues[schema.getNumVars()];
        for (int v = 0; v < schema.getNumVars(); v++) {
            allValues[v] = 0;
        }

        // Iterate through each configuration of input variables
        int config = 0;
        bool keepGoing = false;
        do {
            // Record distribution at this configuration
            for (int output = 0; output < dtree.maxVal; output++) {
                probs[output].push_back(dtree.getProb(output, allValues));
            }

            // Go to next configuration
            config++;
            keepGoing = false;
            list<int>::const_iterator i;
            for (i = dtree.inputs.begin(); i != dtree.inputs.end(); i++) {

                allValues[*i]++;
                if (allValues[*i] >= schema.getRange(*i)) {
                    allValues[*i] = 0.0;
                } else {
                    /* When we find a variable that hasn't wrapped,
                     * we can stop.  If no such variable exists, then
                     * we must have reset our state to the initial one
                     * and we're done.
                     */
                    keepGoing = true;
                    break;
                }
            }
        } while(keepGoing);
    }

    Potential(int rows, int cols)
        :probs(rows, vector<double>(cols, 1.0))
    { /* NOP */ }

    Potential(double** array, int rows, int cols)
        :probs(rows, vector<double>(cols))
    {
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
                probs[r][c] = array[r][c];
            }
        }
    }

    Potential(vector<vector<double> > vals)
        :probs(vals)
    { /* NOP */ }

    Potential(const Potential& other)
        :probs(other.probs) 
    { /* NOP */ }

    Potential& operator=(const Potential& other) 
    { 
        probs = other.probs;
        return *this;
    }

    int dim1() const { return probs.size(); }
    int dim2() const { return probs[0].size(); }

    double get(int i, int j) const 
    { return probs[i][j]; }

    void set(int i, int j, double val) 
    { probs[i][j] = val; }
};

ostream& operator<<(ostream& out, const Potential& p);

#if 1
Distribution operator*(const Distribution& dist, const Potential& potential);
Distribution operator*(const Potential& potential, const Distribution& dist);

#else
inline Distribution operator*(const Distribution& dist, const Potential& potential)
{
    Distribution ret(vector<double>(potential.dim2(), 0.0));

    if (dist.dim() != potential.dim1()) {
        // DEBUG
        cout << "ERROR: dimension " << dist.dim() << " != " 
            << potential.dim1() << endl;

        // TODO: throw an exception?
        return ret;
    }

    for (int j = 0; j < potential.dim2(); j++) {
        for (int i = 0; i < potential.dim1(); i++) {
            ret[j] += dist[i] * potential.get(i,j);
        }
    }

    return ret;
}

inline Distribution operator*(const Potential& potential, const Distribution& dist)
{
    Distribution ret(vector<double>(potential.dim1(), 0.0));

    if (dist.dim() != potential.dim2()) {
        // DEBUG
        cout << "ERROR: dimension " << dist.dim() << " != " 
            << potential.dim2() << endl;

        return ret;
    }

    for (int i = 0; i < potential.dim1(); i++) {
        for (int j = 0; j < potential.dim2(); j++) {
            ret[i] += dist[j] * potential.get(i,j);
        }
    }

    return ret;
}
#endif

#endif // ndef POTENTIAL_H