/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/PhysicsTools/MVATrainer/src/LeastSquares.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <cstring>
#include <vector>
#include <cmath>

#include <TMatrixD.h>
#include <TVectorD.h>
#include <TDecompSVD.h>

#include "FWCore/Utilities/interface/Exception.h"

#include "PhysicsTools/MVATrainer/interface/XMLDocument.h"
#include "PhysicsTools/MVATrainer/interface/XMLSimpleStr.h"
#include "PhysicsTools/MVATrainer/interface/XMLUniStr.h"
#include "PhysicsTools/MVATrainer/interface/LeastSquares.h"

XERCES_CPP_NAMESPACE_USE

namespace PhysicsTools {

LeastSquares::LeastSquares(unsigned int n) :
        coeffs(n + 2), covar(n + 1), corr(n + 1), rotation(n, n),
        weights(n + 1), variance(n + 1), trace(n), n(n)
{
}

LeastSquares::~LeastSquares()
{
}

void LeastSquares::add(const std::vector<double> &values,
                       double dest, double weight)
{
        if (values.size() != n)
                throw cms::Exception("LeastSquares")
                        << "add(): invalid array size!" << std::endl;

        for(unsigned int i = 0; i < n; i++) {
                for(unsigned int j = 0; j < n; j++)
                        coeffs(i, j) += values[i] * values[j] * weight;

                coeffs(n, i) += values[i] * dest * weight;
                coeffs(i, n) += values[i] * dest * weight;
                coeffs(n + 1, i) += values[i] * weight;
                coeffs(i, n + 1) += values[i] * weight;
        } 

        coeffs(n, n) += dest * dest * weight;
        coeffs(n + 1, n) += dest * weight;
        coeffs(n, n + 1) += dest * weight;
        coeffs(n + 1, n + 1) += weight;
}

void LeastSquares::add(const LeastSquares &other, double weight)
{
        if (other.getSize() != n)
                throw cms::Exception("LeastSquares")
                        << "add(): invalid array size!" << std::endl;

        coeffs += weight * other.coeffs;
}

TVectorD LeastSquares::solveFisher(const TMatrixDSym &coeffs)
{
        unsigned int n = coeffs.GetNrows() - 2;

        TMatrixDSym tmp;
        coeffs.GetSub(0, n, tmp);
        tmp[n] = TVectorD(n + 1, coeffs[n + 1].GetPtr());
        tmp(n, n) = coeffs(n + 1, n + 1);

        TDecompSVD decCoeffs(tmp);
        bool ok;
        return decCoeffs.Solve(TVectorD(n + 1, coeffs[n].GetPtr()), ok);
}

TMatrixD LeastSquares::solveRotation(const TMatrixDSym &covar, TVectorD &trace)
{
        TMatrixDSym tmp;
        covar.GetSub(0, covar.GetNrows() - 2, tmp);
        TDecompSVD decCovar(tmp);
        trace = decCovar.GetSig();
        return decCovar.GetU();
}

void LeastSquares::calculate()
{
        double N = coeffs(n + 1, n + 1);

        for(unsigned int i = 0; i <= n; i++) {
                double M = coeffs(n + 1, i);
                for(unsigned int j = 0; j <= n; j++)
                        covar(i, j) = coeffs(i, j) * N - M * coeffs(n + 1, j);
        }

        for(unsigned int i = 0; i <= n; i++) {
                double c = covar(i, i);
                variance[i] = c > 0.0 ? std::sqrt(c) : 0.0;
        }

        for(unsigned int i = 0; i <= n; i++) {
                double M = variance[i];
                for(unsigned int j = 0; j <= n; j++) {
                        double v = M * variance[j];
                        double w = covar(i, j);

                        corr(i, j) = (v >= 1.0e-9) ? (w / v) : (i == j);
                }
        }

        weights = solveFisher(coeffs);
        rotation = solveRotation(covar, trace);
}

std::vector<double> LeastSquares::getWeights() const
{
        std::vector<double> results;
        results.reserve(n);
        
        for(unsigned int i = 0; i < n; i++)
                results.push_back(weights[i]);
        
        return results;
}

std::vector<double> LeastSquares::getMeans() const
{
        std::vector<double> results;
        results.reserve(n);

        double N = coeffs(n + 1, n + 1);
        for(unsigned int i = 0; i < n; i++)
                results.push_back(coeffs(n + 1, i) / N);

        return results;
}

double LeastSquares::getConstant() const
{
        return weights[n];
}

static void loadMatrix(DOMElement *elem, unsigned int n, TMatrixDBase &matrix)
{
        if (std::strcmp(XMLSimpleStr(elem->getNodeName()),
                                     "matrix") != 0)
                throw cms::Exception("LeastSquares")
                        << "Expected matrix in data file."
                        << std::endl;

        unsigned int row = 0;
        for(DOMNode *node = elem->getFirstChild();
            node; node = node->getNextSibling()) {
                if (node->getNodeType() != DOMNode::ELEMENT_NODE)
                        continue;

                if (std::strcmp(XMLSimpleStr(node->getNodeName()),
                                             "row") != 0)
                        throw cms::Exception("LeastSquares")
                                << "Expected row tag in data file."
                                << std::endl;

                if (row >= n)
                        throw cms::Exception("LeastSquares")
                                << "Too many rows in data file." << std::endl;

                elem = static_cast<DOMElement*>(node);

                unsigned int col = 0;
                for(DOMNode *subNode = elem->getFirstChild();
                    subNode; subNode = subNode->getNextSibling()) {
                        if (subNode->getNodeType() != DOMNode::ELEMENT_NODE)
                                continue;

                        if (std::strcmp(XMLSimpleStr(subNode->getNodeName()),
                                                     "value") != 0)
                        throw cms::Exception("LeastSquares")
                                << "Expected value tag in data file."
                                << std::endl;

                        if (col >= n)
                                throw cms::Exception("LeastSquares")
                                        << "Too many columns in data file."
                                        << std::endl;

                        matrix(row, col) =
                                XMLDocument::readContent<double>(subNode);
                        col++;
                }

                if (col != n)
                        throw cms::Exception("LeastSquares")
                                << "Missing columns in data file."
                                << std::endl;
                row++;
        }

        if (row != n)
                throw cms::Exception("LeastSquares")
                        << "Missing rows in data file."
                        << std::endl;
}

static void loadVector(DOMElement *elem, unsigned int n, TVectorD &vector)
{
        if (std::strcmp(XMLSimpleStr(elem->getNodeName()),
                                     "vector") != 0)
                throw cms::Exception("LeastSquares")
                        << "Expected matrix in data file."
                        << std::endl;

        unsigned int col = 0;
        for(DOMNode *node = elem->getFirstChild();
            node; node = node->getNextSibling()) {
                if (node->getNodeType() != DOMNode::ELEMENT_NODE)
                        continue;

                if (std::strcmp(XMLSimpleStr(node->getNodeName()),
                                             "value") != 0)
                throw cms::Exception("LeastSquares")
                        << "Expected value tag in data file."
                        << std::endl;

                if (col >= n)
                        throw cms::Exception("LeastSquares")
                                << "Too many columns in data file."
                                << std::endl;

                vector(col) = XMLDocument::readContent<double>(node);
                col++;
        }

        if (col != n)
                throw cms::Exception("LeastSquares")
                        << "Missing columns in data file."
                        << std::endl;
}

static DOMElement *saveMatrix(DOMDocument *doc, unsigned int n,
                              const TMatrixDBase &matrix)
{
        DOMElement *root = doc->createElement(XMLUniStr("matrix"));
        XMLDocument::writeAttribute<unsigned int>(root, "size", n);

        for(unsigned int i = 0; i < n; i++) {
                DOMElement *row = doc->createElement(XMLUniStr("row"));
                root->appendChild(row);

                for(unsigned int j = 0; j < n; j++) {
                        DOMElement *value =
                                doc->createElement(XMLUniStr("value"));
                        row->appendChild(value);

                        XMLDocument::writeContent<double>(value, doc,
                                                          matrix(i, j));
                }
        }

        return root;
}

static DOMElement *saveVector(DOMDocument *doc, unsigned int n,
                              const TVectorD &vector)
{
        DOMElement *root = doc->createElement(XMLUniStr("vector"));
        XMLDocument::writeAttribute<unsigned int>(root, "size", n);

        for(unsigned int i = 0; i < n; i++) {
                DOMElement *value =
                        doc->createElement(XMLUniStr("value"));
                root->appendChild(value);

                XMLDocument::writeContent<double>(value, doc, vector(i));
        }

        return root;
}

void LeastSquares::load(DOMElement *elem)
{
        if (std::strcmp(XMLSimpleStr(elem->getNodeName()),
                                     "LinearAnalysis") != 0)
                throw cms::Exception("LeastSquares")
                        << "Expected LinearAnalysis in data file."
                        << std::endl;

        unsigned int version = XMLDocument::readAttribute<unsigned int>(
                                                        elem, "version", 1);

        enum Position {
                POS_COEFFS, POS_COVAR, POS_CORR, POS_ROTATION,
                POS_SUMS, POS_WEIGHTS, POS_VARIANCE, POS_TRACE, POS_DONE
        } pos = POS_COEFFS;

        for(DOMNode *node = elem->getFirstChild();
            node; node = node->getNextSibling()) {
                if (node->getNodeType() != DOMNode::ELEMENT_NODE)
                        continue;

                DOMElement *elem = static_cast<DOMElement*>(node);

                switch(pos) {
                    case POS_COEFFS:
                        if (version < 2) {
                                loadMatrix(elem, n + 1, coeffs);
                                coeffs.ResizeTo(n + 2, n + 2);
                                for(unsigned int i = 0; i <= n; i++) {
                                        coeffs(n + 1, i) = coeffs(n, i);
                                        coeffs(i, n + 1) = coeffs(i, n);
                                }
                                coeffs(n + 1, n + 1) = coeffs(n + 1, n);
                        } else
                                loadMatrix(elem, n + 2, coeffs);
                        break;
                    case POS_COVAR:
                        if (version < 2)
                                loadMatrix(elem, n, covar);
                        else
                                loadMatrix(elem, n + 1, covar);
                        break;
                    case POS_CORR:
                        if (version < 2)
                                loadMatrix(elem, n, corr);
                        else
                                loadMatrix(elem, n + 1, corr);
                        break;
                    case POS_ROTATION:
                        loadMatrix(elem, n, rotation);
                        break;
                    case POS_SUMS:
                        if (version < 2) {
                                TVectorD tmp(n + 1);
                                loadVector(elem, n + 1, tmp);

                                double M = coeffs(n + 1, n);
                                double N = coeffs(n + 1, n + 1);

                                for(unsigned int i = 0; i <= n; i++) {
                                        double v = coeffs(n, i) * N -
                                                   M * coeffs(n + 1, i);
                                        double w = coeffs(n, i) * N - v;

                                        covar(n, i) = w;
                                        covar(i, n) = w;
                                }

                                break;
                        } else
                                pos = (Position)(pos + 1);
                    case POS_WEIGHTS:
                        loadVector(elem, n + 1, weights);
                        break;
                    case POS_VARIANCE:
                        if (version < 2) {
                                loadVector(elem, n, variance);

                                double M = covar(n, n);
                                M = M > 0.0 ? std::sqrt(M) : 0.0;
                                variance[n] = M;

                                for(unsigned int i = 0; i <= n; i++) {
                                        double v = M * variance[i];
                                        double w = covar(n, i);
                                        double c = (v >= 1.0e-9)
                                                        ? (w / v) : (i == n);

                                        corr(n, i) = c;
                                        corr(i, n) = c;
                                }
                        } else
                                loadVector(elem, n + 1, variance);
                        break;
                    case POS_TRACE:
                        loadVector(elem, n, trace);
                        break;
                    default:
                        throw cms::Exception("LeastSquares")
                                << "Superfluous content in data file."
                                << std::endl;
                }

                pos = (Position)(pos + 1);
        }

        if (pos != POS_DONE)
                throw cms::Exception("LeastSquares")
                        << "Missing objects in data file."
                        << std::endl;
}

DOMElement *LeastSquares::save(DOMDocument *doc) const
{
        DOMElement *root = doc->createElement(XMLUniStr("LinearAnalysis"));
        XMLDocument::writeAttribute<unsigned int>(root, "version", 2);
        XMLDocument::writeAttribute<unsigned int>(root, "size", n);

        root->appendChild(saveMatrix(doc, n + 2, coeffs));
        root->appendChild(saveMatrix(doc, n + 1, covar));
        root->appendChild(saveMatrix(doc, n + 1, corr));
        root->appendChild(saveMatrix(doc, n, rotation));
        root->appendChild(saveVector(doc, n + 1, weights));
        root->appendChild(saveVector(doc, n + 1, variance));
        root->appendChild(saveVector(doc, n, trace));

        return root;
}

} // namespace PhysicsTools