Here are the examples of the python api numpy.multiply taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
115 Examples
3
Example 1
def _eval(self, v):
"""Evaluate the function on v (ignoring parameters).
"""
# Square
vsum = v.copy()
np.multiply(v, v, vsum)
# Sum along dimensions and keep dimensions
for d in self.group_dims:
vsum = np.sum(vsum, axis=d, keepdims=True)
# Sqrt
np.sqrt(vsum, vsum)
# L1 norm is then sum of norms
return np.sum(vsum)
3
Example 2
def __array_wrap__(self, obj, context=None):
obj = super(ComplicatedSubArray, self).__array_wrap__(obj, context)
if context is not None and context[0] is np.multiply:
obj.info['multiplied'] = obj.info.get('multiplied', 0) + 1
return obj
3
Example 3
Project: iris Source File: test_Cube__operators.py
def test_lazy_biggus_mul_scalar(self):
c1 = self.build_lazy_cube([1, 2])
cube = c1 * 5
self.assertEqual(c1 * 5, 5 * c1)
result = cube.lazy_data()
self.assertTrue(cube.has_lazy_data())
self.assert_elementwise(c1, None, result, np.multiply)
3
Example 4
Project: automl_gpu Source File: libscores.py
def acc_stat (solution, prediction):
''' Return accuracy statistics TN, FP, TP, FN
Assumes that solution and prediction are binary 0/1 vectors.'''
# This uses floats so the results are floats
TN = sum(np.multiply((1-solution), (1-prediction)))
FN = sum(np.multiply(solution, (1-prediction)))
TP = sum(np.multiply(solution, prediction))
FP = sum(np.multiply((1-solution), prediction))
#print "TN =",TN
#print "FP =",FP
#print "TP =",TP
#print "FN =",FN
return (TN, FP, TP, FN)
3
Example 5
def times(A, b, offset=0):
"""
Times the view of A with b in place (!).
Returns modified A
Broadcasting is allowed, thus b can be scalar.
if offset is not zero, make sure b is of right shape!
:param ndarray A: 2 dimensional array
:param ndarray-like b: either one dimensional or scalar
:param int offset: same as in view.
:rtype: view of A, which is adjusted inplace
"""
return _diag_ufunc(A, b, offset, np.multiply)
3
Example 6
Project: pylon Source File: ac_pf.py
def _evaluate_function(self, Ybus, V, Sbus, pv, pq):
""" Evaluates F(x).
"""
mis = multiply(V, conj(Ybus * V)) - Sbus
F = r_[mis[pv].real, mis[pq].real, mis[pq].imag]
return F
3
Example 7
@property
def bin_sizes(self):
sizes1 = 0.5 * (self.get_bin_right_edges(0) ** 2 - self.get_bin_left_edges(0) ** 2)
sizes2 = self.get_bin_widths(1)
sizes3 = self.get_bin_widths(2)
return reduce(np.multiply, np.ix_(sizes1, sizes2, sizes3))
3
Example 8
Project: pylon Source File: ac_pf.py
def _evaluate_mismatch(self, Ybus, V, Sbus, pq, pvpq):
""" Evaluates the mismatch.
"""
mis = (multiply(V, conj(Ybus * V)) - Sbus) / abs(V)
P = mis[pvpq].real
Q = mis[pq].imag
return P, Q
3
Example 9
Project: openpathsampling Source File: pes.py
def dVdx(self, sys):
"""Derivative of potential energy (-force)
Parameters
----------
sys : :class:`.ToyEngine`
engine contains its state, including velocities and masses
Returns
-------
np.array
the derivatives of the potential at this point
"""
dx = sys.positions - self.x0
exp_part = self.A*np.exp(-np.dot(self.alpha, np.multiply(dx, dx)))
for i in range(len(dx)):
self._local_dVdx[i] = -2*self.alpha[i]*dx[i]*exp_part
return self._local_dVdx
3
Example 10
Project: RLScore Source File: sqerror_measure.py
def sqerror_multitask(Y, Y_predicted):
Y = np.mat(Y)
Y_predicted = np.mat(Y_predicted)
sqerror = Y - Y_predicted
multiply(sqerror, sqerror, sqerror)
performances = mean(sqerror, axis = 0)
performances = np.array(performances)[0]
return performances
3
Example 11
Project: deep_nets_iclr04 Source File: tools.py
def blend_mask(im, color, mask):
mask = mask * (mask > 80) # 120
for i in range(0, 3):
im[:, :, i] = np.multiply((255 - mask) / 255.0, im[:, :, i])
for i in range(0, 3):
im[:, :, i] = im[:, :, i] + np.multiply((mask / 255.0), color[i])
return im
3
Example 12
def _prep_data(self):
if self.basemap:
self.prep_data_for_basemap()
else:
self.plot_data = self.data[0]
if self.mult_factor[0]:
self.plot_data = np.multiply(float(self.mult_factor[0]),
self.plot_data)
3
Example 13
@Operation.factory
def Complex(a, b):
"""
Complex number op.
"""
return np.add(a, np.multiply(b, 1j)),
3
Example 14
Project: ldsc Source File: regressions.py
def _prop(self, jknife, M, Nbar, cat, tot):
'''Convert total h2 and per-category h2 to per-category proportion h2 or gencov.'''
n_annot = self.n_annot
n_blocks = jknife.delete_values.shape[0]
numer_delete_vals = np.multiply(
M, jknife.delete_values[:, 0:n_annot]) / Nbar # (n_blocks, n_annot)
denom_delete_vals = np.sum(
numer_delete_vals, axis=1).reshape((n_blocks, 1))
denom_delete_vals = np.dot(denom_delete_vals, np.ones((1, n_annot)))
prop = jk.RatioJackknife(
cat / tot, numer_delete_vals, denom_delete_vals)
return prop.est, prop.jknife_cov, prop.jknife_se
3
Example 15
Project: nimfa Source File: linalg.py
def multiply(X, Y):
"""
Compute element-wise multiplication of matrices :param:`X` and :param:`Y`.
:param X: First input matrix.
:type X: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil, dia
or :class:`numpy.matrix`
:param Y: Second input matrix.
:type Y: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil, dia
or :class:`numpy.matrix`
"""
if sp.isspmatrix(X) and sp.isspmatrix(Y):
return X.multiply(Y)
elif sp.isspmatrix(X) or sp.isspmatrix(Y):
return _op_spmatrix(X, Y, np.multiply)
else:
with warnings.catch_warnings():
warnings.simplefilter('ignore')
return np.multiply(np.mat(X), np.mat(Y))
3
Example 16
def computeGradients(self, pixelBlock, props):
# pixel size in input raster SR...
p = props['cellSize'] if self.sr is None else computeCellSize(props, self.sr, self.proj)
if p is not None and len(p) == 2:
p = np.multiply(p, 1.11e5 if isGeographic(props['spatialReference']) else 1.) # conditional degrees to meters conversion
xs, ys = (self.zf + (np.power(p, self.ce) * self.cf)) / (8*p)
else:
xs, ys = 1., 1. # degenerate case. shouldn't happen.
return (ndimage.convolve(pixelBlock, self.xKernel)*xs, ndimage.convolve(pixelBlock, self.yKernel)*ys)
3
Example 17
Project: iris Source File: test_Cube__operators.py
def test_lazy_biggus_mul_cubes(self):
c1 = self.build_lazy_cube([1, 2])
cube = c1 * c1
result = cube.lazy_data()
self.assertTrue(cube.has_lazy_data())
self.assert_elementwise(c1, c1.lazy_data(), result, np.multiply)
3
Example 18
Project: AST-text-analysis Source File: relevance.py
def relevance(self, keyphrase, text, synonimizer=None):
# Based on: https://janav.wordpress.com/2013/10/27/tf-idf-and-cosine-similarity/,
# but query vectors are defined here in the same vector space as docuement vectors
# (not in the reduced one as in the article).
# TF-IDF for query tokens
query_tokens = self._preprocess_tokens([utils.tokenize_and_filter(keyphrase)])
query_tf, query_idf = self._tf_idf(query_tokens)
query_tf = query_tf[0]
# Weighting for both text and query (either TF or TF-IDF)
if self.term_weighting == consts.TermWeighting.TF:
text_vector = self.tf[text]
query_vector = query_tf
elif self.term_weighting == consts.TermWeighting.TF_IDF:
text_vector = np.multiply(self.tf[text], self.idf)
query_vector = np.multiply(query_tf, query_idf)
return self._cosine_similarity(text_vector, query_vector)
3
Example 19
Project: polar2grid Source File: rescale.py
def linear_scale(img, m, b, **kwargs):
LOG.debug("Running 'linear_scale' with (m: %f, b: %f)..." % (m, b))
if m != 1:
numpy.multiply(img, m, img)
if b != 0:
numpy.add(img, b, img)
return img
3
Example 20
def V(self, sys):
"""Potential energy
Parameters
----------
sys : :class:`.ToyEngine`
engine contains its state, including velocities and masses
Returns
-------
float
the potential energy
"""
dx = sys.positions - self.x0
return self.A*np.exp(-np.dot(self.alpha, np.multiply(dx, dx)))
3
Example 21
def update_gradients_full(self, dL_dK, X, X2=None):
if len(self.parts)==2:
self.parts[0].update_gradients_full(dL_dK*self.parts[1].K(X,X2), X, X2)
self.parts[1].update_gradients_full(dL_dK*self.parts[0].K(X,X2), X, X2)
else:
for combination in itertools.combinations(self.parts, len(self.parts) - 1):
prod = reduce(np.multiply, [p.K(X, X2) for p in combination])
to_update = list(set(self.parts) - set(combination))[0]
to_update.update_gradients_full(dL_dK * prod, X, X2)
3
Example 22
def _update_inernal(self, varp_list):
"""Make an update of the internal status by gathering the variational posteriors for all the individual models."""
# The probability for the binary variable for the same latent dimension of any of the models is on.
if self.group_spike:
self._b_prob_all = 1.-param_to_array(varp_list[0].gamma_group)
[np.multiply(self._b_prob_all, 1.-vp.gamma_group, self._b_prob_all) for vp in varp_list[1:]]
else:
self._b_prob_all = 1.-param_to_array(varp_list[0].binary_prob)
[np.multiply(self._b_prob_all, 1.-vp.binary_prob, self._b_prob_all) for vp in varp_list[1:]]
3
Example 23
def kinetic_energy(self, sys):
"""Default kinetic energy implementation.
Parameters
----------
sys : :class:`.ToyEngine`
engine contains its state, including velocities and masses
"""
v = sys.velocities
m = sys.mass
return 0.5*np.dot(m, np.multiply(v,v))
3
Example 24
Project: pyvolve Source File: evolver.py
def _exponentiate_matrix(self, Q, t):
'''
Perform exponentiation on instantaneous matrix to produce transition matrix, from a given instantaneous matrix Q and a given branch length t.
Assert that all rows sum to 1.
Return P
'''
P = linalg.expm( np.multiply(Q, float(t) ) )
assert( np.allclose( np.sum(P, axis = 1), np.ones(len(self._code))) ), "Rows in transition matrix do not each sum to 1."
return P
3
Example 25
Project: RLScore Source File: linalg.py
def svd_economy_sized(X):
"""Returns the reduced singular value decomposition of the data matrix X so that only the singular vectors corresponding to the nonzero singular values are returned.
@param X: data matrix whose rows and columns correspond to the data and features, respectively.
@type X: numpy matrix of floats
@return: the nonzero singular values and the corresponding left and right singular vectors of X. The singular vectors are contained in a r*1-matrix, where r is the number of nonzero singular values.
@rtype: a tuple of three numpy matrices"""
evecs, svals, U = la.svd(X, full_matrices=0)
evals = np.multiply(svals, svals)
evecs = evecs[:, svals > 0]
svals = svals[svals > 0]
U = U[svals > 0]
return evecs, svals, U
3
Example 26
Project: statsmodels Source File: kernels.py
def smooth(self, xs, ys, x):
"""Returns the kernel smoothing estimate for point x based on x-values
xs and y-values ys.
Not expected to be called by the user.
Special implementation optimised for Gaussian.
"""
w = np.sum(exp(multiply(square(divide(subtract(xs, x),
self.h)),-0.5)))
v = np.sum(multiply(ys, exp(multiply(square(divide(subtract(xs, x),
self.h)), -0.5))))
return v/w
3
Example 27
Project: statsmodels Source File: kernels.py
def smooth(self, xs, ys, x):
"""Returns the kernel smoothing estimate for point x based on x-values
xs and y-values ys.
Not expected to be called by the user.
Special implementation optimised for Biweight.
"""
xs, ys = self.in_domain(xs, ys, x)
if len(xs) > 0:
w = np.sum(square(subtract(1, square(divide(subtract(xs, x),
self.h)))))
v = np.sum(multiply(ys, square(subtract(1, square(divide(
subtract(xs, x), self.h))))))
return v / w
else:
return np.nan
3
Example 28
@Operation.factory
def Mul(a, b):
"""
Multiplication op.
"""
return np.multiply(a, b),
3
Example 29
Project: TensorLog Source File: mutil.py
def multiplyByBroadcastRowVec(m,v):
(dm,dv,i) = codensify(m,v)
if dm!=None:
dp = NP.multiply(dm,dv)
return undensify(dp, i)
else:
bv = repeat(v, numRows(m))
return m.multiply(bv)
3
Example 30
def callback(self, learner):
print
print 'LOOCV mean squared error', learner.bestlooperf
print 'The indices of selected features', learner.selected
if not self.test_features is None:
mod = learner.predictor
tpreds = mod.predict(self.test_features)
if not self.test_measure is None:
test_perf = self.test_measure(self.test_labels, tpreds)
else:
testdiff = self.test_labels - tpreds
test_perf = np.mean(np.multiply(testdiff, testdiff))
print 'Test performance', test_perf
3
Example 31
def _cat(self, jknife, M, Nbar, coef, coef_cov):
'''Convert coefficients to per-category h2 or gencov.'''
cat = np.multiply(M, coef)
cat_cov = np.multiply(np.dot(M.T, M), coef_cov)
cat_se = np.sqrt(np.diag(cat_cov))
return cat, cat_cov, cat_se
3
Example 32
Project: tfdeploy Source File: tfdeploy.py
@Operation.factory
def SegmentProd(a, ids):
"""
Segmented prod op.
"""
func = lambda idxs: reduce(np.multiply, a[idxs])
return seg_map(func, a, ids),
3
Example 33
def forward(self, inputs, outputs):
"""The forward operator.
Reads from inputs and writes to outputs.
"""
np.multiply(inputs[0], self.scalar, outputs[0])
3
Example 34
Project: broca Source File: wikipedia.py
def compute_bridge_similarity(self, vec1, vec2):
EWP = 1 - np.multiply(vec1, vec2)
# not sure exactly how to sort the EWP vector
#EWP = sorted(EWP, reverse=True)
EWP = sorted(EWP, reverse=True)
k = 10
EWP = EWP[:k]
# The paper does not mention using logs but start to get into underflow
# issues multiplying so many decimal values
lEWP = -1 * np.log(EWP)
return 1/np.sum(lEWP)
3
Example 35
@Cache_this(limit=3, force_kwargs=['which_parts'])
def K(self, X, X2=None, which_parts=None):
if which_parts is None:
which_parts = self.parts
elif not isinstance(which_parts, (list, tuple)):
# if only one part is given
which_parts = [which_parts]
return reduce(np.multiply, (p.K(X, X2) for p in which_parts))
3
Example 36
def update_gradients_diag(self, dL_dKdiag, X):
if len(self.parts)==2:
self.parts[0].update_gradients_diag(dL_dKdiag*self.parts[1].Kdiag(X), X)
self.parts[1].update_gradients_diag(dL_dKdiag*self.parts[0].Kdiag(X), X)
else:
for combination in itertools.combinations(self.parts, len(self.parts) - 1):
prod = reduce(np.multiply, [p.Kdiag(X) for p in combination])
to_update = list(set(self.parts) - set(combination))[0]
to_update.update_gradients_diag(dL_dKdiag * prod, X)
3
Example 37
def prepMessages(self):
""" Multiplies together incoming messages to make new outgoing
"""
# compute new messages if no observation has been made
if self.enabled and self.observed < 0 and len(self.nbrs) > 1:
# switch reference for old messages
self.nextStep()
for i in range(0, len(self.incoming)):
# multiply together all excluding message at current index
curr = self.incoming[:]
del curr[i]
self.outgoing[i] = reduce(np.multiply, curr)
# normalize once finished with all messages
self.normalizeMessages()
3
Example 38
Project: statsmodels Source File: kernels.py
def smoothvar(self, xs, ys, x):
"""
Returns the kernel smoothing estimate of the variance at point x.
"""
xs, ys = self.in_domain(xs, ys, x)
if len(xs) > 0:
fittedvals = np.array([self.smooth(xs, ys, xx) for xx in xs])
rs = square(subtract(ys, fittedvals))
w = np.sum(square(subtract(1.0, square(divide(subtract(xs, x),
self.h)))))
v = np.sum(multiply(rs, square(subtract(1, square(divide(
subtract(xs, x), self.h))))))
return v / w
else:
return np.nan
2
Example 39
Project: CrisisMappingToolkit Source File: histogram.py
def __computeJT(relHist, binVals, T):
'''As part of the Kittler/Illingworth method, compute J(T) for a given T'''
FAIL_VAL = 999999
# Split the hostogram at the threshold T.
histogram1 = relHist[0:T]
histogram2 = relHist[T+1:]
# Compute the number of pixels in the two classes.
P1 = sum(histogram1)
P2 = sum(histogram2)
# Only continue if both classes contain at least one pixel.
if (P1 <= 0) or (P2 <= 0):
return FAIL_VAL
# Compute the standard deviations of the classes.
weightedBins1 = numpy.multiply(histogram1, binVals[0:T])
weightedBins2 = numpy.multiply(histogram2, binVals[T+1:])
mean1 = sum(weightedBins1) / P1;
mean2 = sum(weightedBins2) / P2;
diffs1 = numpy.subtract(binVals[0:T], mean1)
diffs2 = numpy.subtract(binVals[T+1:], mean2)
diffsSq1 = [d*d for d in diffs1]
diffsSq2 = [d*d for d in diffs2]
weightedBins1 = numpy.multiply(histogram1, diffsSq1)
weightedBins2 = numpy.multiply(histogram2, diffsSq2)
sigma1 = math.sqrt(sum(weightedBins1) / P1)
sigma2 = math.sqrt(sum(weightedBins2) / P2)
# Make sure both classes contain at least two intensity values.
if (sigma1 <= 0) or (sigma2 <= 0):
return FAIL_VAL
# Compute J(T).
J = 1 + 2*(P1*math.log(sigma1) + P2*math.log(sigma2)) - 2*(P1*math.log(P1) + P2*math.log(P2))
return J
0
Example 40
Project: RLScore Source File: global_rankrls.py
def holdout(self, indices):
"""Computes hold-out predictions for a trained RankRLS
Parameters
----------
indices : list of indices, shape = [n_hsamples]
list of indices of training examples belonging to the set for which the
hold-out predictions are calculated. The list can not be empty.
Returns
-------
F : array, shape = [n_hsamples, n_labels]
holdout predictions
Notes
-----
The algorithm is a modification of the ones published in [1,2] for the regular RLS method.
References
----------
[1] Tapio Pahikkala, Jorma Boberg, and Tapio Salakoski.
Fast n-Fold Cross-Validation for Regularized Least-Squares.
Proceedings of the Ninth Scandinavian Conference on Artificial Intelligence,
83-90, Otamedia Oy, 2006.
[2] Tapio Pahikkala, Hanna Suominen, and Jorma Boberg.
Efficient cross-validation for kernelized least-squares regression with sparse basis expansions.
Machine Learning, 87(3):381--407, June 2012.
"""
indices = array_tools.as_index_list(indices, self.Y.shape[0])
if len(indices) != len(np.unique(indices)):
raise IndexError('Hold-out can have each index only once.')
Y = self.Y
m = self.size
evals, V = self.evals, self.svecs
#results = []
C = np.mat(np.zeros((self.size, 3), dtype = np.float64))
onevec = np.mat(np.ones((self.size, 1), dtype = np.float64))
for i in range(self.size):
C[i, 0] = 1.
VTY = V.T * Y
VTC = V.T * onevec
CTY = onevec.T * Y
holen = len(indices)
newevals = multiply(evals, 1. / ((m - holen) * evals + self.regparam))
R = np.mat(np.zeros((holen, holen + 1), dtype = np.float64))
for i in range(len(indices)):
R[i, 0] = -1.
R[i, i + 1] = sqrt(self.size - float(holen))
Vho = V[indices]
Vhov = multiply(Vho, newevals)
Ghoho = Vhov * Vho.T
GCho = Vhov * VTC
GBho = Ghoho * R
for i in range(len(indices)):
GBho[i, 0] += GCho[i, 0]
CTGC = multiply(VTC.T, newevals) * VTC
RTGCho = R.T * GCho
BTGB = R.T * Ghoho * R
for i in range(len(indices) + 1):
BTGB[i, 0] += RTGCho[i, 0]
BTGB[0, i] += RTGCho[i, 0]
BTGB[0, 0] += CTGC[0, 0]
BTY = R.T * Y[indices]
BTY[0] = BTY[0] + CTY[0]
GDYho = Vhov * (self.size - holen) * VTY
GLYho = GDYho - GBho * BTY
CTGDY = multiply(VTC.T, newevals) * (self.size - holen) * VTY
BTGLY = R.T * GDYho - BTGB * BTY
BTGLY[0] = BTGLY[0] + CTGDY[0]
F = GLYho - GBho * la.inv(-mat(eye(holen + 1)) + BTGB) * BTGLY
#results.append(F)
#return results
F = np.squeeze(np.array(F))
return F
0
Example 41
Project: RoBO Source File: stochastic_local_search.py
def maximize(self):
"""
Maximizes the given acquisition function.
Returns
-------
np.ndarray(N,D)
Point with highest acquisition value.
"""
if hasattr(self.objective_func, "_get_most_probable_minimum"):
xx = self.objective_func._get_most_probable_minimum()
else:
xx = np.add(
np.multiply(
(self.X_lower - self.X_upper),
self.rng.uniform(
size=(
1,
self.X_lower.shape[0]))),
self.X_lower)
def fun_p(x):
acq_v = self.objective_func(np.array([x]))[0]
log_acq_v = np.log(acq_v) if acq_v > 0 else -np.inf
return log_acq_v
sc_fun = self._scipy_optimizer_fkt_wrapper(self.objective_func, False)
S0 = 0.5 * np.linalg.norm(self.X_upper - self.X_lower)
D = self.X_lower.shape[0]
Xstart = np.zeros((self.Ne, D))
restarts = np.zeros((self.Ne, D))
if self.starts is None and hasattr(self.objective_func, "BestGuesses"):
self.starts = self.objective_func.BestGuesses
if self.starts is not None and self.Ne > self.starts.shape[0]:
restarts[self.starts.shape[0]:self.Ne, ] = self.X_lower + \
(self.X_upper - self.X_lower) * self.rng.uniform(
size=(self.Ne - self.starts.shape[0], D))
elif self.starts is not None:
restarts[0:self.Ne] = self.starts[0:self.Ne]
else:
restarts = self.X_lower + (self.X_upper - self.X_lower) * \
self.rng.uniform(size=(self.Ne, D))
sampler = emcee.EnsembleSampler(self.Ne, D, fun_p)
Xstart, logYstart, _ = sampler.run_mcmc(restarts, 20)
search_cons = []
for i in range(0, self.X_lower.shape[0]):
xmin = self.X_lower[i]
xmax = self.X_upper[i]
search_cons.append({'type': 'ineq',
'fun': lambda x: x - xmin})
search_cons.append({'type': 'ineq',
'fun': lambda x: xmax - x})
search_cons = tuple(search_cons)
minima = []
jacobian = False
i = 0
while i < self.Ne:
# try:
minima.append(scipy.optimize.minimize(fun=sc_fun,
x0=Xstart[i, np.newaxis],
jac=jacobian,
method='L-BFGS-B',
constraints=search_cons,
options={'ftol': np.spacing(1),
'maxiter': 20}))
i += 1
# no derivatives
# except BayesianOptimizationError, e:
# if e.errno == BayesianOptimizationError.NO_DERIVATIVE:
# jacobian = False
# sc_fun = self._scipy_optimizer_fkt_wrapper(self.objective_func, False)
# else:
# raise e
# X points:
Xend = np.array([res.x for res in minima])
# Objective function values:
Xdh = np.array([res.fun for res in minima])
new_x = Xend[np.nanargmin(Xdh)]
if len(new_x.shape):
new_x = np.array([new_x])
return new_x
0
Example 42
Project: RLScore Source File: global_rankrls.py
def _reference(self, pairs):
evals, evecs = self.evals, self.svecs
Y = self.Y
m = self.size
results = []
D = mat(zeros((self.size, 1), dtype=float64))
C = mat(zeros((self.size, 3), dtype=float64))
for i in range(self.size):
D[i, 0] = self.size - 2.
C[i, 0] = 1.
lamb = self.regparam
G = evecs * multiply(multiply(evals, 1. / ((m - 2.) * evals + lamb)).T, evecs.T)
DY = multiply(D, Y)
GDY = G * DY
GC = G * C
CTG = GC.T
CT = C.T
CTGC = CT * GC
CTY = CT * Y
CTGDY = CT * GDY
minusI3 = -mat(eye(3))
for i, j in pairs:
hoinds = [i, j]
R = mat(zeros((2, 3), dtype=float64))
R[0, 0] = -1.
R[1, 0] = -1.
R[0, 1] = sqrt(self.size - 2.)
R[1, 2] = sqrt(self.size - 2.)
GBho = GC[hoinds] + G[np.ix_(hoinds, hoinds)] * R
BTGB = CTGC \
+ R.T * GC[hoinds] \
+ CTG[:, hoinds] * R \
+ R.T * G[np.ix_(hoinds, hoinds)] * R
BTY = CTY + R.T * Y[hoinds]
GLYho = GDY[hoinds] - GBho * BTY
BTGLY = CTGDY + R.T * GDY[hoinds] - BTGB * BTY
F = GLYho - GBho * la.inv(minusI3 + BTGB) * BTGLY
results.append(F)
return results
0
Example 43
Project: auto-sklearn Source File: classification_metrics.py
def bac_metric(solution, prediction, task=BINARY_CLASSIFICATION):
"""
Compute the normalized balanced accuracy.
The binarization and
the normalization differ for the multi-label and multi-class case.
:param solution:
:param prediction:
:param task:
:return:
"""
if task == BINARY_CLASSIFICATION:
if len(solution.shape) == 1:
# Solution won't be touched - no copy
solution = solution.reshape((-1, 1))
elif len(solution.shape) == 2:
if solution.shape[1] > 1:
raise ValueError('Solution array must only contain one class '
'label, but contains %d' % solution.shape[1])
else:
raise ValueError('Solution.shape %s' % solution.shape)
if len(prediction.shape) == 2:
if prediction.shape[1] > 2:
raise ValueError('A prediction array with probability values '
'for %d classes is not a binary '
'classification problem' % prediction.shape[1])
# Prediction will be copied into a new binary array - no copy
prediction = prediction[:, 1].reshape((-1, 1))
else:
raise ValueError('Invalid prediction shape %s' % prediction.shape)
elif task == MULTICLASS_CLASSIFICATION:
if len(solution.shape) == 1:
solution = create_multiclass_solution(solution, prediction)
elif len(solution.shape) == 2:
if solution.shape[1] > 1:
raise ValueError('Solution array must only contain one class '
'label, but contains %d' % solution.shape[1])
else:
solution = create_multiclass_solution(solution.reshape((-1, 1)),
prediction)
else:
raise ValueError('Solution.shape %s' % solution.shape)
elif task == MULTILABEL_CLASSIFICATION:
pass
else:
raise NotImplementedError('bac_metric does not support task type %s'
% task)
bin_prediction = binarize_predictions(prediction, task)
fn = np.sum(np.multiply(solution, (1 - bin_prediction)), axis=0,
dtype=float)
tp = np.sum(np.multiply(solution, bin_prediction), axis=0, dtype=float)
# Bounding to avoid division by 0
eps = 1e-15
tp = sp.maximum(eps, tp)
pos_num = sp.maximum(eps, tp + fn)
tpr = tp / pos_num # true positive rate (sensitivity)
if task in (BINARY_CLASSIFICATION, MULTILABEL_CLASSIFICATION):
tn = np.sum(np.multiply((1 - solution), (1 - bin_prediction)),
axis=0, dtype=float)
fp = np.sum(np.multiply((1 - solution), bin_prediction), axis=0,
dtype=float)
tn = sp.maximum(eps, tn)
neg_num = sp.maximum(eps, tn + fp)
tnr = tn / neg_num # true negative rate (specificity)
bac = 0.5 * (tpr + tnr)
base_bac = 0.5 # random predictions for binary case
elif task == MULTICLASS_CLASSIFICATION:
label_num = solution.shape[1]
bac = tpr
base_bac = 1. / label_num # random predictions for multiclass case
bac = np.mean(bac) # average over all classes
# Normalize: 0 for random, 1 for perfect
score = (bac - base_bac) / sp.maximum(eps, (1 - base_bac))
return score
0
Example 44
Project: RLScore Source File: global_rankrls.py
def solve(self, regparam=1.0):
"""Re-trains RankRLS for the given regparam
Parameters
----------
regparam : float, optional
regularization parameter, regparam > 0 (default=1.0)
Notes
-----
Computational complexity of re-training:
m = n_samples, d = n_features, l = n_labels, b = n_bvectors
O(lm^2): basic case
O(lmd): Linear Kernel, d < m
O(lmb): Sparse approximation with basis vectors
"""
if not hasattr(self, "multiplyright"):
#Eigenvalues of the kernel matrix
self.evals = multiply(self.svals, self.svals)
#Temporary variables
ssvecs = multiply(self.svecs, self.svals)
J = mat(ones((self.size, 1), dtype=float64))
#These are cached for later use in solve function
self.lowrankchange = ssvecs.T * J[range(ssvecs.shape[0])]
self.multipleright = ssvecs.T * (self.size * self.Y - J * (J.T * self.Y))
self.regparam = regparam
#Compute the eigenvalues determined by the given regularization parameter
neweigvals = 1. / (self.size * self.evals + regparam)
P = self.lowrankchange
nP = multiply(neweigvals.T, P)
fastinv = 1. / (-1. + P.T * nP)
self.A = self.svecs * multiply(1. / self.svals.T, \
(multiply(neweigvals.T, self.multipleright) \
- nP * (fastinv * (nP.T * self.multipleright))))
self.predictor = self.svdad.createModel(self)
0
Example 45
def train(self):
#Cached results
self.evals = np.multiply(self.svals, self.svals)
self.newevals = 1. / (self.evals + self.regparam)
newevalslamtilde = np.multiply(self.evals, self.newevals)
self.D = np.sqrt(newevalslamtilde)
#self.D = -newevalslamtilde
self.VTY = self.svecs.T * self.Y
self.DVTY = np.multiply(self.D.T, self.svecs.T * self.Y)
self.sqrtR = np.multiply(np.sqrt(newevalslamtilde), self.svecs)
#self.R = self.svecs * multiply(newevalslamtilde.T, self.svecs.T)
'''
#Global variation
self.R = self.sqrtR * self.sqrtR.T
self.minus_diagRx2 = - 2 * np.diag(self.R)
'''
#Space efficient variation
self.R = None
self.minus_diagRx2 = - 2 * np.array(np.sum(np.multiply(self.sqrtR, self.sqrtR), axis = 1)).reshape((self.size))
#'''
self.RY = self.sqrtR * (self.sqrtR.T * self.Y)
self.Y_Schur_RY = np.multiply(self.Y, self.RY)
self.classFitnessList = []
for i in range(self.labelcount):
#DVTY_i = DVTY[:,i]
#self.DVTY_list.append(DVTY_i)
YTRY_i = self.Y[:,i].T * self.RY[:,i]
#self.YTRY_list.append(YTRY_i)
fitness_i = self.size - YTRY_i
self.classFitnessList.append(fitness_i[0, 0])
self.classFitnessRowVec = np.array(self.classFitnessList)
if not self.callbackfun is None:
self.callbackfun.callback(self)
0
Example 46
def leave_pair_out(self, pairs_start_inds, pairs_end_inds):
"""Computes leave-pair-out predictions for a trained RankRLS.
Parameters
----------
pairs_start_inds : list of indices, shape = [n_pairs]
list of indices from range [0, n_samples-1]
pairs_end_inds : list of indices, shape = [n_pairs]
list of indices from range [0, n_samples-1]
Returns
-------
P1 : array, shape = [n_pairs]
holdout predictions for pairs_start_inds
P2 : array, shape = [n_pairs]
holdout predictions for pairs_end_inds
Notes
-----
Computes the leave-pair-out cross-validation predictions, where each (i,j) pair with
i= pair_start_inds[k] and j = pairs_end_inds[k] is left out in turn.
When estimating area under ROC curve with leave-pair-out, one should leave out all
positive-negative pairs, while for estimating the general ranking error one should
leave out all pairs with different labels.
Computational complexity of leave-pair-out with most pairs left out:
m = n_samples, d = n_features, l = n_labels, b = n_bvectors
O(lm^2+m^3): basic case
O(lm^2+dm^2): Linear Kernel, d < m
O(lm^2+bm^2): Sparse approximation with basis vectors
The leave-pair-out cross-validation algorithm is described in [1,2]. The use of
leave-pair-out cross-validation for AUC estimation has been analyzed in [3]
[1] Tapio Pahikkala, Evgeni Tsivtsivadze, Antti Airola, Jouni Jarvinen, and Jorma Boberg.
An efficient algorithm for learning to rank from preference graphs.
Machine Learning, 75(1):129-165, 2009.
[2] Tapio Pahikkala, Antti Airola, Jorma Boberg, and Tapio Salakoski.
Exact and efficient leave-pair-out cross-validation for ranking RLS.
In Proceedings of the 2nd International and Interdisciplinary Conference
on Adaptive Knowledge Representation and Reasoning (AKRR'08), pages 1-8,
Espoo, Finland, 2008.
[3] Antti Airola, Tapio Pahikkala, Willem Waegeman, Bernard De Baets, Tapio Salakoski.
An Experimental Comparison of Cross-Validation Techniques for Estimating the Area Under the ROC Curve.
Computational Statistics & Data Analysis 55(4), 1828-1844, 2011.
"""
pairs_start_inds = array_tools.as_index_list(pairs_start_inds, self.Y.shape[0])
pairs_end_inds = array_tools.as_index_list(pairs_end_inds, self.Y.shape[0])
evals, svecs = self.evals, self.svecs
m = self.size
Y = self.Y
modevals = np.squeeze(np.array(np.multiply(evals, 1. / ((m - 2.) * evals + self.regparam))))
GDY = (self.size - 2.) * (svecs * np.multiply(np.mat(modevals).T, (svecs.T * Y)))
GC = np.squeeze(np.array(svecs * np.multiply(np.mat(modevals).T, np.sum(svecs.T, axis = 1))))
CTGC = np.sum(GC)
pairslen = len(pairs_start_inds)
sm2Gdiag = np.zeros((self.Y.shape[0]))
BTY = np.zeros((self.Y.shape))
sqrtsm2GDY = np.zeros((self.Y.shape))
BTGBBTY = np.zeros((self.Y.shape))
results_first = np.zeros((pairslen, self.Y.shape[1]))
results_second = np.zeros((pairslen, self.Y.shape[1]))
_global_rankrls.leave_pair_out(pairslen,
self.Y.shape[0],
pairs_start_inds,
pairs_end_inds,
self.Y.shape[1],
Y,
svecs,
modevals,
svecs.shape[1],
np.zeros((self.Y.shape[0])),
np.squeeze(np.array(GC)),
sm2Gdiag,
CTGC,
GDY,
BTY,
sqrtsm2GDY,
BTGBBTY,
np.array(np.sum(Y, axis = 0))[0], #CTY
np.array(np.sum(GDY, axis = 0))[0], #CTGDY
results_first,
results_second)
return np.squeeze(results_first), np.squeeze(results_second)
0
Example 47
def updateA(self):
self.A = self.svecs * np.multiply(self.newevals.T, self.VTY)
0
Example 48
Project: auto-sklearn Source File: classification_metrics.py
def acc_metric(solution, prediction, task=BINARY_CLASSIFICATION):
"""
Compute the accuracy.
Get the accuracy stats
acc = (tpr + fpr) / (tn + fp + tp + fn)
Normalize, so 1 is the best and zero mean random...
:param solution:
:param prediction:
:param task:
:return:
"""
if task == BINARY_CLASSIFICATION:
if len(solution.shape) == 1:
# Solution won't be touched - no copy
solution = solution.reshape((-1, 1))
elif len(solution.shape) == 2:
if solution.shape[1] > 1:
raise ValueError('Solution array must only contain one class '
'label, but contains %d' % solution.shape[1])
else:
raise ValueError('Solution.shape %s' % solution.shape)
if len(prediction.shape) == 2:
if prediction.shape[1] > 2:
raise ValueError('A prediction array with probability values '
'for %d classes is not a binary '
'classification problem' % prediction.shape[1])
# Prediction will be copied into a new binary array - no copy
prediction = prediction[:, 1].reshape((-1, 1))
else:
raise ValueError('Invalid prediction shape %s' % prediction.shape)
elif task == MULTICLASS_CLASSIFICATION:
if len(solution.shape) == 1:
solution = create_multiclass_solution(solution, prediction)
elif len(solution.shape ) == 2:
if solution.shape[1] > 1:
raise ValueError('Solution array must only contain one class '
'label, but contains %d' % solution.shape[1])
else:
solution = create_multiclass_solution(solution.reshape((-1, 1)),
prediction)
else:
raise ValueError('Solution.shape %s' % solution.shape)
elif task == MULTILABEL_CLASSIFICATION:
pass
else:
raise NotImplementedError('acc_metric does not support task type %s'
% task)
bin_predictions = binarize_predictions(prediction, task)
tn = np.sum(np.multiply((1 - solution), (1 - bin_predictions)), axis=0,
dtype=float)
fn = np.sum(np.multiply(solution, (1 - bin_predictions)), axis=0,
dtype=float)
tp = np.sum(np.multiply(solution, bin_predictions), axis=0,
dtype=float)
fp = np.sum(np.multiply((1 - solution), bin_predictions), axis=0,
dtype=float)
# Bounding to avoid division by 0, 1e-7 because of float32
eps = np.float(1e-7)
"""
tp = np.sum(tp)
fp = np.sum(fp)
tn = np.sum(tn)
fn = np.sum(fn)
"""
if task in (BINARY_CLASSIFICATION, MULTILABEL_CLASSIFICATION):
accuracy = (np.sum(tp) + np.sum(tn)) / (
np.sum(tp) + np.sum(fp) + np.sum(tn) + np.sum(fn)
)
elif task == MULTICLASS_CLASSIFICATION:
accuracy = np.sum(tp) / (np.sum(tp) + np.sum(fp))
if task in (BINARY_CLASSIFICATION, MULTILABEL_CLASSIFICATION):
base_accuracy = 0.5 # random predictions for binary case
elif task == MULTICLASS_CLASSIFICATION:
label_num = solution.shape[1]
base_accuracy = 1. / label_num
# Normalize: 0 for random, 1 for perfect
score = (accuracy - base_accuracy) / sp.maximum(eps, (1 - base_accuracy))
return score
0
Example 49
Project: auto-sklearn Source File: classification_metrics.py
def f1_metric(solution, prediction, task=BINARY_CLASSIFICATION):
"""
Compute the normalized f1 measure.
The binarization differs
for the multi-label and multi-class case.
A non-weighted average over classes is taken.
The score is normalized.
:param solution:
:param prediction:
:param task:
:return:
"""
if task == BINARY_CLASSIFICATION:
if len(solution.shape) == 1:
# Solution won't be touched - no copy
solution = solution.reshape((-1, 1))
elif len(solution.shape) == 2:
if solution.shape[1] > 1:
raise ValueError('Solution array must only contain one class '
'label, but contains %d' % solution.shape[1])
else:
raise ValueError('Solution.shape %s' % solution.shape)
if len(prediction.shape) == 2:
if prediction.shape[1] > 2:
raise ValueError('A prediction array with probability values '
'for %d classes is not a binary '
'classification problem' % prediction.shape[1])
# Prediction will be copied into a new binary array - no copy
prediction = prediction[:, 1].reshape((-1, 1))
else:
raise ValueError('Invalid prediction shape %s' % prediction.shape)
elif task == MULTICLASS_CLASSIFICATION:
if len(solution.shape) == 1:
solution = create_multiclass_solution(solution, prediction)
elif len(solution.shape) == 2:
if solution.shape[1] > 1:
raise ValueError('Solution array must only contain one class '
'label, but contains %d' % solution.shape[1])
else:
solution = create_multiclass_solution(solution.reshape((-1, 1)),
prediction)
else:
raise ValueError('Solution.shape %s' % solution.shape)
elif task == MULTILABEL_CLASSIFICATION:
pass
else:
raise NotImplementedError('f1_metric does not support task type %s'
% task)
bin_prediction = binarize_predictions(prediction, task)
# Bounding to avoid division by 0
eps = 1e-15
fn = np.sum(np.multiply(solution, (1 - bin_prediction)), axis=0, dtype=float)
tp = np.sum(np.multiply(solution, bin_prediction), axis=0, dtype=float)
fp = np.sum(np.multiply((1 - solution), bin_prediction), axis=0, dtype=float)
true_pos_num = sp.maximum(eps, tp + fn)
found_pos_num = sp.maximum(eps, tp + fp)
tp = sp.maximum(eps, tp)
tpr = tp / true_pos_num # true positive rate (recall)
ppv = tp / found_pos_num # positive predictive value (precision)
arithmetic_mean = 0.5 * sp.maximum(eps, tpr + ppv)
# Harmonic mean:
f1 = tpr * ppv / arithmetic_mean
# Average over all classes
f1 = np.mean(f1)
# Normalize: 0 for random, 1 for perfect
if task in (BINARY_CLASSIFICATION, MULTILABEL_CLASSIFICATION):
# How to choose the "base_f1"?
# For the binary/multilabel classification case, one may want to predict all 1.
# In that case tpr = 1 and ppv = frac_pos. f1 = 2 * frac_pos / (1+frac_pos)
# frac_pos = mvmean(solution.ravel())
# base_f1 = 2 * frac_pos / (1+frac_pos)
# or predict random values with probability 0.5, in which case
# base_f1 = 0.5
# the first solution is better only if frac_pos > 1/3.
# The solution in which we predict according to the class prior frac_pos gives
# f1 = tpr = ppv = frac_pos, which is worse than 0.5 if frac_pos<0.5
# So, because the f1 score is used if frac_pos is small (typically <0.1)
# the best is to assume that base_f1=0.5
base_f1 = 0.5
# For the multiclass case, this is not possible (though it does not make much sense to
# use f1 for multiclass problems), so the best would be to assign values at random to get
# tpr=ppv=frac_pos, where frac_pos=1/label_num
elif task == MULTICLASS_CLASSIFICATION:
label_num = solution.shape[1]
base_f1 = 1. / label_num
score = (f1 - base_f1) / sp.maximum(eps, (1 - base_f1))
return score
0
Example 50
Project: RLScore Source File: global_rankrls.py
def _leave_pair_out_python(self, pairs_start_inds, pairs_end_inds, oind=0):
"""Computes leave-pair-out predictions for a trained RankRLS.
Parameters
----------
pairs_start_inds : list of indices, shape = [n_pairs]
list of indices from range [0, n_samples-1]
pairs_end_inds : list of indices, shape = [n_pairs]
list of indices from range [0, n_samples-1]
oind : index from range [0, n_samples-1]
column of Y, for which pairwise cv is computed
Returns
-------
P1 : array, shape = [n_pairs]
holdout predictions for pairs_start_inds
P2 : array, shape = [n_pairs]
holdout predictions for pairs_end_inds
Notes
-----
Computes the leave-pair-out cross-validation predicitons, where each (i,j) pair with
i= pair_start_inds[k] and j = pairs_end_inds[k] is left out in turn.
When estimating area under ROC curve with leave-pair-out, one should leave out all
positive-negative pairs, while for estimating the general ranking error one should
leave out all pairs with different labels.
Computational complexity of holdout with most pairs left out: m = n_samples
O(TODO)
The leave-pair-out cross-validation algorithm is described in [1,2]. The use of
leave-pair-out cross-validation for AUC estimation has been analyzed in [3]
References
----------
[1] Tapio Pahikkala, Evgeni Tsivtsivadze, Antti Airola, Jouni Jarvinen, and Jorma Boberg.
An efficient algorithm for learning to rank from preference graphs.
Machine Learning, 75(1):129-165, 2009.
[2] Tapio Pahikkala, Antti Airola, Jorma Boberg, and Tapio Salakoski.
Exact and efficient leave-pair-out cross-validation for ranking RLS.
In Proceedings of the 2nd International and Interdisciplinary Conference
on Adaptive Knowledge Representation and Reasoning (AKRR'08), pages 1-8,
Espoo, Finland, 2008.
[3] Antti Airola, Tapio Pahikkala, Willem Waegeman, Bernard De Baets, Tapio Salakoski.
An Experimental Comparison of Cross-Validation Techniques for Estimating the Area Under the ROC Curve.
Computational Statistics & Data Analysis 55(4), 1828-1844, 2011.
"""
evals, svecs = self.evals, self.svecs
m = self.size
Y = self.Y
#This is, in the worst case, a cubic operation.
#If there are multiple outputs,
#this operation should be common for them all. THIS IS TO BE FIXED!
def computeG():
regparam = self.regparam
G = svecs * multiply(multiply(evals, 1. / ((m - 2.) * evals + regparam)).T, svecs.T)
return G
G = computeG()
GDY = (self.size - 2.) * G * Y
GC = sum(G, axis=1)
CTGC = sum(GC)
CTY = sum(Y, axis=0)[0, oind]
CTGDY = sum(GDY, axis=0)[0, oind]
sm2 = self.size - 2.
sqrtsm2 = sqrt(sm2)
#Array is faster to access than matrix
G = array(G)
#Lists are faster to access than matrices or arrays
def hack():
GDY_ = []
sqrtsm2GDY_ = []
GC_ = []
Y_ = []
BTY_ = []
Gdiag_ = []
sm2Gdiag_ = []
BTGBBTY_ = []
for i in range(m):
GDYi = GDY[i, oind]
GDY_.append(GDYi)
sqrtsm2GDY_.append(sqrtsm2 * GDYi)
GC_.append(GC[i, 0])
Yi = Y[i, oind]
Y_.append(Yi)
BTY_.append(sqrtsm2 * Yi)
Gii = G[i, i]
Gdiag_.append(Gii)
sm2Gdiag_.append(sm2 * Gii - 1.)
BTGBBTY_.append(sm2 * Gii * sqrtsm2 * Yi)
return GDY_, sqrtsm2GDY_, GC_, Y_, BTY_, Gdiag_, sm2Gdiag_, BTGBBTY_
GDY_, sqrtsm2GDY_, GC_, Y_, BTY_, Gdiag_, sm2Gdiag_, BTGBBTY_ = hack()
results_start, results_end = [], []
#This loops through the list of hold-out pairs.
#Each pair is handled in a constant time.
def looppairs(results_start, results_end):
for pairind in range(len(pairs_start_inds)):
i, j = pairs_start_inds[pairind], pairs_end_inds[pairind]
Gii = Gdiag_[i]
Gij = G[i, j]
Gjj = Gdiag_[j]
GCi = GC_[i]
GCj = GC_[j]
Yi = Y_[i]
Yj = Y_[j]
GDYi = GDY_[i]
GDYj = GDY_[j]
BTY0 = CTY - Yi - Yj
BTY1 = BTY_[i]
BTY2 = BTY_[j]
GiipGij = Gii + Gij
GijpGjj = Gij + Gjj
GCipGCj = GCi + GCj
BTGB00 = GiipGij + GijpGjj + CTGC - GCipGCj - GCipGCj
BTGB01 = sqrtsm2 * (GCi - GiipGij)
BTGB02 = sqrtsm2 * (GCj - GijpGjj)
BTGB12 = sm2 * Gij
BTGLY0 = CTGDY - (GDYi + GDYj + BTGB00 * BTY0 + BTGB01 * BTY1 + BTGB02 * BTY2)
BTGLY1 = sqrtsm2GDY_[i] - (BTGB01 * BTY0 + BTGBBTY_[i] + BTGB12 * BTY2)
BTGLY2 = sqrtsm2GDY_[j] - (BTGB02 * BTY0 + BTGB12 * BTY1 + BTGBBTY_[j])
print CTGDY, BTGLY0
BTGB00m1 = BTGB00 - 1.
BTGB11m1 = sm2Gdiag_[i]
BTGB22m1 = sm2Gdiag_[j]
CF00 = BTGB11m1 * BTGB22m1 - BTGB12 * BTGB12
CF01 = -BTGB01 * BTGB22m1 + BTGB12 * BTGB02
CF02 = BTGB01 * BTGB12 - BTGB11m1 * BTGB02
CF11 = BTGB00m1 * BTGB22m1 - BTGB02 * BTGB02
CF12 = -BTGB00m1 * BTGB12 + BTGB01 * BTGB02
CF22 = BTGB00m1 * BTGB11m1 - BTGB01 * BTGB01
invdeter = 1. / (BTGB00m1 * CF00 + BTGB01 * CF01 + BTGB02 * CF02)
b0 = invdeter * (CF00 * BTGLY0 + CF01 * BTGLY1 + CF02 * BTGLY2) + BTY0
b1 = invdeter * (CF01 * BTGLY0 + CF11 * BTGLY1 + CF12 * BTGLY2) + BTY1
b2 = invdeter * (CF02 * BTGLY0 + CF12 * BTGLY1 + CF22 * BTGLY2) + BTY2
t1 = -b0 + sqrtsm2 * b1
t2 = -b0 + sqrtsm2 * b2
F0 = GDYi - (Gii * t1 + Gij * t2 + GCi * b0)
F1 = GDYj - (Gij * t1 + Gjj * t2 + GCj * b0)
results_start.append(F0), results_end.append(F1)
looppairs(results_start, results_end)
return np.array(results_start), np.array(results_end)