Regularized Upper/Complementary Incomplete Gamma Function
Q(s,x) = 1 - P(s,x) = 1 - <i><big>&#8747;</big><sub><small>0</small></sub><sup><small>x</small></sup> e<sup>-t</sup> t<sup>(s-1)</sup> dt</i>
487231
public static double regularizedUpperIncompleteGamma(double s, double x) {
        if (s < 0.0) {
            throw new IllegalArgumentException("Invalid s: " + s);
        }

        if (x < 0.0) {
            throw new IllegalArgumentException("Invalid x: " + x);
        }

        double igf = 0.0;

        if (x != 0.0) {
            if (x == 1.0 / 0.0) {
                igf = 1.0;
            } else {
                if (x < s + 1.0) {
                    // Series representation
                    igf = 1.0 - regularizedIncompleteGammaSeries(s, x);
                } else {
                    // Continued fraction representation
                    igf = 1.0 - regularizedIncompleteGammaFraction(s, x);
                }
            }
        }
        return igf;
    }
1171203
def primal_dual_solve(func, x0, Dfunc, A, b, G, h, args=(), tol=1e-10,
                      maxiter=100, show_progress=True, full_output=False):
    """"""
    fcalls, func = wrap_function(func, args)
    Dfcalls, Dfunc = wrap_function(Dfunc, args)

    M, N = G.shape
    P, N = A.shape

    """Total number of inequality constraints"""
    m = M

    def gap(z):
        r"""Gap-function"""
        l = z[N+P:N+P+M]
        s = z[N+P+M:]
        return mydot(l, s)/m

    def centrality(z):
        r"""Centrality function"""
        l = z[N+P:N+P+M]
        s = z[N+P+M:]
        return np.min(l*s)

    def KKT(z, sigma=0.0):
        r"""KKT system (possible perturbed).Primal variable"""
        x = z[0:N]

        """Multiplier for equality constraints"""
        nu = z[N:N+P]

        """Multiplier for inequality constraints"""
        l = z[N+P:N+P+M]

        """Slacks"""
        s = z[N+P+M:]

        """Evaluate objective function"""
        F = func(x)

        """Dual infeasibility"""
        rdual = F+mydot(A.transpose(), nu)+mydot(G.transpose(), l)

        """Primal infeasibilities"""
        rprim1 = mydot(A, x)-b
        rprim2 = mydot(G, x)-h+s

        """Complementary slackness (perturbed)"""
        mu = gap(z)
        rslack = l*s-sigma*mu

        return np.hstack((rdual, rprim1, rprim2, rslack))

    def step_fast(z, KKTval, LU, G, A, mu, beta, gamma, alpha0):
        r"""Affine scaling step."""
        dz = solve_factorized(z, KKTval, LU, G, A)

        """Reduce step length until slacks s and multipliers l are positive"""
        alpha = 1.0*alpha_0
        k = 0
        for k in range(10):
            z_new = z + alpha*dz
            if np.all( z_new[N+P:] > 0.0 ):
                break
            alpha *= 0.5
            k += 1
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")

        """Reduce step length until iterates lie in correct neighborhood"""
        for k in range(10):
            z_new = z + alpha*dz
            KKTval_new = KKT(z_new)
            dual = mynorm(KKTval_new[0:N])
            prim = mynorm(KKTval_new[N:N+P+M])
            mu_new = gap(z_new)
            cent_new = centrality(z_new)

            if (dual <= beta*mu_new and prim <= beta*mu_new and
                cent_new >= gamma*mu_new):
                break
            alpha *= 0.5
            # alpha *= 0.95
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")
        return z_new, mu_new

    def step_safe(z, KKTval, LU, G, A, mu, beta, gamma):
        r"""Centering step."""
        dz = solve_factorized(z, KKTval, LU, G, A)

        """Reduce step length until slacks s and multipliers l are positive"""
        alpha = 1.0
        k = 0
        for k in range(10):
            z_new = z + alpha*dz
            if np.all( z_new[N+P:] > 0.0 ):
                break
            alpha *= 0.5
            k += 1
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction (pos.) reached")

        """Reduce step length until iterates lie in correct neighborhood
        and mu fulfills Armijo condition"""
        k = 0
        for k in range(10):
            z_new = z+alpha*dz
            KKTval_new = KKT(z_new, sigma=SIGMA)
            dual = mynorm(KKTval_new[0:N])
            prim = mynorm(KKTval_new[N:N+P+M])
            mu_new = gap(z_new)
            cent_new = centrality(z_new)

            if (dual <= beta*mu_new and prim <= beta*mu_new and
                cent_new >= gamma*mu_new and
                mu_new<=(1.0-KAPPA*alpha*(1.0-SIGMA))*mu):
                break
            alpha *= 0.5
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")
        return z_new, mu_new

    """INITIALIZATIONInitial Slacks for inequality constraints"""
    s0 = -1.0*(mydot(G, x0)-h)

    """Initial multipliers for inequality constraints"""
    l0 = 1.0*np.ones(M)

    """Initial multipliers for equality constraints"""
    nu0 = np.zeros(P)

    """Initial point"""
    z0 = np.hstack((x0, nu0, l0, s0))

    """Initial KKT-values"""
    KKTval0 = KKT(z0, sigma=0.0)
    mu0 = gap(z0)
    dual0 = mynorm(KKTval0[0:N])
    prim0 = mynorm(KKTval0[N:N+P+M])

    """Initial neighborhood"""
    beta = BETA * np.sqrt(dual0**2 + prim0**2)/mu0
    gamma = 1.0 * GAMMA_MAX

    """Number of fast steps"""
    t = 0

    """Number of iterations"""
    n = 0

    """Dummy variable for step type"""
    step_type = " "

    if show_progress:
        print("%s %s %s %s %s %s %s" %("iter", "gap", "dual", "primal",
                                       "min", "max", "step"))
    """MAIN LOOP"""
    z = z0
    x = z0[0:N]
    KKTval = KKTval0
    dual = dual0
    prim = prim0
    mu = mu0
    Dfunc_val = Dfunc(x)
    LU = factor(z, Dfunc_val, G, A)

    if full_output:
        info = {'z': []}
        info['z'].append(z)

    for n in range(maxiter):
        if show_progress:
            l=z[N+P:N+P+M]
            s=z[N+P+M:]
            print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                     np.min(l*s), np.max(l*s),
                                                     step_type))
        """Attempt fast step"""
        beta_new = (1.0 + GAMMA_BAR**(t+1)) * beta
        gamma_new = GAMMA_MIN + GAMMA_BAR**(t+1)*(GAMMA_MAX - GAMMA_MIN)
        alpha_0 = 1.0 - np.sqrt(mu)/GAMMA_BAR**t

        if alpha_0 > 0.0:
            z_new, mu_new = step_fast(z, KKTval, LU, G, A, mu,
                                      beta_new, gamma_new, alpha_0)
            if mu_new < RHO * mu:
                """Fast successful"""
                z = z_new
                mu = mu_new
                beta = beta_new
                gamma = gamma_new
                t += 1
                step_type = "f"
            else:
                """Perturbed right-had side"""
                KKTval_pert = 1.0*KKTval
                KKTval_pert[N+P+M:] -= SIGMA * mu
                z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
                step_type = "s"
        else:
            """Perturbed right-hand side"""
            KKTval_pert = 1.0*KKTval
            KKTval_pert[N+P+M:] -= SIGMA * mu
            z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
            step_type = "s"

        """Compute new iterates"""
        KKTval = KKT(z, sigma=0.0)
        dual = mynorm(KKTval[0:N])
        prim = mynorm(KKTval[N:N+P+M])
        x = z[0:N]
        Dfunc_val = Dfunc(x)
        LU = factor(z, Dfunc_val, G, A)
        if full_output:
            info['z'].append(z)
        if mu < tol and dual < tol and prim < tol:
            break
    if n == maxiter - 1:
        raise RuntimeError("Maximum number of iterations reached")

    if show_progress:
        l=z[N+P:N+P+M]
        s=z[N+P+M:]
        print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                 np.min(l*s), np.max(l*s),
                                                 step_type))
    if full_output:
        return z[0:N], info
    else:
        return z[0:N]
1190905
def lnprior(pars):
    
    # Limit norm and B to be positive
    logprob = (
        naima.uniform_prior(pars[0], 0.0, np.inf)
        + naima.uniform_prior(pars[1], -1, 5)
        + naima.uniform_prior(pars[3], 0, np.inf)
    )

    return logprob
270377
def prox_max_entropy(X, step, gamma=1):
    
    from scipy.special import lambertw
    gamma_ = _step_gamma(step, gamma)
    # minimize entropy: return gamma_ * np.real(lambertw(np.exp((X - gamma_) / gamma_) / gamma_))
    above = X > 0
    X[above] = gamma_ * np.real(lambertw(np.exp(X[above]/gamma_ - 1) / gamma_))
    return X
950058
def _msge_with_gradient_underdetermined(data, delta, xvschema, skipstep, p):
    
    t, m, l = data.shape
    d = None
    j, k = 0, 0
    nt = np.ceil(t / skipstep)
    for trainset, testset in xvschema(t, skipstep):

        a, b = _construct_var_eqns(atleast_3d(data[trainset, :, :]), p)
        c, d = _construct_var_eqns(atleast_3d(data[testset, :, :]), p)

        e = sp.linalg.inv(np.eye(a.shape[0]) * delta ** 2 + a.dot(a.T))

        cc = c.transpose().dot(c)

        be = b.transpose().dot(e)
        bee = be.dot(e)
        bea = be.dot(a)
        beea = bee.dot(a)
        beacc = bea.dot(cc)
        dc = d.transpose().dot(c)

        j += np.sum(beacc * bea - 2 * bea * dc) + np.sum(d ** 2)
        k += np.sum(beea * dc - beacc * beea) * 4 * delta

    return j / (nt * d.size), k / (nt * d.size)
878478
def _protected_inverse(x1):
    """"""
    with np.errstate(divide='ignore', invalid='ignore'):
        return np.where(np.abs(x1) > 0.001, 1. / x1, 0.)
1248388
function norminv(p) {

  function rateval(a, b, x) {
    var u = a[a.length - 1];
    for (var i = a.length - 1; i > 0; --i) {
      u = x * u + a[i - 1];
    }
    var v = b[b.length - 1];
    for (var j = b.length - 1; j > 0; --j) {
      v = x * v + b[j - 1];
    }
    r = u / v;
    return r;
  };

  var a = [
    3.3871328727963666080,     1.3314166789178437745e+2,
    1.9715909503065514427e+3,  1.3731693765509461125e+4,
    4.5921953931549871457e+4,  6.7265770927008700853e+4,
    3.3430575583588128105e+4,  2.5090809287301226727e+3];
  var b = [
    1.0,                       4.2313330701600911252e+1,
    6.8718700749205790830e+2,  5.3941960214247511077e+3,
    2.1213794301586595867e+4,  3.9307895800092710610e+4,
    2.8729085735721942674e+4,  5.2264952788528545610e+3];
  var c = [
    1.42343711074968357734,     4.63033784615654529590,
    5.76949722146069140550,     3.64784832476320460504,
    1.27045825245236838258,     2.41780725177450611770e-1,
    2.27238449892691845833e-2,  7.74545014278341407640e-4];
  var d = [
    1.0,                        2.05319162663775882187,
    1.67638483018380384940,     6.89767334985100004550e-1,
    1.48103976427480074590e-1,  1.51986665636164571966e-2,
    5.47593808499534494600e-4,  1.05075007164441684324e-9];
  var e = [
    6.65790464350110377720,     5.46378491116411436990,
    1.78482653991729133580,     2.96560571828504891230e-1,
    2.65321895265761230930e-2,  1.24266094738807843860e-3,
    2.71155556874348757815e-5,  2.01033439929228813265e-7];
  var f = [
    1.0,                        5.99832206555887937690e-1,
    1.36929880922735805310e-1,  1.48753612908506148525e-2,
    7.86869131145613259100e-4,  1.84631831751005468180e-5,
    1.42151175831644588870e-7,  2.04426310338993978564e-15];

  if (p <= 0.0) {
    return Number.NEGATIVE_INFINITY;
  }
  if (p >= 1.0) {
    return Number.POSITIVE_INFINITY;
  }

  var q = p - 0.5;
  if (Math.abs(q) <= 0.425) {
    return q * rateval(a, b, 0.180625 - q*q);
  }

  var r = Math.sqrt(-Math.log(q < 0.0 ? p : 1.0 - p));
  var x;
  if (r <= 5.0) {
    x = rateval(c, d, r - 1.6);
  } else {
    x = rateval(e, f, r - 0.5);
  }
  if (p < 0.5) {
    return -x;
  }
  return x;
}
848550
def GaP_Eg_Gamma(self, **kwargs):
    '''
    
    '''
    T = kwargs.get('T', 300.)
    if T < 1e-4:
        return self.Eg_Gamma_0()
    return self.Eg_Gamma_0() + 0.1081 * (1 - 1. / tanh(164. / T))  # eV
113064
function pmf( x, p ) {
	var q;
	if ( x !== x) {
		return NaN;
	}
	if ( isNonNegativeInteger( x ) ) {
		q = 1 - p;
		return p * pow( q, x );
	}
	return 0;
}
88943
function partial( lambda ) {

	/**
	* FUNCTION: cdf( x )
	*	Evaluates the cumulative distribution function (CDF) for a Poisson distribution.
	*
	* @private
	* @param {Number} x - input value
	* @returns {Number} evaluated CDF
	*/
	return function cdf( x ) {
		if ( x < 0 ) {
			return 0;
		}
		return gammainc( lambda, floor( x ) + 1, {
			'tail': 'upper',
			'regularized': true
		});
	};
}
1316581
def reg_incomplete_beta(a, b, x):
  

  if (x == 0):
    return 0
  elif (x == 1):
    return 1
  else:
    lbeta = (math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) +
             a * math.log(x) + b * math.log(1 - x))
    if (x < (a + 1) / (a + b + 2)):
      return math.exp(lbeta) * __contfractbeta(a, b, x) / a
    else:
      return 1 - math.exp(lbeta) * __contfractbeta(b, a, 1 - x) / b
780009
private void applyL1Reg(final double eta_t, Vec x)
    {
        //apply l1 regularization
        if(lambda1 > 0)
        {
            l1U += eta_t*lambda1;//line 6: in Tsuruoka et al paper, figure 2
            for(int k = 0; k < ws.length; k++)
            {
                final Vec w_k = ws[k];
                final double[] l1Q_k = l1Q[k];
                for(IndexValue iv : x)
                {
                    final int i = iv.getIndex();
                    //see "APPLYPENALTY(i)" on line 15: from Figure 2 in Tsuruoka et al paper
                    final double z = w_k.get(i);
                    double newW_i = 0;
                    if (z > 0)
                        newW_i = Math.max(0, z - (l1U + l1Q_k[i]));
                    else if(z < 0)
                        newW_i = Math.min(0, z + (l1U - l1Q_k[i]));
                    l1Q_k[i] += (newW_i - z);
                    w_k.set(i, newW_i);
                }
            }
        }
    }
1363597
def forward(Q, p, G, h, A, b, Q_LU, S_LU, R, eps=1e-12, verbose=0, notImprovedLim=3,
            maxIter=20, solver=KKTSolvers.LU_PARTIAL):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    # Find initial values
    if solver == KKTSolvers.LU_FULL:
        D = torch.eye(nineq).repeat(nBatch, 1, 1).type_as(Q)
        x, s, z, y = factor_solve_kkt(
            Q, D, G, A, p,
            torch.zeros(nBatch, nineq).type_as(Q),
            -h, -b if b is not None else None)
    elif solver == KKTSolvers.LU_PARTIAL:
        d = torch.ones(nBatch, nineq).type_as(Q)
        factor_kkt(S_LU, R, d)
        x, s, z, y = solve_kkt(
            Q_LU, d, G, A, S_LU,
            p, torch.zeros(nBatch, nineq).type_as(Q),
            -h, -b if neq > 0 else None)
    elif solver == KKTSolvers.IR_UNOPT:
        D = torch.eye(nineq).repeat(nBatch, 1, 1).type_as(Q)
        x, s, z, y = solve_kkt_ir(
            Q, D, G, A, p,
            torch.zeros(nBatch, nineq).type_as(Q),
            -h, -b if b is not None else None)
    else:
        assert False

    # Make all of the slack variables >= 1.
    M = torch.min(s, 1)[0]
    M = M.view(M.size(0), 1).repeat(1, nineq)
    I = M < 0
    s[I] -= M[I] - 1

    # Make all of the inequality dual variables >= 1.
    M = torch.min(z, 1)[0]
    M = M.view(M.size(0), 1).repeat(1, nineq)
    I = M < 0
    z[I] -= M[I] - 1

    best = {'resids': None, 'x': None, 'z': None, 's': None, 'y': None}
    nNotImproved = 0

    for i in range(maxIter):
        # affine scaling direction
        rx = (torch.bmm(y.unsqueeze(1), A).squeeze(1) if neq > 0 else 0.) + \
            torch.bmm(z.unsqueeze(1), G).squeeze(1) + \
            torch.bmm(x.unsqueeze(1), Q.transpose(1, 2)).squeeze(1) + \
            p
        rs = z
        rz = torch.bmm(x.unsqueeze(1), G.transpose(1, 2)).squeeze(1) + s - h
        ry = torch.bmm(x.unsqueeze(1), A.transpose(
            1, 2)).squeeze(1) - b if neq > 0 else 0.0
        mu = torch.abs((s * z).sum(1).squeeze() / nineq)
        z_resid = torch.norm(rz, 2, 1).squeeze()
        y_resid = torch.norm(ry, 2, 1).squeeze() if neq > 0 else 0
        pri_resid = y_resid + z_resid
        dual_resid = torch.norm(rx, 2, 1).squeeze()
        resids = pri_resid + dual_resid + nineq * mu

        d = z / s
        try:
            factor_kkt(S_LU, R, d)
        except:
            return best['x'], best['y'], best['z'], best['s']

        if verbose == 1:
            print('iter: {}, pri_resid: {:.5e}, dual_resid: {:.5e}, mu: {:.5e}'.format(
                i, pri_resid.mean(), dual_resid.mean(), mu.mean()))
        if best['resids'] is None:
            best['resids'] = resids
            best['x'] = x.clone()
            best['z'] = z.clone()
            best['s'] = s.clone()
            best['y'] = y.clone() if y is not None else None
            nNotImproved = 0
        else:
            I = resids < best['resids']
            if I.sum() > 0:
                nNotImproved = 0
            else:
                nNotImproved += 1
            I_nz = I.repeat(nz, 1).t()
            I_nineq = I.repeat(nineq, 1).t()
            best['resids'][I] = resids[I]
            best['x'][I_nz] = x[I_nz]
            best['z'][I_nineq] = z[I_nineq]
            best['s'][I_nineq] = s[I_nineq]
            if neq > 0:
                I_neq = I.repeat(neq, 1).t()
                best['y'][I_neq] = y[I_neq]
        if nNotImproved == notImprovedLim or best['resids'].max() < eps or mu.min() > 1e32:
            if best['resids'].max() > 1. and verbose >= 0:
                print(INACC_ERR)
            return best['x'], best['y'], best['z'], best['s']

        if solver == KKTSolvers.LU_FULL:
            D = bdiag(d)
            dx_aff, ds_aff, dz_aff, dy_aff = factor_solve_kkt(
                Q, D, G, A, rx, rs, rz, ry)
        elif solver == KKTSolvers.LU_PARTIAL:
            dx_aff, ds_aff, dz_aff, dy_aff = solve_kkt(
                Q_LU, d, G, A, S_LU, rx, rs, rz, ry)
        elif solver == KKTSolvers.IR_UNOPT:
            D = bdiag(d)
            dx_aff, ds_aff, dz_aff, dy_aff = solve_kkt_ir(
                Q, D, G, A, rx, rs, rz, ry)
        else:
            assert False

        # compute centering directions
        alpha = torch.min(torch.min(get_step(z, dz_aff),
                                    get_step(s, ds_aff)),
                          torch.ones(nBatch).type_as(Q))
        alpha_nineq = alpha.repeat(nineq, 1).t()
        t1 = s + alpha_nineq * ds_aff
        t2 = z + alpha_nineq * dz_aff
        t3 = torch.sum(t1 * t2, 1).squeeze()
        t4 = torch.sum(s * z, 1).squeeze()
        sig = (t3 / t4)**3

        rx = torch.zeros(nBatch, nz).type_as(Q)
        rs = ((-mu * sig).repeat(nineq, 1).t() + ds_aff * dz_aff) / s
        rz = torch.zeros(nBatch, nineq).type_as(Q)
        ry = torch.zeros(nBatch, neq).type_as(Q) if neq > 0 else torch.Tensor()

        if solver == KKTSolvers.LU_FULL:
            D = bdiag(d)
            dx_cor, ds_cor, dz_cor, dy_cor = factor_solve_kkt(
                Q, D, G, A, rx, rs, rz, ry)
        elif solver == KKTSolvers.LU_PARTIAL:
            dx_cor, ds_cor, dz_cor, dy_cor = solve_kkt(
                Q_LU, d, G, A, S_LU, rx, rs, rz, ry)
        elif solver == KKTSolvers.IR_UNOPT:
            D = bdiag(d)
            dx_cor, ds_cor, dz_cor, dy_cor = solve_kkt_ir(
                Q, D, G, A, rx, rs, rz, ry)
        else:
            assert False

        dx = dx_aff + dx_cor
        ds = ds_aff + ds_cor
        dz = dz_aff + dz_cor
        dy = dy_aff + dy_cor if neq > 0 else None
        alpha = torch.min(0.999 * torch.min(get_step(z, dz),
                                            get_step(s, ds)),
                          torch.ones(nBatch).type_as(Q))
        alpha_nineq = alpha.repeat(nineq, 1).t()
        alpha_neq = alpha.repeat(neq, 1).t() if neq > 0 else None
        alpha_nz = alpha.repeat(nz, 1).t()

        x += alpha_nz * dx
        s += alpha_nineq * ds
        z += alpha_nineq * dz
        y = y + alpha_neq * dy if neq > 0 else None

    if best['resids'].max() > 1. and verbose >= 0:
        print(INACC_ERR)
    return best['x'], best['y'], best['z'], best['s']
397481
def dBinaryRochedx (r, D, q, F):
    
    return -r[0]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*(r[0]-D)*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5 -q/D/D + F*F*(1+q)*r[0]
352337
private static function _logBeta($p, $q) {
		if ($p != self::$_logBetaCache_p || $q != self::$_logBetaCache_q) {
			self::$_logBetaCache_p = $p;
			self::$_logBetaCache_q = $q;
			if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > LOG_GAMMA_X_MAX_VALUE)) {
				self::$_logBetaCache_result = 0.0;
			} else {
				self::$_logBetaCache_result = self::_logGamma($p) + self::_logGamma($q) - self::_logGamma($p + $q);
			}
		}
		return self::$_logBetaCache_result;
	}
199118
def solve_kkt(Q_LU, d, G, A, S_LU, rx, rs, rz, ry):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    invQ_rx = rx.btrisolve(*Q_LU)
    if neq > 0:
        h = torch.cat((invQ_rx.unsqueeze(1).bmm(A.transpose(1, 2)).squeeze(1) - ry,
                       invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz), 1)
    else:
        h = invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz

    w = -(h.btrisolve(*S_LU))

    g1 = -rx - w[:, neq:].unsqueeze(1).bmm(G).squeeze(1)
    if neq > 0:
        g1 -= w[:, :neq].unsqueeze(1).bmm(A).squeeze(1)
    g2 = -rs - w[:, neq:]

    dx = g1.btrisolve(*Q_LU)
    ds = g2 / d
    dz = w[:, neq:]
    dy = w[:, :neq] if neq > 0 else None

    return dx, ds, dz, dy
757602
function getSqSegDist(p, p1, p2) {

    var x = p1.x,
        y = p1.y,
        dx = p2.x - x,
        dy = p2.y - y;

    if (dx !== 0 || dy !== 0) {

        var t = ((p.x - x) * dx + (p.y - y) * dy) / (dx * dx + dy * dy);

        if (t > 1) {
            x = p2.x;
            y = p2.y;

        } else if (t > 0) {
            x += dx * t;
            y += dy * t;
        }
    }

    dx = p.x - x;
    dy = p.y - y;

    return dx * dx + dy * dy;
}