static private void runActiveConjugate(int ntax,double[] d,double[] W,double[] x){
  final boolean collapse_many_negs=true;
  int npairs=d.length;
  if (W.length != npairs || x.length != npairs)   throw new IllegalArgumentException("Vectors d,W,x have different dimensions");
  CircularSplitWeights.runUnconstrainedLS(ntax,d,x);
  boolean all_positive=true;
  for (int k=0; k < npairs && all_positive; k++)   if (x[k] < 0.0)   all_positive=false;
  if (all_positive)   return;
  double[] r=new double[npairs];
  double[] w=new double[npairs];
  double[] p=new double[npairs];
  double[] y=new double[npairs];
  double[] old_x=new double[npairs];
  Arrays.fill(old_x,1.0);
  boolean[] active=new boolean[npairs];
  Arrays.fill(active,false);
  double[] AtWd=new double[npairs];
  for (int k=0; k < npairs; k++)   y[k]=W[k] * d[k];
  CircularSplitWeights.calculateAtx(ntax,y,AtWd);
  boolean first_pass=true;
  while (true) {
    while (true) {
      if (!first_pass)       CircularSplitWeights.circularConjugateGrads(ntax,npairs,r,w,p,y,W,AtWd,active,x);
      first_pass=false;
      if (collapse_many_negs) {
        int[] entriesToContract=worstIndices(x,0.6);
        if (entriesToContract != null) {
          int numToContract=entriesToContract.length;
          for (          int index : entriesToContract) {
            x[index]=0.0;
            active[index]=true;
          }
          CircularSplitWeights.circularConjugateGrads(ntax,npairs,r,w,p,y,W,AtWd,active,x);
        }
      }
      int min_i=-1;
      double min_xi=-1.0;
      for (int i=0; i < npairs; i++) {
        if (x[i] < 0.0) {
          double xi=(old_x[i]) / (old_x[i] - x[i]);
          if ((min_i == -1) || (xi < min_xi)) {
            min_i=i;
            min_xi=xi;
          }
        }
      }
      if (min_i == -1)       break;
 else {
        for (int i=0; i < npairs; i++)         if (!active[i])         old_x[i]+=min_xi * (x[i] - old_x[i]);
        active[min_i]=true;
        x[min_i]=0.0;
      }
    }
    calculateAb(ntax,x,y);
    for (int i=0; i < npairs; i++)     y[i]*=W[i];
    calculateAtx(ntax,y,r);
    int min_i=-1;
    double min_grad=1.0;
    for (int i=0; i < npairs; i++) {
      r[i]-=AtWd[i];
      r[i]*=2.0;
      if (active[i]) {
        double grad_ij=r[i];
        if ((min_i == -1) || (grad_ij < min_grad)) {
          min_i=i;
          min_grad=grad_ij;
        }
      }
    }
    if ((min_i == -1) || (min_grad > -0.0001))     return;
 else     active[min_i]=false;
  }
}

static private void runActiveConjugate ( int ntax , double [ ] d , double [ ] W , double [ ] x ) { final boolean collapse_many_negs = true ; int npairs = d . length ; if ( W . length != npairs || x . length != npairs ) throw new IllegalArgumentException ( "Vectors d,W,x have different dimensions" ) ; CircularSplitWeights . runUnconstrainedLS ( ntax , d , x ) ; boolean all_positive = true ; for ( int k = 0 ; k < npairs && all_positive ; k ++ ) if ( x [ k ] < 0.0 ) all_positive = false ; if ( all_positive ) return ; double [ ] r = new double [ npairs ] ; double [ ] w = new double [ npairs ] ; double [ ] p = new double [ npairs ] ; double [ ] y = new double [ npairs ] ; double [ ] old_x = new double [ npairs ] ; Arrays . fill ( old_x , 1.0 ) ; boolean [ ] active = new boolean [ npairs ] ; Arrays . fill ( active , false ) ; double [ ] AtWd = new double [ npairs ] ; for ( int k = 0 ; k < npairs ; k ++ ) y [ k ] = W [ k ] * d [ k ] ; CircularSplitWeights . calculateAtx ( ntax , y , AtWd ) ; boolean first_pass = true ; while ( true ) { while ( true ) { if ( ! first_pass ) CircularSplitWeights . circularConjugateGrads ( ntax , npairs , r , w , p , y , W , AtWd , active , x ) ; first_pass = false ; if ( collapse_many_negs ) { int [ ] entriesToContract = worstIndices ( x , 0.6 ) ; if ( entriesToContract != null ) { int numToContract = entriesToContract . length ; for ( int index : entriesToContract ) { x [ index ] = 0.0 ; active [ index ] = true ; } CircularSplitWeights . circularConjugateGrads ( ntax , npairs , r , w , p , y , W , AtWd , active , x ) ; } } int min_i = - 1 ; double min_xi = - 1.0 ; for ( int i = 0 ; i < npairs ; i ++ ) { if ( x [ i ] < 0.0 ) { double xi = ( old_x [ i ] ) / ( old_x [ i ] - x [ i ] ) ; if ( ( min_i == - 1 ) || ( xi < min_xi ) ) { min_i = i ; min_xi = xi ; } } } if ( min_i == - 1 ) break ; else { for ( int i = 0 ; i < npairs ; i ++ ) if ( ! active [ i ] ) old_x [ i ] += min_xi * ( x [ i ] - old_x [ i ] ) ; active [ min_i ] = true ; x [ min_i ] = 0.0 ; } } calculateAb ( ntax , x , y ) ; for ( int i = 0 ; i < npairs ; i ++ ) y [ i ] *= W [ i ] ; calculateAtx ( ntax , y , r ) ; int min_i = - 1 ; double min_grad = 1.0 ; for ( int i = 0 ; i < npairs ; i ++ ) { r [ i ] -= AtWd [ i ] ; r [ i ] *= 2.0 ; if ( active [ i ] ) { double grad_ij = r [ i ] ; if ( ( min_i == - 1 ) || ( grad_ij < min_grad ) ) { min_i = i ; min_grad = grad_ij ; } } } if ( ( min_i == - 1 ) || ( min_grad > - 0.0001 ) ) return ; else active [ min_i ] = false ; } }
uses an active set method with the conjugate gradient algorithm to find x that minimises < p / > ( ax - d ) w ( ax - d ) < p / > here , a is the design matrix for the set of cyclic splits with ordering 0 , 1 , 2 , . . . , n - 1 d is the distance vector , with pairs in order ( 0 , 1 ) , ( 0 , 2 ) , . . . , ( 0 , n - 1 ) , ( 1 , 2 ) , ( 1 , 3 ) , . . . , ( 1 , n - 1 ) , . . . , ( n - 1 , n ) w is a vector of variances for d , with pairs in same order as d . x is a vector of split weights , with pairs in same order as d . the split ( i , j ) , for i < j , is { i , i + 1 , . . . , j - 1 } | rest
