Run QuicGraphicalLasso with mode='path' and gamma; use EBIC criteria for model
    selection.

    The EBIC criteria is built into InverseCovarianceEstimator base class
    so we demonstrate those utilities here.
481524
def quic_graph_lasso_ebic_manual(X, gamma=0):
    
    print("QuicGraphicalLasso (manual EBIC) with:")
    print("   mode: path")
    print("   gamma: {}".format(gamma))
    model = QuicGraphicalLasso(
        lam=1.0,
        mode="path",
        init_method="cov",
        path=np.logspace(np.log10(0.01), np.log10(1.0), num=100, endpoint=True),
    )
    model.fit(X)
    ebic_index = model.ebic_select(gamma=gamma)
    covariance_ = model.covariance_[ebic_index]
    precision_ = model.precision_[ebic_index]
    lam_ = model.lam_at_index(ebic_index)
    print("   len(path lams): {}".format(len(model.path_)))
    print("   lam_scale_: {}".format(model.lam_scale_))
    print("   lam_: {}".format(lam_))
    print("   ebic_index: {}".format(ebic_index))
    return covariance_, precision_, lam_
975103
def maxcut_qaoa(
        graph,
        steps=DEFAULT_STEPS,
        learning_rate=LEARNING_RATE,
        verbose=False):
    """"""

    if not isinstance(graph, nx.Graph):
        graph = nx.from_edgelist(graph)

    init_scale = 0.01
    init_bias = 0.5

    init_beta = normal(loc=init_bias, scale=init_scale, size=[steps])
    init_gamma = normal(loc=init_bias, scale=init_scale, size=[steps])

    beta = tf.get_variable('beta', initializer=init_beta)
    gamma = tf.get_variable('gamma', initializer=init_gamma)

    circ = qubo_circuit(graph, steps, beta, gamma)
    cuts = graph_cuts(graph)
    maxcut = cuts.max()
    expect = circ.run().expectation(cuts)
    loss = - expect

    # === Optimization ===
    opt = tf.train.GradientDescentOptimizer(learning_rate=learning_rate)
    train = opt.minimize(loss, var_list=[beta, gamma])

    with tf.Session() as sess:
        init_op = tf.global_variables_initializer()
        sess.run(init_op)

        block = 10
        min_difference = 0.0001
        last_ratio = -1
        for step in range(0, MAX_OPT_STEPS, block):
            for _ in range(block):
                sess.run(train)
            ratio = sess.run(expect) / maxcut
            if ratio - last_ratio < min_difference:
                break
            last_ratio = ratio

            if verbose:
                print("# step: {}  ratio: {:.4f}%".format(step, ratio))

        opt_beta = sess.run(beta)
        opt_gamma = sess.run(gamma)

    return ratio, opt_beta, opt_gamma
646328
@Override
	public boolean process(PairwiseImageGraph pairwiseGraph ) {
		this.graph = new MetricSceneGraph(pairwiseGraph);
		for (int i = 0; i < graph.edges.size(); i++) {
			decomposeEssential(graph.edges.get(i));
		}

		declareModelFitting();

		for (int i = 0; i < graph.edges.size(); i++) {
			Motion e = graph.edges.get(i);
			e.triangulationAngle = medianTriangulationAngle(e);
		}

		if( verbose != null )
			verbose.println("Selecting root");
		// Select the view which will act as the origin
		View origin = selectOriginNode();
		// Select the motion which will define the coordinate system
		Motion baseMotion = selectCoordinateBase( origin );

		this.graph.sanityCheck();

		if( verbose != null )
			verbose.println("Stereo triangulation");
		// Triangulate features in all motions which exceed a certain angle
		for (int i = 0; i < graph.edges.size() && !stopRequested ; i++) {
			Motion e = graph.edges.get(i);
			if( e.triangulationAngle > Math.PI/10 || e == baseMotion) {
				triangulateStereoEdges(e);
				if( verbose != null ) {
					int a = e.viewSrc.index;
					int b = e.viewDst.index;
					verbose.println("   Edge[" + i + "] "+a+"->"+b+"  feat3D="+e.stereoTriangulations.size());
				}
			}
		}
		if( stopRequested )
			return false;

		if( verbose != null )
			verbose.println("Defining the coordinate system");

		// Using the selecting coordinate frames and triangulated points define the coordinate system
		defineCoordinateSystem(origin, baseMotion);
		if( stopRequested )
			return false;

		if( verbose != null )
			verbose.println("Estimate all features");

		// Now estimate all the other view locations and 3D features
		estimateAllFeatures(origin, baseMotion.destination(origin));
		if( stopRequested )
			return false;

		// Convert the graph into the output format
		convertToOutput(origin);

		return viewsAdded.size() >= 2;
	}
1109785
def discriminator1(ndf, no_bias=True, fix_gamma=True, eps=1e-5 + 1e-12):
    ''''''
    BatchNorm = mx.sym.BatchNorm

    data = mx.sym.Variable('data')

    d1 = mx.sym.Convolution(data, name='d1', kernel=(5,5), stride=(2,2), pad=(2,2), num_filter=ndf, no_bias=no_bias)
    dact1 = mx.sym.LeakyReLU(d1, name='dact1', act_type='leaky', slope=0.2)

    d2 = mx.sym.Convolution(dact1, name='d2', kernel=(5,5), stride=(2,2), pad=(2,2), num_filter=ndf*2, no_bias=no_bias)
    dbn2 = BatchNorm(d2, name='dbn2', fix_gamma=fix_gamma, eps=eps)
    dact2 = mx.sym.LeakyReLU(dbn2, name='dact2', act_type='leaky', slope=0.2)

    d3 = mx.sym.Convolution(dact2, name='d3', kernel=(5,5), stride=(2,2), pad=(2,2), num_filter=ndf*4, no_bias=no_bias)
    dbn3 = BatchNorm(d3, name='dbn3', fix_gamma=fix_gamma, eps=eps)
    dact3 = mx.sym.LeakyReLU(dbn3, name='dact3', act_type='leaky', slope=0.2)

    return dact3
77313
def visualize_model(X, y, estimator, path, **kwargs):
     
    y = LabelEncoder().fit_transform(y)
    model = Pipeline([
         ('one_hot_encoder', OneHotEncoder()), 
         ('estimator', estimator)
    ])

    _, ax = plt.subplots()

    # Instantiate the classification model and visualizer
    visualizer = ClassificationReport(
        model, classes=['edible', 'poisonous'], 
        cmap="YlGn", size=(600, 360), ax=ax, **kwargs
    )
    visualizer.fit(X, y)  
    visualizer.score(X, y)
    visualizer.poof(outpath=path)
947423
def isomap(self, num_dims=None, directed=None):
    '''
    '''
    W = -0.5 * self.shortest_path(directed=directed) ** 2
    kpca = KernelPCA(n_components=num_dims, kernel='precomputed')
    return kpca.fit_transform(W)
1298652
def run_false_positive_experiment_dim(
                              numActive = 128,
                              dim = 500,
                              numSamples = 1000,
                              numDendrites = 500,
                              synapses = 24,
                              numTrials = 10000,
                              seed = 42,
                              nonlinearity = sigmoid_nonlinearity(11.5, 5)):
  
  numpy.random.seed(seed)
  fps = []
  fns = []
  totalUnclassified = 0

  for trial in range(numTrials):

    # data = generate_evenly_distributed_data_sparse(dim = dim,
    #                                                num_active = numActive,
    #                                                num_samples = numSamples)
    # labels = numpy.asarray([1 for i in range(numSamples / 2)] +
    #                        [-1 for i in range(numSamples / 2)])
    # flipped_labels = labels * -1

    negData = generate_evenly_distributed_data_sparse(dim = dim,
                                                   num_active = numActive,
                                                   num_samples = numSamples/2)
    posData = generate_evenly_distributed_data_sparse(dim = dim,
                                                   num_active = numActive,
                                                   num_samples = numSamples/2)
    halfLabels = numpy.asarray([1 for _ in range(numSamples / 2)])
    flippedHalfLabels = halfLabels * -1

    neuron = Neuron(size =synapses * numDendrites,
                    num_dendrites = numDendrites,
                    dendrite_length = synapses,
                    dim = dim, nonlinearity = nonlinearity)
    neg_neuron = Neuron(size =synapses * numDendrites,
                        num_dendrites = numDendrites,
                        dendrite_length = synapses,
                        dim = dim, nonlinearity = nonlinearity)

    neuron.HTM_style_initialize_on_positive_data(posData)
    neg_neuron.HTM_style_initialize_on_positive_data(negData)

    # Get error for positively labeled data
    fp, fn, uc = get_error(posData, halfLabels, [neuron], [neg_neuron])
    totalUnclassified += uc
    fps.append(fp)
    fns.append(fn)

    # Get error for negatively labeled data
    fp, fn, uc = get_error(negData, flippedHalfLabels, [neuron], [neg_neuron])
    totalUnclassified += uc
    fps.append(fp)
    fns.append(fn)


  print "Error with n = {} : {} FP, {} FN, {} unclassified".format(
    dim, sum(fps), sum(fns), totalUnclassified)

  result = {
    "dim": dim,
    "totalFP": sum(fps),
    "totalFN": sum(fns),
    "total mistakes": sum(fns + fps) + totalUnclassified,
    "error": float(sum(fns + fps) + totalUnclassified) / (numTrials * numSamples),
    "totalSamples": numTrials * numSamples,
    "a": numActive,
    "num_dendrites": numDendrites,
    "totalUnclassified": totalUnclassified,
    "synapses": 24,
    "seed": seed,
  }

  return result
401926
def _init_parameters_jaro(self, X_bin):
        

        _, n_features = X_bin.shape

        class_prior = [.9, .1]
        feature_prob = np.zeros((2, n_features))

        for i, bin in enumerate(self._binarizers):

            if bin.classes_.shape[0] != 2:
                raise ValueError("Only binary labels are allowed for "
                                 "'jaro'method. "
                                 "Column {} has {} different labels.".format(
                                     i, bin.classes_.shape[0]))

            # TODO: ensure classes are [0, 1] (not [1, 0])
            # TODO: check with bin.y_type_

            feature_prob[0, :] = np.tile([.9, .1], int(n_features / 2))
            feature_prob[1, :] = np.tile([.1, .9], int(n_features / 2))

        return np.log(class_prior), np.log(feature_prob)
795569
def setup_thresholds_panel(self, classification):
        
        # Set text in the label
        layer_purpose = self.parent.step_kw_purpose.selected_purpose()
        layer_subcategory = self.parent.step_kw_subcategory.\
            selected_subcategory()

        if is_raster_layer(self.parent.layer):
            active_band = self.parent.step_kw_band_selector.selected_band()
            layer_extent = self.parent.layer.extent()
            statistics = self.parent.layer.dataProvider().bandStatistics(
                active_band, QgsRasterBandStats.All, layer_extent, 0)
            description_text = continuous_raster_question % (
                layer_purpose['name'],
                layer_subcategory['name'],
                classification['name'],
                statistics.minimumValue,
                statistics.maximumValue)
        else:
            field_name = self.parent.step_kw_field.selected_fields()
            field_index = self.parent.layer.fields().lookupField(field_name)
            min_value_layer = self.parent.layer.minimumValue(field_index)
            max_value_layer = self.parent.layer.maximumValue(field_index)
            description_text = continuous_vector_question % (
                layer_purpose['name'],
                layer_subcategory['name'],
                field_name,
                classification['name'],
                min_value_layer,
                max_value_layer)

        # Set description
        description_label = QLabel(description_text)
        description_label.setWordWrap(True)
        self.right_layout.addWidget(description_label)

        if self.thresholds:
            thresholds = self.thresholds
        else:
            thresholds = self.parent.get_existing_keyword('thresholds')
        selected_unit = self.parent.step_kw_unit.selected_unit()['key']

        self.threshold_classes = OrderedDict()
        classes = classification.get('classes')
        # Sort by value, put the lowest first
        classes = sorted(classes, key=lambda the_key: the_key['value'])

        grid_layout_thresholds = QGridLayout()

        for i, the_class in enumerate(classes):
            class_layout = QHBoxLayout()

            # Class label
            class_label = QLabel(the_class['name'])

            # Min label
            min_label = QLabel(tr('Min >'))

            # Min value as double spin
            min_value_input = QDoubleSpinBox()
            # TODO(IS) We can set the min and max depends on the unit, later
            min_value_input.setMinimum(0)
            min_value_input.setMaximum(999999)

            if thresholds.get(self.active_exposure['key']):
                exposure_thresholds = thresholds.get(
                    self.active_exposure['key'])
                if exposure_thresholds.get(classification['key']):
                    exposure_thresholds_classifications = exposure_thresholds\
                        .get(classification['key'])
                    min_value_input.setValue(
                        exposure_thresholds_classifications['classes'][
                            the_class['key']][0])
                else:
                    default_min = the_class['numeric_default_min']
                    if isinstance(default_min, dict):
                        default_min = the_class[
                            'numeric_default_min'][selected_unit]
                    min_value_input.setValue(default_min)
            else:
                default_min = the_class['numeric_default_min']
                if isinstance(default_min, dict):
                    default_min = the_class[
                        'numeric_default_min'][selected_unit]
                min_value_input.setValue(default_min)
            min_value_input.setSingleStep(0.1)

            # Max label
            max_label = QLabel(tr('Max <='))

            # Max value as double spin
            max_value_input = QDoubleSpinBox()
            # TODO(IS) We can set the min and max depends on the unit, later
            max_value_input.setMinimum(0)
            max_value_input.setMaximum(999999)
            if thresholds.get(self.active_exposure['key']):
                exposure_thresholds = thresholds.get(
                    self.active_exposure['key'])
                if exposure_thresholds.get(classification['key']):
                    exposure_thresholds_classifications = exposure_thresholds \
                        .get(classification['key'])
                    max_value_input.setValue(
                        exposure_thresholds_classifications['classes'][
                            the_class['key']][1])
                else:
                    default_max = the_class['numeric_default_max']
                    if isinstance(default_max, dict):
                        default_max = the_class[
                            'numeric_default_max'][selected_unit]
                    max_value_input.setValue(default_max)
            else:
                default_max = the_class['numeric_default_max']
                if isinstance(default_max, dict):
                    default_max = the_class[
                        'numeric_default_max'][selected_unit]
                max_value_input.setValue(default_max)
            max_value_input.setSingleStep(0.1)

            # Add to class_layout
            class_layout.addWidget(min_label)
            class_layout.addWidget(min_value_input)
            class_layout.addWidget(max_label)
            class_layout.addWidget(max_value_input)

            class_layout.setStretch(0, 1)
            class_layout.setStretch(1, 2)
            class_layout.setStretch(2, 1)
            class_layout.setStretch(3, 2)

            # Add to grid_layout
            grid_layout_thresholds.addWidget(class_label, i, 0)
            grid_layout_thresholds.addLayout(class_layout, i, 1)

            self.threshold_classes[the_class['key']] = [
                min_value_input, max_value_input]

        grid_layout_thresholds.setColumnStretch(0, 1)
        grid_layout_thresholds.setColumnStretch(0, 2)

        def min_max_changed(double_spin_index, mode):
            """Slot when min or max value change.

            :param double_spin_index: The index of the double spin.
            :type double_spin_index: int

            :param mode: The flag to indicate the min or max value.
            :type mode: int
            """
            if mode == MAX_VALUE_MODE:
                current_max_value = list(self.threshold_classes.values())[
                    double_spin_index][1]
                target_min_value = list(self.threshold_classes.values())[
                    double_spin_index + 1][0]
                if current_max_value.value() != target_min_value.value():
                    target_min_value.setValue(current_max_value.value())
            elif mode == MIN_VALUE_MODE:
                current_min_value = list(self.threshold_classes.values())[
                    double_spin_index][0]
                target_max_value = list(self.threshold_classes.values())[
                    double_spin_index - 1][1]
                if current_min_value.value() != target_max_value.value():
                    target_max_value.setValue(current_min_value.value())

        # Set behaviour
        for k, v in list(self.threshold_classes.items()):
            index = list(self.threshold_classes.keys()).index(k)
            if index < len(self.threshold_classes) - 1:
                # Max value changed
                v[1].valueChanged.connect(partial(
                    min_max_changed,
                    double_spin_index=index,
                    mode=MAX_VALUE_MODE))
            if index > 0:
                # Min value
                v[0].valueChanged.connect(partial(
                    min_max_changed,
                    double_spin_index=index,
                    mode=MIN_VALUE_MODE))

        grid_layout_thresholds.setSpacing(0)

        self.right_layout.addLayout(grid_layout_thresholds)
225582
function decision2DStumpTrain(data, labels, ix, options) {
    
    options = options || {};
    var numtries = options.numTries || 10;
    
    // choose a dimension at random and pick a best split
    var N = ix.length;
    
    var ri1= 0;
    var ri2= 1;
    if(data[0].length > 2) {
      // more than 2D data. Pick 2 random dimensions
      ri1= randi(0, data[0].length);
      ri2= randi(0, data[0].length);
      while(ri2 == ri1) ri2= randi(0, data[0].length); // must be distinct!
    }
    
    // evaluate class entropy of incoming data
    var H= entropy(labels, ix);
    var bestGain=0; 
    var bestw1, bestw2, bestthr;
    var dots= new Array(ix.length);
    for(var i=0;i<numtries;i++) {
        
        // pick random line parameters
        var alpha= randf(0, 2*Math.PI);
        var w1= Math.cos(alpha);
        var w2= Math.sin(alpha);
        
        // project data on this line and get the dot products
        for(var j=0;j<ix.length;j++) {
          dots[j]= w1*data[ix[j]][ri1] + w2*data[ix[j]][ri2];
        }
        
        // we are in a tricky situation because data dot product distribution
        // can be skewed. So we don't want to select just randomly between
        // min and max. But we also don't want to sort as that is too expensive
        // let's pick two random points and make the threshold be somewhere between them.
        // for skewed datasets, the selected points will with relatively high likelihood
        // be in the high-desnity regions, so the thresholds will make sense
        var ix1= ix[randi(0, N)];
        var ix2= ix[randi(0, N)];
        while(ix2==ix1) ix2= ix[randi(0, N)]; // enforce distinctness of ix2
        var a= Math.random();
        var dotthr= dots[ix1]*a + dots[ix2]*(1-a);
        
        // measure information gain we'd get from split with thr
        var l1=1, r1=1, lm1=1, rm1=1; //counts for Left and label 1, right and label 1, left and minus 1, right and minus 1
        for(var j=0;j<ix.length;j++) {
            if(dots[j] < dotthr) {
              if(labels[ix[j]]==1) l1++;
              else lm1++;
            } else {
              if(labels[ix[j]]==1) r1++;
              else rm1++;
            }
        }
        var t= l1+lm1; 
        l1=l1/t;
        lm1=lm1/t;
        t= r1+rm1;
        r1=r1/t;
        rm1= rm1/t;
        
        var LH= -l1*Math.log(l1) -lm1*Math.log(lm1); // left and right entropy
        var RH= -r1*Math.log(r1) -rm1*Math.log(rm1);
        
        var informationGain= H - LH - RH;
        //console.log("Considering split %f, entropy %f -> %f, %f. Gain %f", thr, H, LH, RH, informationGain);
        if(informationGain > bestGain || i === 0) {
            bestGain= informationGain;
            bestw1= w1;
            bestw2= w2;
            bestthr= dotthr;
        }
    }
    
    model= {};
    model.w1= bestw1;
    model.w2= bestw2;
    model.dotthr= bestthr;
    return model;
}
342452
private void qcrit( RenderedImage slope, RenderedImage ab, RandomIter trasmissivityRI, RandomIter frictionRI,
            RandomIter cohesionRI, RandomIter hsIter, RandomIter effectiveRI, RandomIter densityRI ) {
        HashMap<String, Double> regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inSlope);
        int cols = regionMap.get(CoverageUtilities.COLS).intValue();
        int rows = regionMap.get(CoverageUtilities.ROWS).intValue();

        RandomIter slopeRI = RandomIterFactory.create(slope, null);
        RandomIter abRI = RandomIterFactory.create(ab, null);

        WritableRaster qcritWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, null);
        WritableRandomIter qcritIter = RandomIterFactory.createWritable(qcritWR, null);
        WritableRaster classiWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, null);
        WritableRandomIter classiIter = RandomIterFactory.createWritable(classiWR, null);

        pm.beginTask("Creating qcrit map...", rows);
        for( int j = 0; j < rows; j++ ) {
            pm.worked(1);
            for( int i = 0; i < cols; i++ ) {
                double slopeValue = slopeRI.getSampleDouble(i, j, 0);
                double tanPhiValue = frictionRI.getSampleDouble(i, j, 0);
                double cohValue = cohesionRI.getSampleDouble(i, j, 0);
                double rhoValue = densityRI.getSampleDouble(i, j, 0);
                double hsValue = hsIter.getSampleDouble(i, j, 0);

                if (!isNovalue(slopeValue) && !isNovalue(tanPhiValue) && !isNovalue(cohValue) && !isNovalue(rhoValue)
                        && !isNovalue(hsValue)) {
                    if (hsValue <= EPS || slopeValue > pRock) {
                        qcritIter.setSample(i, j, 0, ROCK);
                    } else {
                        double checkUnstable = tanPhiValue + cohValue / (9810.0 * rhoValue * hsValue) * (1 + pow(slopeValue, 2));
                        if (slopeValue >= checkUnstable) {
                            /*
                             * uncond unstable
                             */
                            qcritIter.setSample(i, j, 0, 5);
                        } else {
                            double checkStable = tanPhiValue * (1 - 1 / rhoValue) + cohValue / (9810 * rhoValue * hsValue)
                                    * (1 + pow(slopeValue, 2));
                            if (slopeValue < checkStable) {
                                /*
                                 * uncond. stable
                                 */
                                qcritIter.setSample(i, j, 0, 0);
                            } else {
                                double qCrit = trasmissivityRI.getSampleDouble(i, j, 0)
                                        * sin(atan(slopeValue))
                                        / abRI.getSampleDouble(i, j, 0)
                                        * rhoValue
                                        * (1 - slopeValue / tanPhiValue + cohValue / (9810 * rhoValue * hsValue * tanPhiValue)
                                                * (1 + pow(slopeValue, 2))) * 1000;
                                qcritIter.setSample(i, j, 0, qCrit);
                                /*
                                 * see the Qcrit (critical effective
                                 * precipitation) that leads the slope to
                                 * instability (see article of Montgomery et Al,
                                 * Hydrological Processes, 12, 943-955, 1998)
                                 */
                                double value = qcritIter.getSampleDouble(i, j, 0);
                                if (value > 0 && value < 50)
                                    qcritIter.setSample(i, j, 0, 1);
                                if (value >= 50 && value < 100)
                                    qcritIter.setSample(i, j, 0, 2);
                                if (value >= 100 && value < 200)
                                    qcritIter.setSample(i, j, 0, 3);
                                if (value >= 200)
                                    qcritIter.setSample(i, j, 0, 4);
                            }
                        }
                    }
                } else {
                    qcritIter.setSample(i, j, 0, doubleNovalue);
                }
            }
        }
        pm.done();

        /*
         * build the class matrix 1=inc inst 2=inc stab 3=stab 4=instab
         * rock=presence of rock
         */
        pm.beginTask("Creating stability map...", rows);
        double Tq = 0;
        for( int j = 0; j < rows; j++ ) {
            pm.worked(1);
            for( int i = 0; i < cols; i++ ) {
                Tq = trasmissivityRI.getSampleDouble(i, j, 0) / (effectiveRI.getSampleDouble(i, j, 0) / 1000.0);
                double slopeValue = slopeRI.getSampleDouble(i, j, 0);
                double abValue = abRI.getSampleDouble(i, j, 0);
                double tangPhiValue = frictionRI.getSampleDouble(i, j, 0);
                double cohValue = cohesionRI.getSampleDouble(i, j, 0);
                double rhoValue = densityRI.getSampleDouble(i, j, 0);
                double hsValue = hsIter.getSampleDouble(i, j, 0);

                if (!isNovalue(slopeValue) && !isNovalue(abValue) && !isNovalue(tangPhiValue) && !isNovalue(cohValue)
                        && !isNovalue(rhoValue) && !isNovalue(hsValue)) {
                    if (hsValue <= EPS || slopeValue > pRock) {
                        classiIter.setSample(i, j, 0, ROCK);
                    } else {
                        double checkUncondUnstable = tangPhiValue + cohValue / (9810 * rhoValue * hsValue)
                                * (1 + pow(slopeValue, 2));
                        double checkUncondStable = tangPhiValue * (1 - 1 / rhoValue) + cohValue / (9810 * rhoValue * hsValue)
                                * (1 + pow(slopeValue, 2));
                        double checkStable = Tq
                                * sin(atan(slopeValue))
                                * rhoValue
                                * (1 - slopeValue / tangPhiValue + cohValue / (9810 * rhoValue * hsValue * tangPhiValue)
                                        * (1 + pow(slopeValue, 2)));
                        if (slopeValue >= checkUncondUnstable) {
                            classiIter.setSample(i, j, 0, 1);
                        } else if (slopeValue < checkUncondStable) {
                            classiIter.setSample(i, j, 0, 2);
                        } else if (abValue < checkStable && classiIter.getSampleDouble(i, j, 0) != 1
                                && classiIter.getSampleDouble(i, j, 0) != 2) {
                            classiIter.setSample(i, j, 0, 3);
                        } else {
                            classiIter.setSample(i, j, 0, 4);
                        }
                    }
                } else {
                    classiIter.setSample(i, j, 0, doubleNovalue);
                }
            }
        }
        pm.done();

        outQcrit = CoverageUtilities.buildCoverage("qcrit", qcritWR, regionMap, inSlope.getCoordinateReferenceSystem());
        outShalstab = CoverageUtilities.buildCoverage("classi", classiWR, regionMap, inSlope.getCoordinateReferenceSystem());

    }
1285837
def estimate_parameters_rippe(
        self, max_dist_kb, size_bin_kb, display_graph
    ):
        
        logger.info("estimation of the parameters of the model")
        self.bins = np.arange(
            size_bin_kb, max_dist_kb + size_bin_kb, size_bin_kb
        )

        self.mean_contacts = np.zeros_like(self.bins, dtype=np.float32)
        self.dict_collect = dict()
        self.gpu_vect_frags.copy_from_gpu()
        epsi = self.mean_value_trans
        for k in self.bins:
            self.dict_collect[k] = []
        for i in range(0, self.n_frags // 10):
            # print "frag i = ", i
            start = self.sparse_matrix.indptr[i]
            end = self.sparse_matrix.indptr[i + 1]
            id_j = self.sparse_matrix.indices[start:end]
            data = self.sparse_matrix.data[start:end]

            info_i = self.np_sub_frags_2_frags[i]
            init_id_fi = int(info_i[0])
            # pos_i = self.S_o_A_frags["pos"][init_id_fi]
            id_c_i = self.S_o_A_frags["id_c"][init_id_fi]
            s_i = (
                self.S_o_A_frags["start_bp"][init_id_fi] / 1000.0
                + self.np_sub_frags_2_frags[i][1]
            )
            len_kb_c_i = self.S_o_A_frags["l_cont_bp"][init_id_fi] / 1000
            # local_bins = np.arange(size_bin_kb, min(len_kb_c_i, max_dist_kb)
            # + size_bin_kb, size_bin_kb) local_storage =
            # np.zeros_like(local_bins, dtype=np.int32)
            #
            # for fj, dj in zip(id_j, data):
            #     info_j = self.np_sub_frags_2_frags[fj]
            #     init_id_fj = info_j[0]
            #     id_c_j = self.S_o_A_frags['id_c'][init_id_fj]
            #     if id_c_i == id_c_j:
            #         pos_j = self.S_o_A_frags['pos'][init_id_fj]
            #         s_j = self.S_o_A_frags['start_bp'][init_id_fj]/1000.0 +
            #         self.np_sub_frags_2_frags[fj][1]
            #         d = np.abs(s_i - s_j)
            #         if d < max_dist_kb:
            #             id_bin = d / size_bin_kb
            #             local_storage[id_bin] += dj
            #             # self.dict_collect[self.bins[id_bin]].append(dj)
            # # we have to add also the zeros
            # for bin, val in zip(local_bins, local_storage):
            #     # print "bin = ", bin
            #     self.dict_collect[bin].append(val)

            if size_bin_kb < len_kb_c_i:
                # local_bins = np.arange(size_bin_kb, min(len_kb_c_i,
                # max_dist_kb) + size_bin_kb, size_bin_kb)
                local_bins = np.arange(
                    size_bin_kb, max_dist_kb + size_bin_kb, size_bin_kb
                )
                local_storage = np.zeros_like(local_bins, dtype=np.int32)

                for fj, dj in zip(id_j, data):
                    info_j = self.np_sub_frags_2_frags[fj]
                    init_id_fj = int(info_j[0])
                    id_c_j = self.S_o_A_frags["id_c"][init_id_fj]
                    if id_c_i == id_c_j:
                        # pos_j = self.S_o_A_frags["pos"][init_id_fj]
                        s_j = (
                            self.S_o_A_frags["start_bp"][init_id_fj] / 1000.0
                            + self.np_sub_frags_2_frags[fj][1]
                        )
                        d = np.abs(s_i - s_j)
                        if d < max_dist_kb:
                            id_bin = int(d / size_bin_kb)
                            local_storage[id_bin] += dj
                            # self.dict_collect[self.bins[id_bin]].append(dj)
                # we have to add also the zeros
                for bin, val in zip(local_bins, local_storage):
                    # print "bin = ", bin
                    self.dict_collect[bin].append(val)

        for id_bin in range(0, len(self.bins)):
            k = self.bins[id_bin]
            self.mean_contacts[id_bin] = np.mean(self.dict_collect[k])

        for id_bin in range(0, len(self.bins)):
            k = self.bins[id_bin]
            tmp = np.mean(self.dict_collect[k])
            if np.isnan(tmp) or tmp == 0:
                # if np.isnan(tmp):
                # if np.isnan(tmp):
                #     print "removing nan"
                self.mean_contacts[id_bin] = np.nan
            else:
                self.mean_contacts[id_bin] = tmp + epsi

        self.mean_contacts_upd = []
        self.bins_upd = []

        for count, ele in zip(self.mean_contacts, self.bins):
            if not np.isnan(count):
                self.bins_upd.append(ele)
                self.mean_contacts_upd.append(count)

        self.bins_upd = np.array(self.bins_upd)
        # self.mean_contacts_upd.sort()
        # self.mean_contacts_upd.reverse()
        self.mean_contacts_upd = np.array(self.mean_contacts_upd)

        # self.mean_contacts_upd =
        # ndi.filters.gaussian_filter1d(self.mean_contacts_upd,
        # sigma=len(self.bins_upd) / 5.)

        # print "size mean contacts vector = ", self.mean_contacts_upd.shape
        # print "mean contacts vector = ", self.mean_contacts_upd

        p, self.y_estim = opti.estimate_param_rippe(
            self.mean_contacts_upd, self.bins_upd
        )
        ##########################################
        logger.info("p from estimate parameters  = {}".format(p))
        # p = list(p[0])
        # p[3] = 2
        # p = tuple(p)
        ##########################################
        fit_param = p
        logger.info("mean value trans = {}".format(self.mean_value_trans))
        ##########################################
        logger.info("BEWARE!!! : I will lower mean value trans  !!!")
        self.mean_value_trans = self.mean_value_trans / 10.0
        ##########################################

        estim_max_dist = opti.estimate_max_dist_intra(
            fit_param, self.mean_value_trans
        )
        logger.info("estimate max dist cis trans = {}".format(estim_max_dist))
        self.param_simu = self.setup_rippe_parameters(
            fit_param, estim_max_dist
        )
        self.param_simu_test = self.param_simu
        self.gpu_param_simu = cuda.mem_alloc(self.param_simu.nbytes)
        self.gpu_param_simu_test = cuda.mem_alloc(self.param_simu.nbytes)

        cuda.memcpy_htod(self.gpu_param_simu, self.param_simu)
        cuda.memcpy_htod(self.gpu_param_simu_test, self.param_simu_test)
        display_graph = False
        if display_graph:
            plt.figure()
            plt.loglog(self.bins_upd, self.mean_contacts_upd, "-*b")
            plt.loglog(self.bins_upd, self.y_estim, "-*r")
            plt.xlabel("genomic distance (kb)")
            plt.ylabel("frequency of contact")
            plt.title(
                r"$\mathrm{Frequency\ of\ contact\ versus\ genomic\ distance"
                " \(data\ tricho test):}\ slope=%.3f,\ max\ cis\ distance(kb)"
                "=%.3f\ d=%.3f\ scale\ factor=%.3f\ $"
                % (
                    self.param_simu["slope"],
                    estim_max_dist,
                    self.param_simu["d"],
                    self.param_simu["fact"],
                )
            )
            plt.legend(["obs", "fit"])
            plt.savefig()
            plt.close()

        self.eval_likelihood_init()
1281248
public void convolveDerivOrder() {
		T blur = GeneralizedImageOps.createSingleBand(imageType, width, height);
		T blurDeriv = GeneralizedImageOps.createSingleBand(imageType, width, height);
		T deriv = GeneralizedImageOps.createSingleBand(imageType, width, height);
		T derivBlur = GeneralizedImageOps.createSingleBand(imageType, width, height);

		BlurStorageFilter<T> funcBlur = FactoryBlurFilter.gaussian(ImageType.single(imageType),sigma,radius);
		ImageGradient<T,T> funcDeriv = FactoryDerivative.three(imageType,imageType);

		funcBlur.process(input,blur);
		funcDeriv.process(blur,blurDeriv,derivY);

		funcDeriv.process(input,deriv,derivY);
		funcBlur.process(deriv,derivBlur);

		printIntensity("Blur->Deriv",blurDeriv);
		printIntensity("Deriv->Blur",blurDeriv);
	}
970884
def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids):
        """"""

        npt = bd.shape[0]
        n = self.X_ADJUSTED.shape[0]
        zero_index = None
        zero_value = False

        a_inv = scipy.linalg.inv(a)

        if np.any(np.absolute(bd) <= self.eps):
            zero_value = True
            zero_index = np.where(np.absolute(bd) <= self.eps)

        if self.UNBIAS:
            b = np.zeros((npt, n_withdrifts+1, 1))
        else:
            b = np.zeros((npt, n_withdrifts, 1))
        b[:, :n, 0] = - self.variogram_function(self.variogram_model_parameters, bd)
        if zero_value:
            b[zero_index[0], zero_index[1], 0] = 0.0

        i = n
        if self.regional_linear_drift:
            b[:, i, 0] = xyz[:, 2]
            i += 1
            b[:, i, 0] = xyz[:, 1]
            i += 1
            b[:, i, 0] = xyz[:, 0]
            i += 1
        if self.specified_drift:
            for spec_vals in spec_drift_grids:
                b[:, i, 0] = spec_vals.flatten()
                i += 1
        if self.functional_drift:
            for func in self.functional_drift_terms:
                b[:, i, 0] = func(xyz[:, 2], xyz[:, 1], xyz[:, 0])
                i += 1
        if i != n_withdrifts:
            warnings.warn("Error in setting up kriging system. "
                          "Kriging may fail.", RuntimeWarning)
        if self.UNBIAS:
            b[:, n_withdrifts, 0] = 1.0

        if (~mask).any():
            mask_b = np.repeat(mask[:, np.newaxis, np.newaxis],
                               n_withdrifts+1, axis=1)
            b = np.ma.array(b, mask=mask_b)

        if self.UNBIAS:
            x = np.dot(a_inv, b.reshape((npt, n_withdrifts+1)).T).reshape((1, n_withdrifts+1, npt)).T
        else:
            x = np.dot(a_inv, b.reshape((npt, n_withdrifts)).T).reshape((1, n_withdrifts, npt)).T
        kvalues = np.sum(x[:, :n, 0] * self.VALUES, axis=1)
        sigmasq = np.sum(x[:, :, 0] * -b[:, :, 0], axis=1)

        return kvalues, sigmasq
61829
public static int cudnnBatchNormalizationBackward(
        cudnnHandle handle, 
        int mode, 
        Pointer alphaDataDiff, 
        Pointer betaDataDiff, 
        Pointer alphaParamDiff, 
        Pointer betaParamDiff, 
        cudnnTensorDescriptor xDesc, /** same desc for x, dx, dy */
        Pointer x, 
        cudnnTensorDescriptor dyDesc, 
        Pointer dy, 
        cudnnTensorDescriptor dxDesc, 
        Pointer dx, 
        /** Shared tensor desc for the 4 tensors below */
        cudnnTensorDescriptor dBnScaleBiasDesc, 
        Pointer bnScale, /** bnBias doesn't affect backpropagation */
        /** scale and bias diff are not backpropagated below this layer */
        Pointer dBnScaleResult, 
        Pointer dBnBiasResult, 
        /** Same epsilon as forward pass */
        double epsilon, 
        /** Optionally cached intermediate results from
                                           forward pass */
        Pointer savedMean, 
        Pointer savedInvVariance)
    {
        return checkResult(cudnnBatchNormalizationBackwardNative(handle, mode, alphaDataDiff, betaDataDiff, alphaParamDiff, betaParamDiff, xDesc, x, dyDesc, dy, dxDesc, dx, dBnScaleBiasDesc, bnScale, dBnScaleResult, dBnBiasResult, epsilon, savedMean, savedInvVariance));
    }
619129
def inference(self, kern, X, likelihood, Y, Y_metadata=None):

        
        N = X.shape[0] #number of training points
        D = X.shape[1] #number of dimensions

        Kds = np.zeros(D, dtype=object) #vector for holding covariance per dimension
        Qs = np.zeros(D, dtype=object) #vector for holding eigenvectors of covariance per dimension
        QTs = np.zeros(D, dtype=object) #vector for holding transposed eigenvectors of covariance per dimension
        V_kron = 1 # kronecker product of eigenvalues

        # retrieve the one-dimensional variation of the designated kernel
        oneDkernel = kern.get_one_dimensional_kernel(D)

        for d in range(D):
            xg = list(set(X[:,d])) #extract unique values for a dimension
            xg = np.reshape(xg, (len(xg), 1))
            oneDkernel.lengthscale = kern.lengthscale[d]
            Kds[d] = oneDkernel.K(xg)
            [V, Q] = np.linalg.eig(Kds[d])
            V_kron = np.kron(V_kron, V)
            Qs[d] = Q
            QTs[d] = Q.T

        noise = likelihood.variance + 1e-8

        alpha_kron = self.kron_mvprod(QTs, Y)
        V_kron = V_kron.reshape(-1, 1)
        alpha_kron = alpha_kron / (V_kron + noise)
        alpha_kron = self.kron_mvprod(Qs, alpha_kron)

        log_likelihood = -0.5 * (np.dot(Y.T, alpha_kron) + np.sum((np.log(V_kron + noise))) + N*log_2_pi)

        # compute derivatives wrt parameters Thete
        derivs = np.zeros(D+2, dtype='object')
        for t in range(len(derivs)):
            dKd_dTheta = np.zeros(D, dtype='object')
            gamma = np.zeros(D, dtype='object')
            gam = 1
            for d in range(D):
                xg = list(set(X[:,d]))
                xg = np.reshape(xg, (len(xg), 1))
                oneDkernel.lengthscale = kern.lengthscale[d]
                if t < D:
                    dKd_dTheta[d] = oneDkernel.dKd_dLen(xg, (t==d), lengthscale=kern.lengthscale[t]) #derivative wrt lengthscale
                elif (t == D):
                    dKd_dTheta[d] = oneDkernel.dKd_dVar(xg) #derivative wrt variance
                else:
                    dKd_dTheta[d] = np.identity(len(xg)) #derivative wrt noise
                gamma[d] = np.diag(np.dot(np.dot(QTs[d], dKd_dTheta[d].T), Qs[d]))
                gam = np.kron(gam, gamma[d])
            
            gam = gam.reshape(-1,1)
            kappa = self.kron_mvprod(dKd_dTheta, alpha_kron)
            derivs[t] = 0.5*np.dot(alpha_kron.T,kappa) - 0.5*np.sum(gam / (V_kron + noise))

        # separate derivatives
        dL_dLen = derivs[:D]
        dL_dVar = derivs[D]
        dL_dThetaL = derivs[D+1]

        return GridPosterior(alpha_kron=alpha_kron, QTs=QTs, Qs=Qs, V_kron=V_kron), \
                log_likelihood, {'dL_dLen':dL_dLen, 'dL_dVar':dL_dVar, 'dL_dthetaL':dL_dThetaL}
1066380
def solution_path(self, min_lambda, max_lambda, lambda_bins, verbose=0):
        ''''''
        lambda_grid = np.exp(np.linspace(np.log(max_lambda), np.log(min_lambda), lambda_bins))
        aic_trace = np.zeros(lambda_grid.shape) # The AIC score for each lambda value
        aicc_trace = np.zeros(lambda_grid.shape) # The AICc score for each lambda value (correcting for finite sample size)
        bic_trace = np.zeros(lambda_grid.shape) # The BIC score for each lambda value
        dof_trace = np.zeros(lambda_grid.shape) # The degrees of freedom of each final solution
        log_likelihood_trace = np.zeros(lambda_grid.shape)
        beta_trace = []
        best_idx = None
        best_plateaus = None

        # Solve the series of lambda values with warm starts at each point
        for i, lam in enumerate(lambda_grid):
            if verbose:
                print('#{0} Lambda = {1}'.format(i, lam))

            # Fit to the final values
            beta = self.solve(lam)

            if verbose:
                print('Calculating degrees of freedom')

            # Count the number of free parameters in the grid (dof)
            plateaus = calc_plateaus(beta, self.edges)
            dof_trace[i] = len(plateaus)

            if verbose:
                print('Calculating AIC')

            # Get the negative log-likelihood
            log_likelihood_trace[i] = self.log_likelihood(beta)

            # Calculate AIC = 2k - 2ln(L)
            aic_trace[i] = 2. * dof_trace[i] - 2. * log_likelihood_trace[i]
            
            # Calculate AICc = AIC + 2k * (k+1) / (n - k - 1)
            aicc_trace[i] = aic_trace[i] + 2 * dof_trace[i] * (dof_trace[i]+1) / (len(beta) - dof_trace[i] - 1.)

            # Calculate BIC = -2ln(L) + k * (ln(n) - ln(2pi))
            bic_trace[i] = -2 * log_likelihood_trace[i] + dof_trace[i] * (np.log(len(beta)) - np.log(2 * np.pi))

            # Track the best model thus far
            if best_idx is None or bic_trace[i] < bic_trace[best_idx]:
                best_idx = i
                best_plateaus = plateaus

            # Save the trace of all the resulting parameters
            beta_trace.append(np.array(beta))
            
            if verbose:
                print('DoF: {0} AIC: {1} AICc: {2} BIC: {3}'.format(dof_trace[i], aic_trace[i], aicc_trace[i], bic_trace[i]))

        if verbose:
            print('Best setting (by BIC): lambda={0} [DoF: {1}, AIC: {2}, AICc: {3} BIC: {4}]'.format(lambda_grid[best_idx], dof_trace[best_idx], aic_trace[best_idx], aicc_trace[best_idx], bic_trace[best_idx]))

        return {'aic': aic_trace,
                'aicc': aicc_trace,
                'bic': bic_trace,
                'dof': dof_trace,
                'loglikelihood': log_likelihood_trace,
                'beta': np.array(beta_trace),
                'lambda': lambda_grid,
                'best_idx': best_idx,
                'best': beta_trace[best_idx],
                'plateaus': best_plateaus}