d/dx  (y - lmbda)^T C = d/dx -exp(Cx + d)^T C
            = -C^T exp(Cx + d)^T C
557120
def hessian_local_log_likelihood(self, x):
        
        # Observation likelihoods
        lmbda = np.exp(np.dot(x, self.C.T) + np.dot(self.inputs, self.D.T))
        return np.einsum('tn, ni, nj ->tij', -lmbda, self.C, self.C)
56082
def vap(x, a, b, c):
    
    return np.exp(a+b/x+c*np.log(x))
966088
def d2Sbus_dV2(Ybus, V, lam):
    
    n = len(V)
    Ibus = Ybus * V
    diaglam = spdiag(lam)
    diagV = spdiag(V)

    A = spmatrix(mul(lam, V), range(n), range(n))
    B = Ybus * diagV
    C = A * conj(B)
    D = Ybus.H * diagV
    E = conj(diagV) * (D * diaglam - spmatrix(D*lam, range(n), range(n)))
    F = C - A * spmatrix(conj(Ibus), range(n), range(n))
    G = spmatrix(div(matrix(1.0, (n, 1)), abs(V)), range(n), range(n))

    Gaa = E + F
    Gva = 1j * G * (E - F)
    Gav = Gva.T
    Gvv = G * (C + C.T) * G

    return Gaa, Gav, Gva, Gvv
397481
def dBinaryRochedx (r, D, q, F):
    
    return -r[0]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*(r[0]-D)*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5 -q/D/D + F*F*(1+q)*r[0]
1186747
def _derX(self,x,y):
        '''
        
        '''
        x_pos, y_pos = self.findSector(x,y)
        alpha, beta = self.findCoords(x,y,x_pos,y_pos)

        # Get four corners data for each point
        xA = self.x_values[x_pos,y_pos]
        xB = self.x_values[x_pos+1,y_pos]
        xC = self.x_values[x_pos,y_pos+1]
        xD = self.x_values[x_pos+1,y_pos+1]
        yA = self.y_values[x_pos,y_pos]
        yB = self.y_values[x_pos+1,y_pos]
        yC = self.y_values[x_pos,y_pos+1]
        yD = self.y_values[x_pos+1,y_pos+1]
        fA = self.f_values[x_pos,y_pos]
        fB = self.f_values[x_pos+1,y_pos]
        fC = self.f_values[x_pos,y_pos+1]
        fD = self.f_values[x_pos+1,y_pos+1]

        # Calculate components of the alpha,beta --> x,y delta translation matrix
        alpha_x = (1-beta)*(xB-xA) + beta*(xD-xC)
        alpha_y = (1-beta)*(yB-yA) + beta*(yD-yC)
        beta_x  = (1-alpha)*(xC-xA) + alpha*(xD-xB)
        beta_y  = (1-alpha)*(yC-yA) + alpha*(yD-yB)

        # Invert the delta translation matrix into x,y --> alpha,beta
        det = alpha_x*beta_y - beta_x*alpha_y
        x_alpha = beta_y/det
        x_beta  = -alpha_y/det

        # Calculate the derivative of f w.r.t. alpha and beta
        dfda = (1-beta)*(fB-fA) + beta*(fD-fC)
        dfdb = (1-alpha)*(fC-fA) + alpha*(fD-fB)

        # Calculate the derivative with respect to x (and return it)
        dfdx = x_alpha*dfda + x_beta*dfdb
        return dfdx
1131953
def ystep(self):
        
        """

        self.Y[..., 0:-1] = sp.prox_l2(
            self.AX[..., 0:-1] + self.U[..., 0:-1],
            (self.lmbda/self.rho)*self.Wtvna, axis=self.saxes)
        self.Y[..., -1] = sp.prox_l1(
            self.AX[..., -1] + self.U[..., -1] - self.S,
            (1.0/self.rho)*self.Wdf)
317955
def ystep(self):
        """

        self.Y = np.asarray(sp.prox_l1(self.AX + self.U,
                                       (self.lmbda / self.rho) * self.wl1),
                            dtype=self.dtype)
        super(BPDN, self).ystep()
459067
def _Aterm(N, h, m, k, dW):
    """"""
    sqrt2h = np.sqrt(2.0/h)
    Xk = np.random.normal(0.0, 1.0, (N, m, 1))
    Yk = np.random.normal(0.0, 1.0, (N, m, 1))
    term1 = _dot(Xk, _t(Yk + sqrt2h*dW))
    term2 = _dot(Yk + sqrt2h*dW, _t(Xk))
    return (term1 - term2)/k
480154
def derivatives(self, x, y, amp, sigma, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_**2 + y_**2)
        sigma_x, sigma_y = sigma, sigma
        if isinstance(R, int) or isinstance(R, float):
            R = max(R, self.ds)
        else:
            R[R <= self.ds] = self.ds
        alpha = self.alpha_abs(R, amp, sigma)
        return alpha / R * x_, alpha / R * y_
707982
def F_t(X, Y, S, M_E, E, m0, rho):
    '''  '''
    r = X.shape[1]
    out1 = (((np.dot(np.dot(X, S), Y.T) - M_E) * E)**2).sum() / 2
    out2 = rho * G(Y, m0, r)
    out3 = rho * G(X, m0, r)

    return out1 + out2 + out3
233193
def _derY(self,x,y):
        '''
        
        '''
        m = len(x)
        temp = np.zeros((m,self.funcCount))
        for j in range(self.funcCount):
            temp[:,j] = self.functions[j](x,y)
        temp[np.isnan(temp)] = np.inf
        i = np.argmin(temp,axis=1)
        y = temp[np.arange(m),i]
        dfdy = np.zeros_like(x)
        for j in range(self.funcCount):
            c = i == j
            dfdy[c] = self.functions[j].derivativeY(x[c],y[c])
        return dfdy
1327673
def ystep(self):
        
        """

        self.Y = np.asarray(sp.prox_l2(
            self.AX + self.U, (self.lmbda/self.rho)*self.Wtvna,
            axis=self.saxes), dtype=self.dtype)
628553
def _dK_ode_dtheta(self, target):
        """"""
        t_ode = self._t[self._index>0]
        dL_dK_ode = self._dL_dK[self._index>0, :]
        index_ode = self._index[self._index>0]-1
        if self._t2 is None:
            if t_ode.size==0:
                return        
            t2_ode = t_ode
            dL_dK_ode = dL_dK_ode[:, self._index>0]
            index2_ode = index_ode
        else:
            t2_ode = self._t2[self._index2>0]
            dL_dK_ode = dL_dK_ode[:, self._index2>0]
            if t_ode.size==0 or t2_ode.size==0:
                return
            index2_ode = self._index2[self._index2>0]-1

        h1 = self._compute_H(t_ode, index_ode, t2_ode, index2_ode, stationary=self.is_stationary, update_derivatives=True)
        #self._dK_ddelay = self._dh_ddelay
        self._dK_dsigma = self._dh_dsigma

        if self._t2 is None:
            h2 = h1
        else:
            h2 = self._compute_H(t2_ode, index2_ode, t_ode, index_ode, stationary=self.is_stationary, update_derivatives=True)

        #self._dK_ddelay += self._dh_ddelay.T
        self._dK_dsigma += self._dh_dsigma.T
        # C1 = self.sensitivity
        # C2 = self.sensitivity

        # K = 0.5 * (h1 + h2.T)
        # var2 = C1*C2
        # if self.is_normalized:
        #     dk_dD1 = (sum(sum(dL_dK.*dh1_dD1)) + sum(sum(dL_dK.*dh2_dD1.T)))*0.5*var2
        #     dk_dD2 = (sum(sum(dL_dK.*dh1_dD2)) + sum(sum(dL_dK.*dh2_dD2.T)))*0.5*var2
        #     dk_dsigma = 0.5 * var2 * sum(sum(dL_dK.*dK_dsigma))
        #     dk_dC1 = C2 * sum(sum(dL_dK.*K))
        #     dk_dC2 = C1 * sum(sum(dL_dK.*K))
        # else:
        #     K = np.sqrt(np.pi) * K
        #     dk_dD1 = (sum(sum(dL_dK.*dh1_dD1)) + * sum(sum(dL_dK.*K))
        #     dk_dC2 = self.sigma * C1 * sum(sum(dL_dK.*K))


        # dk_dSim1Variance = dk_dC1
        # Last element is the length scale.
        (dL_dK_ode[:, :, None]*self._dh_ddelay[:, None, :]).sum(2)

        target[-1] += (dL_dK_ode*self._dK_dsigma/np.sqrt(2)).sum()
493127
def ystep(self):
        """

        self.Y = sp.prox_l1(self.AX + self.U,
                            (self.lmbda / self.rho) * self.wl1)
        super(ConvBPDN, self).ystep()
9314
def f(self, y):
        
        d = self.d
        mpsi = self.psi
        z = d * y.copy()
        for i in range(len(mpsi)):
            a, b, c = mpsi[i]
            z += a * np.tanh(b * (y + c))
        return z
1153747
def Itaax(mt, x, t):
    
    return (Sx(mt, x) - Sx(mt, x + t)) / Dx(mt, x)
1084635
def adapt(self, d, x):
        
        y = np.dot(self.w, x)
        e = d - y
        R1 = np.dot(np.dot(np.dot(self.R,x),x.T),self.R)
        R2 = self.mu + np.dot(np.dot(x,self.R),x.T)
        self.R = 1/self.mu * (self.R - R1/R2)
        dw = np.dot(self.R, x.T) * e
        self.w += dw