Package stallone.stat

Source Code of stallone.stat.EventBinningCorrelator$Window

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

import java.util.Arrays;
import java.util.List;
import static stallone.api.API.*;
import stallone.api.doubles.IDoubleArray;
import stallone.api.doubles.IDoubleList;
import stallone.util.CommandLineParser;

/**
*
* @author noe
*/
public class EventBinningCorrelator
{

    private IDoubleArray lagtimes;
    private IDoubleArray averageWidths;
    private double samplingTime = 1;
    private IDoubleArray correlation;
    private IDoubleArray weights;
    private double sum = 0;
    private double sumOfWeights = 0;

    /**
     *
     * @param _lagtimes
     * @param _averageWidths
     * @param samplingTime the time step at which t is advanced when calculating
     * the average < x * y >_t
     */
    public EventBinningCorrelator(IDoubleArray _lagtimes, IDoubleArray _averageWidths, double _samplingTime)
    {
        if (_lagtimes.size() != _averageWidths.size())
        {
            throw (new IllegalArgumentException("Number of lagtimes does not match number of average widths in correlate."));
        }

        lagtimes = _lagtimes;
        averageWidths = _averageWidths;
        samplingTime = _samplingTime;
        correlation = doublesNew.array(lagtimes.size());
        weights = doublesNew.array(lagtimes.size());
    }

    public static String getUsageString()
    {
        String res = "EventCorrelator" + "\n"
                + "Correlates the given columns of all files passed."
                + "-i <trajectorie(s)>" + "\n"
                + "\n"
                + "-columns <time> <D> <A>" + "\n"
                + "-samplingtime <dt>" + "\n"
                + "-maxtime <maxy>" + "\n"
                + "-latmult <mult>" + "\n"
                + "-windowfraction <fw>" + "\n"
                + "-subtractmean" + "\n";

        return res;
    }

    public static CommandLineParser parseArguments(String[] args)
    {
        CommandLineParser parser = new CommandLineParser();

        // input
        parser.addStringArrayCommand("i", true);

        parser.addCommand("columns", true);
        parser.addIntArgument("columns", true);
        parser.addIntArgument("columns", true);
        parser.addIntArgument("columns", true);

        parser.addIntCommand("maxtime", true);
        parser.addDoubleCommand("lagmult", true);

        parser.addIntCommand("samplingtime", true);

        parser.addDoubleCommand("windowfraction", true);

        parser.addCommand("subtractmean", true);

        // parse
        if (!parser.parse(args))
        {
            throw (new IllegalArgumentException("Parsing error!"));
        }

        return parser;
    }

    private static IDoubleArray cumulate(IDoubleArray data)
    {
        IDoubleArray res = doublesNew.array(data.size() + 1);
        res.set(0, 0);
        for (int i = 1; i < res.size(); i++)
        {
            res.set(i, res.get(i - 1) + data.get(i - 1));
        }
        return res;
    }

    class Window
    {
        public IDoubleArray times;
        public IDoubleArray data;
        public double width;
        public double t, tl, tr;
        public int l=0, r=0; // first event index included and first event index after window [t-width, t+width]

        public Window(IDoubleArray _times, IDoubleArray _data, double _width)
        {
            this.times = _times;
            this.data = _data;
            this.width = _width;
        }

        /**
         * Positions the window to time t.
         */
        public void init(double _t)
        {
            this.t = _t;
            this.tl = t-width/2.0;
            this.tr = t+width/2.0;

            l=0;
            r=0;
            advance(0);
        }

        public void advance(double dt)
        {
            this.t += dt;
            this.tl = t-width/2.0;
            this.tr = t+width/2.0;

            // move left index
            while (l<times.size())
            {
                if (times.get(l) < tl)
                    l++;
                else
                    break;
            }

            // move right index
            while (r<times.size())
            {
                if (times.get(r) < tr)
                    r++;
                else
                    break;
            }
        }

        public boolean isEmpty()
        {
            return (l==r);
        }

        public boolean endOfData()
        {
            return (l>=times.size());
        }
    }


    public double correlate(IDoubleArray time, IDoubleArray data, double tau, double windowSize)
    {
        if (time.size() != data.size())
        {
            throw (new IllegalArgumentException("Number of time points does not match number of data points in correlate."));
        }

        // cumulate data
        IDoubleArray cumdata = cumulate(data);

        Window w1 = new Window(time, data, windowSize);
        w1.init(windowSize/2.0);

        Window w2 = new Window(time, data, windowSize);
        w2.init(windowSize/2.0 + tau);

        // estimate
        double sum = 0;
        double count = 0;

        while(!w2.endOfData())
        {
            /*
                if (tau == 0)
                {
                    System.out.println("t = "+w1.t+"\t times ["+w1.tl+","+w1.tr+"] -> ["+w2.tl+","+w2.tr+"]");
                    System.out.println("  w1=("+w1.l+", "+w1.r+") w2="+w2.l+", "+w2.r);
                }*/

            // update convolution sums
            if ((!w1.isEmpty()) && (!w2.isEmpty()))
            {
                double x1 = (cumdata.get(w1.r) - cumdata.get(w1.l))/((double)(w1.r-w1.l));
                double x2 = (cumdata.get(w2.r) - cumdata.get(w2.l))/((double)(w2.r-w2.l));
                sum += x1*x2;
                count += 1.0;

                /*if (tau == 0)
                {
                    System.out.println("t = "+w1.t+"\t times ["+w1.tl+","+w1.tr+"] -> ["+w2.tl+","+w2.tr+"]");
                    System.out.println("  w1=("+w1.l+", "+w1.r+") w2="+w2.l+", "+w2.r);
                    System.out.println("  x1= "+x1+"  x2= "+x2);
                }*/
            }

            // advance window by dt.
            w1.advance(samplingTime);
            w2.advance(samplingTime);
        }

        return sum / count;
    }

    /**
     * Correlates the given set of events on the tau-grid specified
     *
     * @param time
     * @param data
     * @return the autocorrelations of the data on the tau grid specified
     */
    public IDoubleArray correlate(IDoubleArray time, IDoubleArray data)
    {
        IDoubleArray res = doublesNew.array(lagtimes.size());

        for (int i = 0; i < res.size(); i++)
        {
            res.set(i, correlate(time, data, lagtimes.get(i), averageWidths.get(i)));
        }

        return res;
    }

    /**
     * Adds data set to a cumulative correlation
     *
     * @param time
     * @param data
     */
    public void add(IDoubleArray time, IDoubleArray data)
    {
        IDoubleArray c = correlate(time, data);
        for (int i = 0; i < c.size(); i++)
        {
            if (!Double.isNaN(c.get(i)))
            {
                correlation.set(i, correlation.get(i) + data.size() * c.get(i));
                weights.set(i, weights.get(i) + data.size());
            }
        }

        sum += doubles.sum(data);
        sumOfWeights += time.size();
    }

    public IDoubleArray getCorrelation()
    {
        IDoubleArray res = doublesNew.array(correlation.size());
        for (int i = 0; i < res.size(); i++)
        {
            res.set(i, correlation.get(i) / weights.get(i));
        }
        return res;
    }

    public IDoubleArray getCorrelationMeanFree()
    {
        IDoubleArray res = doublesNew.array(correlation.size());
        double shift = (sum / sumOfWeights) * (sum / sumOfWeights);
        System.out.println("# square mean = "+shift);
        for (int i = 0; i < res.size(); i++)
        {
            res.set(i, (correlation.get(i) / weights.get(i)) - shift);
        }
        return res;
    }

    public void reset()
    {
        correlation.zero();
        weights.zero();
        sum = 0;
        sumOfWeights = 0;
    }

    public static void main(String[] args)
    {
        if (args.length == 0)
        {
            System.out.println(getUsageString());
            System.exit(0);
        }

        CommandLineParser parser = parseArguments(args);

        // load photon trajectories
        List<String> inputFiles = Arrays.asList(parser.getStringArray("i"));

        int timecol = parser.getInt("columns", 0);
        int col1 = parser.getInt("columns", 1);
        int col2 = parser.getInt("columns", 2);
        double lagtmult = parser.getDouble("lagmult",0);
        int maxtime = parser.getInt("maxtime", 0);
        int samplingtime = parser.getInt("samplingtime",0);
        //double windowfraction = parser.getDouble("windowfraction",0);
        boolean subtractMean = parser.hasCommand("subtractmean");


        IDoubleList lagtimes = doublesNew.listFrom(0, 1);
        while (lagtimes.get(lagtimes.size() - 1) < maxtime)
        {
            lagtimes.append((int) ((lagtimes.get(lagtimes.size() - 1) + 1) * lagtmult));
        }
        //IDoubleArray lagtimes = doublesNew.arrayFrom(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000, 100000, 200000, 500000, 1000000);
        IDoubleArray averageWidths = doublesNew.array(lagtimes.size());
        for (int i = 0; i < averageWidths.size(); i++)
        {
            averageWidths.set(i, 1.0);
            //averageWidths.set(i, Math.max(1.0, windowfraction * lagtimes.get(i)));
            //System.out.println(lagtimes.get(i)+"\t"+averageWidths.get(i));
        }
        //System.out.println();
        //System.exit(0);
        EventBinningCorrelator correlator = new EventBinningCorrelator(lagtimes, averageWidths, samplingtime);


        for (int i = 0; i < inputFiles.size(); i++)
        {
            IDoubleArray time = data.readColumn(inputFiles.get(i), timecol);
            IDoubleArray don = data.readColumn(inputFiles.get(i), col1);
            IDoubleArray acc = data.readColumn(inputFiles.get(i), col2);
            IDoubleArray E = doublesNew.array(don.size());
            for (int j = 0; j < E.size(); j++)
            {
                E.set(j, acc.get(j) / (don.get(j) + acc.get(j)));
            }

            correlator.add(time, E);
        }

        IDoubleArray corr = null;
        if (subtractMean)
        {
            corr = correlator.getCorrelationMeanFree();
        }
        else
        {
            corr = correlator.getCorrelation();
        }


        //IDoubleArray corr = correlator.correlate(time, E);

        for (int i = 0; i < lagtimes.size(); i++)
        {
            System.out.println(lagtimes.get(i) + "\t" + corr.get(i));
        }
    }
}
TOP

Related Classes of stallone.stat.EventBinningCorrelator$Window

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.