edu.stanford.cfuller.imageanalysistools.clustering.GaussianLikelihoodObjectiveFunction.java Source code

Java tutorial

Introduction

Here is the source code for edu.stanford.cfuller.imageanalysistools.clustering.GaussianLikelihoodObjectiveFunction.java

Source

/* ***** BEGIN LICENSE BLOCK *****
 * 
 * Copyright (c) 2011 Colin J. Fuller
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 * 
 * ***** END LICENSE BLOCK ***** */

package edu.stanford.cfuller.imageanalysistools.clustering;

import edu.stanford.cfuller.imageanalysistools.fitting.ObjectiveFunction;

import org.apache.commons.math3.analysis.function.Exp;
import org.apache.commons.math3.linear.*;

/**
 * A function that calculates the likelihood of observing a set of points given a set of Gaussian sources generating them.
 *
 * This is used for clustering to determine the likelihood of having observed the ClusterObjects at their locations given the
 * set of Clusters currently under consideration.
 *
 * @author Colin J. Fuller
 *
 */
public class GaussianLikelihoodObjectiveFunction implements ObjectiveFunction {

    final double negLog2PI = -1.0 * Math.log(Math.PI * 2.0);
    final int nParametersEach = 5;
    final int numDim = 2;
    final double detCutoff = 1.0e-30;
    RealVector mean;
    RealVector x;
    RealVector pk;
    RealMatrix clusterProbs;
    java.util.Vector<RealMatrix> abdMatrices;
    RealVector det;
    java.util.Vector<ClusterObject> objects;

    /**
     * Constructs a new GaussianLikelihoodObjectiveFunction that will determine the likelihood of having observed the
     * supplied ClusterObjects at their locations.
     *
     * @param objects   A Vector containing the observed ClusterObjects (with locations already determined and assigned).
     */
    public GaussianLikelihoodObjectiveFunction(java.util.Vector<ClusterObject> objects) {
        mean = new ArrayRealVector(numDim);
        x = new ArrayRealVector(numDim);
        pk = new ArrayRealVector();
        clusterProbs = new Array2DRowRealMatrix();
        abdMatrices = new java.util.Vector<RealMatrix>();
        det = new ArrayRealVector();
        this.objects = objects;

    }

    /**
     * Evaluates the function with the specified parameters.
     *
     * The parameters describe a set of gaussian generators (which are the Clusters).
     *
     * @param parameters    A RealVector containing the values of all the parameters of each Gaussian, ordered so that all the parameters of a single gaussian are together, then the next gaussian, etc.
     * @return              The negative log likelihood of having observed the ClusterObjects at their locations, given the parameters describing the Gaussian clusters.
     */
    public double evaluate(RealVector parameters) {

        int nClusters = parameters.getDimension() / nParametersEach;

        //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("nClusters: " + nClusters + "  abdMatrices_size: " + abdMatrices.size() + "  det_dim: " + det.getDimension());

        if (det.getDimension() != nClusters) {

            clusterProbs = new Array2DRowRealMatrix(this.objects.size(), nClusters);
            det = new ArrayRealVector(nClusters);
            pk = new ArrayRealVector(nClusters);

            if (abdMatrices.size() < nClusters) {
                int originalSize = abdMatrices.size();
                for (int i = 0; i < nClusters - originalSize; i++) {
                    abdMatrices.add(new Array2DRowRealMatrix(numDim, numDim));
                }
            } else {
                abdMatrices.setSize(nClusters);
            }

        }

        pk.mapMultiplyToSelf(0.0);

        //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("nClusters: " + nClusters + "  abdMatrices_size: " + abdMatrices.size() + "  det_dim: " + det.getDimension());

        for (int i = 0; i < nClusters; i++) {
            /*
            double ct = Math.cos(parameters.getEntry(nParametersEach*i+3));
            double st = Math.sin(parameters.getEntry(nParametersEach*i+3));
            double sin2t = Math.sin(2*parameters.getEntry(nParametersEach*i+3));
            double a = (ct*ct/(2*parameters.getEntry(nParametersEach*i+2)) + st*st/(2*parameters.getEntry(nParametersEach*i+4)));
            double b = (sin2t/(4*parameters.getEntry(nParametersEach*i+4)) - sin2t/(4*parameters.getEntry(nParametersEach*i+2)));
            double d = (st*st/(2*parameters.getEntry(nParametersEach*i+2)) + ct*ct/(2*parameters.getEntry(nParametersEach*i+4)));
            */

            double a = parameters.getEntry(nParametersEach * i + 2);
            double d = parameters.getEntry(nParametersEach * i + 4);
            double b = Math.sqrt(a * d) * parameters.getEntry(nParametersEach * i + 3);

            abdMatrices.get(i).setEntry(0, 0, a);
            abdMatrices.get(i).setEntry(1, 0, b);
            abdMatrices.get(i).setEntry(0, 1, b);
            abdMatrices.get(i).setEntry(1, 1, d);

            LUDecomposition abdLU = (new LUDecomposition(abdMatrices.get(i)));

            det.setEntry(i, (abdLU).getDeterminant());
            //det.setEntry(i, a*d-b*b);
            try {
                abdMatrices.set(i, abdLU.getSolver().getInverse());
            } catch (org.apache.commons.math3.linear.SingularMatrixException e) {
                return Double.MAX_VALUE;
            }

        }

        for (int n = 0; n < this.objects.size(); n++) {

            ClusterObject c = this.objects.get(n);

            double max = -1.0 * Double.MAX_VALUE;
            int maxIndex = 0;

            for (int k = 0; k < nClusters; k++) {

                mean.setEntry(0, c.getCentroid().getX() - parameters.getEntry(nParametersEach * k));
                mean.setEntry(1, c.getCentroid().getY() - parameters.getEntry(nParametersEach * k + 1));

                x = abdMatrices.get(k).operate(mean);

                double dot = x.dotProduct(mean);

                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("k, n: " + k + ", " + this.objects.size());
                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("parameters: " + parameters.toString());
                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("abd matrix: " + abdMatrices.get(k).toString());
                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("det: " + det.getEntry(k));
                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("mean: " + mean.toString());
                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("dot: " + dot);

                double logN = negLog2PI - 0.5 * Math.log(det.getEntry(k)) - 0.5 * dot;

                //                java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("logN: " + logN);

                clusterProbs.setEntry(n, k, logN);
                if (logN > max) {
                    max = logN;
                    maxIndex = k;
                }

                if (Double.isInfinite(logN) || Double.isNaN(logN)) {
                    return Double.MAX_VALUE;
                }

            }

            c.setMostProbableCluster(maxIndex);

        }

        for (int k = 0; k < nClusters; k++) {

            double tempMax = -1.0 * Double.MAX_VALUE;

            for (int n = 0; n < this.objects.size(); n++) {
                if (clusterProbs.getEntry(n, k) > tempMax)
                    tempMax = clusterProbs.getEntry(n, k);
            }

            pk.setEntry(k,
                    tempMax + Math.log(
                            clusterProbs.getColumnVector(k).mapSubtract(tempMax).mapToSelf(new Exp()).getL1Norm())
                            - Math.log(this.objects.size()));

        }

        double pkMax = -1.0 * Double.MAX_VALUE;

        for (int k = 0; k < nClusters; k++) {
            if (pk.getEntry(k) > pkMax)
                pkMax = pk.getEntry(k);
        }

        double logSumPk = pkMax + Math.log(pk.mapSubtract(pkMax).mapToSelf(new Exp()).getL1Norm());

        pk.mapSubtractToSelf(logSumPk);

        //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("pk: " + pk.toString());

        double L = 0;

        for (int n = 0; n < this.objects.size(); n++) {

            RealVector toSum = clusterProbs.getRowVector(n).add(pk);

            double tempMax = -1.0 * Double.MAX_VALUE;

            for (int k = 0; k < nClusters; k++) {
                if (toSum.getEntry(k) > tempMax)
                    tempMax = toSum.getEntry(k);
            }

            double pn = tempMax + Math.log(toSum.mapSubtract(tempMax).mapToSelf(new Exp()).getL1Norm());

            //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("pn: " + pn);

            L += pn;

        }

        return -1.0 * L;
    }

}