/*
 * SubstitutionModelDelegate.java
 *
 * Copyright (c) 2002-2016 Alexei Drummond, Andrew Rambaut and Marc Suchard
 *
 * This file is part of BEAST.
 * See the NOTICE file distributed with this work for additional
 * information regarding copyright ownership and licensing.
 *
 * BEAST is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as
 * published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 *  BEAST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with BEAST; if not, write to the
 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
 * Boston, MA  02110-1301  USA
 */

package dr.evomodel.treedatalikelihood;

import beagle.Beagle;
import dr.evolution.tree.Tree;
import dr.evomodel.branchmodel.BranchModel;
import dr.evomodel.substmodel.EigenDecomposition;
import dr.evomodel.substmodel.SubstitutionModel;
import dr.util.Timer;

import java.io.Serializable;
import java.util.ArrayDeque;
import java.util.ArrayList;
import java.util.Deque;
import java.util.List;

/**
 * A simple substitution model delegate with the same substitution model over the whole tree
 * @author Andrew Rambaut
 * @version $Id$
 */
public final class HomogenousSubstitutionModelDelegate implements EvolutionaryProcessDelegate, Serializable {

    private final SubstitutionModel substitutionModel;

    private final int eigenCount = 1;
    private final int nodeCount;

    private final BufferIndexHelper eigenBufferHelper;
    private final BufferIndexHelper matrixBufferHelper;

    /**
     * A class which handles substitution models including epoch models where multiple
     * substitution models on a branch are convolved.
     * @param tree
     * @param branchModel Describes which substitution models use on each branch
     */
    public HomogenousSubstitutionModelDelegate(Tree tree, BranchModel branchModel) {
        this(tree, branchModel, 0);
    }

    /**
     * A class which handles substitution models including epoch models where multiple
     * substitution models on a branch are convolved.
     * @param tree
     * @param branchModel Describes which substitution models use on each branch
     * @param partitionNumber which data partition is this (used to offset eigen and matrix buffer numbers)
     */
    public HomogenousSubstitutionModelDelegate(Tree tree, BranchModel branchModel, int partitionNumber) {

        assert(branchModel.getSubstitutionModels().size() == 1) : "this delegate should only be used with simple branch models";

        this.substitutionModel = branchModel.getRootSubstitutionModel();


        nodeCount = tree.getNodeCount();

        // two eigen buffers for each decomposition for store and restore.
        eigenBufferHelper = new BufferIndexHelper(eigenCount, 0, partitionNumber);

        // two matrices for each node less the root
        matrixBufferHelper = new BufferIndexHelper(nodeCount, 0, partitionNumber);

    }// END: Constructor

    @Override
    public boolean canReturnComplexDiagonalization() {
        return substitutionModel.canReturnComplexDiagonalization();
    }

    @Override
    public int getEigenBufferCount() {
        return eigenBufferHelper.getBufferCount();
    }

    @Override
    public int getMatrixBufferCount() {
        return matrixBufferHelper.getBufferCount();
    }

    @Override
    public int getSubstitutionModelCount() {
        return 1;
    }

    @Override
    public SubstitutionModel getSubstitutionModel(int index) {
        assert(index == 0);
        return substitutionModel;
    }

    @Override
    public int getEigenIndex(int bufferIndex) {
        return eigenBufferHelper.getOffsetIndex(bufferIndex);
    }

    @Override
    public int getMatrixIndex(int branchIndex) {
        return matrixBufferHelper.getOffsetIndex(branchIndex);
    }

    @Override
    public double[] getRootStateFrequencies() {
        return substitutionModel.getFrequencyModel().getFrequencies();
    }

    @Override
    public void updateSubstitutionModels(Beagle beagle, boolean flip) {
        if (flip) {
            eigenBufferHelper.flipOffset(0);
        }
        EigenDecomposition ed = substitutionModel.getEigenDecomposition();

        beagle.setEigenDecomposition(
                eigenBufferHelper.getOffsetIndex(0),
                ed.getEigenVectors(),
                ed.getInverseEigenVectors(),
                ed.getEigenValues());
    }

    @Override
    public void updateTransitionMatrices(Beagle beagle, int[] branchIndices, double[] edgeLengths, int updateCount, boolean flip) {

        int[] probabilityIndices = new int[updateCount];

        for (int i = 0; i < updateCount; i++) {
            if (flip) {
                matrixBufferHelper.flipOffset(branchIndices[i]);
            }
            probabilityIndices[i] = matrixBufferHelper.getOffsetIndex(branchIndices[i]);
        }// END: i loop

        beagle.updateTransitionMatrices(eigenBufferHelper.getOffsetIndex(0),
                probabilityIndices,
                null, // firstDerivativeIndices
                null, // secondDerivativeIndices
                edgeLengths,
                updateCount);

    }

    @Override
    public void flipTransitionMatrices(int[] branchIndices, int updateCount) {
        for (int i = 0; i < updateCount; i++) {
            matrixBufferHelper.flipOffset(branchIndices[i]);
        }
    }


    @Override
    public void storeState() {
        eigenBufferHelper.storeState();
        matrixBufferHelper.storeState();
    }

    @Override
    public void restoreState() {
        eigenBufferHelper.restoreState();
        matrixBufferHelper.restoreState();
    }

}// END: class

@@ -1884,16 +1884,88 @@
 rHelper;
+%0A    private final BranchModel branchModel;%0A    private final Tree tree;
 %0A%0A    /*
@@ -3019,16 +3019,82 @@
 del();%0A%0A
+        this.branchModel = branchModel;%0A        this.tree = tree;%0A
 %0A       
@@ -4899,32 +4899,427 @@
 lues());%0A    %7D%0A%0A
+    @Override%0A    public SubstitutionModel getSubstitutionModelForBranch(int branchIndex) %7B%0A        BranchModel.Mapping mapping = branchModel.getBranchModelMapping(tree.getNode(branchIndex));%0A        int%5B%5D order = mapping.getOrder();%0A%0A        if (order.length %3E 1) %7B%0A            throw new RuntimeException(%22Not yet implemented%22);%0A        %7D%0A%0A        return getSubstitutionModel(order%5B0%5D);%0A    %7D%0A%0A
     @Override%0A  

copy new function over