calculate Debye energy using old fortran routine

    :params x: Debye x value
    :return: Debye energy
    :note: this is a wraper function (vectorization) for debye_E_single
177100
def debye_E(x):
    
    if isuncertainties([x]):
        f_u = uct.wrap(debye_E_single)
        f_v = np.vectorize(f_u)
    else:
        f_v = np.vectorize(debye_E_single)
    return f_v(x)
985715
def _f(self, x, user_data=None):
        
        p_gen = x[self._Pg.i1:self._Pg.iN + 1] # Active generation in p.u.
        q_gen = x[self._Qg.i1:self._Qg.iN + 1] # Reactive generation in p.u.

        # Polynomial cost of P and Q.
        xx = r_[p_gen, q_gen] * self._base_mva
        if len(self._ipol) > 0:
            f = sum([g.total_cost(xx[i]) for i,g in enumerate(self._gn)])
        else:
            f = 0

        # Piecewise linear cost of P and Q.
        if self._ny:
            y = self.om.get_var("y")
            self._ccost = csr_matrix((ones(self._ny),
                (range(y.i1, y.iN + 1), zeros(self._ny))),
                shape=(self._nxyz, 1)).T
            f = f + self._ccost * x
        else:
            self._ccost = zeros((1, self._nxyz))
        # TODO: Generalised cost term.

        return f
903302
def _hammer_function_precompute(self,x0, L, Min, model):
        
        if x0 is None: return None, None
        if len(x0.shape)==1: x0 = x0[None,:]
        m = model.predict(x0)[0]
        pred = model.predict(x0)[1].copy()
        pred[pred<1e-16] = 1e-16
        s = np.sqrt(pred)
        r_x0 = (m-Min)/L
        s_x0 = s/L
        r_x0 = r_x0.flatten()
        s_x0 = s_x0.flatten()
        return r_x0, s_x0
555525
def prox_xline(x, step):
    """"""
    if not np.isscalar(x):
        x= x[0]
    if x > 0.5:
        return np.array([0.5])
    else:
        return np.array([x])
1206852
def _generate_energy_edges_single(ene):
    """"""
    midene = np.sqrt((ene[1:] * ene[:-1]))
    elo, ehi = np.zeros(len(ene)) * ene.unit, np.zeros(len(ene)) * ene.unit
    elo[1:] = ene[1:] - midene
    ehi[:-1] = midene - ene[:-1]
    elo[0] = ene[0] * (1 - ene[0] / (ene[0] + ehi[0]))
    ehi[-1] = elo[-1]
    return u.Quantity([elo, ehi])
333466
def _call(self, x):
        """"""
        # The algorithm takes finite differences in one dimension at a time
        # reusing the dx vector to improve efficiency.
        dfdx = self.domain.zero()
        dx = self.domain.zero()

        if self.method == 'backward':
            fx = self.functional(x)
            for i in range(self.domain.size):
                dx[i - 1] = 0  # reset step from last iteration
                dx[i] = self.step
                dfdx[i] = fx - self.functional(x - dx)
        elif self.method == 'forward':
            fx = self.functional(x)
            for i in range(self.domain.size):
                dx[i - 1] = 0  # reset step from last iteration
                dx[i] = self.step
                dfdx[i] = self.functional(x + dx) - fx
        elif self.method == 'central':
            for i in range(self.domain.size):
                dx[i - 1] = 0  # reset step from last iteration
                dx[i] = self.step / 2
                dfdx[i] = self.functional(x + dx) - self.functional(x - dx)
        else:
            raise RuntimeError('unknown method')

        dfdx /= self.step
        return dfdx
30103
def pre_process(self, x0, params=()):
        
        # Should be used by all methods matching "solve_*"
        if self.x_by_name and isinstance(x0, dict):
            x0 = [x0[k] for k in self.names]
        if self.par_by_name and isinstance(params, dict):
            params = [params[k] for k in self.param_names]
        for pre_processor in self.pre_processors:
            x0, params = pre_processor(x0, params)
        return x0, np.atleast_1d(params)
481188
def get_no_mixing_energy(self):
        
        assert self.grand == 1, \
            'Please provide grand potential phase diagram ' \
            'for computing no_mixing_energy!'

        energy1 = self.pd.get_hull_energy(self.comp1) - \
            self._get_grand_potential(self.c1)
        energy2 = self.pd.get_hull_energy(self.comp2) - \
            self._get_grand_potential(self.c2)
        unit = 'eV/f.u.'
        if self.norm:
            unit = 'eV/atom'
        return [(self.c1_original.reduced_formula +
                 ' ({0})'.format(unit), energy1),
                (self.c2_original.reduced_formula +
                 ' ({0})'.format(unit), energy2)]
1203508
def bethe_lattice(energy, hopping):
    """"""
    energy = np.asarray(energy).clip(-2*hopping, 2*hopping)
    return np.sqrt(4*hopping**2 - energy**2) / (2*np.pi*hopping**2)
578972
def _F(self, X, tau):
        
        t2 = tau ** 2
        #Fx = self.F(X)

        _F = self.F(X)
        a = t2*(t2+1)**-2
        if isinstance(X, np.ndarray):
            #b = (t2 + 1) * (X ** 2 - 1) ** -1 * (1 - _F)
            b = np.ones_like(X)
            b[X == 1] = (t2+1) * 1./3
            b[X != 1] = (t2 + 1) * (X[X != 1] ** 2 - 1) ** -1 * (1 - _F[X != 1])

        elif isinstance(X, float) or isinstance(X, int):
            if X == 1:
                b = (t2+1)* 1./3
            else:
                b = (t2+1)*(X**2-1)**-1*(1-_F)
        else:
            raise ValueError("The variable type is not compatible with the function, please use float, int or ndarray's.")

        c = 2*_F
        d = -np.pi*(t2+X**2)**-0.5
        e = (t2-1)*(tau*(t2+X**2)**0.5)**-1*self.L(X,tau)
        result = a * (b + c + d + e)

        return result
1285746
def sum(x, weights=None):
    '''
    
    '''
    x = to_potential(x)
    if is_const_potential(x): return PotentialConstant(np.sum(x.c))
    else:                     return PotentialSum(x, weights=weights)
1369396
def get_energy(self, y_true, input_energy, mask):
        
        input_energy = K.sum(input_energy * y_true, 2)  # (B, T)
        chain_energy = K.sum(K.dot(y_true[:, :-1, :], self.chain_kernel) * y_true[:, 1:, :], 2)  # (B, T-1)

        if mask is not None:
            mask = K.cast(mask, K.floatx())
            chain_mask = mask[:, :-1] * mask[:, 1:]  # (B, T-1), mask[:,:-1]*mask[:,1:] makes it work with any padding
            input_energy = input_energy * mask
            chain_energy = chain_energy * chain_mask
        total_energy = K.sum(input_energy, -1) + K.sum(chain_energy, -1)  # (B, )

        return total_energy
793811
private static double getEpsilon(IntDoubleVector x, IntDoubleVector d) {
        double machineEpsilon = 2.2204460492503131e-16;
        double xInfNorm = DoubleArrays.infinityNorm(x.toNativeArray());
        double dInfNorm = DoubleArrays.infinityNorm(d.toNativeArray());
        return machineEpsilon * (1.0 + xInfNorm) / dInfNorm;
    }
699384
double getBetheFreeEnergy() {
        // 
        // G_{Bethe} = \sum_a \sum_{x_a} - b(x_a) ln \chi(x_a)
        //              + \sum_a \sum_{x_a} b(x_a) ln b(x_a)
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //           = \sum_a \sum_{x_a} b(x_a) ln (b(x_a) / \chi(x_a))
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //
        //     where n_i is the number of neighbors of the variable x_i,
        //     b(x_a) and b(x_i) are normalized distributions and x_a is 
        //     the set of variables participating in factor a. 
        //
                
        double bethe = 0.0;
        for (int a=0; a<fg.getFactors().size(); a++) {
            Factor f = fg.getFactors().get(a);
            if (!(f instanceof GlobalFactor)) {
                int numConfigs = f.getVars().calcNumConfigs();
                VarTensor beliefs = getFactorBeliefs(a);
                for (int c=0; c<numConfigs; c++) {                
                    // Since we want multiplication by 0 to always give 0 (not the case for Double.POSITIVE_INFINITY or Double.NaN.
                    double b_c = beliefs.getValue(c);
                    if (b_c != s.zero()) {
                        double r_b_c = s.toReal(b_c);
                        double log_b_c = s.toLogProb(b_c);
                        double log_chi_c = f.getLogUnormalizedScore(c);
                        bethe += r_b_c * (log_b_c - log_chi_c);
                    }
                }
            } else {
                VarTensor[] inMsgs = getMsgs(f.getId(), msgs, IN_MSG);
                bethe += ((GlobalFactor) f).getExpectedLogBelief(inMsgs);
            }
        }
        for (int v=0; v<fg.getVars().size(); v++) {
            Var var = fg.getVars().get(v);
            int numNeighbors = bg.numNbsT1(v);
            VarTensor beliefs = getVarBeliefs(v);
            double sum = 0.0;
            for (int c=0; c<var.getNumStates(); c++) {
                double b_c = beliefs.getValue(c);
                if (b_c != s.zero()) {
                    double r_b_c = s.toReal(b_c);
                    double log_b_c = s.toLogProb(b_c);
                    sum += r_b_c * log_b_c;
                }
            }
            bethe -= (numNeighbors - 1) * sum;
        }
        
        assert !Double.isNaN(bethe);        
        return bethe;
    }
1298051
def energy_minimum(self, x):
    
    E      = self.energy
    y_min  = self.find_energy_minimum(E, x)
    return y_min
1097008
def _ttm_me_compute(self, V, edims, sdims, transp):
        
        shapeY = np.copy(self.shape)

        # Determine size of Y
        for n in np.union1d(edims, sdims):
            shapeY[n] = V[n].shape[1] if transp else V[n].shape[0]

        # Allocate Y (final result) and v (vectors for elementwise computations)
        Y = zeros(shapeY)
        shapeY = array(shapeY)
        v = [None for _ in range(len(edims))]

        for i in range(np.prod(shapeY[edims])):
            rsubs = unravel_index(shapeY[edims], i)
948225
def DebyeChain(q, Rg):
    
    a = (q * Rg) ** 2
    return 2 * (np.exp(-a) - 1 + a) / a ** 2