TString MillePedeTrees::PhiSwaps(double swapAround, const TString &tree1, const TString &tree2) const
{
  // 'true'/1 if one phi is above, one below 'swapAround'
  return 
    Parenth((Phi(tree1) += Form(">%f", swapAround)) += 
            AndL() += Phi(tree2) += Form("<%f", swapAround))
    += OrL()
    += Parenth((Phi(tree1) += Form("<%f", swapAround)) += 
               AndL() += Phi(tree2) += Form(">%f", swapAround));
}

double chisq(double z,int n)
{double s=0.,t=1.,q,h;
 int i;
 if(z<=0.) return (0.);
 if(n>3000) return Phi((exp(log(z/n)/3.)-1.+2./(9*n))/sqrt(2./(9*n)));
 h=.5*z;
 if(n&1){q=sqrt(z); t=2*exp(-.5*z-.918938533204673)/q;
    for(i=1;i<=(n-2);i+=2){ t=t*z/i; s+=t;}
    return(2*Phi(q)-1-s);
        }
 for(i=1;i<n/2;i++) { t=t*h/i; s+=t;}
 return (1.-(1+s)*exp(-h));
}

void gspmodel(void)
{
    real beta_a, mcut, vcut, vfac;
    static real *sig2 = NULL;
    real r, vmax2, sigr, sig, x, y, vr, v1, v2;
    bodyptr bp;
    vector rhat, vec1, vec2, vtmp;

    beta_a = getdparam("beta_a");
    assert(beta_a <= 1.0);
    nbody = getiparam("nbody");
    assert(nbody > 0);
    mcut = getdparam("mcut");
    assert(0.0 < mcut && mcut <= 1.0);
    vcut = getdparam("vcut");
    assert(vcut > 0.0);
    if (sig2 == NULL)
        sig2 = calc_sig2_gsp(gsp, ggsp, beta_a);
    if (btab == NULL)
	btab = (bodyptr) allocate(nbody * SizeofBody);
    vfac = rsqrt(1 - beta_a);
    for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
	Mass(bp) = gsp->mtot / nbody;
	r = r_mass_gsp(gsp, xrandom(0.0, mcut * gsp->mtot));
	vmax2 = -2 * rsqr(vcut) * phi_gsp(ggsp, r);
	if (vfac > 1.0)
	    vmax2 = vmax2 / rsqr(vfac);
	sigr = rsqrt(sig2_gsp(gsp, ggsp, beta_a, sig2, r));
	sig = fixsigma(sigr, rsqrt(vmax2));
	do {
	    vr = grandom(0.0, sig);
	    v1 = grandom(0.0, sig);
	    v2 = grandom(0.0, sig);
	} while (vr*vr + v1*v1 + v2*v2 > vmax2);
	picktriad(rhat, vec1, vec2);
	MULVS(Pos(bp), rhat, r);
	MULVS(Vel(bp), rhat, vr);
	MULVS(vtmp, vec1, v1 * vfac);
	ADDV(Vel(bp), Vel(bp), vtmp);
	MULVS(vtmp, vec2, v2 * vfac);
	ADDV(Vel(bp), Vel(bp), vtmp);
	Phi(bp) = phi_gsp(ggsp, r);
	Aux(bp) = Phi(bp) + 0.5 * dotvp(Vel(bp), Vel(bp));
    }
    if (getbparam("besort"))
	qsort(btab, nbody, SizeofBody, berank);
    if (getbparam("zerocm"))
	snapcenter(btab, nbody, MassField.offset);
    if (! strnull(getparam("auxvar")))
	setauxvar(btab, nbody);
}

local void hqmforces(bodyptr btab, int nbody, real M, real a, real b,
		     real tol)
{
  bodyptr bp;
  double r, mr3i, params[4], phi0, aR0, az0, abserr[3];
  static gsl_integration_workspace *wksp = NULL;
  gsl_function FPhi, F_aR, F_az;
  static double maxerr = 0.0;
  int stat[3];

  if (a == b) {					// spherical case is easy!
    for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
      r = absv(Pos(bp));
      Phi(bp) = - M / (a + r);
      mr3i = M * rsqr(r / (a + r)) / rqbe(r);
      MULVS(Acc(bp), Pos(bp), - mr3i);
    }
  } else {					// flattened case is harder
    if (wksp == NULL) {				// on first call, initialze
      wksp = gsl_integration_workspace_alloc(1000);
      gsl_set_error_handler_off();		// handle errors below
    }
    FPhi.function = &intPhi;
    F_aR.function = &int_aR;
    F_az.function = &int_az;
    FPhi.params = F_aR.params = F_az.params = params;
    a2(params) = rsqr(a);
    b2(params) = rsqr(b);
    for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
      R(params) = rsqrt(rsqr(Pos(bp)[0]) + rsqr(Pos(bp)[1]));
      z(params) = Pos(bp)[2];
      stat[0] = gsl_integration_qagiu(&FPhi, 0.0, tol, 0.0,
				      1000, wksp, &phi0, &abserr[0]);
      stat[1] = gsl_integration_qagiu(&F_aR, 0.0, tol, 0.0,
				      1000, wksp, &aR0, &abserr[1]);
      stat[2] = gsl_integration_qagiu(&F_az, 0.0, tol, 0.0,
				      1000, wksp, &az0, &abserr[2]);
      if (stat[0] || stat[1] || stat[2])	// any errors reported?
	for (int i = 0; i < 3; i++)
	  if (stat[i] != 0 && abserr[i] > maxerr) {
	    eprintf("[%s.hqmforces: warning: %s  abserr[%d] = %g]\n",
		    getprog(), gsl_strerror(stat[i]), i+1, abserr[i]);
	    maxerr = abserr[i];			// adjust reporting threshold
	  }
      Phi(bp) = - M * phi0;
      Acc(bp)[0] = - M * (Pos(bp)[0] / R(params)) * aR0;
      Acc(bp)[1] = - M * (Pos(bp)[1] / R(params)) * aR0;
      Acc(bp)[2] = - M * az0;
    }
  }
}

    void RBFInterpolation::interpolate(
        const matrix & values,
        matrix & valuesInterpolation
        )
    {
        if ( cpu && not computed )
            compute( positions, positionsInterpolation );

        assert( computed );

        if ( cpu )
        {
            matrix B, valuesLU( n_A, values.cols() ), Phi( n_B, n_A );

            if ( polynomialTerm )
            {
                Phi.resize( n_B, n_A + dimGrid + 1 );
                valuesLU.resize( n_A + dimGrid + 1, values.cols() );
            }

            valuesLU.setZero();
            valuesLU.topLeftCorner( values.rows(), values.cols() ) = values;

            B = lu.solve( valuesLU );

            evaluatePhi( positions, positionsInterpolation, Phi );

            if ( polynomialTerm )
            {
                // Include polynomial contributions in matrix Phi

                for ( int i = 0; i < Phi.rows(); i++ )
                    Phi( i, n_A ) = 1;

                Phi.topRightCorner( n_B, dimGrid ) = positionsInterpolation.block( 0, 0, n_B, dimGrid );
            }

            valuesInterpolation.noalias() = Phi * B;
        }

        if ( not cpu )
        {
            valuesInterpolation.noalias() = Hhat * values;
        }

        assert( valuesInterpolation.rows() == n_B );
        assert( values.cols() == valuesInterpolation.cols() );
    }

void polymodel(void)
{
  gsl_interp_accel *pmsplacc = gsl_interp_accel_alloc();
  bodyptr p;
  real rad, phi, vel, psi, vr, vp, a, E, J;
  vector rhat, vtmp, vper;

  for (p = btab; p < NthBody(btab, nbody); p = NextBody(p)) {
    rad = rad_m(xrandom(0.0, mtot));
    phi = gsl_spline_eval(pmspline, (double) rad, pmsplacc);
    vel = pick_v(phi);
    psi = pick_psi();
    vr = vel * rcos(psi);
    vp = vel * rsin(psi);
    Mass(p) = mtot / nbody;
    pickshell(rhat, NDIM, 1.0);
    MULVS(Pos(p), rhat, rad);
    pickshell(vtmp, NDIM, 1.0);
    a = dotvp(vtmp, rhat);
    MULVS(vper, rhat, - a);
    ADDV(vper, vper, vtmp);
    a = absv(vper);
    MULVS(vper, vper, vp / a);
    MULVS(Vel(p), rhat, vr);
    ADDV(Vel(p), Vel(p), vper);
    Phi(p) = phi;
    E = phi + 0.5 * rsqr(vel);
    J = rad * ABS(vp);
    Aux(p) = Kprime * rpow(phi1 - E, npol - 1.5) * rpow(J, 2 * mpol);
  }
  gsl_interp_accel_free(pmsplacc);
}

local void sum1force(bodyptr btab, int nbody, int i0, real eps2)
{
    bodyptr bp0, bp;
    double phi0, acc0[NDIM];
    vector pos0, dr;
    real dr2, drab, mri, mr3i;
    int j;

    bp0 = NthBody(btab, i0);
    phi0 = 0.0;
    CLRV(acc0);
    SETV(pos0, Pos(bp0));
    for (j = 0; j < nbody; j++)
	if (j != i0) {
	    bp = NthBody(btab, j);
	    DOTPSUBV(dr2, dr, Pos(bp), pos0);
	    dr2 += eps2;
	    drab = rsqrt(dr2);
	    mri = Mass(bp) / drab;
	    phi0 -= mri;
	    mr3i = mri / dr2;
	    ADDMULVS(acc0, dr, mr3i);
	}
    Phi(bp0) = phi0;
    SETV(Acc(bp0), acc0);
}

/********************************************************************
* The executable "e" runs as:  ./e <File> <Size>, where:
* File: is the text file with the network (graph) info
* Size: is the number of Monte Carlo trials.
* Seed: is the seed for the random numbers generator
********************************************************************/
int main(int argc, char *argv[]){
  int i,j,k,S,size;
  double X,V,Q,t;
  pt_net n;
  clock_t  ti;


  if(argc<4){
    printf("\n Some input data is missing! Run as:\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }
  if(argc>4){
    printf("\n You've entered more data than necessary! Run as:\n\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }
  if((size=atoi(argv[2]))<1){
    printf("\n The number of trials can not be less than 1! Run as:\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }

  if((seed=atoi(argv[3]))<0){
    printf("\n The seed can not be negative! Run as:\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }

  n=Initialize(argv[1]);

  //---------------- THE CRUDE MONTE CARLO ALGORITHM ----------------
  ti = clock();

  S=0;
  X=0.0;
  V=0.0;
  for(k=0;k<size;k++){
    Fail(n);
    X=(1-Phi(n));
    S+=X;
    V+=X*X;
  }
  Q=(double)S/size;
  V=(V/size-Q*Q)/(size-1);

  t=(double)(clock()-ti)/CLOCKS_PER_SEC;
 //------------------ PRINT OF OUTPUT AND RESULTS ------------------
 printf("\n  CRUDE MONTE CARLO:");
 printf("\n  Network: %s   Replications: %d   ExecTime=%f",argv[1],size,t);
 printf("\n\n  Q =%1.16f = %1.2e  (Unreliability)",Q,Q);
 printf("\n  V =%1.16f = %1.2e  (Variance)",V,V);
 printf("\n  SD=%1.16f = %1.2e  (Standard Deviation)",sqrt(V),sqrt(V));
 printf("\n  RE=%1.16f = %1.2f%%     (Relative Error)",sqrt(V)/Q,100*sqrt(V)/Q);
 printf("\n  ------------------------------------------------------------\n\n");
 //-----------------------------------------------------------------

  return 0;

}

EXTERN_ENV
#define global extern

#include "stdinc.h"

/*
 * HACKGRAV: evaluate grav field at a given particle.
 */

void hackgrav(bodyptr p, long ProcessId)
{
   Local[ProcessId].pskip = p;
   SETV(Local[ProcessId].pos0, Pos(p));
   Local[ProcessId].phi0 = 0.0;
   CLRV(Local[ProcessId].acc0);
   Local[ProcessId].myn2bterm = 0;
   Local[ProcessId].mynbcterm = 0;
   Local[ProcessId].skipself = FALSE;
   hackwalk(ProcessId);
   Phi(p) = Local[ProcessId].phi0;
   SETV(Acc(p), Local[ProcessId].acc0);
#ifdef QUADPOLE
   Cost(p) = Local[ProcessId].myn2bterm + NDIM * Local[ProcessId].mynbcterm;
#else
   Cost(p) = Local[ProcessId].myn2bterm + Local[ProcessId].mynbcterm;
#endif
}

//function reads in the real and imaginary data, and rotates them to output signal into R and noise into I.
void rotate_data(double * R, double * I, const int size)
{

double phi;
double sumR=0, sumI=0;
std::complex<double> m(0,-1);

for(unsigned int h=0; h<30;h++)
{
sumR=sumR+R[h];
sumI=sumI+I[h];
}
phi = atan2(sumI,sumR);

std::complex<double> Phi(phi,0);

for(unsigned int ii=0;ii<size;ii++)
{
std::complex<double> M(R[ii],I[ii]);
M=M*exp(m*Phi);
R[ii]=real(M);
I[ii]=imag(M);
}


}

inline void trajectory(double* act) {
    /*trajectory tau, t, d, p*/
    int k, j = 1;
    p_str[1].t = p_str[1].p = 0.;
    while ( p_str[j].t <= Tmax ) {
        p_str[j].tau = log( 1/(RAND() + 0.0000001 ) ) / kappa;
        p_str[j].d = randd();
        j++;
        p_str[j].t = p_str[j-1].t + p_str[j-1].tau;
        p_str[j].p = p_str[j-1].p + p_str[j-1].d;
    }


    if (j>1) p_str[j-1].tau = Tmax - p_str[j-1].t;
    else p_str[j].t = Tmax;
    /*end trajectory*/
    j--;

    if (j == 1) act[1] = act[2] = act[3] = act[4] = 0.;
    else {
        act[1] = j-1;
        act[2] = p_str[j].p;
        act[3] = act[4] = 0.;
        for (k=1; k<=j-1; k++)  {
            act[3] += p_str[k].tau * (p_str[k].p + E0*(p_str[k+1].t-p_str[k].tau/2.));
            act[4] += Phi(p_str[k].d) +
                      ( M_PI + p_str[k].tau*( p_str[k].p + E0*( p_str[k].t-p_str[k].tau/2. )*( p_str[k].t-p_str[k].tau/2. )
                                              + E0*E0 * p_str[k].tau * p_str[k].tau * p_str[k].tau / 12. )
                      ) / 2.;
        }
    }
}

void XZHydrostatic_GeoPotential<EvalT, Traits>::
evaluateFields(typename Traits::EvalData workset)
{
#ifndef ALBANY_KOKKOS_UNDER_DEVELOPMENT
  for (int cell=0; cell < workset.numCells; ++cell) {
    for (int node=0; node < numNodes; ++node) {
      for (int level=0; level < numLevels; ++level) {
        ScalarT sum =
        PhiSurf(cell,node) +
        0.5 * Pi(cell,node,level) * E.delta(level) / density(cell,node,level);
        for (int j=level+1; j < numLevels; ++j) sum += Pi(cell,node,j)     * E.delta(j)     / density(cell,node,j);

        Phi(cell,node,level) = sum;
        
        //std::cout <<"Inside GeoP, cell, node, PhiSurf(cell,node)="<<cell<<
        //", "<<node<<", "<<PhiSurf(cell,node) <<std::endl;
      }
    }
  }

  /* OG Debugging statements
  std::cout << "Printing PHI at level 0 ----------------------------------------- \n";
  //for(int level=0; level < numLevels; ++level){
  for (int node=0; node < numNodes; ++node) {
	 //std::cout << "lev= " << level <<  ", phi = " << Phi(23,0,level) <<"\n";
	std::cout << "node = " << node <<  ", phi = " << Phi(23,node,0) <<"\n";
  }
  //}*/

#else
  Kokkos::parallel_for(XZHydrostatic_GeoPotential_Policy(0,workset.numCells),*this);
  cudaCheckError();

#endif
}

void setauxvar(bodyptr btab, int nbody)
{
    bodyptr bp;
    vector jvec;
    real jtot, etot, r0, r1, r;

    if (streq(getparam("auxvar"), "mass"))
	for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp))
	    Aux(bp) = mass_gsp(ggsp, absv(Pos(bp)));
    else if (streq(getparam("auxvar"), "rperi"))
	for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
	    CROSSVP(jvec, Pos(bp), Vel(bp));
	    jtot = absv(jvec);
	    etot = 0.5 * dotvp(Vel(bp), Vel(bp)) + Phi(bp);
	    r0 = 0.0;
	    r1 = absv(Pos(bp));
	    r = 0.5 * (r0 + r1);
	    while ((r1 - r0) > TOL * r) {
		if (rsqrt(2 * (etot - phi_gsp(ggsp, r))) > jtot/r)
		    r1 = r;
		else
		    r0 = r;
		r = 0.5 * (r0 + r1);
	    }
	    Aux(bp) = r;
	}
    else
	error("%s: unknown auxvar option %s\n",
	      getargv0(), getparam("auxvar"));
}

local void sumforces(bodyptr btab, int nbody, bodyptr gtab, int ngrav,
		     real eps2)
{
  bodyptr bp, gp;
  double phi0, acc0[NDIM];
  vector pos0, dr;
  real dr2, dr2i, dr1i, mdr1i, mdr3i;

  for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
    phi0 = 0.0;
    CLRV(acc0);
    SETV(pos0, Pos(bp));
    for (gp = gtab; gp < NthBody(gtab, ngrav); gp = NextBody(gp)) {
      DOTPSUBV(dr2, dr, Pos(gp), pos0);
      dr2i = ((real) 1.0) / (dr2 + eps2);
      dr1i = rsqrt(dr2i);
      mdr1i = Mass(gp) * dr1i;
      mdr3i = mdr1i * dr2i;
      phi0 -= mdr1i;
      ADDMULVS(acc0, dr, mdr3i);
    }
    Phi(bp) = phi0;
    SETV(Acc(bp), acc0);
  }
}

    void RBFInterpolation::buildPhi(
        const matrix & positions,
        const matrix & positionsInterpolation
        )
    {
        n_A = positions.rows();
        n_B = positionsInterpolation.rows();
        int phiColsOld = Phi.cols();

        if ( polynomialTerm )
            Phi.conservativeResize( n_B, n_A + dimGrid + 1 );
        else
            Phi.conservativeResize( n_B, n_A );

        int nNewPoints = Phi.cols() - phiColsOld;

        if ( nNewPoints == Phi.cols() )
            nNewPoints = n_A;

        scalar r = 0;

        for ( int i = 0; i < nNewPoints; i++ )
        {
            int index = Phi.cols() - (i + 1);

            if ( polynomialTerm )
                index = Phi.cols() - 1 - dimGrid - (i + 1);

            for ( int j = 0; j < n_B; j++ )
            {
                r = ( positions.row( index ) - positionsInterpolation.row( j ) ).norm();
                Phi( j, index ) = rbfFunction->evaluate( r );
            }
        }
    }

nemo_main()
{
    stream instr, quadstr, outstr;
    Body *btab = NULL, *bp;
    int nbody, bits;
    real eps_r, eps_t, tsnap = 0.0;

    if (!hasvalue("out"))
        warning("No output supplied (out=)");

    quadstr = stropen(getparam("quad"), "r");
    get_history(quadstr);
    instr = stropen(getparam("in"), "r");
    get_history(instr);
    get_quadfield(quadstr, &eps_r, &eps_t);
    get_snap(instr, &btab, &nbody, &tsnap, &bits);
    if (bits & PhaseSpaceBit == 0)
	error("not enuf info: bits = %o", bits);
    for (bp = btab; bp < btab+nbody; bp++) {
	CLRV(Acc(bp));
	Phi(bp) = 0.0;
    }
    quadinter(btab, nbody, eps_r, eps_t);
    if (hasvalue("out")) {
	outstr = stropen(getparam("out"), "w");
	put_history(outstr);
	bits = bits | PotentialBit | AccelerationBit;
	put_snap(outstr, &btab, &nbody, &tsnap, &bits);
    }
}

Vec3f HairBcsdf::NpIntegrand(float beta, float cosThetaD, float phi, int p, float h) const
{
    float iorPrime = std::sqrt(Eta*Eta - (1.0f - cosThetaD*cosThetaD))/cosThetaD;
    float cosThetaT = std::sqrt(1.0f - (1.0f - cosThetaD*cosThetaD)*sqr(1.0f/Eta));
    Vec3f sigmaAPrime = _sigmaA/cosThetaT;

    float gammaI = std::asin(clamp(h, -1.0f, 1.0f));
    float gammaT = std::asin(clamp(h/iorPrime, -1.0f, 1.0f));
    // The correct internal path length (the one in d'Eon et al.'s paper
    // as well as Marschner et al.'s paper is wrong).
    // The correct factor is also mentioned in "Light Scattering from Filaments", eq. (20)
    float l = 2.0f*std::cos(gammaT);

    float f = Fresnel::dielectricReflectance(1.0f/Eta, cosThetaD*trigInverse(h));
    Vec3f T = std::exp(-sigmaAPrime*l);
    Vec3f Aph = (1.0f - f)*(1.0f - f)*T;
    for (int i = 1; i < p; ++i)
        Aph *= f*T;

    float deltaPhi = phi - Phi(gammaI, gammaT, p);
    deltaPhi = std::fmod(deltaPhi, TWO_PI);
    if (deltaPhi < 0.0f)
        deltaPhi += TWO_PI;

    return Aph*D(beta, deltaPhi);
}

void CParam::S8_Phi(double f_Sigma, double a_Phi, double b_Phi) {
  for (int i_var=1; i_var<=n_var_independent; i_var++) {
    ColumnVector Sigma_inv_temp = Sigma_k_inv_ll.column(i_var) ;
    double a_Phi_tilde = a_Phi + 0.5 * K * f_Sigma ;			// (n_var-#. of balance edit+1)
    double b_Phi_tilde = b_Phi + 0.5 * Sigma_inv_temp.sum();
    Phi(i_var,i_var) = rgamma_fn( a_Phi_tilde, b_Phi_tilde );
  }       // end: for (int i_var)
}

/*
   SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
                                          the semismooth jacobian.
*/
PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db)
{
  PetscErrorCode ierr;
  PetscScalar    *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
  PetscInt       i,nlocal;

  PetscFunctionBegin;
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
  ierr = VecGetArray(snes->xl,&l);CHKERRQ(ierr);
  ierr = VecGetArray(snes->xu,&u);CHKERRQ(ierr);
  ierr = VecGetArray(Da,&da);CHKERRQ(ierr);
  ierr = VecGetArray(Db,&db);CHKERRQ(ierr);
  ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr);

  for (i=0; i< nlocal; i++) {
    if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
      da[i] = 0;
      db[i] = 1;
    } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {                     /* upper bound on variable only */
      da[i] = DPhi(u[i] - x[i], -f[i]);
      db[i] = DPhi(-f[i],u[i] - x[i]);
    } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {                      /* lower bound on variable only */
      da[i] = DPhi(x[i] - l[i], f[i]);
      db[i] = DPhi(f[i],x[i] - l[i]);
    } else if (l[i] == u[i]) {                              /* fixed variable */
      da[i] = 1;
      db[i] = 0;
    } else {                                                /* upper and lower bounds on variable */
      da1   = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
      db1   = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]);
      da2   = DPhi(u[i] - x[i], -f[i]);
      db2   = DPhi(-f[i],u[i] - x[i]);
      da[i] = da1 + db1*da2;
      db[i] = db1*db2;
    }
  }

  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  ierr = VecRestoreArray(snes->xl,&l);CHKERRQ(ierr);
  ierr = VecRestoreArray(snes->xu,&u);CHKERRQ(ierr);
  ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr);
  ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

void Hammersley::GenerateSamples() {
	for (int i = 0; i < m_NumSets; i++)
		for (int j = 0; j < m_NumSamples; j++) {
			glm::vec2 point((float)j / (float)m_NumSamples, Phi(j));  //Pi = (xi, yi) = [1/n,	Phi2(i)]

			m_Samples.push_back(point);
		}
}