Efficient conversion from scipy sparse matrix to cvxopt sparse matrix
486847
def scipy_sparse_to_spmatrix(A):
    """"""
    coo = A.tocoo()
    SP = spmatrix(coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=A.shape)
    return SP
900762
def _create_objective(self, m, n):
        
        # S is the completed matrix
        shape = (m, n)
        S = cvxpy.Variable(shape, name="S")
        norm = cvxpy.norm(S, "nuc")
        objective = cvxpy.Minimize(norm)
        return S, objective
1273455
def _ncc_c_3dim(x, y):
    
    den = norm(x, axis=1)[:, None] * norm(y, axis=1)
    den[den == 0] = np.Inf
    x_len = x.shape[-1]
    fft_size = 1 << (2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size))[:, None])
    cc = np.concatenate((cc[:,:,-(x_len-1):], cc[:,:,:x_len]), axis=2)
    return np.real(cc) / den.T[:, :, None]
611029
def decompose_covariance(c):
    

    # make it a kickass copy of the original
    c = _n.array(c)

    # first get the error vector
    e = []
    for n in range(0, len(c[0])): e.append(_n.sqrt(c[n][n]))

    # now cycle through the matrix, dividing by e[1]*e[2]
    for n in range(0, len(c[0])):
        for m in range(0, len(c[0])):
            c[n][m] = c[n][m] / (e[n]*e[m])

    return [_n.array(e), _n.array(c)]
1343145
def _row_sparse_ndarray_from_definition(data, indices, shape=None, ctx=None,
                                        dtype=None, indices_type=None):
    """"""
    storage_type = 'row_sparse'
    # context
    ctx = current_context() if ctx is None else ctx
    # types
    dtype = _prepare_default_dtype(data, dtype)
    indices_type = _STORAGE_AUX_TYPES[storage_type][0] if indices_type is None else indices_type
    # prepare src array and types
    data = _prepare_src_array(data, dtype)
    indices = _prepare_src_array(indices, indices_type)

    # TODO(junwu): Convert data, indptr, and indices to mxnet NDArrays
    # if they are not for now. In the future, we should provide a c-api
    # to accept np.ndarray types to copy from to result.data and aux_data
    if not isinstance(data, NDArray):
        data = _array(data, ctx, dtype)
    if not isinstance(indices, NDArray):
        indices = _array(indices, ctx, indices_type)
    if shape is None:
        num_indices = indices.shape[0]
        if num_indices == 0:
            raise ValueError('invalid shape')
        dim0 = indices[num_indices - 1].asscalar() + 1
        shape = (dim0, ) + data.shape[1:]
    # verify shapes
    if data.ndim != len(shape) or indices.ndim != 1 or np.prod(shape[1:]) == 0:
        raise ValueError("invalid shape")
    result = RowSparseNDArray(_new_alloc_handle(storage_type, shape, ctx, False, dtype,
                                                [indices_type], [indices.shape]))
    check_call(_LIB.MXNDArraySyncCopyFromNDArray(result.handle, data.handle, ctypes.c_int(-1)))
    check_call(_LIB.MXNDArraySyncCopyFromNDArray(result.handle, indices.handle, ctypes.c_int(0)))
    return result
1218098
def cublasDtpsv(handle, uplo, trans, diag, n, AP, x, incx):
    

    status = _libcublas.cublasDtpsv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, int(AP), int(x), incx)
    cublasCheckStatus(status)
17299
def spm(x, model, y=None, n_samples=None, dx_min=-0.1,
        dx_max=0.1, n_dxs=5, dy_min=-0.1, dy_max=0.1, n_dys=5,
        angle_min=-30, angle_max=30, n_angles=31, black_border_size=0):
  
  if y is None:
    preds = model.get_probs(x)
    # Using model predictions as ground truth to avoid label leaking
    preds_max = reduce_max(preds, 1, keepdims=True)
    y = tf.to_float(tf.equal(preds, preds_max))
    y = tf.stop_gradient(y)
    del preds
  y = y / reduce_sum(y, 1, keepdims=True)

  # Define the range of transformations
  dxs = np.linspace(dx_min, dx_max, n_dxs)
  dys = np.linspace(dy_min, dy_max, n_dys)
  angles = np.linspace(angle_min, angle_max, n_angles)

  if n_samples is None:
    import itertools
    transforms = list(itertools.product(*[dxs, dys, angles]))
  else:
    sampled_dxs = np.random.choice(dxs, n_samples)
    sampled_dys = np.random.choice(dys, n_samples)
    sampled_angles = np.random.choice(angles, n_samples)
    transforms = zip(sampled_dxs, sampled_dys, sampled_angles)
  transformed_ims = parallel_apply_transformations(
      x, transforms, black_border_size)

  def _compute_xent(x):
    preds = model.get_logits(x)
    return tf.nn.softmax_cross_entropy_with_logits_v2(
        labels=y, logits=preds)

  all_xents = tf.map_fn(
      _compute_xent,
      transformed_ims,
      parallel_iterations=1)  # Must be 1 to avoid keras race conditions

  # Return the adv_x with worst accuracy

  # all_xents is n_total_samples x batch_size (SB)
  all_xents = tf.stack(all_xents)  # SB

  # We want the worst case sample, with the largest xent_loss
  worst_sample_idx = tf.argmax(all_xents, axis=0)  # B

  batch_size = tf.shape(x)[0]
  keys = tf.stack([
      tf.range(batch_size, dtype=tf.int32),
      tf.cast(worst_sample_idx, tf.int32)
  ], axis=1)
  transformed_ims_bshwc = tf.einsum('sbhwc->bshwc', transformed_ims)
  after_lookup = tf.gather_nd(transformed_ims_bshwc, keys)  # BHWC
  return after_lookup
1107735
def get_scipy_eig_vec(self):
    
    if not self.params['has_conv']:
      matrix_m = self.sess.run(self.dual_object.matrix_m)
      min_eig_vec_val, estimated_eigen_vector = eigs(matrix_m, k=1, which='SR',
                                                     tol=1E-4)
      min_eig_vec_val = np.reshape(np.real(min_eig_vec_val), [1, 1])
      return np.reshape(estimated_eigen_vector, [-1, 1]), min_eig_vec_val
    else:
      dim = self.dual_object.matrix_m_dimension
      input_vector = tf.placeholder(tf.float32, shape=(dim, 1))
      output_vector = self.dual_object.get_psd_product(input_vector)

      def np_vector_prod_fn(np_vector):
        np_vector = np.reshape(np_vector, [-1, 1])
        output_np_vector = self.sess.run(output_vector, feed_dict={input_vector:np_vector})
        return output_np_vector
      linear_operator = LinearOperator((dim, dim), matvec=np_vector_prod_fn)
      # Performing shift invert scipy operation when eig val estimate is available
      min_eig_vec_val, estimated_eigen_vector = eigs(linear_operator,
                                                     k=1, which='SR', tol=1E-4)
      min_eig_vec_val = np.reshape(np.real(min_eig_vec_val), [1, 1])
      return np.reshape(estimated_eigen_vector, [-1, 1]), min_eig_vec_val
255886
def getoptT(X, W, Y, Z, S, M_E, E, m0, rho):
    ''' 
    '''
    iter_max = 20
    norm2WZ = np.linalg.norm(W, ord='fro')**2 + np.linalg.norm(Z, ord='fro')**2
    f = np.zeros(iter_max + 1)
    f[0] = F_t(X, Y, S, M_E, E, m0, rho)

    t = -1e-1
    for i in range(iter_max):
        f[i + 1] = F_t(X + t * W, Y + t * Z, S, M_E, E, m0, rho)

        if f[i + 1] - f[0] <= 0.5 * t * norm2WZ:
            return t
        t /= 2
    return t
529293
def cross_validate(self, ax):
        '''
        

        '''

        # The CDPP to beat
        cdpp_opt = self.get_cdpp_arr()

        # Loop over all chunks
        for b, brkpt in enumerate(self.breakpoints):

            log.info("Cross-validating chunk %d/%d..." %
                     (b + 1, len(self.breakpoints)))

            # Mask for current chunk
            m = self.get_masked_chunk(b)

            # Mask transits and outliers
            time = self.time[m]
            flux = self.fraw[m]
            ferr = self.fraw_err[m]
            med = np.nanmedian(self.fraw)

            # Setup the GP
            gp = GP(self.kernel, self.kernel_params, white=False)
            gp.compute(time, ferr)

            # The masks
            masks = list(Chunks(np.arange(0, len(time)),
                                len(time) // self.cdivs))

            # The pre-computed matrices
            pre_v = [self.cv_precompute(mask, b) for mask in masks]

            # Initialize with the nPLD solution
            log_lam_opt = np.log10(self.lam[b])
            scatter_opt = self.validation_scatter(
                log_lam_opt, b, masks, pre_v, gp, flux, time, med)
            log.info("Iter 0/%d: " % (self.piter) +
                     "logL = (%s), s = %.3f" %
                     (", ".join(["%.3f" % l for l in log_lam_opt]),
                      scatter_opt))

            # Do `piter` iterations
            for p in range(self.piter):

                # Perturb the initial condition a bit
                log_lam = np.array(
                    np.log10(self.lam[b])) * \
                    (1 + self.ppert * np.random.randn(len(self.lam[b])))
                scatter = self.validation_scatter(
                    log_lam, b, masks, pre_v, gp, flux, time, med)
                log.info("Initializing at: " +
                         "logL = (%s), s = %.3f" %
                         (", ".join(["%.3f" % l for l in log_lam]), scatter))

                # Call the minimizer
                log_lam, scatter, _, _, _, _ = \
                    fmin_powell(self.validation_scatter, log_lam,
                                args=(b, masks, pre_v, gp, flux, time, med),
                                maxfun=self.pmaxf, disp=False,
                                full_output=True)

                # Did it improve the CDPP?
                tmp = np.array(self.lam[b])
                self.lam[b] = 10 ** log_lam
                self.compute()
                cdpp = self.get_cdpp_arr()[b]
                self.lam[b] = tmp
                if cdpp < cdpp_opt[b]:
                    cdpp_opt[b] = cdpp
                    log_lam_opt = log_lam

                # Log it
                log.info("Iter %d/%d: " % (p + 1, self.piter) +
                         "logL = (%s), s = %.3f" %
                         (", ".join(["%.3f" % l for l in log_lam]), scatter))

            # The best solution
            log.info("Found minimum: logL = (%s), s = %.3f" %
                     (", ".join(["%.3f" % l for l in log_lam_opt]),
                      scatter_opt))
            self.lam[b] = 10 ** log_lam_opt

        # We're just going to plot lambda as a function of chunk number
        bs = np.arange(len(self.breakpoints))
        color = ['k', 'b', 'r', 'g', 'y']
        for n in range(self.pld_order):
            ax[0].plot(bs + 1, [np.log10(self.lam[b][n])
                                for b in bs], '.', color=color[n])
            ax[0].plot(bs + 1, [np.log10(self.lam[b][n])
                                for b in bs], '-', color=color[n], alpha=0.25)
        ax[0].set_ylabel(r'$\log\Lambda$', fontsize=5)
        ax[0].margins(0.1, 0.1)
        ax[0].set_xticks(np.arange(1, len(self.breakpoints) + 1))
        ax[0].set_xticklabels([])

        # Now plot the CDPP
        cdpp_arr = self.get_cdpp_arr()
        ax[1].plot(bs + 1, cdpp_arr, 'b.')
        ax[1].plot(bs + 1, cdpp_arr, 'b-', alpha=0.25)
        ax[1].margins(0.1, 0.1)
        ax[1].set_ylabel(r'Scatter (ppm)', fontsize=5)
        ax[1].set_xlabel(r'Chunk', fontsize=5)
        ax[1].set_xticks(np.arange(1, len(self.breakpoints) + 1))
433757
def build_constraint_matrix(constrs)
      cbeg = []
      cind = []
      cval = []
      constrs.each.map do |constr|
        cbeg << cind.length
        constr.expression.terms.each do |var, coeff|
          cind << var.index
          cval << coeff
        end
      end

      [cbeg, cind, cval]
    end
1287251
def _run_opf(self, P, q, AA, ll, uu, xmin, xmax, x0, opt):
        
        AAcvx = tocvx(AA)
        nx = x0.shape[0] # number of variables
        # add var limits to linear constraints
        eyex = spmatrix(1.0, range(nx), range(nx))
        A = sparse([eyex, AAcvx])
        l = matrix([-xmin, ll])
        u = matrix([ xmax, uu])

        Ae, be, Ai, bi = split_linear_constraints(A, l, u)


        if len(P.V) > 0:
            cvx_sol = solvers.qp(P, q, Ai, bi, Ae, be, self.solver, {"x": x0})
        else:
            cvx_sol = solvers.lp(q, Ai, bi, Ae, be, self.solver, {"x": x0})

        return cvx_sol
916154
def choi(self) -> bk.BKTensor:
        """"""
        # Put superop axes in [ok, ib, ob, ik] and reshape to matrix
        N = self.qubit_nb
        return bk.reshape(self.sharp.tensor, [2**(N*2)] * 2)
657239
def xstep(self):
        
        """

        self.YU[:] = self.Y - self.U
        self.block_sep0(self.YU)[:] += self.S
        Zf = sl.rfftn(self.YU, None, self.cri.axisN)
        Z0f = self.block_sep0(Zf)
        Z1f = self.block_sep1(Zf)

        DZ0f = np.conj(self.Df) * Z0f
        DZ0fBQ = sl.dot(self.B.dot(self.Q).T, DZ0f, axis=self.cri.axisC)
        Z1fQ = sl.dot(self.Q.T, Z1f, axis=self.cri.axisC)
        b = DZ0fBQ + Z1fQ

        Xh = sl.solvedbd_sm(self.gDf, (self.mu / self.rho) * self.GHGf + 1.0,
                            b, self.c, axis=self.cri.axisM)
        self.Xf[:] = sl.dot(self.Q, Xh, axis=self.cri.axisC)
        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            DDXf = np.conj(self.Df) *  sl.inner(self.Df, self.Xf,
                                                axis=self.cri.axisM)
            DDXfBB = sl.dot(self.B.T.dot(self.B), DDXf, axis=self.cri.axisC)
            ax = self.rho * (DDXfBB + self.Xf) + \
                 self.mu * self.GHGf * self.Xf
            b = self.rho * (sl.dot(self.B.T, DZ0f, axis=self.cri.axisC)
                            + Z1f)
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
203750
def __push_facvar_sparse(self, polynomial, block_index, row_offset, i, j):
        
        width = self.block_struct[block_index - 1]
        # Preprocess the polynomial for uniform handling later
        # DO NOT EXPAND THE POLYNOMIAL HERE!!!!!!!!!!!!!!!!!!!
        # The simplify_polynomial bypasses the problem.
        # Simplifying here will trigger a bug in SymPy related to
        # the powers of daggered variables.
        # polynomial = polynomial.expand()
        if is_number_type(polynomial) or polynomial.is_Mul:
            elements = [polynomial]
        else:
            elements = polynomial.as_coeff_mul()[1][0].as_coeff_add()[1]
        # Identify its constituent monomials
        for element in elements:
            results = self._get_index_of_monomial(element)
            # k identifies the mapped value of a word (monomial) w
            for (k, coeff) in results:
                if k > -1 and coeff != 0:
                    self.F[row_offset + i * width + j, k] += coeff
14790
def getoptS(X, Y, M_E, E):
    ''' 
    '''
    n, r = X.shape
    C = np.dot(np.dot(X.T, M_E), Y)
    C = C.flatten()

    A = np.zeros((r * r, r * r))
    for i in range(r):
        for j in range(r):
            ind = j * r + i
            temp = np.dot(
                np.dot(X.T, np.dot(X[:, i, None], Y[:, j, None].T) * E), Y)
            A[:, ind] = temp.flatten()

    S = np.linalg.solve(A, C)

    return np.reshape(S, (r, r)).T
1361461
def eighCovarianceMatrix(cfile):
    
    # precompute eigenvalue decomposition
    K = np.loadtxt(cfile+'.cov')
    K+= 1e-4*sp.eye(K.shape[0])
    S,U = la.eigh(K); S=S[::-1]; U=U[:,::-1]
    np.savetxt(cfile+'.cov.eval',S,fmt='%.6f')
    np.savetxt(cfile+'.cov.evec',U,fmt='%.6f')