/* Copyright (c) 2012,2013 Genome Research Ltd.
 *
 * Author: Stephan Schiffels <stephan.schiffels@sanger.ac.uk>
 *
 * This file is part of msmc.
 * msmc is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 3 of the License, or (at your option) any later
 * version.
 * 
 * This program 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 General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License along with
 * this program.  If not, see <http://www.gnu.org/licenses/>.
 */

module model.transition_rate;
import std.math;
import std.conv;
import std.stdio;
import std.exception;
import std.parallelism;
import std.range;
import model.time_intervals;
import model.triple_index_marginal;
import model.coalescence_rate;
import model.rate_integrator;

class TransitionRate {
  double rho;
  const TimeIntervals timeIntervals;
  const PiecewiseConstantCoalescenceRate coal;
  const MarginalTripleIndex marginalIndex;
  private CoalescenceRateIntegrator coalIntegrator;

  private double[] transitionProbabilitiesQ1;
  private double[][] transitionProbabilitiesQ2;
  
  this(in MarginalTripleIndex marginalIndex, in PiecewiseConstantCoalescenceRate coal, in TimeIntervals timeIntervals, double rho) {
    enforce(rho > 0.0, "need positive recombination rate");
    this.coal = coal;
    this.timeIntervals = timeIntervals;
    this.rho = rho;
    this.marginalIndex = marginalIndex;
    coalIntegrator = new CoalescenceRateIntegrator(timeIntervals, coal);
    fillTransitionProbabilitiesParallel();
    // fillTransitionProbabilitiesSingleThread();
  }
  
  private void fillTransitionProbabilitiesSingleThread() {
    transitionProbabilitiesQ2 = new double[][](marginalIndex.nrStates, marginalIndex.nrStates);
    transitionProbabilitiesQ1 = new double[](marginalIndex.nrStates);
    foreach(bv; 0 .. marginalIndex.nrMarginals) {
      auto sum = 0.0;
      foreach(au; 0 .. marginalIndex.nrMarginals) {
        auto aij = marginalIndex.getIndexFromMarginalIndex(au);
        auto bkl = marginalIndex.getIndexFromMarginalIndex(bv);
        transitionProbabilitiesQ2[au][bv] = transitionProbabilityOffDiagonal(aij, bkl);
        sum += transitionProbabilitiesQ2[au][bv] * marginalIndex.getDegeneracyForMarginalIndex(au);
      }
      transitionProbabilitiesQ1[bv] = 1.0 - sum;
    }
  }

  private void fillTransitionProbabilitiesParallel() {
    transitionProbabilitiesQ2 = new double[][](marginalIndex.nrStates, marginalIndex.nrStates);
    transitionProbabilitiesQ1 = new double[](marginalIndex.nrStates);
    foreach(bv; taskPool.parallel(iota(marginalIndex.nrMarginals))) {
      auto sum = 0.0;
      foreach(au; 0 .. marginalIndex.nrMarginals) {
        auto aij = marginalIndex.getIndexFromMarginalIndex(au);
        auto bkl = marginalIndex.getIndexFromMarginalIndex(bv);
        transitionProbabilitiesQ2[au][bv] = transitionProbabilityOffDiagonal(aij, bkl);
        sum += transitionProbabilitiesQ2[au][bv] * marginalIndex.getDegeneracyForMarginalIndex(au);
      }
      transitionProbabilitiesQ1[bv] = 1.0 - sum;
    }
  }
  
  private double transitionProbabilityOffDiagonal(size_t aij, size_t bkl) const {
    auto aij_triple = marginalIndex.getTripleFromIndex(aij);
    auto bkl_triple = marginalIndex.getTripleFromIndex(bkl);
    auto a = aij_triple.time;
    auto i = aij_triple.ind1;
    auto j = aij_triple.ind2;
    auto b = bkl_triple.time;
    auto k = bkl_triple.ind1;
    auto l = bkl_triple.ind2;
    
    double ret;
    if(a < b) {
      ret = q2IntegralSmaller(timeIntervals.leftBoundary(a), i, j,
                                   timeIntervals.rightBoundary(a), a, b);
    }
    if(a > b) {
      ret = q2IntegralGreater(timeIntervals.leftBoundary(a), i, j,
                                   timeIntervals.rightBoundary(a), k, l, a, b);
    }
    if(a == b) {
      ret = q2IntegralSmaller(timeIntervals.leftBoundary(a), i, j, 
                                   timeIntervals.meanTimeWithLambda(b, coal.getTotalMarginalLambda(b)), a, b) +
                 q2IntegralGreater(timeIntervals.meanTimeWithLambda(b, coal.getTotalMarginalLambda(b)),i, j, 
                                   timeIntervals.rightBoundary(a), k, l, a, b);
      
    }
    return ret;
  }
  
  private double q2IntegralSmaller(double t1, size_t i, size_t j, double t2, size_t a, size_t b) const
    in {
      assert(t1 >= timeIntervals.leftBoundary(a) && t1 <= timeIntervals.rightBoundary(a));
      assert(t2 >= timeIntervals.leftBoundary(a) && t2 <= timeIntervals.rightBoundary(a));
      assert(t1 < t2 && t2 <= timeIntervals.meanTimeWithLambda(b, coal.getTotalMarginalLambda(b)));
    }
  body {
    auto meanTime = timeIntervals.meanTimeWithLambda(b, coal.getTotalMarginalLambda(b));
    double sum = 0.0;
    foreach(g; 0 .. a) {
      foreach(pair_index; 0 .. 2) {
        auto m = pair_index == 0 ? i : j;
        double integ = exp(-(t1 - timeIntervals.leftBoundary(a)) * 
                           coal.getSemiMarginalLambda(a, m)) -
                       exp(-(t2 - timeIntervals.leftBoundary(a)) * 
                           coal.getSemiMarginalLambda(a, m));
        integ /= coal.getSemiMarginalLambda(a, m);
        sum += (1.0 - exp(-timeIntervals.delta(g) * coal.getSemiMarginalLambda(g, m))) *
          coalIntegrator.integrateSemiMarginalLambda(m, timeIntervals.rightBoundary(g),
              timeIntervals.leftBoundary(a), g + 1, a) / coal.getSemiMarginalLambda(g, m) * integ;
      }
    }
    foreach(pair_index; 0 .. 2) {
      auto m = pair_index == 0 ? i : j;
      double integ = exp(-(t1 - timeIntervals.leftBoundary(a)) * 
                         coal.getSemiMarginalLambda(a, m)) -
                     exp(-(t2 - timeIntervals.leftBoundary(a)) * 
                         coal.getSemiMarginalLambda(a, m));
      integ /= coal.getSemiMarginalLambda(a, m);
      sum += ((t2 - t1) - integ) / coal.getSemiMarginalLambda(a, m);
    }
    
    double ret = (1.0 - exp(-rho * marginalIndex.nrIndividuals * meanTime)) /
        (meanTime * marginalIndex.nrIndividuals) * coal.getLambda(a, i, j) * sum;
    return ret;
  }

  private double q2IntegralGreater(double t1, size_t i, size_t j, double t2, size_t k, size_t l, size_t a, size_t b) const
    in {
      assert(t1 >= timeIntervals.leftBoundary(a) && t1 <= timeIntervals.rightBoundary(a));
      assert(t2 >= timeIntervals.leftBoundary(a) && t2 <= timeIntervals.rightBoundary(a));
      assert(t1 < t2 && t1 >= timeIntervals.meanTimeWithLambda(b, coal.getTotalMarginalLambda(b)));
    }
  body {
    auto meanTime = timeIntervals.meanTimeWithLambda(b, coal.getTotalMarginalLambda(b));
    double integ = coalIntegrator.integrateTotalMarginalLambda(meanTime, timeIntervals.leftBoundary(a), b, a) / 
                   coal.getTotalMarginalLambda(a) *
                   (exp(-(t1 - timeIntervals.leftBoundary(a)) * coal.getTotalMarginalLambda(a)) -
                    exp(-(t2 - timeIntervals.leftBoundary(a)) * coal.getTotalMarginalLambda(a)));
    double sum = 0.0;
    foreach(g; 0 .. b) {
      foreach(pair_index; 0 .. 2) {
        auto m = pair_index == 0 ? k : l;
        sum += (1.0 - exp(-coal.getSemiMarginalLambda(g, m) * timeIntervals.delta(g))) / 
               coal.getSemiMarginalLambda(g, m) * 
               coalIntegrator.integrateSemiMarginalLambda(m, timeIntervals.rightBoundary(g), meanTime, g + 1, b);
      }
    }
    foreach(pair_index; 0 .. 2) {
      auto m = pair_index == 0 ? k : l;
      sum += (1.0 - exp(-coal.getSemiMarginalLambda(b, m) * (meanTime - timeIntervals.leftBoundary(b)))) / 
             coal.getSemiMarginalLambda(b, m);
    }
      
    double ret = integ * (1.0 - exp(-rho * marginalIndex.nrIndividuals * meanTime)) /
          (meanTime * marginalIndex.nrIndividuals) * coal.getLambda(a, i, j) * sum;
    
    return ret;
  }
  
  double transitionProbabilityQ1(size_t au) const {
    return transitionProbabilitiesQ1[au];
  }
  
  double transitionProbabilityQ2(size_t au, size_t bv) const {
    return transitionProbabilitiesQ2[au][bv];
  }
  
  double transitionProbability(size_t aij, size_t bkl) const
    in {
      assert(aij < marginalIndex.nrStates && bkl < marginalIndex.nrStates);
    }
  body {
    auto au = marginalIndex.getMarginalIndexFromIndex(aij);
    auto bv = marginalIndex.getMarginalIndexFromIndex(bkl);
    double ret = transitionProbabilitiesQ2[au][bv];
    if(aij == bkl) {
      ret += transitionProbabilitiesQ1[au];
    }
    return ret;
  }
  
  double transitionProbabilityMarginal(size_t au, size_t bv) const {
    double ret = marginalIndex.getDegeneracyForMarginalIndex(au) * transitionProbabilitiesQ2[au][bv];
    if(au == bv)
      ret += transitionProbabilitiesQ1[au];
    return ret;
  }
  
  double equilibriumProbability(size_t aij) const {
    // return 1.0 / marginalIndex.nrStates;
    
    auto triple = marginalIndex.getTripleFromIndex(aij);
    auto a = triple.time;
    auto i = triple.ind1;
    auto j = triple.ind2;
    return coal.getLambda(a, i, j) / coal.getTotalMarginalLambda(a) * 
        coalIntegrator.integrateTotalMarginalLambda(0.0, timeIntervals.leftBoundary(a), 0, a) *
        (1.0 - exp(-timeIntervals.delta(a) * coal.getTotalMarginalLambda(a)));
  }
}

unittest {
  writeln("test transitionRate.transitionProbabilityMarginals");
  auto subpop_labels = [0UL, 0, 1, 1];
  auto T = 4UL;
  auto r = 0.01;
  auto lambda_subpop_rates = [1, 0.1, 1, 2, 0.5, 4, 2, 1.0, 4, 1., 2, 1];
  auto lvl = 1.0e-8;
  
  auto timeIntervals = TimeIntervals.standardIntervals(T, 4UL);
  auto marginalIndex = new MarginalTripleIndex(T, subpop_labels);
  auto coal = new PiecewiseConstantCoalescenceRate(marginalIndex, lambda_subpop_rates);
  auto transitionRate = new TransitionRate(marginalIndex, coal, timeIntervals, r);
  
  foreach(aij; 0 .. marginalIndex.nrStates) {
    auto au = marginalIndex.getMarginalIndexFromIndex(aij);
    foreach(bkl; 0 .. marginalIndex.nrStates) {
      auto bv = marginalIndex.getMarginalIndexFromIndex(bkl);
      auto q = transitionRate.transitionProbability(aij, bkl);
      auto qFast = transitionRate.transitionProbabilitiesQ2[au][bv];
      if(aij == bkl)
        qFast += transitionRate.transitionProbabilitiesQ1[au];
      assert(approxEqual(q, qFast, lvl, 0.0), text(q, " ", qFast));
    }
  }
}

unittest {
  writeln("test transitionRate.equilibriumProbability");
  auto subpop_labels = [0UL, 0, 1, 1];
  auto T = 4UL;
  auto r = 0.01;
  auto lambda_subpop_rates = [1, 0.1, 1, 2, 0.5, 4, 2, 1.0, 4, 1., 2, 1];
  auto lvl = 1.0e-8;
  
  auto boundaries = TimeIntervals.getLiAndDurbinBoundaries(T, 1.0 / 6.0);
  auto timeIntervals = new TimeIntervals(boundaries);
  auto marginalIndex = new MarginalTripleIndex(T, subpop_labels);
  auto coal = new PiecewiseConstantCoalescenceRate(marginalIndex, lambda_subpop_rates);
  auto transitionRate = new TransitionRate(marginalIndex, coal, timeIntervals, r);
  
  auto q0Th = [0.03973331222724081,0.003973331222724081,0.003973331222724081,0.003973331222724081,
   0.003973331222724081,0.03973331222724081,0.1560399350214053,0.03900998375535133,
   0.03900998375535133,0.03900998375535133,0.03900998375535133,0.3120798700428106,0.0557645826740214,
   0.0278822913370107,0.0278822913370107,0.0278822913370107,0.0278822913370107,0.1115291653480428,
   0.00016573971988938216,0.0003314794397787643,0.0003314794397787643,0.0003314794397787643,
   0.0003314794397787643,0.00016573971988938216];
  
  foreach(aij; 0 .. marginalIndex.nrStates) {
    assert(approxEqual(transitionRate.equilibriumProbability(aij), q0Th[aij], lvl, 0.0),
      text(aij, ", ", transitionRate.equilibriumProbability(aij), ", ", q0Th[aij]));
  }
}

unittest {  
  writeln("test transitionRate.transitionProbability");
  auto subpop_labels = [0UL, 0, 1, 1];
  auto T = 4UL;
  auto r = 0.01;
  auto lambda_subpop_rates = [1, 0.1, 1, 2, 0.5, 4, 2, 1.0, 4, 1., 2, 1];
  auto lvl = 1.0e-8;
  
  auto boundaries = TimeIntervals.getLiAndDurbinBoundaries(T, 1.0 / 6.0);
  auto timeIntervals = new TimeIntervals(boundaries);
  auto marginalIndex = new MarginalTripleIndex(T, subpop_labels);
  auto coal = new PiecewiseConstantCoalescenceRate(marginalIndex, lambda_subpop_rates);
  auto transitionRate = new TransitionRate(marginalIndex, coal, timeIntervals, r);

  auto qTh =[
    [0.9996011921363851,0.000012455678497214317,0.000012455678497214317,0.000012455678497214317,0.000012455678497214317,0.000012455678497214317,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687],
    [1.2455678497214319e-6,0.9995899820257377,1.2455678497214319e-6,1.2455678497214319e-6,1.2455678497214319e-6,1.2455678497214319e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6],
    [1.2455678497214319e-6,1.2455678497214319e-6,0.9995899820257377,1.2455678497214319e-6,1.2455678497214319e-6,1.2455678497214319e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6],
    [1.2455678497214319e-6,1.2455678497214319e-6,1.2455678497214319e-6,0.9995899820257377,1.2455678497214319e-6,1.2455678497214319e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6],
    [1.2455678497214319e-6,1.2455678497214319e-6,1.2455678497214319e-6,1.2455678497214319e-6,0.9995899820257377,1.2455678497214319e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6880908166079487e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6818944004689653e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6,1.6646828868169687e-6],
    [0.000012455678497214317,0.000012455678497214317,0.000012455678497214317,0.000012455678497214317,0.000012455678497214317,0.9996011921363851,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.000016880908166079486,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.00001681894400468965,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687,0.000016646828868169687],
    [0.00006578189470246562,0.00006578189470246562,0.00006578189470246562,0.00006578189470246562,0.00006578189470246562,0.00006578189470246562,0.9981438612501496,0.00033341450658784576,0.00033341450658784576,0.00033341450658784576,0.00033341450658784576,0.0003188708985366713,0.000503404910386131,0.000503404910386131,0.000503404910386131,0.000503404910386131,0.000503404910386131,0.000503404910386131,0.0004982533619386321,0.0004982533619386321,0.0004982533619386321,0.0004982533619386321,0.0004982533619386321,0.0004982533619386321],
    [0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.00008535958825037931,0.9979934996319003,0.00008172368623758571,0.00008172368623758571,0.00008172368623758571,0.00007808778422479208,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984],
    [0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.00008535958825037931,0.00008172368623758571,0.9979934996319003,0.00008172368623758571,0.00008172368623758571,0.00007808778422479208,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984],
    [0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.00008535958825037931,0.00008172368623758571,0.00008172368623758571,0.9979934996319003,0.00008172368623758571,0.00007808778422479208,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984],
    [0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.000016445473675616406,0.00008535958825037931,0.00008172368623758571,0.00008172368623758571,0.00008172368623758571,0.9979934996319003,0.00007808778422479208,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011498273827118802,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984,0.00011380607285809984],
    [0.00013156378940493125,0.00013156378940493125,0.00013156378940493125,0.00013156378940493125,0.00013156378940493125,0.00013156378940493125,0.0006698371827280286,0.0006407499666256798,0.0006407499666256798,0.0006407499666256798,0.0006407499666256798,0.9986393115063381,0.0008329139915667465,0.0008329139915667465,0.0008329139915667465,0.0008329139915667465,0.0008329139915667465,0.0008329139915667465,0.0008243904418523335,0.0008243904418523335,0.0008243904418523335,0.0008243904418523335,0.0008243904418523335,0.0008243904418523335],
    [0.000023508724898444774,0.000023508724898444774,0.000023508724898444774,0.000023508724898444774,0.000023508724898444774,0.000023508724898444774,0.0001599192048218725,0.000148447718781806,0.000148447718781806,0.000148447718781806,0.000148447718781806,0.00013697623274173948,0.9942786407732027,0.0009833078915005644,0.0009833078915005644,0.0009833078915005644,0.0009833078915005644,0.0008946768863262026,0.0030932299624244542,0.0030932299624244542,0.0030932299624244542,0.0030932299624244542,0.0030932299624244542,0.0030932299624244542],
    [0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.00007995960241093625,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.00006848811637086974,0.0004987122153383788,0.9941112673182289,0.0004543967127511977,0.0004543967127511977,0.0004543967127511977,0.0004100812101640168,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968],
    [0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.00007995960241093625,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.00006848811637086974,0.0004987122153383788,0.0004543967127511977,0.9941112673182289,0.0004543967127511977,0.0004543967127511977,0.0004100812101640168,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968],
    [0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.00007995960241093625,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.00006848811637086974,0.0004987122153383788,0.0004543967127511977,0.0004543967127511977,0.9941112673182289,0.0004543967127511977,0.0004100812101640168,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968],
    [0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.000011754362449222387,0.00007995960241093625,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.000074223859390903,0.00006848811637086974,0.0004987122153383788,0.0004543967127511977,0.0004543967127511977,0.0004543967127511977,0.9941112673182289,0.0004100812101640168,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968,0.0012705077347738968],
    [0.00004701744979688955,0.00004701744979688955,0.00004701744979688955,0.00004701744979688955,0.00004701744979688955,0.00004701744979688955,0.000319838409643745,0.000296895437563612,0.000296895437563612,0.000296895437563612,0.000296895437563612,0.00027395246548347896,0.0018458199293571766,0.0016685579190084529,0.0016685579190084529,0.0016685579190084529,0.0016685579190084529,0.9955983352430873,0.003977601953342267,0.003977601953342267,0.003977601953342267,0.003977601953342267,0.003977601953342267,0.003977601953342267],
    [6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,4.753010412908107e-7,4.412062665821081e-7,4.412062665821081e-7,4.412062665821081e-7,4.412062665821081e-7,4.071114918734055e-7,4.407511543766052e-6,3.7061412359603403e-6,3.7061412359603403e-6,3.7061412359603403e-6,3.7061412359603403e-6,3.004770928154629e-6,0.9806496072786828,0.0006160033952568582,0.0006160033952568582,0.0006160033952568582,0.0006160033952568582,0.0005810561655501418],
    [1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,9.506020825816214e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.14222983746811e-7,8.815023087532104e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,6.009541856309258e-6,0.0012061320013822918,0.9814843664927553,0.001136237541968859,0.001136237541968859,0.001136237541968859,0.001066343082555426],
    [1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,9.506020825816214e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.14222983746811e-7,8.815023087532104e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,6.009541856309258e-6,0.0012061320013822918,0.001136237541968859,0.9814843664927553,0.001136237541968859,0.001136237541968859,0.001066343082555426],
    [1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,9.506020825816214e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.14222983746811e-7,8.815023087532104e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,6.009541856309258e-6,0.0012061320013822918,0.001136237541968859,0.001136237541968859,0.9814843664927553,0.001136237541968859,0.001066343082555426],
    [1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,1.397420833363443e-7,9.506020825816214e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.824125331642162e-7,8.14222983746811e-7,8.815023087532104e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,7.4122824719206806e-6,6.009541856309258e-6,0.0012061320013822918,0.001136237541968859,0.001136237541968859,0.001136237541968859,0.9814843664927553,0.001066343082555426],
    [6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,6.987104166817215e-8,4.753010412908107e-7,4.412062665821081e-7,4.412062665821081e-7,4.412062665821081e-7,4.412062665821081e-7,4.071114918734055e-7,4.407511543766052e-6,3.7061412359603403e-6,3.7061412359603403e-6,3.7061412359603403e-6,3.7061412359603403e-6,3.004770928154629e-6,0.0005551813764187171,0.0005202341467120008,0.0005202341467120008,0.0005202341467120008,0.0005202341467120008,0.9811828881648589]
  ];
  
  foreach(aij; 0 .. marginalIndex.nrStates) {
    foreach(bkl; 0 .. marginalIndex.nrStates) {
      assert(approxEqual(transitionRate.transitionProbability(aij, bkl), qTh[aij][bkl], lvl, 0.0),
             text([aij, bkl, transitionRate.transitionProbability(aij, bkl), qTh[aij][bkl]]));
    }
  }
}