gt:         List<Integer> iv = new ArrayList<Integer>(2);
context: package dr.evomodel.coalescent;
  
  import dr.evolution.coalescent.DemographicFunction;
  import dr.evolution.coalescent.IntervalType;
  import dr.evolution.coalescent.TreeIntervals;
  import dr.evolution.tree.Tree;
  
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.List;
  
  /**
   * @author Joseph Heled
   */
  public class VDdemographicFunction extends DemographicFunction.Abstract {
      private double[] values;
      private double[] times;
      private double[] intervals;
      private double[][] ttimes;
      private double[] alltimes;
      private boolean[] dirtyTrees;
      boolean dirty;
  
      private final VariableDemographicModel.Type type;
  
      TreeIntervals[] ti;
  
      public VDdemographicFunction(Tree[] trees, VariableDemographicModel.Type type,
                                   double[] indicatorParameter, double[] popSizeParameter, boolean logSpace,
                                   boolean mid) {
          super(trees[0].getUnits());
          this.type = type;
  
          ti = new TreeIntervals[trees.length];
          dirtyTrees = new boolean[trees.length];
          Arrays.fill(dirtyTrees, true);
          ttimes = new double[ti.length][];
          int tot = 0;
          for (int k = 0; k < ti.length; ++k) {
              ttimes[k] = new double[trees[k].getTaxonCount() - 1];
              tot += ttimes[k].length;
          }
          alltimes = new double[tot];
  
          setDirty();
  
          assert !(type == VariableDemographicModel.Type.EXPONENTIAL && !logSpace);
  
          setup(trees, indicatorParameter, popSizeParameter, logSpace, mid);
      }
  
      /**
       * Reduce memory footprint of object. After a call to freeze only population/intensity
       * are allowed.
       */
      public void freeze() {
          ttimes = null;
          alltimes = null;
          dirtyTrees = null;
          ti = null;
      }
  
      public VDdemographicFunction(VDdemographicFunction demoFunction) {
          super(demoFunction.getUnits());
          type = demoFunction.type;
  
          this.ti = demoFunction.ti.clone();
          this.values = demoFunction.values.clone();
          this.times = demoFunction.times.clone();
          this.intervals = demoFunction.intervals.clone();
          this.ttimes = demoFunction.ttimes.clone();
          for (int k = 0; k < ttimes.length; ++k) {
              ttimes[k] = ttimes[k].clone();
          }
  
          this.alltimes = demoFunction.alltimes.clone();
          this.dirtyTrees = demoFunction.dirtyTrees.clone();
          this.dirty = demoFunction.dirty;
      }
  
      // Hack so that VDdemo can be used as just a linear piecewise demography (the BEAST one is broken)
      // Alexei fixed PiecewiseLinearPopulation, but did not say yet if it is tested or not.
  
      public VDdemographicFunction(double[] t, double[] p, Type units) {
          this(t, p, units, VariableDemographicModel.Type.LINEAR);
      }
  
      public VDdemographicFunction(double[] t, double[] p, Type units, VariableDemographicModel.Type type) {
          super(units);
  // seem safe to remove them for now
  //        assert t[0] >= 0;
  //        for(int k = 1; k < t.length; ++k) {
  //            assert t[k-1] <= t[k];
  //        }
  //        assert t.length + 1 == p.length;
  //
          this.type = type;
          final int tot = p.length;
          times = new double[tot + 1];
          values = p;
          intervals = new double[tot - 1];
  
          times[0] = 0.0;
          times[tot] = Double.POSITIVE_INFINITY;
  
  //         boolean logSpace = false;
  //        values[0] = logSpace ? Math.exp(p[0]) : p[0];
          System.arraycopy(t, 0, times, 1, t.length);
  
          for (int n = 0; n < intervals.length; ++n) {
              intervals[n] = times[n + 1] - times[n];
          }
          dirty = false;
      }
  
      public int numberOfChanges() {
          return values.length - 1;
      }
  
      public void treeChanged(int nt) {
          dirtyTrees[nt] = true;
          setDirty();
      }
  
      public void setDirty() {
          dirty = true;
      }
  
      private boolean setTreeTimes(int nt, Tree[] trees) {
          if (dirtyTrees[nt]) {
              /*double[] doubles = null;
              if( ! dirtyTrees[nt] ) {
                 doubles = ttimes[nt].clone();
  
              }*/
              ti[nt] = new TreeIntervals(trees[nt]);
  
              TreeIntervals nti = ti[nt];
              // make sure we get each coalescent event individually
              nti.setMultifurcationLimit(0);
              // code probably incorrect for serial samples
              final int nLineages = nti.getIntervalCount();
              assert nLineages >= ttimes[nt].length : nLineages + " " + ttimes[nt].length;
  
              int iCount = 0;
              for (int k = 0; k < ttimes[nt].length; ++k) {
                  double timeToCoal = nti.getInterval(iCount);
                  while (nti.getIntervalType(iCount) != IntervalType.COALESCENT) {
                      ++iCount;
                      timeToCoal += nti.getInterval(iCount);
                  }
  
                  int linAtStart = nti.getLineageCount(iCount);
                  ++iCount;
  
                  assert !(iCount == nLineages && linAtStart != 2);
  
                  int linAtEnd = (iCount == nLineages) ? 1 : nti.getLineageCount(iCount);
  
                  while (linAtStart <= linAtEnd) {
                      ++iCount;
                      timeToCoal += nti.getInterval(iCount);
  
                      linAtStart = linAtEnd;
                      ++iCount;
                      linAtEnd = nti.getLineageCount(iCount);
                  }
                  ttimes[nt][k] = timeToCoal + (k == 0 ? 0 : ttimes[nt][k - 1]);
              }
  
              /*if( doubles != null ) {
                  if( ! Arrays.equals(doubles, ttimes[nt]) ) {
                     System.out.println(Arrays.toString(doubles) + " != " + Arrays.toString(ttimes[nt])
                             + Arrays.toString(dirtyTrees) + " " + dirtyTrees);
                  }
              }*/
              dirtyTrees[nt] = false;
              // System.out.print(nt + " " + Arrays.toString(dirtyTrees) + " " + dirtyTrees);
              return true;
          }
          return false;
      }
  
      void setup(Tree[] trees, double[] indicatorParameter, double[] popSizes, boolean logSpace, boolean mid) {
          // boolean was = dirty;
          if (dirty) {
              // for exponential we do the exp in the code
              if (type == VariableDemographicModel.Type.EXPONENTIAL) logSpace = false;
  
              boolean any = false;
              for (int nt = 0; nt < ti.length; ++nt) {
                  if (setTreeTimes(nt, trees)) {
                      any = true;
                  }
              }
  
              final int nd = indicatorParameter.length;
  
              assert nd == alltimes.length + (type == VariableDemographicModel.Type.STEPWISE ? -1 : 0) :
                      " nd=" + nd + " alltimes.length=" + alltimes.length + " type=" + type;
  
              if (any) {
                  // now we want to merge times together
                  int[] inds = new int[ttimes.length];
  
                  for (int k = 0; k < alltimes.length; ++k) {
                      int j = 0;
                      while (inds[j] == ttimes[j].length) {
                          ++j;
                      }
                      for (int l = j + 1; l < inds.length; ++l) {
                          if (inds[l] < ttimes[l].length) {
                              if (ttimes[l][inds[l]] < ttimes[j][inds[j]]) {
                                  j = l;
                              }
                          }
                      }
                      alltimes[k] = ttimes[j][inds[j]];
                      inds[j]++;
                  }
              }
  
              // assumes lowest node has time 0. this is probably problematic when we come
              // to deal with multiple trees
  
              int tot = 1;
  
              for (int k = 0; k < nd; ++k) {
                  if (indicatorParameter[k] > 0) {
                      ++tot;
                  }
              }
  
              times = new double[tot + 1];
              values = new double[tot];
              intervals = new double[tot - 1];
  
              times[0] = 0.0;
              times[tot] = Double.POSITIVE_INFINITY;
  
              final boolean xx = type == VariableDemographicModel.Type.LINEAR && !logSpace && false;
              if (xx) {
                  double[] a = alltimes;
                  if (mid) {
                      a = new double[alltimes.length];
                      for (int k = 0; k < a.length; ++k) {
                          a[k] = ((alltimes[k] + (k > 0 ? alltimes[k - 1] : 0)) / 2);
                      }
                  }
                  bestLinearFit(a, popSizes, indicatorParameter, times, values);
                  for (int n = 0; n < intervals.length; ++n) {
                      intervals[n] = times[n + 1] - times[n];
                  }
                  for (int n = 0; n < values.length; ++n) {
                      if (values[n] <= 0) {
                          values[n] = 1e-30;
                      }
                  }
              }
  
  
              if (!xx) {
  
                  values[0] = logSpace ? Math.exp(popSizes[0]) : popSizes[0];
  
                  int n = 0;
                  for (int k = 0; k < nd && n + 1 < tot; ++k) {
  
                      if (indicatorParameter[k] > 0) {
                          times[n + 1] = mid ? ((alltimes[k] + (k > 0 ? alltimes[k - 1] : 0)) / 2) : alltimes[k];
  
                          values[n + 1] = logSpace ? Math.exp(popSizes[k + 1]) : popSizes[k + 1];
                          intervals[n] = times[n + 1] - times[n];
                          ++n;
                      }
                  }
              }
              dirty = false;
          }
          //
          /*System.out.println("after setup " + (was ? "(dirty)" : "") + " , alltimes " + Arrays.toString(alltimes)
         + " times " + Arrays.toString(times) + " values " + Arrays.toString(values) +
     " inds " + Arrays.toString(indicatorParameter.getParameterValues())) ;*/
      }
  
      private int ti2f(int i, int j) {
          return (i == 0) ? j : 2 * i + j + 1;
      }
  
      private void
      bestLinearFit(double[] xs, double[] ys, double[] use, double[] ot, double[] oz) {
  
          assert (xs.length + 1) == ys.length;
          assert ys.length == use.length + 2 || ys.length == use.length + 1;
  
          int N = ys.length;
          if (N == 2) {
              // cheaper
              assert xs.length == ot.length;
              assert ys.length == oz.length;
              System.arraycopy(xs, 0, ot, 0, xs.length);
              System.arraycopy(ys, 0, oz, 0, ys.length);
              //return new VDdemographicFunction(xs, ys, getUnits());
          }
  
  