Package stallone.mc.correlations

Source Code of stallone.mc.correlations.DynamicalExpectationsSpectral

/*
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
package stallone.mc.correlations;

import stallone.api.algebra.Algebra;
import stallone.api.algebra.IEigenvalueDecomposition;
import stallone.api.doubles.Doubles;
import stallone.api.doubles.IDoubleArray;
import stallone.api.mc.IDynamicalExpectationsSpectral;
import stallone.api.mc.MarkovModel;
import stallone.mc.StationaryDistribution;


/**
* Calculates expectations and correlations of functions of states
*
* @author noe
*/
public class DynamicalExpectationsSpectral implements IDynamicalExpectationsSpectral
{
    // input
    private IDoubleArray T, K;
    private StationaryDistribution statdist = new StationaryDistribution();

    // stationary distribution
    private IDoubleArray pi;

    // spectral components
    private IDoubleArray timescales;
    private IDoubleArray R; // eigenvectors, normalized with sqrt()

    // result
    private IDoubleArray amplitudes;

    public DynamicalExpectationsSpectral(IDoubleArray M)
    {
        if (MarkovModel.util.isTransitionMatrix(M))
        {
            setT(M);
        }
        else if (MarkovModel.util.isTransitionMatrix(M))
        {
            setK(M);
        }
        else
        {
            throw(new IllegalArgumentException("Trying to construct DynamicalExpectationsSpectral with a Matrix that is neither a transition nor a rate matrix"));
        }
    }

    public DynamicalExpectationsSpectral()
    {
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#setT(stallone.api.doubles.IDoubleArray)
     */
    public final void setT(IDoubleArray _T)
    {
        this.K = null;
        this.T = _T;
        statdist.setT(_T);
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#setK(stallone.api.doubles.IDoubleArray)
     */
    public final void setK(IDoubleArray _K)
    {
        this.K = _K;
        this.T = null;
        statdist.setK(_K);
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#setStationaryDistribution(stallone.api.doubles.IDoubleArray)
     */
    public void setStationaryDistribution(IDoubleArray _pi)
    {
        pi = Doubles.create.array(_pi.getArray());
    }


    private void calculateAlgebra()
    {
        if (K==null && T==null)
            throw(new RuntimeException("Trying to calculate dynamical expectations before setting T or K."));

        // calculates the stationary distribution if not given
        if (pi == null)
        {
            pi = Doubles.create.array(statdist.calculate().getArray());
        }

        // eigenvalues and timescales
        if (T != null)
        {
            IEigenvalueDecomposition evd = Algebra.util.evd(T);
            IDoubleArray evalT = evd.getEval();
            timescales = Doubles.create.array(evalT.size());
            for (int i=0; i<timescales.size(); i++)
                timescales.set(i, -1/Math.log(Math.abs(evalT.get(i))));

            R = evd.getRightEigenvectorMatrix();
        }
        else
        {
            IEigenvalueDecomposition evd = Algebra.util.evd(K);
            IDoubleArray evalK = evd.getEval();
            timescales = Doubles.create.array(evalK.size());
            for (int i=0; i<timescales.size(); i++)
                timescales.set(i, -1/Math.abs(evalK.get(i)));

            R = evd.getRightEigenvectorMatrix();
        }

        // normalize eigenvectors
        for (int i=0; i<R.columns(); i++)
        {
            IDoubleArray ri = R.viewColumn(i);
            double s = Math.sqrt(Algebra.util.dot(ri, ri, pi));

            for (int j=0; j<R.size(); j++)
                R.set(j,i, R.get(j,i) / s);
        }
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#calculatePerturbationExpectation(stallone.api.doubles.IDoubleArray, stallone.api.doubles.IDoubleArray)
     */
    public void calculatePerturbationExpectation(IDoubleArray p0, IDoubleArray a)
    {
        if (pi == null || timescales == null || R == null)
            calculateAlgebra();

        amplitudes = Doubles.create.array(timescales.size());
        for (int i=0; i<amplitudes.size(); i++)
        {
            double a1 = Algebra.util.dot(p0, R.viewColumn(i));
            double a2 = Algebra.util.dot(a, R.viewColumn(i), pi);
            amplitudes.set(i, a1*a2);
        }
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#calculateAutocorrelation(stallone.api.doubles.IDoubleArray)
     */
    public void calculateAutocorrelation(IDoubleArray a)
    {
        if (pi == null || timescales == null || R == null)
            calculateAlgebra();

        amplitudes = Doubles.create.array(timescales.size());
        for (int i=0; i<amplitudes.size(); i++)
        {
            double a1 = Algebra.util.dot(a, R.viewColumn(i), pi);
            amplitudes.set(i, a1*a1);
        }
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#calculateCorrelation(stallone.api.doubles.IDoubleArray, stallone.api.doubles.IDoubleArray)
     */
    public void calculateCorrelation(IDoubleArray a, IDoubleArray b)
    {
        if (pi == null || timescales == null || R == null)
            calculateAlgebra();

        amplitudes = Doubles.create.array(timescales.size());
        for (int i=0; i<amplitudes.size(); i++)
        {
            double a1 = Algebra.util.dot(a, R.viewColumn(i), pi);
            double a2 = Algebra.util.dot(b, R.viewColumn(i), pi);
            amplitudes.set(i, a1*a2);
        }
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#getAmplitudes()
     */
    public IDoubleArray getAmplitudes()
    {
        return(amplitudes);
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#getTimescales()
     */
    public IDoubleArray getTimescales()
    {
        return(timescales);
    }

    /* (non-Javadoc)
     * @see stallone.mc.correlations.IDynamicalExpectationsSpectral#getValue(double)
     */
    public double getValue(double t)
    {
        double res = 0;
        for (int i=0; i<timescales.size(); i++)
        {
            res += amplitudes.get(i) * Math.exp(-t/timescales.get(i));
        }
        return(res);
    }

}
TOP

Related Classes of stallone.mc.correlations.DynamicalExpectationsSpectral

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.