def torque_raise(self,F,fs,dm,lead,dc,fc,alpha=2):
     if self.type == 'Square':
         return F*dm/2.0*((lead+np.pi*fs*dm)/(np.pi*dm-fs*lead)) + F*fc*dc/2.0
     elif self.type == 'ACME':#beep beep
         return F*dm/2.0*((lead+np.pi*fs*dm*np.reciprocal(np.cos(np.radians(14.5))))/(np.pi*dm-fs*lead*np.reciprocal(np.cos(np.radians(14.5))))) + F*fc*dc/2
     else:
         return F*dm/2.0*((lead+np.pi*fs*dm*np.reciprocal(np.cos(np.radians(alpha))))/(np.pi*dm-fs*lead*np.reciprocal(np.cos(np.radians(alpha))))) + F*fc*dc/2

def invert_sart(csc_A,csc_b,max_iterations=50,lam_start=1.0):
    # This code has been checked against Scott Silburn's Matlab code

    shap = csc_A.shape
    lam = lam_start
    colsum = (csc_A.transpose()).dot(sps.csc_matrix(np.ones(shap[0])).transpose())
    lamda = colsum
    #lamda = lamda.multiply(colsum != 0)
    np.reciprocal(lamda.data,out=lamda.data)
    np.multiply(lamda.data,lam,out=lamda.data)
    
    # Initialise output
    sol = sps.csc_matrix(np.zeros((shap[1],1))+np.exp(-1))
    # Create an array to monitor the convergence
    conv = np.zeros(max_iterations)
    
    for i in range(max_iterations):
        # Calculate sol_new = sol+lambda*(x'*(b-Ax))
        tmp = csc_b.transpose()-csc_A.dot(sol)
        tmp2 = csc_A.transpose().dot(tmp)
        #newsol = sol+tmp2*lamda
        newsol = sol+tmp2.multiply(lamda)
        # Eliminate negative values
        newsol = newsol.multiply(newsol > 0.0)
        newsol.eliminate_zeros()
        # Calculate how quickly the code is converging
        conv[i] = (sol.multiply(sol).sum()-newsol.multiply(newsol).sum())/sol.multiply(sol).sum()
        # Set the new solution to be the old solution and repeat
        sol = newsol
        
    return newsol.todense(), conv

 def torque_lower(self,l,k,C):
     if self.type == 'Square':
         return F*dm/2.0*((-lead+np.pi*fs*dm)/(np.pi*dm+fs*lead)) + F*fc*dc/2.0
     elif self.type == 'ACME':#beep beep
         return F*dm/2.0*((-lead+np.pi*fs*dm*np.reciprocal(np.cos(np.radians(14.5))))/(np.pi*dm+fs*lead*np.reciprocal(np.cos(np.radians(14.5))))) + F*fc*dc/2
     else:
         return F*dm/2.0*((-lead+np.pi*fs*dm*np.reciprocal(np.cos(np.radians(alpha))))/(np.pi*dm+fs*lead*np.reciprocal(np.cos(np.radians(alpha))))) + F*fc*dc/2

    def predict(self):

        m, n = self.A.shape # m observations

        convgraph = np.zeros(self.maxiter / 25)
        prevdist = 0.
        converged = False

        eps = 1e-6

        dd = np.array(self.A.sum(1))[:,0]
        D = diags(dd,0, format="csr")

        m, n = self.A.shape


        # random initialization, will initialize with K-means if told to
        H = csr_matrix(np.random.rand(m, self.k))

        EPS = csr_matrix(np.ones(H.shape) * eps)

        if self._embedding:
            # Apply eigenspace embedding K-means for initialization (Ng Weiss Jordan)

            Dz = diags(1 / (np.sqrt(dd) + eps), 0, format="csr")
            DAD = Dz.dot(self.A).dot(Dz)

            V = eigs(DAD, self.k)[1].real
            km_data = V / (np.linalg.norm(V, 2, axis=1).T * np.ones((self.k,1))).T

            km_predict = KMeans(n_clusters=self.k).fit_predict(km_data)

            indices = km_predict
            indptr = range(len(indices)+1)
            data = np.ones(len(indices))
            H = csr_matrix((data, indices, indptr))

        # Run separately for sparse and dense versions

        for i in range(self.maxiter):

            AH = self.A.dot(H)
            alpha = H.T.dot(AH)

            M1 = AH + EPS
            M2 = D.dot(H).dot(alpha) + EPS

            np.reciprocal(M2.data, out=M2.data)
            d1 = M1.multiply(M2).sqrt()

            H = H.multiply(d1)

            if i % 25 == 0:
                dist = sptrace(alpha)
                convgraph[i/25] = dist

                diff = dist / prevdist - 1
                prevdist = dist

        return NMFResult((H.toarray(),), convgraph, pdist)

def get_reconstructed(
        S0, S1, D0, D1, L,
        U0, U1, lam0, lam1, XQ,
        ):
    """
    Return the reconstructed matrix given a spectral form.
    """
    R11 = ndot(
            np.diag(np.reciprocal(np.sqrt(D0))),
            U0,
            np.diag(lam0),
            U0.T,
            np.diag(np.reciprocal(D0)),
            )
    R22 = ndot(
            np.diag(np.reciprocal(np.sqrt(D1))),
            U1,
            np.diag(lam1),
            U1.T,
            np.diag(np.reciprocal(D1)),
            )
    Q_reconstructed = build_block_2x2([
        [R11, ndot(R11, XQ) - ndot(XQ, R22)],
        [np.zeros_like(np.diag(L)), R22],
        ])
    return Q_reconstructed

def directions(degrees):
	if degrees == 0:
		return "You are going East"
	elif degrees == 90:
		return "You are going North"
	elif degrees == 180:
		return "You are going West"
	elif degrees == 270:
		return "You are going South"
	directions = [["North ","East "],["North ","West "],["South ", "West "],["South ", "East "]]
	main_dir = int(degrees/90)
	if degrees<90 or (degrees>180 and degrees<270):
		if degrees%90>45:
			num_say = int(numpy.reciprocal(round((90-(degrees%90)))/90)/2)
			if num_say == 0:
				num_say = 1
			return "You are going " + directions[main_dir][0] * num_say + directions[main_dir][1]
		else:
			num_say = int(numpy.reciprocal(round((degrees%90))/90)/2)
			if num_say == 0:
				num_say = 1
			return "You are going " + directions[main_dir][0] + directions[main_dir][1] * num_say
	else:
		if degrees%90>45:
			num_say = int(numpy.reciprocal(round((90-(degrees%90)))/90)/2)
			if num_say == 0:
				num_say = 1
			return "You are going " + directions[main_dir][0] + directions[main_dir][1] * num_say
		else:
			num_say = int(numpy.reciprocal(round((degrees%90))/90)/2)
			if num_say == 0:
				num_say = 1
			return "You are going " + directions[main_dir][0] * num_say + directions[main_dir][1]

def evaluation(dat , parent , p_order , trial , t_order):
    sum_dist = 0
    sum_dist_2 = 0

    # 親の評価
    for i in range(len(parent) - 1):
        sum_dist += distance(dat.ix[parent[i]] , dat.ix[parent[i + 1]])

    sum_dist += distance(dat.ix[parent[len(parent) - 1]] , dat.ix[parent[0]])

    fitness_parent = np.reciprocal(sum_dist)     # 小さい方が適応度が高いので逆数

    # 子の評価
    for i in range(len(trial) - 1):
        sum_dist_2 += distance(dat.ix[trial[i]] , dat.ix[trial[i + 1]])

    sum_dist_2 += distance(dat.ix[trial[len(trial) - 1]] , dat.ix[trial[0]])

    fitness_trial = np.reciprocal(sum_dist_2)     # 小さい方が適応度が高いので逆数

    # 親と子の比較
    if fitness_parent > fitness_trial:
        return fitness_parent , parent , p_order
    elif fitness_parent < fitness_trial:
        return fitness_trial , trial , t_order

    def predict_proba(self, X):
        """Probability estimates.

        The returned estimates for all classes are ordered by the
        label of classes.

        Parameters
        ----------
        X : array-like, shape = [n_samples, n_features]

        Returns
        -------
        T : array-like, shape = [n_samples, n_classes]
            Returns the probability of the sample for each class in the model,
            where classes are ordered as they are in ``self.classes_``.
        """
        # 1. / (1. + np.exp(-scores)), computed in-place
        prob = self.decision_function(X)
        prob *= -1
        np.exp(prob, prob)
        prob += 1
        np.reciprocal(prob, prob)
        if len(prob.shape) == 1:
            return np.vstack([1 - prob, prob]).T
        else:
            # OvR, not softmax, like Liblinear's predict_probability
            prob /= prob.sum(axis=1).reshape((prob.shape[0], -1))
            return prob

def _normalize_by_index(workspace, index):
    """
    Normalize each spectra of the specified workspace by the
    y-value at the specified index in that spectra.

    @param workspace    The workspace to normalize.
    @param index        The index of the y-value to normalize by.
    """
    number_of_histograms = workspace.getNumberHistograms()

    for idx in range(0, number_of_histograms):
        y_values = workspace.readY(idx)
        y_errors = workspace.readE(idx)

        # Avoid divide by zero
        if y_values[index] == 0.0:
            scale = np.reciprocal(1.0e-8)
        else:
            scale = np.reciprocal(y_values[index])

        # Normalise y values
        y_values_normalised = scale * y_values

        # Propagate y errors: C = A / B ; dC = sqrt( (dA/B)^2 + (A*dB/B^2)^2 )
        a = (y_errors*scale)
        b = (y_values*y_errors[index]*(scale ** 2))
        y_errors_propagated = np.sqrt(a ** 2 + b ** 2)

        workspace.setY(idx, y_values_normalised)
        workspace.setE(idx, y_errors_propagated)

    def forwardTrain(self, layer):
        forwardSink = layer.forward
        '''
        The logistic neuron applies the following transformation to the 
        linear combination of the incoming signals:
        1/1+e^(-x)
        
        This is done in place through the following logic:
        0. x = x + tiny (in place)
        1. x = x * -1 (in place)
        2. x = e^x (in place)
        3. x = 1 + x + tiny(in place)
        5. x = 1/x (in place)

        In this case, x is the "forwardSink" of this layer.
        
        The original data is "lost", but is no longer needed (techincally
        it could be recovered because each of the applied operations has
        an inverse.)
        '''
        np.add(forwardSink, tiny, out=forwardSink)
        np.multiply(forwardSink, -1, out=forwardSink)
        np.exp(forwardSink, out=forwardSink)
        np.add(forwardSink, 1 + tiny, out=forwardSink)
        np.reciprocal(forwardSink,out=forwardSink)

def estimate_dirichlet_param(samples, param):
    """
    Uses a Newton-Raphson scheme to estimating the parameter of a
    K-dimensional Dirichlet distribution

    :param samples: an NxK matrix of K-dimensional vectors drawn from
    a Dirichlet distribution
    :param param: the old value of the paramter. This is overwritten
    :return: a K-dimensional vector which is the new
    """

    N, K = samples.shape
    p = np.sum(np.log(samples), axis=0)

    for _ in range(60):
        g = -N * fns.digamma(param)
        g += N * fns.digamma(param.sum())
        g += p

        q = -N * fns.polygamma(1, param)
        np.reciprocal(q, out=q)

        z = N * fns.polygamma(1, param.sum())

        b = np.sum(g * q)
        b /= 1 / z + q.sum()

        param -= (g - b) * q

        print("%.2f" % param.mean(), end=" --> ")
    print

    return param

 def calMPRC_Fileter_scores(self):
     '''Calculate the weight of every term per document, using PubMed Related Citation (PRC) algorithm, Jimmy Lin and John Wilbur 2007.
        input: idf vector, docLen vector, occurrence count matrix (n documents, all terms in the vocabulary)
        output: a matrix of PRC scores.
     '''
     la = 0.022
     mu = 0.013
     score_threshold = 0.5 # the PRC weight threshold 
     div = mu/la
     ## generate m1
     reciSqrtIdf = np.reciprocal(np.sqrt(np.log(len(self.stemmed_corpus)*2.0/(self.df+1)))) # dim 1*19, conversion verified
     expDoclen = np.exp(self.doclen*(la-mu)) # dim 10*1, conversion verified
     m1 = np.dot(expDoclen,reciSqrtIdf) # dim 10*19, product verified
     ## generate m2: matrix
     matrix = np.power(div,self.doc_term_matrix)/div
     ## Hadamard product
     matrix = np.multiply(matrix,m1)
     ## offset
     offset = np.dot(np.ones((matrix.shape[0],1)),reciSqrtIdf)
     ## matrix+offset
     matrix = matrix+offset
     ## reciprocal of recWt
     raw_prc_matrix = np.reciprocal(matrix)
     ## reset scores for the terms that do not occur
     label = (self.doc_term_matrix>0)
     self.prc_matrix = np.multiply(label, raw_prc_matrix)
     
     ## modify the score matrix, remove terms with low scores
     keyword_index_vec = np.where(self.prc_matrix.A[self.pmidList.index(self.query),:]>score_threshold)[0].tolist()
     self.prc_matrix = self.prc_matrix.A[:,keyword_index_vec]

 def forward_cpu(self, inputs):
     self.retain_inputs((0, 1))
     x, gy = inputs
     gx = utils.force_array(numpy.square(x))
     gx += 1
     numpy.reciprocal(gx, out=gx)
     gx *= gy
     return gx,

 def backward_cpu(self, x, gy):
     gx = utils.force_array(numpy.square(x[0]))
     numpy.negative(gx, out=gx)
     gx += 1
     numpy.sqrt(gx, out=gx)
     numpy.reciprocal(gx, out=gx)
     gx *= gy[0]
     return gx,

 def forward(self, bottom_blobs, top_blobs):
   x = bottom_blobs[0].vals
   y = top_blobs[0].vals
   # Compute y = 1 / (1 + exp(-x))
   np.multiply(x, -1.0, out=y)
   np.exp(y, out=y)
   np.add(y, 1, out=y)
   np.reciprocal(y, out=y)

 def predict_proba(self, X):
     """ estimate probability """
     prob = -(np.dot(X, self.coef_.T) + self.intercept_)
     np.exp(prob, prob)
     prob += 1
     np.reciprocal(prob, prob)
     prob /= prob.sum(axis=1).reshape((prob.shape[0], -1))
     return prob

    def learn(self, features, labels):
        """Train the logistic regression model.
        :param features: The training instances' feature vectors
        :param labels: The training instances' labels
        """

        f = self.add_bias(features)
        w = np.random.normal(size=f.shape[1])

        regularization = self.regularization
        num_epochs = self.num_epochs

        # Optimize F_alpha for a suitable alpha (according to beta)
        alpha = None if self.f_beta_sq is None else np.reciprocal(self.f_beta_sq + 1.0)

        lr_init = self.learning_rate
        lr = lr_init

        batch_size = self.batch_size

        for epoch in xrange(num_epochs):
            if batch_size is None:
                f_batch = f
                labels_batch = labels
            else:
                batch = np.random.randint(0, f.shape[0], size=batch_size)
                f_batch = f[batch]
                labels_batch = labels[batch]

            predictions = np.array(map(sigmoid, f_batch.dot(w.T).T))

            if alpha is None:
                # Optimize Accuracy
                errors = labels_batch - predictions
                dw = f_batch.T.dot(errors.T).T

            else:
                # Optimize F_beta^2 
                dp = (f_batch.T.multiply(predictions * (1.0 - predictions))).T    # Matrix of the same shape as "features"
                sum_dp = dp.sum(axis=0)                                     # Vector (len = num features)
                sum_g_dp = dp.T.dot(labels_batch.T)                         # Vector (len = num features)
                sum_p = predictions.sum()                                   # Scalar
                sum_g = labels_batch.sum()                                  # Scalar
                sum_g_p = labels_batch.dot(predictions)                     # Scalar

                # Update rule for optimizing F_beta^2 - see paper for details
                denominator = np.reciprocal(alpha*sum_p + (1 - alpha)*sum_g)

                dw = denominator * sum_g_dp - alpha * denominator * denominator * sum_g_p * sum_dp
                dw *= len(labels_batch)
                if type(dw) != np.ndarray:
                    dw = np.array(dw)[0]

            w += lr * (dw - regularization*w)

            lr = lr_init * (1.0 - (float(epoch) / num_epochs))

        self.weights = w

    def learn(self, path_sets, labels):
        """Train the weighted edge model
        :param path_sets -- the term-pairs represented as path-sets
        :param labels -- the term-pairs gold standard annotations
        """
        regularization = self.regularization
        num_epochs = self.num_epochs

        # Compute alpha = 1/(1+beta^2) - for F measure
        alpha = None if self.f_beta_sq is None else np.reciprocal(self.f_beta_sq + 1.0)

        lr_init = self.learning_rate
        lr = lr_init

        num_pairs = len(labels)

        pair_to_paths, path_features = prebuild_data_structures(path_sets)

        # Initialize the edge types weights randomly
        w = np.random.normal(size=path_features.shape[1])

        # Choose a maximum-score path for each term-pair randomly
        max_paths = np.zeros(num_pairs, dtype=np.int)
        for i, (path_set, label) in enumerate(zip(pair_to_paths, labels)):
            max_paths[i] = np.random.choice(path_set)

        for epoch in xrange(num_epochs):

            # "M-step": Use F_beta derivative to update w
            f = path_features[max_paths]

            predictions = probabilities(w, f)

            dp = (f.T * (predictions * (1.0 - predictions))).T    # Matrix of the same shape as "features"
            sum_dp = np.sum(dp, axis=0)                           # Vector (len = num features)
            sum_g_dp = labels.dot(dp)                             # Vector (len = num features)
            sum_p = np.sum(predictions)                           # Scalar
            sum_g = np.sum(labels)                                # Scalar
            sum_g_p = labels.dot(predictions)                     # Scalar

            denominator = np.reciprocal(alpha*sum_p + (1 - alpha)*sum_g)

            dw = denominator * sum_g_dp - alpha * denominator * denominator * sum_g_p * sum_dp
            dw *= num_pairs

            w += lr * (dw - regularization*w)

            # Reduce the learning rate
            lr = lr_init * (1.0 - (float(epoch) / num_epochs))

            # "E-step": Give a score to each path according to the current weights,
            # and choose the highest-scored path for each pair
            path_probabilities = probabilities(w, path_features)
            for i, (path_set, label) in enumerate(zip(pair_to_paths, labels)):
                max_paths[i] = path_set[np.argmax(path_probabilities[path_set])]

        self.weights = w

def extract_recipcontact(project, close, stride, far=None):
    A,B,C = triplets(project, close, stride, far)
    m = metrics.ContinuousContact(contacts='all', scheme='CA')
    pA, pB, pC = map(m.prepare_trajectory, [A, B, C])

    # reciprocate the maps
    pA, pB, pC = np.reciprocal(pA), np.reciprocal(pB), np.reciprocal(pC)

    return pA - pB, pA - pC

 def numpy_run(self):
     """Forward propagation from batch on CPU only.
     """
     super(All2AllSigmoid, self).numpy_run()
     self.output.map_write()
     mem = self.output.mem
     # 1 / (1 + numpy.exp(-mem))
     numpy.exp(-mem, mem)
     numpy.reciprocal(mem + 1, mem)