Run ModelAverage in default mode (QuicGraphicalLassoCV) to obtain proportion
    matrix.

    NOTE:  Only method = 'binary' really makes sense in this case.
771556
def adaptive_model_average(X, penalization, method):
    
    n_trials = 100
    print("Adaptive ModelAverage with:")
    print("   estimator: QuicGraphicalLasso (default)")
    print("   n_trials: {}".format(n_trials))
    print("   penalization: {}".format(penalization))
    print("   adaptive-method: {}".format(method))

    # if penalization is random, first find a decent scalar lam_ to build
    # random perturbation matrix around. lam doesn't matter for fully-random.
    lam = 0.5
    if penalization == "random":
        cv_model = QuicGraphicalLassoCV(
            cv=2, n_refinements=6, n_jobs=1, init_method="cov", score_metric=metric
        )
        cv_model.fit(X)
        lam = cv_model.lam_
        print("   lam: {}".format(lam))

    model = AdaptiveGraphicalLasso(
        estimator=ModelAverage(
            n_trials=n_trials, penalization=penalization, lam=lam, n_jobs=1
        ),
        method=method,
    )
    model.fit(X)
    lam_norm_ = np.linalg.norm(model.estimator_.lam_)
    print("   ||lam_||_2: {}".format(lam_norm_))
    return model.estimator_.covariance_, model.estimator_.precision_, lam_norm_
761617
def FitCBVs(model):
    '''
    

    '''

    # Get cbvs?
    if model.XCBV is None:
        GetTargetCBVs(model)

    # The number of CBVs to use
    ncbv = model.cbv_num

    # Need to treat short and long cadences differently
    if model.cadence == 'lc':

        # Loop over all the light curve segments
        m = [None for b in range(len(model.breakpoints))]
        weights = [None for b in range(len(model.breakpoints))]
        for b in range(len(model.breakpoints)):

            # Get the indices for this light curve segment
            inds = model.get_chunk(b, pad=False)
            masked_inds = model.get_masked_chunk(b, pad=False)

            # Regress
            mX = model.XCBV[masked_inds, :ncbv + 1]
            A = np.dot(mX.T, mX)
            B = np.dot(mX.T, model.flux[masked_inds])
            try:
                weights[b] = np.linalg.solve(A, B)
            except np.linalg.linalg.LinAlgError:
                # Singular matrix
                log.warn('Singular matrix!')
                weights[b] = np.zeros(mX.shape[1])

            m[b] = np.dot(model.XCBV[inds, :ncbv + 1], weights[b])

            # Vertical alignment
            if b == 0:
                m[b] -= np.nanmedian(m[b])
            else:
                # Match the first finite model point on either side of the
                # break
                # We could consider something more elaborate in the future
                i0 = -1 - np.argmax([np.isfinite(m[b - 1][-i])
                                     for i in range(1, len(m[b - 1]) - 1)])
                i1 = np.argmax([np.isfinite(m[b][i])
                                for i in range(len(m[b]))])
                m[b] += (m[b - 1][i0] - m[b][i1])

        # Join model and normalize
        m = np.concatenate(m)
        m -= np.nanmedian(m)

    else:

        # Interpolate over outliers so we don't have to worry
        # about masking the arrays below
        flux = Interpolate(model.time, model.mask, model.flux)

        # Get downbinned light curve
        newsize = len(model.time) // 30
        time = Downbin(model.time, newsize, operation='mean')
        flux = Downbin(flux, newsize, operation='mean')

        # Get LC breakpoints
        breakpoints = list(Breakpoints(
            model.ID, season=model.season, cadence='lc'))
        breakpoints += [len(time) - 1]

        # Loop over all the light curve segments
        m = [None for b in range(len(breakpoints))]
        weights = [None for b in range(len(breakpoints))]
        for b in range(len(breakpoints)):

            # Get the indices for this light curve segment
            M = np.arange(len(time))
            if b > 0:
                inds = M[(M > breakpoints[b - 1]) & (M <= breakpoints[b])]
            else:
                inds = M[M <= breakpoints[b]]

            # Regress
            A = np.dot(model.XCBV[inds, :ncbv + 1].T,
                       model.XCBV[inds, :ncbv + 1])
            B = np.dot(model.XCBV[inds, :ncbv + 1].T, flux[inds])
            weights[b] = np.linalg.solve(A, B)
            m[b] = np.dot(model.XCBV[inds, :ncbv + 1], weights[b])

            # Vertical alignment
            if b == 0:
                m[b] -= np.nanmedian(m[b])
            else:
                # Match the first finite model point on either side of the
                # break
                # We could consider something more elaborate in the future
                i0 = -1 - np.argmax([np.isfinite(m[b - 1][-i])
                                     for i in range(1, len(m[b - 1]) - 1)])
                i1 = np.argmax([np.isfinite(m[b][i])
                                for i in range(len(m[b]))])
                m[b] += (m[b - 1][i0] - m[b][i1])

        # Join model and normalize
        m = np.concatenate(m)
        m -= np.nanmedian(m)

        # Finally, interpolate back to short cadence
        m = np.interp(model.time, time, m)

    return m
1321343
def loads(s, cls=BinaryQuadraticModel, vartype=None):
    """"""
    return load(s.split('\n'), cls=cls, vartype=vartype)
505663
def bin_matrix(M, subsampling_factor=3):
    

    try:
        from scipy.sparse import issparse
        if issparse(M):
            return bin_sparse(M, subsampling_factor=subsampling_factor)
        else:
            raise ImportError
    except ImportError:
        return bin_dense(M, subsampling_factor=subsampling_factor)
586205
def get_checkerboard_matrix(kernel_width, kernel_type="default", gaussian_param=0.1):

    

    if kernel_type is "gaussian":
        return get_gaussian_kernel(kernel_width, gaussian_param)

    if kernel_type is "default":
        return np.vstack((
            np.hstack((
                -1 * np.ones((kernel_width, kernel_width)), np.ones((kernel_width, kernel_width))
            )),
            np.hstack((
                np.ones((kernel_width, kernel_width)), -1 * np.ones((kernel_width, kernel_width))
            ))
        ))

    elif kernel_type is "bottom_right":
        return np.vstack((
            np.hstack((
                np.ones((kernel_width, kernel_width)), np.ones((kernel_width, kernel_width))
            )),
            np.hstack((
                np.ones((kernel_width, kernel_width)), -1 * np.ones((kernel_width, kernel_width))
            ))
        ))

    elif kernel_type is "top_left":
        return np.vstack((
            np.hstack((
                -1 * np.ones((kernel_width, kernel_width)), np.ones((kernel_width, kernel_width))
            )),
            np.hstack((
                np.ones((kernel_width, kernel_width)), np.ones((kernel_width, kernel_width))
            ))
        ))
74658
function(ctx, noTransform) {
      ctx.save();
      var m = this.transformMatrix;
      this._resetWidthHeight();
      if (this.group) {
        ctx.translate(-this.group.width/2 + this.width/2, -this.group.height/2 + this.height/2);
      }
      if (m) {
        ctx.transform(m[0], m[1], m[2], m[3], m[4], m[5]);
      }
      if (!noTransform) {
        this.transform(ctx);
      }
      this._render(ctx);
      if (this.active && !noTransform) {
        this.drawBorders(ctx);
        this.hideCorners || this.drawCorners(ctx);
      }
      ctx.restore();
    }
1069770
def model(args):
    
    from scipy.stats import binom, poisson

    p = OptionParser(model.__doc__)
    p.add_option("-k", default=23, type="int", help="Kmer size")
    p.add_option("--cov", default=50, type="int", help="Expected coverage")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    erate, = args
    erate = float(erate)
    cov = opts.cov
    k = opts.k

    xy = []
    # Range include c although it is unclear what it means to have c=0
    for c in xrange(0,  cov * 2 + 1):
        Prob_Yk = 0
        for i in xrange(k + 1):
            # Probability of having exactly i errors
            pi_i = binom.pmf(i, k, erate)
            # Expected coverage of kmer with exactly i errors
            mu_i = cov * (erate / 3) ** i * (1 - erate) ** (k - i)
            # Probability of seeing coverage of c
            Prob_Yk_i = poisson.pmf(c, mu_i)
            # Sum i over 0, 1, ... up to k errors
            Prob_Yk += pi_i * Prob_Yk_i
        xy.append((c, Prob_Yk))

    x, y = zip(*xy)
    asciiplot(x, y, title="Model")
1292215
def load_model():
    
    direc = "/home/annaho/TheCannon/code/lamost/mass_age/cn"
    m = model.CannonModel(2)
    m.coeffs = np.load(direc + "/coeffs.npz")['arr_0'][0:3626,:] # no cols
    m.scatters = np.load(direc + "/scatters.npz")['arr_0'][0:3626] # no cols
    m.chisqs = np.load(direc + "/chisqs.npz")['arr_0'][0:3626] # no cols
    m.pivots = np.load(direc + "/pivots.npz")['arr_0']
    return m
575210
def draw(self):
        

        # Perform display related manipulations on the confusion matrix data
        cm_display = self.confusion_matrix_

        # Convert confusion matrix to percent of each row, i.e. the
        # predicted as a percent of true in each class.
        if self.percent == True:
            # Note: div_safe function returns 0 instead of NAN.
            cm_display = div_safe(self.confusion_matrix_, self.class_counts_.reshape(-1,1))
            cm_display = np.round(cm_display* 100, decimals=0)

        # Y axis should be sorted top to bottom in pcolormesh
        cm_display = cm_display[::-1,::]

        # Set up the dimensions of the pcolormesh
        n_classes = len(self.classes_)
        X, Y = np.arange(n_classes+1), np.arange(n_classes+1)
        self.ax.set_ylim(bottom=0, top=cm_display.shape[0])
        self.ax.set_xlim(left=0, right=cm_display.shape[1])

        # Fetch the grid labels from the classes in correct order; set ticks.
        xticklabels = self.classes_
        yticklabels = self.classes_[::-1]
        ticks = np.arange(n_classes) + 0.5

        self.ax.set(xticks=ticks, yticks=ticks)
        self.ax.set_xticklabels(xticklabels, rotation="vertical", fontsize=self.fontsize)
        self.ax.set_yticklabels(yticklabels, fontsize=self.fontsize)

        # Set data labels in the grid enumerating over all x,y class pairs.
        # NOTE: X and Y are one element longer than the confusion matrix, so
        # skip the last element in the enumeration to label grids.
        for x in X[:-1]:
            for y in Y[:-1]:

                # Extract the value and the text label
                value = cm_display[x,y]
                svalue = "{:0.0f}".format(value)
                if self.percent:
                    svalue += "%"

                # Determine the grid and text colors
                base_color = self.cmap(value / cm_display.max())
                text_color = find_text_color(base_color)

                # Make zero values more subtle
                if cm_display[x,y] == 0:
                    text_color = CMAP_MUTEDCOLOR

                # Add the label to the middle of the grid
                cx, cy = x+0.5, y+0.5
                self.ax.text(
                    cy, cx, svalue, va='center', ha='center',
                    color=text_color, fontsize=self.fontsize,
                )

                # Add a dark line on the grid with the diagonal. Note that the
                # tick labels have already been reversed.
                lc = 'k' if xticklabels[x] == yticklabels[y] else 'w'
                self._edgecolors.append(lc)


        # Draw the heatmap with colors bounded by vmin,vmax
        vmin = 0.00001
        vmax = 99.999 if self.percent == True else cm_display.max()
        self.ax.pcolormesh(
            X, Y, cm_display, vmin=vmin, vmax=vmax,
            edgecolor=self._edgecolors, cmap=self.cmap, linewidth='0.01'
        )

        # Return the axes being drawn on
        return self.ax
904306
def doblob(morphed, blobdet, img, anno=True):
    
    keypoints = blobdet.detect(morphed)
    nkey = len(keypoints)
    kpsize = asarray([k.size for k in keypoints])
    final = img.copy()  # is the .copy necessary?

    final = cv2.drawKeypoints(img, keypoints, outImage=final,
                              flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
# %% plot count of blobs
    if anno:
        cv2.putText(final, text=str(nkey), org=(int(img.shape[1]*.9), 25),
                    fontFace=cv2.FONT_HERSHEY_PLAIN, fontScale=2,
                    color=(0, 255, 0), thickness=2)

    return final, nkey, kpsize
421671
private void ar2gvr(float scale) {
        // Real world scale
        Matrix.scaleM(mPoseMatrix, 0, scale, scale, scale);
        mPoseMatrix[12] = mPoseMatrix[12] * scale;
        mPoseMatrix[13] = mPoseMatrix[13] * scale;
        mPoseMatrix[14] = mPoseMatrix[14] * scale;
    }
551010
def apply_model(refwcs):
    
    # apply distortion model to CD matrix
    if 'ocx10' in refwcs.__dict__ and refwcs.ocx10 is not None:
        linmat = np.array([[refwcs.ocx11,refwcs.ocx10],[refwcs.ocy11,refwcs.ocy10]])/refwcs.idcscale
        refwcs.wcs.cd = np.dot(refwcs.wcs.cd,linmat)

        refwcs.wcs.set()
        refwcs.setOrient()
        refwcs.setPscale()
720996
function getShapeOne(x,y,options,modules,mod_matrix){

  var EOL = '\r\n';
  var width = options.width;
  var height = options.height;
  var length = modules.length;
  var style = options.style;
  var xsize = width / (length + 2 * options.padding);
  var ysize = height / (length + 2 * options.padding);
  var px = (x * xsize + options.padding * xsize);
  var py = (y * ysize + options.padding * ysize);

  var orientation = 0;
  if(mod_matrix.right){orientation=90;}
  else if(mod_matrix.bottom){orientation=180;}
  else if(mod_matrix.left){orientation=270;}

  //Returning a 1b
  switch(options.style){
    case 'square':
      return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')"  x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
    case 'circle':
      return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')"  rx="'+xsize+'" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
    case 'sharp-edge':
      return '<path transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')" d="M'+px+' '+py+' L'+(px+xsize)+' '+py+' L'+(px+(xsize/2))+' '+(py+ysize)+' Z" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
  }
}
1129286
def to_bqm(self, model):
        
        linear = ((v, float(model.get_py_value(bias)))
                  for v, bias in self.linear.items())
        quadratic = ((u, v, float(model.get_py_value(bias)))
                     for (u, v), bias in self.quadratic.items())
        offset = float(model.get_py_value(self.offset))

        return dimod.BinaryQuadraticModel(linear, quadratic, offset, dimod.SPIN)
525376
def _hexplot(matrix, fig, colormap):
    
    umatrix_min = matrix.min()
    umatrix_max = matrix.max()
    n_rows, n_columns = matrix.shape
    cmap = plt.get_cmap(colormap)
    offsets = np.zeros((n_columns * n_rows, 2))
    facecolors = []
    for row in range(n_rows):
        for col in range(n_columns):
            if row % 2 == 0:
                offsets[row * n_columns + col] = [col +
                                                  0.5, 2 * n_rows - 2 * row]
                facecolors.append(cmap((matrix[row, col] - umatrix_min) /
                                       (umatrix_max) * 255))
            else:
                offsets[row * n_columns + col] = [col, 2 * n_rows - 2 * row]
                facecolors.append(cmap((matrix[row, col] - umatrix_min) /
                                       (umatrix_max) * 255))
    polygon = np.zeros((6, 2), float)
    polygon[:, 0] = 1.1 * np.array([0.5, 0.5, 0.0, -0.5, -0.5, 0.0])
    polygon[:, 1] = 1.1 * np.array([-np.sqrt(3) / 6, np.sqrt(3) / 6,
                                    np.sqrt(3) / 2 + np.sqrt(3) / 6,
                                    np.sqrt(3) / 6, -np.sqrt(3) / 6,
                                    -np.sqrt(3) / 2 - np.sqrt(3) / 6])
    polygons = np.expand_dims(polygon, 0) + np.expand_dims(offsets, 1)
    ax = fig.gca()
    collection = mcoll.PolyCollection(
        polygons,
        offsets=offsets,
        facecolors=facecolors,
        edgecolors=facecolors,
        linewidths=1.0,
        offset_position="data")
    ax.add_collection(collection, autolim=False)
    corners = ((-0.5, -0.5), (n_columns + 0.5, 2 * n_rows + 0.5))
    ax.update_datalim(corners)
    ax.autoscale_view(tight=True)
    return offsets
417841
private static void balbak(DenseMatrix V, double[] scale) {
        int n = V.nrows();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                V.mul(i, j, scale[i]);
            }
        }
    }
521521
def cublasSgbmv(handle, trans, m, n, kl, ku, alpha, A, lda,
                x, incx, beta, y, incy):
    

    status = _libcublas.cublasSgbmv_v2(handle,
                                       trans, m, n, kl, ku,
                                       ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda,
                                       int(x), incx,
                                       ctypes.byref(ctypes.c_float(beta)),
                                       int(y), incy)
    cublasCheckStatus(status)