bayesnet.cpp

gibbs

#include <string>
#include <iostream>
#include <sstream>
#include <expat.h>
#include <list>
#include "BayesNet.h"
#include "DecisionTree.h"

// HACK -- leaving this as module data prevents us from
// parsing more than one BayesNet at a time.
string elementCdata = "";
DecisionTree* dtree = NULL;
int currVar = -1;
bool inTypeVariables = false;
int maxVarState = -1;


int getVarIndex(const char* colname)
{
    int ret = 0;
    int pow = 1;

    for (int i = strlen(colname) - 1; isdigit(colname[i]); i--) {
        ret += pow * (colname[i] - '0');
        pow *= 10;
    }

    return ret - 1;
}

const XML_Char* getAtt(const XML_Char* attName, const XML_Char** attList)
{
    for (int i = 0; attList[i]; i += 2) {
        if (!strcmp(attList[i], attName)) {
            return attList[i+1];
        }
    }
    return NULL;
}

void characterDataHandler(void* userData, const XML_Char *cdata, int len)
{
    elementCdata.append(cdata, len);
}

void startElementHandler(void* userData, const XML_Char *name,
        const XML_Char **atts)
{
    // DEBUG
    //cout << "Found start element: " << name << endl;

    BayesNet* model = (BayesNet*)userData;
    if (!strcmp(name, "TypeVariables")) {
        inTypeVariables = true;
    } else if (!strcmp(name, "Variable")) {

        // HACK: this only works if every single variable
        // has a typeVariable with a name of the form:
        // "TypeK n-state" OR "TypeK continuous"
        if (!inTypeVariables) {
            const char* typeVarName = getAtt("typeVariable", atts);
            int numStates = -1;

            // DEBUG / HACK -- Really ought to handle general case,
            // when typeVariable might not be specified.
            if (typeVarName != NULL) {

                if (strstr(typeVarName, "ontinuous")) {
                    model->addVar(-1);
                } else {
                    // We expect the number of states to appear after
                    // a space, in a name of the form: "TypeK n-state"
                    for (unsigned int i = 0; i < strlen(typeVarName); i++) {
                        if (isspace(typeVarName[i])) {
                            numStates = atoi(typeVarName + i + 1);
                            break;
                        }
                    }

                    model->addVar(numStates);
                }
            } else {
                const char* type = getAtt("type", atts);
                if (type == NULL) {
                    cout << "ERROR: no type variable or type specified!\n";
                    return;
                }

                // Add continuous vars now.  Add categorical vars once we've
                // counted their states.
                if (!strcmp(type, "continuous")) {
                    model->addVar(-1);
                }
            }
        }
    } else if (!strcmp(name, "State")) {
        if (!inTypeVariables) {
            maxVarState++;
        }
    } else if (!strcmp(name, "LocalModel")) {
        const char* varName = getAtt("variable", atts);
        if (!varName) {
            cout << "ERROR: no variable specified in LocalModel!\n";
        }
        currVar = getVarIndex(varName);
        // DEBUG
        //cout << "Reading model for variable " << currVar << endl;
    } else if (!strcmp(name, "Input")) {
        const char* varName = getAtt("variable", atts);
        if (!varName) {
            cout << "ERROR: no variable specified in Input!\n";
        }
        model->addChild(getVarIndex(varName), currVar);
    } else if (!strcmp(name, "DecisionTree")) {
        // dtree = new DecisionTree(currVar, model->schema.getRange(currVar));
        dtree = model->getDecisionTree(currVar);
    } else if (!strcmp(name, "Vertex")) {
        const char* splitName = getAtt("split", atts);
        if (!splitName) {
            cout << "ERROR: no split variable specified at Vertex!\n";
        }
        dtree->beginVertex(getVarIndex(splitName));
    } else if (!strcmp(name, "Branch")) {
        dtree->beginBranch();
    } else if (!strcmp(name, "Multinomial")) {
        dtree->beginMultinomial();
    }

    elementCdata.clear();
}

void endElementHandler(void* userData, const XML_Char *name)
{
    // DEBUG
    //cout << "Found end element: " << name << endl;

    BayesNet* model = (BayesNet*)userData;
    if (!strcmp(name, "TypeVariables")) {
        inTypeVariables = false;
    } else if (!strcmp(name, "Variable")) {
        if (!inTypeVariables) {
            if (maxVarState >= 0) {
                model->addVar(maxVarState + 1);
            }
            maxVarState = -1;
        }
    } else if (!strcmp(name, "DecisionTree")) {
        /*
        model->setDecisionTree(currVar, dtree);
        dtree = NULL;
        */
    } else if (!strcmp(name, "Vertex")) {
        dtree->endVertex();
    } else if (!strcmp(name, "Branch")) {
        dtree->endBranch();
    } else if (!strcmp(name, "Values")) {
        list<Range> values;
        istringstream in(elementCdata);
        while (in) {
            Range r;
            in >> r;
            if (in) {
                values.push_back(r);
            }
        }
        dtree->endValues(values);
    } else if (!strcmp(name, "Probs")) {
        unsigned int base = 0;
        int index = 0;
        double *probs = new double[dtree->getRange()];

        // DEBUG
        //cout << "Probs: " << elementCdata << " (";

        for (unsigned int i = 0; i < elementCdata.length(); i++) {
            if (isspace(elementCdata[i])) {
                if (i > base) {
                    probs[index++] 
                        = atof(elementCdata.substr(base, i-base).c_str());
                    // DEBUG
                    //cout << probs[index-1] << " ";
                }
                base = i + 1;
            }
        }
        if (elementCdata.length() > base) {
            probs[index++] = atof(elementCdata.substr(base).c_str());
            // DEBUG
            //cout << probs[index-1] << " ";
        }
        // DEBUG
        //cout << ")\n";

        dtree->endProbs(probs);
    } else if (!strcmp(name, "Mean")) {
        dtree->endMean(atof(elementCdata.c_str()));
    } else if (!strcmp(name, "SD")) {
        dtree->endSD(atof(elementCdata.c_str()));
    } else if (!strcmp(name, "ProbMissing")) {
        dtree->endProbMissing(atof(elementCdata.c_str()));
    } else if (!strcmp(name, "BinGaussian")) {
        dtree->endBinGaussian();
    }

    elementCdata.clear();
}

// DEBUG
#include <stdio.h>

void loadModel(FILE* in, BayesNet& model)
{
    XML_Parser parser = XML_ParserCreate(NULL);
    if (!parser) {
        cout << "ERROR: could not allocate memory for parser!\n";
        return;
    }
    XML_SetUserData(parser, &model);
    XML_SetElementHandler(parser, startElementHandler, endElementHandler);
    XML_SetCharacterDataHandler(parser, characterDataHandler);
    char buf[10240];

    while (!feof(in)) {
        int bytesRead = fread(buf, 1, 1000, in);
        /* DEBUG
        ((char*)buf)[bytesRead] = '\0';
        printf((char*)buf);
        cout << "in.eof()? " << feof(in) << endl;
        */

        if ( !XML_Parse(parser, buf, bytesRead, feof(in)) ) {
            cout << "Parse error at line " << XML_GetCurrentLineNumber(parser)
                << ":\n" << XML_ErrorString(XML_GetErrorCode(parser)) << endl;
            return;
        }
    }

    XML_ParserFree(parser);

// We have no reason at present to get the splits of a continuous variable.
// This was debugging code for functionality that's no longer used.
#if 0
    for (int i = 0; i < model.getNumVars(); i++) {
        if (model.getRange(i) >= 0) {
            continue;
        }

        // Find all splits
        list<double> allsplits;
        for (int v = 0; v < model.getNumVars(); v++) {
            DecisionTree* currTree = model.getDecisionTree(v);
            list<double> currsplits = currTree->getSplits(i);
            allsplits.splice(allsplits.end(), currsplits);
        }

        allsplits.sort();

        // Print out values
        cout << i << ":";
        list<double>::iterator iter;
        for (iter = allsplits.begin(); iter != allsplits.end(); iter++) {
            cout << " " << *iter;
        }
        cout << endl;
    }
#endif

    return;
}

vector<double> BayesNet::MBdist(int var, VarSet& allVars) const
{
    const Leaf* l = decisionTrees[var]->getLeaf(allVars.getArray());
    int numVals = decisionTrees[var]->getRange();
    vector<double> distrib(numVals);
    double totalProb = -1.0/0.0000000000001;

    // Compute probability distribution across all values, 
    // using Markov-blanket inference.
    double origVal = allVars[var];
    for (int i = 0; i < numVals; i++) {
        // Compute p(x | parents(x))
        distrib[i] = l->getLogProb(i);

        // DEBUG
        if (isnan(distrib[i])) {
            cout << "distrib[" << i << "] = 0.0!\n";
        }

        // Scale by p(children(x) | parents(children(x)))
        list<int>::const_iterator c;
        for (c = children[var].begin(); c != children[var].end(); c++) {
            allVars[var] = i;
            distrib[i] += getLogProb(*c, allVars.getArray());
            // DEBUG
            if (isnan(distrib[i])) {
                cout << "distrib[" << i << "] = 0.0 after child " << *c << "\n";
            }
        }

        // Keep track of total probability, so we can normalize
        if (totalProb - distrib[i] < -10) {
            // If distrib[i] is much larger, just use it.
            totalProb = distrib[i];
        } else if (totalProb - distrib[i] < 10) {
            // If they're within the same range, add them
            totalProb = log(exp(totalProb - distrib[i]) + 1) + distrib[i];
        }
    }
    allVars[var] = origVal;

    if (isnan(totalProb)) {
        cout << "totalProb = 0!\n";
    }

    // Normalize
    for (int i = 0; i < numVals; i++) {
        distrib[i] = exp(distrib[i] - totalProb);
    }

    return distrib;
}

double BayesNet::MBsample(int var, VarSet& allVars) const
{
    const Leaf* l = decisionTrees[var]->getLeaf(allVars.getArray());
    int numVals;
    vector<double> actualVals;

    if (decisionTrees[var]->getRange() < 0) {
        DT::BinGaussian* g = (DT::BinGaussian*)l;
        double minVal = g->mean - 4 * g->SD;
        double maxVal = g->mean + 4 * g->SD;
        int numIncs = 1000;
        double dx = (maxVal - minVal)/numIncs;

        numVals = numIncs;

        for (double currX = minVal; currX < maxVal; currX += dx) {
            actualVals.push_back(currX);
        }
        actualVals.push_back(VarSet::UNKNOWN);
    } else {
        numVals = decisionTrees[var]->getRange();
        for (int i = 0; i < numVals; i++) {
            actualVals.push_back(i);
        }
    }

    vector<double> distrib(numVals);
    //double totalProb = 0.0;
    double totalProb = -1.0/0.0000000000001;

    // Compute probability distribution across all values, 
    // using Markov-blanket inference.
    double origVal = allVars[var];
    for (int i = 0; i < numVals; i++) {
        allVars[var] = actualVals[i];

        // Compute p(x | parents(x))
        distrib[i] = l->getLogProb(actualVals[i]);

        // Scale by p(children(x) | parents(children(x)))
        list<int>::const_iterator c;
        for (c = children[var].begin(); c != children[var].end(); c++) {
            distrib[i] += getLogProb(*c, allVars.getArray());
        }

        // Keep track of total probability, so we can normalize
        if (totalProb - distrib[i] < -10) {
            // If distrib[i] is much larger, just use it.
            totalProb = distrib[i];
        } else if (totalProb - distrib[i] < 10) {
            // If they're within the same range, add them
            totalProb = log(exp(totalProb - distrib[i]) + 1) + distrib[i];
        }
    }
    allVars[var] = origVal;

    // Normalize
    for (int i = 0; i < numVals; i++) {
        distrib[i] = exp(distrib[i] - totalProb);
    }

    // Sample
    double p = (double)rand() / RAND_MAX;
    for (int i = 0; i < numVals; i++) {
        p -= distrib[i];
        if (p < 0) {
            return actualVals[i];
        }
    }

    // DEBUG
    cout << "Error: defaulting to last value!\n";
    cout << "p = " << p << endl;
    cout << "var = " << var << endl;
    cout << "numVals = " << numVals <<  " (" 
        << decisionTrees[var]->getRange() << ")\n";
    for (int i = 0; i < numVals; i++) {
        cout << distrib[i] << " ";
    }
    cout << endl;

    return actualVals.back();
}

// Generate a complete sample
void BayesNet::wholeSample(VarSet& allvars) const
{
    /* HACK --
     * If you need a variable to be sampled, it must be
     * marked as UNTESTED in allvars before calling this function.
     * UNKNOWN is simply not sufficient, because UNKNOWN is a legal
     * value for continuous variables.  (For discrete variables,
     * the unknown value is a separate value taken care of when building
     * the dataset.  It must be handled specially for continuous variables.)
     */

    int knownVars = 0;
    for (int i = 0; i < numVars; i++) {
        /*
        if (allvars.isObserved(i) && allvars.isTested(i)) {
            knownVars++;
        }
        */
        if (allvars.isTested(i)) {
            knownVars++;
        }
    }

    // Sample each unknown variable whose parents have been specified.
    // Keep going until all variables have been sampled.
    while (knownVars < numVars) {
        for (int i = 0; i < numVars; i++) {
            if (/*!allvars.isObserved(i) || */ !allvars.isTested(i)) {
                bool unspecifiedParent = false;
                list<int>::const_iterator p;
                for (p = parents[i].begin(); p != parents[i].end(); p++) {
                    if (/* !allvars.isObserved(*p) || */ !allvars.isTested(*p)) {
                        unspecifiedParent = true;
                        break;
                    }
                }

                if (!unspecifiedParent) {
                    allvars[i] = decisionTrees[i]->sample(allvars.getArray());
                    knownVars++;
                }
            }
        }
    }
}