Compute the mean of the 3 opbends
556852
def _opbend_transform_mean(rs, fn_low, deriv=0):
    
    v = 0.0
    d = np.zeros((4,3), float)
    dd = np.zeros((4,3,4,3), float)
    #loop over the 3 cyclic permutations
    for p in np.array([[0,1,2], [2,0,1], [1,2,0]]):
        opbend = _opbend_transform([rs[p[0]], rs[p[1]], rs[p[2]], rs[3]], fn_low, deriv)
        v += opbend[0]/3
        index0 = np.where(p==0)[0][0] #index0 is the index of the 0th atom (rs[0])
        index1 = np.where(p==1)[0][0]
        index2 = np.where(p==2)[0][0]
        index3 = 3
        if deriv>0:
            d[0] += opbend[1][index0]/3
            d[1] += opbend[1][index1]/3
            d[2] += opbend[1][index2]/3
            d[3] += opbend[1][index3]/3
        if deriv>1:
            dd[0, :, 0, :] += opbend[2][index0, :, index0, :]/3
            dd[0, :, 1, :] += opbend[2][index0, :, index1, :]/3
            dd[0, :, 2, :] += opbend[2][index0, :, index2, :]/3
            dd[0, :, 3, :] += opbend[2][index0, :, index3, :]/3

            dd[1, :, 0, :] += opbend[2][index1, :, index0, :]/3
            dd[1, :, 1, :] += opbend[2][index1, :, index1, :]/3
            dd[1, :, 2, :] += opbend[2][index1, :, index2, :]/3
            dd[1, :, 3, :] += opbend[2][index1, :, index3, :]/3

            dd[2, :, 0, :] += opbend[2][index2, :, index0, :]/3
            dd[2, :, 1, :] += opbend[2][index2, :, index1, :]/3
            dd[2, :, 2, :] += opbend[2][index2, :, index2, :]/3
            dd[2, :, 3, :] += opbend[2][index2, :, index3, :]/3

            dd[3, :, 0, :] += opbend[2][index3, :, index0, :]/3
            dd[3, :, 1, :] += opbend[2][index3, :, index1, :]/3
            dd[3, :, 2, :] += opbend[2][index3, :, index2, :]/3
            dd[3, :, 3, :] += opbend[2][index3, :, index3, :]/3
    if deriv==0:
        return v,
    elif deriv==1:
        return v, d
    elif deriv==2:
        return v, d, dd
    else:
        raise ValueError("deriv must be 0, 1 or 2.")
874121
def acq(self, graph):
        ''' 
        '''
        mean, std = self.gpr.predict(np.array([graph.extract_descriptor()]))
        if self.optimizemode is OptimizeMode.Maximize:
            return mean + self.beta * std
        return mean - self.beta * std
547341
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        mean, stddevs = super().get_mean_and_stddevs(
            sites, rup, dists, imt, stddev_types)

        idx_rock = sites.vs30 >= self.ROCK_VS30
        idx_soil = sites.vs30 < self.ROCK_VS30

        mean[idx_rock] += 0.275
        mean[idx_soil] += 0.31

        return mean, stddevs
129479
private void computeDandH( DenseMatrix64F param , DenseMatrix64F x , DenseMatrix64F y )
    {
        func.compute(param,x, tempDH);
        subtractEquals(tempDH, y);

        if (jacobianFactory != null)
        {
            jacobianFactory.computeJacobian(param, x, jacobian);
        }
        else
        {
            computeNumericalJacobian(param,x,jacobian);
        }

        int numParam = param.getNumElements();
        int length = y.getNumElements();

        // d = average{ (f(x_i;p) - y_i) * jacobian(:,i) }
        for( int i = 0; i < numParam; i++ ) {
            double total = 0;
            for( int j = 0; j < length; j++ ) {
                total += tempDH.get(j,0)*jacobian.get(i,j);
            }
            d.set(i,0,total/length);
        }

        // compute the approximation of the hessian
        multTransB(jacobian,jacobian,H);
        scale(1.0/length,H);
    }
1146156
def _compute_mean_on_rock(self, C, mag, rrup, rvol, hypo_depth, CN, CR,
                              f4HW):
        

        lnSA_AB = (
            # line 1 of equation 1
            C['c1'] + C['c4as'] * (mag - 6) +
            # line 2
            C['c3as'] * (8.5 - mag) ** 2 +
            # line 3
            C['c5'] * rrup +
            # line 3 and 4
            (C['c8'] + C['c6as'] * (mag - 6)) *
            np.log((rrup ** 2 + C['c10as'] ** 2) ** 0.5) +
            # line 5
            C['c46'] * rvol +
            # line 6
            C['c32'] * CN + C['c33as'] * CR + f4HW
        )

        return lnSA_AB
450498
function mean( out, x ) {
	var len = x.length,
		i;
	if ( out.length !== len ) {
		throw new Error( 'mean()::invalid input arguments. Input and output matrices must be the same length.' );
	}
	for ( i = 0; i < len; i++ ) {
		out.data[ i ] = MEAN( x.data[ i ] );
	}
	return out;
}
493866
def dot_v3(v, w):
    """"""

    return sum([x * y for x, y in zip(v, w)])
670889
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.
        C = self.COEFFS[imt]
        # Implements mean model (equation 12)
        mean = (self._compute_magnitude(rup, C) +
                self._compute_distance(dists, C) +
                self._get_site_amplification(sites, rup, C) +
                self._get_mechanism(rup, C))

        stddevs = self._get_stddevs(rup, np.exp(mean), stddev_types, sites)

        return mean, stddevs
1367907
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extract dictionaries of coefficients specific to required
        # intensity measure type
        C = self.COEFFS[imt]

        mean = self._compute_mean(C, rup, dists, sites, imt)

        stddevs = self._get_stddevs(C, stddev_types, sites.vs30.shape[0])

        return mean, stddevs
1021
def _compute_mean(self, C, mag, rjb):
        
        # line 1686 in hazgridXnga2.f
        ffc = self._compute_finite_fault_correction(mag)
        d = np.sqrt(rjb ** 2 + (C['c7'] ** 2) * (ffc ** 2))

        # lines 1663, 1694-1696 in hazgridXnga2.f
        mean = (
            C['c1'] + C['c2'] * (mag - 6.) +
            C['c3'] * ((mag - 6.) ** 2) -
            C['c4'] * np.log(d) - C['c6'] * d
        )

        factor = np.log(rjb / 100.)
        idx = factor > 0
        mean[idx] -= (C['c5'] - C['c4']) * factor[idx]

        return mean
1343512
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        C = self.COEFFS[imt]
        mean = (
            self._get_magnitude_term(C, rup.mag) +
            self._get_distance_term(C, rup.mag, dists.rrup) +
            self._get_style_of_faulting_term(C, rup.rake) +
            self._get_site_response_term(C, sites.vs30))

        stddevs = self._get_stddevs(C, stddev_types, len(sites.vs30))
        return mean, stddevs
1261814
function ( obj ) {
		var nth = obj.length,
		    n   = 0,
		    mean;

		if ( nth > 0 ) {
			mean = array.mean( obj );

			array.each( obj, function ( i ) {
				n += math.sqr( i - mean );
			} );

			return n / nth;
		}
		else {
			return n;
		}
	}
567051
def _get_warped_mean(self, mean, std, pred_init=None, deg_gauss_hermite=20):
        
        gh_samples, gh_weights = np.polynomial.hermite.hermgauss(deg_gauss_hermite)
        gh_samples = gh_samples[:, None]
        gh_weights = gh_weights[None, :]
        return gh_weights.dot(self._get_warped_term(mean, std, gh_samples)) / np.sqrt(np.pi)
499694
def _expectation(p, mean, none1, none2, none3, nghp=None):
    
    return mean(p.mu)
936233
def _compute_mean(self, C, rup, dists, sites, imt):
        
        mean = (self._get_magnitude_scaling_term(C, rup.mag) +
                self._get_distance_scaling_term(C, rup.mag, dists.rrup) +
                self._get_site_amplification_term(C, sites.vs30))
        # convert from cm/s**2 to g for SA and from m/s**2 to g for PGA (PGV
        # is already in cm/s) and also convert from base 10 to base e.
        if imt.name == "PGA":
            mean = np.log((10 ** mean) * ((2 * np.pi / 0.01) ** 2) *
                          1e-2 / g)
        elif imt.name == "SA":
            mean = np.log((10 ** mean) * ((2 * np.pi / imt.period) ** 2) *
                          1e-2 / g)
        else:
            mean = np.log(10 ** mean)

        return mean
816928
def _compute_mean(self, C, mag, rjb):
        
        d1 = np.sqrt(50. ** 2 + 6. ** 2)
        d = np.sqrt(rjb ** 2 + 6 ** 2)

        mean = np.zeros_like(rjb)

        mean += (
            C['a1'] + C['a2'] * (mag - 6.4) +
            C['a7'] * (8.5 - mag) ** 2
        )

        idx = rjb < 50.
        mean[idx] += (
            C['a3'] * np.log(d[idx]) +
            C['a4'] * (mag - 6.4) * np.log(d[idx]) +
            C['a5'] * rjb[idx]
        )

        idx = rjb >=50.
        mean[idx] += (
            C['a3'] * np.log(d1) +
            C['a4'] * (mag - 6.4) * np.log(d[idx]) +
            C['a5'] * rjb[idx] + C['a6'] * (np.log(d[idx]) - np.log(d1))
        )

        return mean
1064950
def _stieltjes_analytical(dist, order, normed):
    """"""
    dimensions = len(dist)
    mom_order = numpy.arange(order+1).repeat(dimensions)
    mom_order = mom_order.reshape(order+1, dimensions).T
    coeff1, coeff2 = dist.ttr(mom_order)
    coeff2[:, 0] = 1.

    poly = chaospy.poly.collection.core.variable(dimensions)
    if normed:
        orth = [
            poly**0*numpy.ones(dimensions),
            (poly-coeff1[:, 0])/numpy.sqrt(coeff2[:, 1]),
        ]
        for order_ in range(1, order):
            orth.append(
                (orth[-1]*(poly-coeff1[:, order_])
                 -orth[-2]*numpy.sqrt(coeff2[:, order_]))
                /numpy.sqrt(coeff2[:, order_+1])
            )
        norms = numpy.ones(coeff2.shape)
    else:
        orth = [poly-poly, poly**0*numpy.ones(dimensions)]
        for order_ in range(order):
            orth.append(
                orth[-1]*(poly-coeff1[:, order_])
                - orth[-2]*coeff2[:, order_]
            )
        orth = orth[1:]
        norms = numpy.cumprod(coeff2, 1)

    return orth, norms, coeff1, coeff2