Trial function to solve 2 eqns (17 and 18) from Correa et al. (2015c)
        for 1 unknown, i.e. concentration, returned by a minimisation call
771400
def _minimize_c(c, z=0, a_tilde=1, b_tilde=-1,
                Ascaling=900, omega_M_0=0.25, omega_lambda_0=0.75):
    

    # Fn 1 (LHS of Eqn 18)

    Y1 = np.log(2) - 0.5
    Yc = np.log(1+c) - c/(1+c)
    f1 = Y1/Yc

    # Fn 2 (RHS of Eqn 18)

    # Eqn 14 - Define the mean inner density
    rho_2 = 200 * c**3 * Y1 / Yc

    # Eqn 17 rearranged to solve for Formation Redshift
    # essentially when universe had rho_2 density
    zf = (((1 + z)**3 + omega_lambda_0/omega_M_0) *
          (rho_2/Ascaling) - omega_lambda_0/omega_M_0)**(1/3) - 1

    # RHS of Eqn 19
    f2 = ((1 + zf - z)**a_tilde) * np.exp((zf - z) * b_tilde)

    # LHS - RHS should be zero for the correct concentration
    return(f1-f2)
221584
def _g_a (self, x, a, b, s):
        
        return 2*b/(1.0+np.exp(s*(x-a)))
992908
def kinetic(alpha1,lmn1,A,alpha2,lmn2,B):
    
    l1,m1,n1 = lmn1
    l2,m2,n2 = lmn2
    term0 = alpha2*(2*(l2+m2+n2)+3)*\
            overlap(alpha1,(l1,m1,n1),A,\
                           alpha2,(l2,m2,n2),B)
    term1 = -2*pow(alpha2,2)*\
            (overlap(alpha1,(l1,m1,n1),A,
                            alpha2,(l2+2,m2,n2),B)
             + overlap(alpha1,(l1,m1,n1),A,
                              alpha2,(l2,m2+2,n2),B)
             + overlap(alpha1,(l1,m1,n1),A,
                              alpha2,(l2,m2,n2+2),B))
    term2 = -0.5*(l2*(l2-1)*overlap(alpha1,(l1,m1,n1),A,
                                           alpha2,(l2-2,m2,n2),B) +
                  m2*(m2-1)*overlap(alpha1,(l1,m1,n1),A,
                                           alpha2,(l2,m2-2,n2),B) +
                  n2*(n2-1)*overlap(alpha1,(l1,m1,n1),A,
                                           alpha2,(l2,m2,n2-2),B))
    return term0+term1+term2
95755
def _expon_solve_lam_from_mu(mu, b):
    

    def lam_eq(lam, mu, b):
        # Small offset added to denominator to avoid 0/0 erors
        lam, mu, b = Decimal(lam), Decimal(mu), Decimal(b)
        return ( (1 - (lam*b + 1) * np.exp(-lam*b)) /
                 (lam - lam * np.exp(-lam*b) + Decimal(1e-32)) - mu )

    return optim.brentq(lam_eq, -100, 100, args=(mu, b), disp=True)
304942
function eq (y, c) {
  const J0 = base.horner(y * 0.001, c)
  const T = base.J2000Century(J0)
  const W = 35999.373 * D2R * T - 2.47 * D2R
  const Δλ = 1 + 0.0334 * cos(W) + 0.0007 * cos(2 * W)
  let S = 0
  for (let i = terms.length - 1; i >= 0; i--) {
    const t = terms[i]
    S += t.a * cos((t.b + t.c * T) * D2R)
  }
  return J0 + 0.00001 * S / Δλ
}
71392
function q_chi2(df, chi2) {
        var chi, t = 0, k = 0, s = 0;
        if ((df & 1) !== 0) {
            chi = Math.sqrt(chi2);

            if (df === 1) { return 2 * normal__x(chi); }

            t = chi * Math.exp(-0.5 * chi2) / SQ2PI;
            s = t;
            k = 3;
            while (k < df) {
                t = t * (chi2 / k);
                s = s + t;
                k = k + 2;
            }

            return 2 * (normal__x(chi) + s);
        }

        t = Math.exp(-0.5 * chi2);
        s = t;
        k = 2;
        while (k < df) {
            t = t * (chi2 / k);
            s = s + t;
            k = k + 2;
        }

        return s;
    }
845291
def optimize(self, x0, target):
        """"""
        x = x0
        for _ in range(self.maxiter):
            delta = np.linalg.solve(self.h(x, target), -self.g(x, target))
            x = x + delta
            if np.linalg.norm(delta) < self.tol:
                break
        return x
1244593
def _optimize(self, period, tmpid, use_gradient=True):
        """"""
        theta_0 = [self.y.min(), self.y.max() - self.y.min(), 0]
        result = minimize(self._chi2, theta_0, jac=bool(use_gradient),
                          bounds=[(None, None), (0, None), (None, None)],
                          args=(period, tmpid, use_gradient))
        return result.x
898384
def _cho_solve_AATI(A, rho, b, c, lwr, check_finite=True):
    """"""

    N, M = A.shape
    if N >= M:
        x = (b - _cho_solve((c, lwr), b.dot(A).T,
            check_finite=check_finite).T.dot(A.T)) / rho
    else:
        x = _cho_solve((c, lwr), b.T, check_finite=check_finite).T
    return x
360705
def _minimize_constraints_fun_summation(x):
    '''
    
    '''
    summation = sum([x[i] for i in CONSTRAINT_PARAMS_IDX])
    return CONSTRAINT_UPPERBOUND >= summation >= CONSTRAINT_LOWERBOUND
687027
function(poly1, poly2)
{
    var mini1 = findMSA(poly2, poly1.tPlanes);
    if(mini1 == -1) return NONE;
    var min1 = last_MSA_min;
    
    var mini2 = findMSA(poly1, poly2.tPlanes);
    if(mini2 == -1) return NONE;
    var min2 = last_MSA_min;
    
    // There is overlap, find the penetrating verts
    if(min1 > min2)
        return findVerts(poly1, poly2, poly1.tPlanes[mini1].n, min1);
    else
        return findVerts(poly1, poly2, vneg(poly2.tPlanes[mini2].n), min2);
}
430885
@Override
    public void evaluate(BinarySolution solution) {
    int counterOnes;
    int counterZeroes;

    counterOnes = 0;
    counterZeroes = 0;

    BitSet bitset = solution.getVariableValue(0) ;

    for (int i = 0; i < bitset.length(); i++) {
      if (bitset.get(i)) {
        counterOnes++;
      } else {
        counterZeroes++;
      }
    }

    // OneZeroMax is a maximization problem: multiply by -1 to minimize
    solution.setObjective(0, -1.0 * counterOnes);
    solution.setObjective(1, -1.0 * counterZeroes);
  }
64910
def _update_solution_data(self, s, HH, CC, C0):
        
        x = s["x"]
        Va_v = self.om.get_var("Va")
        Pg_v = self.om.get_var("Pg")

        Va = x[Va_v.i1:Va_v.iN + 1]
        Pg = x[Pg_v.i1:Pg_v.iN + 1]
#        f = 0.5 * dot(x.T * HH, x) + dot(CC.T, x)

        s["f"] = s["f"] + C0

        # Put the objective function value in the solution.
#        solution["f"] = f

        return Va, Pg
1179931
function solveLinearIntersect(a,b, c,  d,e, f) {
	// ax + by = c;
	// dx + ey = d;
	console.log('WORKED: ', arguments);
	if(a === 0) {
		var y = c/b;
		var x = (d-e*y)/x;
		return [x, y];
	}
	var denom = e-d/a * b;
	if(denom === 0) {
		// Either infinite solutions or no solutions
		return [undefined, undefined];
	}
	var y = (f-d*c/a)/denom;
	var x = (c-b*y) / a;
	return [x, y];
}
901111
function mpn_cmp(n1, n2, l) {
	while (l > 0) {
		var c = n1[--l] - n2[l];
		if (c != 0) {
			return c;
		}
	}
	return 0;
}
948340
def _initial_guess(self):
        
        a, b, c = np.polyfit(self.volumes, self.energies, 2)
        self.eos_params = [a, b, c]

        v0 = -b/(2*a)
        e0 = a*(v0**2) + b*v0 + c
        b0 = 2 * a * v0
        b1 = 4  # b1 is usually a small number like 4

        vmin, vmax = min(self.volumes), max(self.volumes)

        if not vmin < v0 and v0 < vmax:
            raise EOSError('The minimum volume of a fitted parabola is '
                           'not in the input volumes\n.')

        return e0, b0, b1, v0
1164955
def ccor2(alt, r, h1, zh, h2):
    '''

 '''
    e1 = (alt - zh) / h1;
    e2 = (alt - zh) / h2;
    if ((e1 > 70.0) or (e2 > 70)): # pragma: no cover
        return 1.0 # exp(0)
    if ((e1 < -70) and (e2 < -70)): # pragma: no cover
        return exp(r)
    ex1 = exp(e1);
    ex2 = exp(e2);
    ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
    return exp(ccor2v);