Differentiation function.

        Numerical approximation of a Rosenblatt transformation created from
        copula formulation.
174649
def _diff(self, x, th, eps):
        
        foo = lambda y: self.igen(numpy.sum(self.gen(y, th), 0), th)

        out1 = out2 = 0.
        sign = 1 - 2*(x > .5).T
        for I in numpy.ndindex(*((2,)*(len(x)-1)+(1,))):

            eps_ = numpy.array(I)*eps
            x_ = (x.T + sign*eps_).T
            out1 += (-1)**sum(I)*foo(x_)

            x_[-1] = 1
            out2 += (-1)**sum(I)*foo(x_)

        out = out1/out2
        return out
196837
def prepare_next_tree(self):
        """"""
        for edge in self.edges:
            copula_theta = edge.theta

            if self.level == 1:
                left_u = self.u_matrix[:, edge.L]
                right_u = self.u_matrix[:, edge.R]

            else:
                left_parent, right_parent = edge.parents
                left_u, right_u = Edge.get_conditional_uni(left_parent, right_parent)

            # compute conditional cdfs C(i|j) = dC(i,j)/duj and dC(i,j)/du
            left_u = [x for x in left_u if x is not None]
            right_u = [x for x in right_u if x is not None]
            X_left_right = np.array([[x, y] for x, y in zip(left_u, right_u)])
            X_right_left = np.array([[x, y] for x, y in zip(right_u, left_u)])

            copula = Bivariate(edge.name)
            copula.theta = copula_theta
            left_given_right = copula.partial_derivative(X_left_right)
            right_given_left = copula.partial_derivative(X_right_left)

            # correction of 0 or 1
            left_given_right[left_given_right == 0] = EPSILON
            right_given_left[right_given_left == 0] = EPSILON
            left_given_right[left_given_right == 1] = 1 - EPSILON
            right_given_left[right_given_left == 1] = 1 - EPSILON
            edge.U = np.array([left_given_right, right_given_left])
1326400
def from_dict(cls, copula_dict):
        """"""
        instance = cls()

        instance.fitted = copula_dict['fitted']
        instance.constant_value = copula_dict['constant_value']

        if instance.fitted and not instance.constant_value:
            instance.model = scipy.stats.gaussian_kde([-1, 0, 0])

            for key in ['dataset', 'covariance', 'inv_cov']:
                copula_dict[key] = np.array(copula_dict[key])

            attributes = ['d', 'n', 'dataset', 'covariance', 'factor', 'inv_cov']
            for name in attributes:
                setattr(instance.model, name, copula_dict[name])

        return instance
141532
def _fD(self, neg_pairs, A):
    
    diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
    return np.log(np.sum(np.sqrt(np.sum(np.dot(diff, A) * diff, axis=1))) + 1e-6)
1129582
def _lmder1_rosenbrock():
    """"""

    def func(params, vec):
        vec[0] = 10 * (params[1] - params[0]**2)
        vec[1] = 1 - params[0]

    def jac(params, jac):
        jac[0,0] = -20 * params[0]
        jac[0,1] = -1
        jac[1,0] = 10
        jac[1,1] = 0

    guess = np.asfarray([-1.2, 1])
    norm1s = [0.491934955050e+01, 0.134006305822e+04, 0.1430000511923e+06]

    for i in range(3):
        _lmder1_driver(2, func, jac, guess * 10**i,
                       norm1s[i], 0, [1, 1])
747232
def _deriv_hypot(x, y):
    """"""
    r = np.hypot(x, y)
    df_dx = x / r
    df_dy = y / r
    return np.hstack([df_dx, df_dy])
1076524
def _hessfcn(self, x, lmbda):
        
        Pgen = x[self._Pg.i1:self._Pg.iN + 1] # Active generation in p.u.
        Qgen = x[self._Qg.i1:self._Qg.iN + 1] # Reactive generation in p.u.

        for i, g in enumerate(self._gn):
            g.p = Pgen[i] * self._base_mva # active generation in MW
            g.q = Qgen[i] * self._base_mva # reactive generation in MVAr

        Vang = x[self._Va.i1:self._Va.iN + 1]
        Vmag = x[self._Vm.i1:self._Vm.iN + 1]
        V = Vmag * exp(1j * Vang)
        nxtra = self._nxyz - 2 * self._nb

        #------------------------------------------------------------------
        #  Evaluate d2f.
        #------------------------------------------------------------------

        d2f = self._d2f(x) * self.opt["cost_mult"]
        # TODO: Generalised cost model.

        #------------------------------------------------------------------
        #  Evaluate Hessian of power balance constraints.
        #------------------------------------------------------------------

        nlam = len(lmbda["eqnonlin"]) / 2
        lamP = lmbda["eqnonlin"][:nlam]
        lamQ = lmbda["eqnonlin"][nlam:nlam + nlam]
        Gpaa, Gpav, Gpva, Gpvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamP)
        Gqaa, Gqav, Gqva, Gqvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamQ)

        d2G = vstack([
            hstack([
                vstack([hstack([Gpaa, Gpav]),
                        hstack([Gpva, Gpvv])]).real +
                vstack([hstack([Gqaa, Gqav]),
                        hstack([Gqva, Gqvv])]).imag,
                csr_matrix((2 * self._nb, nxtra))]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        #------------------------------------------------------------------
        #  Evaluate Hessian of flow constraints.
        #------------------------------------------------------------------

        nmu = len(lmbda["ineqnonlin"]) / 2
        muF = lmbda["ineqnonlin"][:nmu]
        muT = lmbda["ineqnonlin"][nmu:nmu + nmu]
        if self.flow_lim == "I":
            dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = \
                self.om.case.dIbr_dV(self._Yf, self._Yt, V)
            Hfaa, Hfav, Hfva, Hfvv = \
                self.om.case.d2AIbr_dV2(dIf_dVa, dIf_dVm, If, self._Yf, V, muF)
            Htaa, Htav, Htva, Htvv = \
                self.om.case.d2AIbr_dV2(dIt_dVa, dIt_dVm, It, self._Yt, V, muT)
        else:
            f = [e.from_bus._i for e in self._ln]
            t = [e.to_bus._i for e in self._ln]
            # Line-bus connection matrices.
            Cf = csr_matrix((ones(self._nl), (range(self._nl), f)), (self._nl, self._nb))
            Ct = csr_matrix((ones(self._nl), (range(self._nl), t)), (self._nl, self._nb))
            dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
                self.om.case.dSbr_dV(self._Yf, self._Yt, V)
            if self.flow_lim == PFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(dSf_dVa.real(), dSf_dVm.real(),
                                            Sf.real(), Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(dSt_dVa.real(), dSt_dVm.real(),
                                            St.real(), Ct, self._Yt, V, muT)
            elif self.flow_lim == SFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(
                        dSf_dVa, dSf_dVm, Sf, Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(
                        dSt_dVa, dSt_dVm, St, Ct, self._Yt, V, muT)
            else:
                raise ValueError

        d2H = vstack([
            hstack([
                vstack([hstack([Hfaa, Hfav]),
                        hstack([Hfva, Hfvv])]) +
                vstack([hstack([Htaa, Htav]),
                        hstack([Htva, Htvv])]),
                csr_matrix((2 * self._nb, nxtra))
            ]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        return d2f + d2G + d2H
1251012
def rosen(self, x, alpha=1e2):
        """"""
        x = [x] if isscalar(x[0]) else x  # scalar into list
        f = [sum(alpha * (x[:-1]**2 - x[1:])**2 + (1. - x[:-1])**2) for x in x]
        return f if len(f) > 1 else f[0]  # 1-element-list into scalar
1038395
def _sample_row(self):
        
        unis = np.random.uniform(0, 1, self.n_var)
        # randomly select a node to start with
        first_ind = np.random.randint(0, self.n_var)
        adj = self.trees[0].get_adjacent_matrix()
        visited = []
        explore = [first_ind]

        sampled = np.zeros(self.n_var)
        itr = 0
        while explore:
            current = explore.pop(0)
            neighbors = np.where(adj[current, :] == 1)[0].tolist()
            if itr == 0:
                new_x = self.ppfs[current](unis[current])

            else:
                for i in range(itr - 1, -1, -1):
                    current_ind = -1

                    if i >= self.truncated:
                        continue

                    current_tree = self.trees[i].edges
                    # get index of edge to retrieve
                    for edge in current_tree:
                        if i == 0:
                            if (edge.L == current and edge.R == visited[0]) or\
                               (edge.R == current and edge.L == visited[0]):
                                current_ind = edge.index
                                break
                        else:
                            if edge.L == current or edge.R == current:
                                condition = set(edge.D)
                                condition.add(edge.L)
                                condition.add(edge.R)

                                visit_set = set(visited)
                                visit_set.add(current)

                                if condition.issubset(visit_set):
                                    current_ind = edge.index
                                break

                    if current_ind != -1:
                        # the node is not indepedent contional on visited node
                        copula_type = current_tree[current_ind].name
                        copula = Bivariate(CopulaTypes(copula_type))
                        copula.theta = current_tree[current_ind].theta
                        derivative = copula.partial_derivative_scalar

                        if i == itr - 1:
                            tmp = optimize.fminbound(
                                derivative, EPSILON, 1.0,
                                args=(unis[visited[0]], unis[current])
                            )
                        else:
                            tmp = optimize.fminbound(
                                derivative, EPSILON, 1.0,
                                args=(unis[visited[0]], tmp)
                            )

                        tmp = min(max(tmp, EPSILON), 0.99)

                new_x = self.ppfs[current](tmp)

            sampled[current] = new_x

            for s in neighbors:
                if s not in visited:
                    explore.insert(0, s)

            itr += 1
            visited.insert(0, current)

        return sampled
463897
def gradF_t(X, Y, S, M_E, E, m0, rho):
    ''' 
    '''
    n, r = X.shape
    m, r = Y.shape

    XS = np.dot(X, S)
    YS = np.dot(Y, S.T)
    XSY = np.dot(XS, Y.T)

    Qx = np.dot(np.dot(X.T, ((M_E - XSY) * E)), YS) / n
    Qy = np.dot(np.dot(Y.T, ((M_E - XSY) * E).T), XS) / m

    W = np.dot((XSY - M_E) * E, YS) + np.dot(X, Qx) + rho * Gp(X, m0, r)
    Z = np.dot(((XSY - M_E) * E).T, XS) + np.dot(Y, Qy) + rho * Gp(Y, m0, r)

    return W, Z
1101510
def _df(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

        iPg = range(self._Pg.i1, self._Pg.iN + 1)
        iQg = range(self._Qg.i1, self._Qg.iN + 1)

        # Polynomial cost of P and Q.
        df_dPgQg = zeros((2 * self._ng, 1))        # w.r.t p.u. Pg and Qg
#            df_dPgQg[ipol] = matrix([g.poly_cost(xx[i], 1) for g in gpol])
#            for i, g in enumerate(gn):
#                der = polyder(list(g.p_cost))
#                df_dPgQg[i] = polyval(der, xx[i]) * base_mva
        for i in self._ipol:
            p_cost = list(self._gn[i].p_cost)
            df_dPgQg[i] = \
                self._base_mva * polyval(polyder(p_cost), xx[i])

        df = zeros((self._nxyz, 1))
        df[iPg] = df_dPgQg[:self._ng]
        df[iQg] = df_dPgQg[self._ng:self._ng + self._ng]

        # Piecewise linear cost of P and Q.
        df = df + self._ccost.T
        # TODO: Generalised cost term.

        return asarray(df).flatten()
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
745665
def _penalized_acquisition(self, x,  model, X_batch, r_x0, s_x0):
        '''
        
        '''
        fval = -self.acq.acquisition_function(x)[:,0]

        if self.transform=='softplus':
            fval_org = fval.copy()
            fval[fval_org>=40.] = np.log(fval_org[fval_org>=40.])
            fval[fval_org<40.] = np.log(np.log1p(np.exp(fval_org[fval_org<40.])))
        elif self.transform=='none':
            fval = np.log(fval+1e-50)

        fval = -fval
        if X_batch is not None:
            h_vals = self._hammer_function(x, X_batch, r_x0, s_x0)
            fval += -h_vals.sum(axis=-1)
        return fval
337141
def Kdiag(self, X, target):
        """"""
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU
        
        k1 = (2*lu+ly)/(lu+ly)**2
        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
642700
def derivatives(self, x, y, Rs, theta_Rs, r_core, center_x=0, center_y=0):
        

        rho0 = self._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs, r_core=r_core)

        if Rs < 0.0000001:
            Rs = 0.0000001
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_ ** 2 + y_ ** 2)

        dx, dy = self.coreBurkAlpha(R, Rs, rho0, r_core, x_, y_)

        return dx, dy
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.")
1209908
private void computeDir(double[] dir, double[] fg, QNInfo qn)
      throws SurpriseConvergence {
    System.arraycopy(fg, 0, dir, 0, fg.length);

    int mmm = qn.size();
    double[] as = new double[mmm];

    for (int i = mmm - 1; i >= 0; i--) {
      as[i] = qn.getRho(i) * ArrayMath.innerProduct(qn.getS(i), dir);
      plusAndConstMult(dir, qn.getY(i), -as[i], dir);
    }

    // multiply by hessian approximation
    qn.applyInitialHessian(dir);

    for (int i = 0; i < mmm; i++) {
      double b = qn.getRho(i) * ArrayMath.innerProduct(qn.getY(i), dir);
      plusAndConstMult(dir, qn.getS(i), as[i] - b, dir);
    }

    ArrayMath.multiplyInPlace(dir, -1);
  }