Incomple gamma function \gamma; computed from NumRec routine gammp.
    >>> np.isclose(gamm_inc(0.5,1),1.49365)
    True
    >>> np.isclose(gamm_inc(1.5,2),0.6545103)
    True
    >>> np.isclose(gamm_inc(2.5,1e-12),0)
    True
961229
def gamm_inc(a,x):
    
    assert (x > 0 and a >= 0), "Invalid arguments in routine gamm_inc: %s,%s" % (x,a)

    if x < (a+1.0): #Use the series representation
        gam,gln = _gser(a,x)
    else: #Use continued fractions
        gamc,gln = _gcf(a,x)
        gam = 1-gamc
    return np.exp(gln)*gam
15179
def projR(gamma, p):
    
    return np.multiply(gamma.T, p / np.maximum(np.sum(gamma, axis=1), 1e-10)).T
40423
function invIncBeta(p, a, b) {
  if(x < 0 || x > 1) {
    throw new RangeError('First argument must be between 0 and 1.');
	}

  // Special cases, there can make trouble otherwise
  else if (a === 1 && b === 1) return p;
  else if (p === 1) return 1;
  else if (p === 0) return 0;
  else if (a === 0) return 0;
  else if (b === 0) return 1;

  else {
    var EPS = 1e-8,
        a1 = a - 1,
        b1 = b - 1,
        j = 0,
        lna, lnb, pp, t, u, err, x, al, h, w, afac;

	if(a >= 1 && b >= 1) {
    pp = (p < 0.5) ? p : 1 - p;
    t = Math.sqrt(-2 * Math.log(pp));

		x = (2.30753 + t * 0.27061) / (1 + t* (0.99229 + t * 0.04481)) - t;
		if(p < 0.5) x = -x;
		al = (x * x - 3) / 6;
		h = 2 / (1 / (2 * a - 1)  + 1 / (2 * b - 1));
		w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) * (al + 5 / 6 - 2 / (3 * h));
		x = a / (a + b * Math.exp(2 * w));
	} else {
		lna = Math.log(a / (a + b));
		lnb = Math.log(b / (a + b));
		t = Math.exp(a * lna) / a;
		u = Math.exp(b * lnb) / b;
		w = t + u;
		if (p < t / w) x = Math.pow(a * w * p, 1 / a);
		else x = 1 - Math.pow(b * w * (1 - p), 1 / b);
	}

	afac = -gammaCollection.logGamma(a) - gammaCollection.logGamma(b) + gammaCollection.logGamma(a + b);

  for(; j < 10; j++) {
		if(x === 0 || x === 1) return x;
		err = regualizedBeta(x, a, b) - p;

    t = Math.exp(a1 * Math.log(x) + b1 * log1p(-x) + afac);
		u = err / t;
		x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x)))));

    if (x <= 0) x = 0.5 * (x + t);
		if (x >= 1) x = 0.5 * (x + t + 1);

		if (Math.abs(t) < EPS * x && j > 0) break;
	}

	return x;
  }
}
392283
def funcGauss1D(x, mu, sig):
    
    arrOut = np.exp(-np.power((x - mu)/sig, 2.)/2)
    # normalize
#    arrOut = arrOut/(np.sqrt(2.*np.pi)*sig)
    # normalize (laternative)
    arrOut = arrOut/np.sum(arrOut)
    return arrOut
790638
def _crps_cdf_single(x, cdf_or_dist, xmin=None, xmax=None, tol=1e-6):
    
    # TODO: this function is pretty slow.  Look for clever ways to speed it up.

    # allow for directly passing in scipy.stats distribution objects.
    cdf = getattr(cdf_or_dist, 'cdf', cdf_or_dist)
    assert callable(cdf)

    # if bounds aren't given, discover them
    if xmin is None or xmax is None:
        # Note that infinite values for xmin and xmax are valid, but
        # it slows down the resulting quadrature significantly.
        xmin, xmax = _discover_bounds(cdf)

    # make sure the bounds haven't clipped the cdf.
    if (tol is not None) and (cdf(xmin) >= tol) or (cdf(xmax) <= (1. - tol)):
        raise ValueError('CDF does not meet tolerance requirements at %s '
                         'extreme(s)! Consider using function defaults '
                         'or using infinities at the bounds. '
                         % ('lower' if cdf(xmin) >= tol else 'upper'))

    # CRPS = int_-inf^inf (F(y) - H(x))**2 dy
    #      = int_-inf^x F(y)**2 dy + int_x^inf (1 - F(y))**2 dy
    def lhs(y):
        # left hand side of CRPS integral
        return np.square(cdf(y))
    # use quadrature to integrate the lhs
    lhs_int, lhs_tol = integrate.quad(lhs, xmin, x)
    # make sure the resulting CRPS will be with tolerance
    if (tol is not None) and (lhs_tol >= 0.5 * tol):
        raise ValueError('Lower integral did not evaluate to within tolerance! '
                         'Tolerance achieved: %f , Value of integral: %f \n'
                         'Consider setting the lower bound to -np.inf.' %
                         (lhs_tol, lhs_int))

    def rhs(y):
        # right hand side of CRPS integral
        return np.square(1. - cdf(y))
    rhs_int, rhs_tol = integrate.quad(rhs, x, xmax)
    # make sure the resulting CRPS will be with tolerance
    if (tol is not None) and (rhs_tol >= 0.5 * tol):
        raise ValueError('Upper integral did not evaluate to within tolerance! \n'
                         'Tolerance achieved: %f , Value of integral: %f \n'
                         'Consider setting the upper bound to np.inf or if '
                         'you already have, set warn_level to `ignore`.' %
                         (rhs_tol, rhs_int))

    return lhs_int + rhs_int
281096
function gamma(x) {
      var expx;

      expx = Math.exp(lgamma(x));

      return x > 0 && x === Math.round(x) ? Math.round(expx) : expx;
   }
731752
public static final double gammaln(double x) {
		double ser = 1.00000000090015;
		double y = x;
		double tmp = x + 5.5;
		tmp -= (x + 0.5) * Math.log(tmp);
		for (int j = 0; j < 5; j++) {
			ser += COFGAMMALN[j] / ++y;
		}
		return -tmp + Math.log(ser * 2.5066282751072975 / x);
	}
765130
private function gamma($x)
    {
        $gamma = 0.577215664901532860606512090; # Euler's gamma constant
        if ($x < 0.001) {
            return 1.0/($x*(1.0 + $gamma*$x));
        }
        if ($x < 12.0) {
            # The algorithm directly approximates gamma over (1,2) and uses
            # reduction identities to reduce other arguments to this interval.
            $y = $x;
            $n = 0;
            $arg_was_less_than_one = ($y < 1.0);
            # Add or subtract integers as necessary to bring y into (1,2)
            # Will correct for this below
            if ($arg_was_less_than_one) {
                $y += 1.0;
            } else {
                $n = floor($y) - 1;  # will use n later
                $y -= $n;
            }
            # numerator coefficients for approximation over the interval (1,2)
            $p =
            array(
                -1.71618513886549492533811E+0,
                 2.47656508055759199108314E+1,
                -3.79804256470945635097577E+2,
                 6.29331155312818442661052E+2,
                 8.66966202790413211295064E+2,
                -3.14512729688483675254357E+4,
                -3.61444134186911729807069E+4,
                 6.64561438202405440627855E+4
            );

            # denominator coefficients for approximation over the interval (1,2)
            $q =
            array(
                -3.08402300119738975254353E+1,
                 3.15350626979604161529144E+2,
                -1.01515636749021914166146E+3,
                -3.10777167157231109440444E+3,
                 2.25381184209801510330112E+4,
                 4.75584627752788110767815E+3,
                -1.34659959864969306392456E+5,
                -1.15132259675553483497211E+5
            );

            $num = 0.0;
            $den = 1.0;

            $z = $y - 1;
            for ($i = 0; $i < 8; $i++) {
                $num = ($num + $p[$i])*$z;
                $den = $den*$z + $q[$i];
            }
            $result = $num/$den + 1.0;

            # Apply correction if argument was not initially in (1,2)
            if ($arg_was_less_than_one) {
                # Use identity gamma(z) = gamma(z+1)/z
                # The variable "result" now holds gamma of the original y + 1
                # Thus we use y-1 to get back the orginal y.
                $result /= ($y-1.0);
            } else {
                # Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
                for ($i = 0; $i < $n; $i++) {
                    $result *= $y++;
                }
            }

            return $result;
        }

        ###########################################################################
        # Third interval: [12, infinity)

        if ($x > 171.624) {
            # Correct answer too large to display.

            return Double.POSITIVE_INFINITY;
        }

        return exp($this->log_gamma($x));
    }
98839
def obfn_reg(self):
        

        rl1 = np.linalg.norm((self.wl1 * self.obfn_gvar()).ravel(), 1)
        rl2 = 0.5*np.linalg.norm(self.obfn_gvar())**2
        return (self.lmbda*rl1 + self.mu*rl2, rl1, rl2)
1084462
def funcGauss1D(x, mu, sig):
    
    arrOut = np.exp(-np.power((x - mu)/sig, 2.)/2)
    # normalize
    arrOut = arrOut/(np.sqrt(2.*np.pi)*sig)

    return arrOut
1096009
def F(self, x):
        
        if isinstance(x, np.ndarray):
            nfwvals = np.ones_like(x)
            inds1 = np.where(x < 1)
            inds2 = np.where(x > 1)
            nfwvals[inds1] = (1 - x[inds1] ** 2) ** -.5 * np.arctanh((1 - x[inds1] ** 2) ** .5)
            nfwvals[inds2] = (x[inds2] ** 2 - 1) ** -.5 * np.arctan((x[inds2] ** 2 - 1) ** .5)
            return nfwvals

        elif isinstance(x, float) or isinstance(x, int):
            if x == 1:
                return 1
            if x < 1:
                return (1 - x ** 2) ** -.5 * np.arctanh((1 - x ** 2) ** .5)
            else:
                return (x ** 2 - 1) ** -.5 * np.arctan((x ** 2 - 1) ** .5)
1167800
def _param_bounds(self, gamma, q):
        
        if gamma < 1.4:
            gamma = 1.4
        if gamma > 2.9:
            gamma = 2.9
        if q < 0.01:
            q = 0.01
        return float(gamma), q
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
1017459
def has_approx_support(m, m_hat, prob=0.01):
    
    m_nz = np.flatnonzero(np.triu(m, 1))
    m_hat_nz = np.flatnonzero(np.triu(m_hat, 1))

    upper_diagonal_mask = np.flatnonzero(np.triu(np.ones(m.shape), 1))
    not_m_nz = np.setdiff1d(upper_diagonal_mask, m_nz)

    intersection = np.in1d(m_hat_nz, m_nz)  # true positives
    not_intersection = np.in1d(m_hat_nz, not_m_nz)  # false positives

    true_positive_rate = 0.0
    if len(m_nz):
        true_positive_rate = 1. * np.sum(intersection) / len(m_nz)
        true_negative_rate = 1. - true_positive_rate

    false_positive_rate = 0.0
    if len(not_m_nz):
        false_positive_rate = 1. * np.sum(not_intersection) / len(not_m_nz)

    return int(np.less_equal(true_negative_rate + false_positive_rate, prob))
322023
private static double regularizedIncompleteGammaSeries(double a, double x) {
        if (a < 0.0 || x < 0.0 || x >= a + 1) {
            throw new IllegalArgumentException(String.format("Invalid a = %f, x = %f", a, x));
        }

        int i = 0;
        double igf = 0.0;
        boolean check = true;

        double acopy = a;
        double sum = 1.0 / a;
        double incr = sum;
        double loggamma = lgamma(a);

        while (check) {
            ++i;
            ++a;
            incr *= x / a;
            sum += incr;
            if (Math.abs(incr) < Math.abs(sum) * INCOMPLETE_GAMMA_EPSILON) {
                igf = sum * Math.exp(-x + acopy * Math.log(x) - loggamma);
                check = false;
            }
            if (i >= INCOMPLETE_GAMMA_MAX_ITERATIONS) {
                check = false;
                igf = sum * Math.exp(-x + acopy * Math.log(x) - loggamma);
                logger.error("Gamma.regularizedIncompleteGammaSeries: Maximum number of iterations wes exceeded");
            }
        }
        return igf;
    }
469687
def _gser(a,x):
    ""
    ITMAX=100
    EPS=3.e-7

    gln=lgamma(a)
    assert(x>=0),'x < 0 in gser'
    if x == 0 : return 0,gln

    ap = a
    delt = sum = 1./a
    for i in range(ITMAX):
        ap=ap+1.
        delt=delt*x/ap
        sum=sum+delt
        if abs(delt) < abs(sum)*EPS: break
    else:
        print('a too large, ITMAX too small in gser')
    gamser=sum*np.exp(-x+a*np.log(x)-gln)
    return gamser,gln
571880
def from_EV(E, V):
        
        a = np.square(E) / V
        b = E / V
        return Gamma(a, b)