Train the posteriori multiple probe algorithm.
@param range the neighborhood search data structure.
@param radius the radius for range search.
@param Nz the number of quantized values.
@param sigma the Parzen window width.
161818
public void learn(RNNSearch<double[], double[]> range, double[][] samples, double radius, int Nz, double sigma) {
        TrainSample[] training = new TrainSample[samples.length];
        for (int i = 0; i < samples.length; i++) {
            training[i] = new TrainSample();
            training[i].query = samples[i];
            training[i].neighbors = new ArrayList<>();
            ArrayList<Neighbor<double[], double[]>> neighbors = new ArrayList<>();
            range.range(training[i].query, radius, neighbors);
            for (Neighbor<double[], double[]> n : neighbors) {
                training[i].neighbors.add(keys.get(n.index));
            }
        }

        model = new ArrayList<>(hash.size());
        for (int i = 0; i < hash.size(); i++) {
            model.add(new PosterioriModel(hash.get(i), training, Nz, sigma));
        }
    }
166578
public static void setPreferredRadius(double radius) {
		assert !Double.isNaN(radius);
		final Preferences prefs = Preferences.userNodeForPackage(MapElementConstants.class);
		if (prefs != null) {
			prefs.putDouble("RADIUS", radius); //$NON-NLS-1$
			try {
				prefs.flush();
			} catch (BackingStoreException exception) {
				//
			}
		}
	}
1074451
def sigma_s2(self, r, R, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        beta = self.anisotropy.beta_r(r, kwargs_anisotropy)
        return (1 - beta * R**2/r**2) * self.sigma_r2(r, kwargs_profile, kwargs_anisotropy, kwargs_light)
144133
def sigma_s2(self, r, R, r_ani, a, gamma, rho0_r0_gamma):
        
        beta = self._beta_ani(r, r_ani)
        return (1 - beta * R**2/r**2) * self.sigma_r2(r, a, gamma, rho0_r0_gamma, r_ani)
611939
@Override
    public Map<Integer, Integer> pGetRanges(int partition) {
        Map<Integer, Integer> ranges = new TreeMap<Integer, Integer>();
        Integer first = null; // start of the very first token on the ring
        Integer start = null; // start of a range
        UnmodifiableIterator<Map.Entry<Integer,Integer>> iter = m_tokensMap.get().entrySet().iterator();

        // Iterate through the token map to find the ranges assigned to
        // the given partition
        while (iter.hasNext()) {
            Map.Entry<Integer, Integer> next = iter.next();
            int token = next.getKey();
            int pid = next.getValue();

            if (first == null) {
                first = token;
            }

            // if start is not null, there's an open range, now is
            // the time to close it.
            // else there is no open range, keep on going.
            if (start != null) {
                //Range end is inclusive so do token - 1
                ranges.put(start, token - 1);
                start = null;
            }

            if (pid == partition) {
                // if start is null, there's no open range, start one.
                start = token;
            }
        }

        // if there is an open range when we get here
        // It is the last token which implicity ends at the next max value
        if (start != null) {
            assert first != null;
            ranges.put(start, Integer.MAX_VALUE);
        }

        return ranges;
    }
1291950
public void orizzonte7( double delta, int quadrata, double beta, double alfa, RandomIter elevImageIterator,
            WritableRaster curvatureImage, int[][] shadow ) {
        int rows = curvatureImage.getHeight();
        int cols = curvatureImage.getWidth();
        /*=====================*/

        int y, I, J;
        double zenith;

        /*======================*/

        for( int i = quadrata - 1; i >= 0; i-- ) {
            I = -1;
            J = -1;
            y = 0;
            for( int ii = i; ii < quadrata - 1; ii++ ) {
                for( int j = (int) floor(1 / tan(beta) * (ii - i)); j <= (int) floor(1 / tan(beta) * (ii - i + 1)) - 1
                        && j < cols; j++ ) {
                    if (ii >= (rows + 2 * cols) && !isNovalue(elevImageIterator.getSampleDouble(j, ii - (rows + 2 * cols), 0))) {
                        /*shadow.element[ii-(Z0.nrh+2*Z0.nch)][j]=i;}}}}*/
                        if (curvatureImage.getSampleDouble(j, ii - (rows + 2 * cols), 0) == 1 && I == -1) {
                            I = ii - (rows + 2 * cols);
                            J = j;
                            y = 1;
                        } else if (curvatureImage.getSampleDouble(j, ii - (rows + 2 * cols), 0) == 1 && I != -1) {
                            zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(j, ii
                                    - (rows + 2 * cols), 0))
                                    / sqrt(pow((double) (I - (ii - (rows + 2 * cols))) * (double) delta, (double) 2)
                                            + pow((double) (J - j) * (double) delta, (double) 2));
                            if (zenith <= tan(alfa)) {
                                shadow[ii - (rows + 2 * cols)][j] = 0;
                                I = ii - (rows + 2 * cols);
                                J = j;
                            } else {
                                shadow[ii - (rows + 2 * cols)][j] = 1;
                            }
                        } else if (curvatureImage.getSampleDouble(j, ii - (rows + 2 * cols), 0) == 0 && y == 1) {
                            zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(j, ii
                                    - (rows + 2 * cols), 0))
                                    / sqrt(pow((double) (I - (ii - (rows + 2 * cols))) * (double) delta, (double) 2)
                                            + pow((double) (J - j) * (double) delta, (double) 2));
                            if (zenith <= tan(alfa)) {
                                shadow[ii - (rows + 2 * cols)][j] = 0;
                                y = 0;
                            } else {
                                shadow[ii - (rows + 2 * cols)][j] = 1;
                            }
                        }
                    }
                }
            }
        }

    }
675129
private void computeDerivatives(double sigma) {
		this.Lx2 = new float[width][height];
		this.Ly2 = new float[width][height];
		this.Lxy = new float[width][height];

		float[][][] grad = new float[width][height][];
		for (int y = 0; y < height; y++)
			for (int x = 0; x < width; x++)
				grad[x][y] = sobel(x, y);

		int radius = (int) (2 * sigma);
		int window = 1 + 2 * radius;
		float[][] gaussian = new float[window][window];
		for (int j = -radius; j <= radius; j++)
			for (int i = -radius; i <= radius; i++)
				gaussian[i + radius][j + radius] = (float) gaussian(i, j, sigma);

		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {

				for (int dy = -radius; dy <= radius; dy++) {
					for (int dx = -radius; dx <= radius; dx++) {
						int xk = x + dx;
						int yk = y + dy;
						if (xk < 0 || xk >= width)
							continue;
						if (yk < 0 || yk >= height)
							continue;

						double gw = gaussian[dx + radius][dy + radius];

						this.Lx2[x][y] += gw * grad[xk][yk][0]
								* grad[xk][yk][0];
						this.Ly2[x][y] += gw * grad[xk][yk][1]
								* grad[xk][yk][1];
						this.Lxy[x][y] += gw * grad[xk][yk][0]
								* grad[xk][yk][1];
					}
				}
			}
		}
	}
1001134
private void depthFirstSearch(final Range<C> query, final Node node, final List<Range<C>> result, final Set<Node> visited) {
        if (node == null || visited.contains(node) || query.isEmpty()) {
            return;
        }
        if (node.left() != null && Ranges.isLessThan(query, node.center())) {
            depthFirstSearch(query, node.left(), result, visited);
        }
        else if (node.right() != null && Ranges.isGreaterThan(query, node.center())) {
            depthFirstSearch(query, node.right(), result, visited);
        }
        if (Ranges.isGreaterThan(query, node.center())) {
            for (Range<C> range : node.overlapByUpperEndpoint()) {
                if (Ranges.intersect(range, query)) {
                    result.add(range);
                }
                if (Ranges.isGreaterThan(query, range.upperEndpoint())) {
                    break;
                }
            }
        }
        else if (Ranges.isLessThan(query, node.center())) {
            for (Range<C> range : node.overlapByLowerEndpoint()) {
                if (Ranges.intersect(range, query)) {
                    result.add(range);
                }
                if (Ranges.isLessThan(query, range.lowerEndpoint())) {
                    break;
                }
            }
        }
        else {
            result.addAll(node.overlapByLowerEndpoint());
        }
        visited.add(node);
    }
1173132
def search(self, coords=None, distance_modulus=None, tolerance=1.e-2):
        
        nmoduli = len(self.distance_modulus_array)
        npixels    = len(self.roi.pixels_target)
        self.log_likelihood_sparse_array       = np.zeros([nmoduli, npixels])
        self.richness_sparse_array             = np.zeros([nmoduli, npixels])
        self.richness_lower_sparse_array       = np.zeros([nmoduli, npixels])
        self.richness_upper_sparse_array       = np.zeros([nmoduli, npixels])
        self.richness_upper_limit_sparse_array = np.zeros([nmoduli, npixels])
        self.stellar_mass_sparse_array         = np.zeros([nmoduli, npixels])
        self.fraction_observable_sparse_array  = np.zeros([nmoduli, npixels])

        # Specific pixel/distance_modulus
        coord_idx, distance_modulus_idx = None, None
        if coords is not None:
            # Match to nearest grid coordinate index
            coord_idx = self.roi.indexTarget(coords[0],coords[1])
        if distance_modulus is not None:
            # Match to nearest distance modulus index
            distance_modulus_idx=np.fabs(self.distance_modulus_array-distance_modulus).argmin()

        lon, lat = self.roi.pixels_target.lon, self.roi.pixels_target.lat
            
        logger.info('Looping over distance moduli in grid search ...')
        for ii, distance_modulus in enumerate(self.distance_modulus_array):

            # Specific pixel
            if distance_modulus_idx is not None:
                if ii != distance_modulus_idx: continue

            logger.info('  (%-2i/%i) Distance Modulus=%.1f ...'%(ii+1,nmoduli,distance_modulus))

            # Set distance_modulus once to save time
            self.loglike.set_params(distance_modulus=distance_modulus)

            for jj in range(0, npixels):
                # Specific pixel
                if coord_idx is not None:
                    if jj != coord_idx: continue

                # Set kernel location
                self.loglike.set_params(lon=lon[jj],lat=lat[jj])
                # Doesn't re-sync distance_modulus each time
                self.loglike.sync_params()
                                         
                args = (jj+1, npixels, self.loglike.source.lon, self.loglike.source.lat)
                msg = '    (%-3i/%i) Candidate at (%.2f, %.2f) ... '%(args)

                self.log_likelihood_sparse_array[ii][jj], self.richness_sparse_array[ii][jj], parabola = self.loglike.fit_richness()
                self.stellar_mass_sparse_array[ii][jj] = self.stellar_mass_conversion * self.richness_sparse_array[ii][jj]
                self.fraction_observable_sparse_array[ii][jj] = self.loglike.f
                if self.config['scan']['full_pdf']:
                    #n_pdf_points = 100
                    #richness_range = parabola.profileUpperLimit(delta=25.) - self.richness_sparse_array[ii][jj]
                    #richness = np.linspace(max(0., self.richness_sparse_array[ii][jj] - richness_range),
                    #                          self.richness_sparse_array[ii][jj] + richness_range,
                    #                          n_pdf_points)
                    #if richness[0] > 0.:
                    #    richness = np.insert(richness, 0, 0.)
                    #    n_pdf_points += 1
                    # 
                    #log_likelihood = np.zeros(n_pdf_points)
                    #for kk in range(0, n_pdf_points):
                    #    log_likelihood[kk] = self.loglike.value(richness=richness[kk])
                    #parabola = ugali.utils.parabola.Parabola(richness, 2.*log_likelihood)
                    #self.richness_lower_sparse_array[ii][jj], self.richness_upper_sparse_array[ii][jj] = parabola.confidenceInterval(0.6827)
                    self.richness_lower_sparse_array[ii][jj], self.richness_upper_sparse_array[ii][jj] = self.loglike.richness_interval(0.6827)
                    
                    self.richness_upper_limit_sparse_array[ii][jj] = parabola.bayesianUpperLimit(0.95)

                    args = (
                        2. * self.log_likelihood_sparse_array[ii][jj],
                        self.stellar_mass_conversion*self.richness_sparse_array[ii][jj],
                        self.stellar_mass_conversion*self.richness_lower_sparse_array[ii][jj],
                        self.stellar_mass_conversion*self.richness_upper_sparse_array[ii][jj],
                        self.stellar_mass_conversion*self.richness_upper_limit_sparse_array[ii][jj]
                    )
                    msg += 'TS=%.1f, Stellar Mass=%.1f (%.1f -- %.1f @ 0.68 CL, < %.1f @ 0.95 CL)'%(args)
                else:
                    args = (
                        2. * self.log_likelihood_sparse_array[ii][jj], 
                        self.stellar_mass_conversion * self.richness_sparse_array[ii][jj],
                        self.fraction_observable_sparse_array[ii][jj]
                    )
                    msg += 'TS=%.1f, Stellar Mass=%.1f, Fraction=%.2g'%(args)
                logger.debug(msg)
                
                #if coords is not None and distance_modulus is not None:
                #    results = [self.richness_sparse_array[ii][jj],
                #               self.log_likelihood_sparse_array[ii][jj],
                #               self.richness_lower_sparse_array[ii][jj],
                #               self.richness_upper_sparse_array[ii][jj],
                #               self.richness_upper_limit_sparse_array[ii][jj],
                #               richness, log_likelihood, self.loglike.p, self.loglike.f]
                #    return results

            jj_max = self.log_likelihood_sparse_array[ii].argmax()
            args = (
                jj_max+1, npixels, lon[jj_max], lat[jj_max],
                2. * self.log_likelihood_sparse_array[ii][jj_max], 
                self.stellar_mass_conversion * self.richness_sparse_array[ii][jj_max]
            )
            msg = '  (%-3i/%i) Maximum at (%.2f, %.2f) ... TS=%.1f, Stellar Mass=%.1f'%(args)
            logger.info(msg)
505414
public void setRangeList(int i, Annotation v) { 
    if (OntRelationMention_Type.featOkTst && ((OntRelationMention_Type)jcasType).casFeat_rangeList == null)
      jcasType.jcas.throwFeatMissing("rangeList", "de.julielab.jules.types.OntRelationMention");
    jcasType.jcas.checkArrayBounds(jcasType.ll_cas.ll_getRefValue(addr, ((OntRelationMention_Type)jcasType).casFeatCode_rangeList), i);
    jcasType.ll_cas.ll_setRefArrayValue(jcasType.ll_cas.ll_getRefValue(addr, ((OntRelationMention_Type)jcasType).casFeatCode_rangeList), i, jcasType.ll_cas.ll_getFSRef(v));}
105064
private void gaussNewton( double betas[] ) {

		A_temp.reshape(L_full.numRows, numControl);
		v_temp.reshape(L_full.numRows, 1);
		x.reshape(numControl,1,false);

		// don't check numControl inside in hope that the JVM can optimize the code better
		if( numControl == 4 ) {
			for( int i = 0; i < numIterations; i++ ) {
				UtilLepetitEPnP.jacobian_Control4(L_full,betas, A_temp);
				UtilLepetitEPnP.residuals_Control4(L_full,y,betas,v_temp.data);

				if( !solver.setA(A_temp) )
					return;

				solver.solve(v_temp,x);

				for( int j = 0; j < numControl; j++ ) {
					betas[j] -= x.data[j];
				}
			}
		} else {
			for( int i = 0; i < numIterations; i++ ) {
				UtilLepetitEPnP.jacobian_Control3(L_full,betas, A_temp);
				UtilLepetitEPnP.residuals_Control3(L_full,y,betas,v_temp.data);

				if( !solver.setA(A_temp) )
					return;

				solver.solve(v_temp,x);

				for( int j = 0; j < numControl; j++ ) {
					betas[j] -= x.data[j];
				}
			}
		}
	}
827249
@Override
    public IDocumentQuery<T> withinRadiusOf(String fieldName, double radius, double latitude, double longitude) {
        return withinRadiusOf(fieldName, radius, latitude, longitude, null, Constants.Documents.Indexing.Spatial.DEFAULT_DISTANCE_ERROR_PCT);
    }
1236354
public static double quantile(double val, double mu, double sigma, double xi) {
    if(val < 0.0 || val > 1.0) {
      return Double.NaN;
    }
    if(xi == 0.) {
      return mu - sigma * FastMath.log(1 - val);
    }
    return mu - sigma / xi * (1 - FastMath.pow(1 - val, -xi));
  }
91112
protected static double calcPosterior(double f, double alpha, double mu, double sigma, double lambda) {
    final double pi = calcP_i(f, mu, sigma);
    final double qi = calcQ_i(f, lambda);
    return (alpha * pi) / (alpha * pi + (1.0 - alpha) * qi);
  }
45017
def _set_mode_kappa_Chaput(self, i_sigma, i_temp, weights):
        

        X = self._get_X(i_temp, weights, self._gv).ravel()
        num_ir_grid_points = len(self._ir_grid_points)
        num_band = self._primitive.get_number_of_atoms() * 3
        size = num_ir_grid_points * num_band * 3
        v = self._collision_matrix[i_sigma, i_temp].reshape(size, size)
        solver = _select_solver(self._pinv_solver)
        if solver in [1, 2, 4, 5]:
            v = v.T
        e = self._get_eigvals_pinv(i_sigma, i_temp)
        t = self._temperatures[i_temp]

        omega_inv = np.empty(v.shape, dtype='double', order='C')
        np.dot(v, (e * v).T, out=omega_inv)
        Y = np.dot(omega_inv, X)
        self._set_f_vectors(Y, num_ir_grid_points, weights)
        elems = ((0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1))
        for i, vxf in enumerate(elems):
            mat = self._get_I(vxf[0], vxf[1], num_ir_grid_points * num_band)
            self._mode_kappa[i_sigma, i_temp, :, :, i] = 0
            if mat is not None:
                np.dot(mat, omega_inv, out=mat)
                # vals = (X ** 2 * np.diag(mat)).reshape(-1, 3).sum(axis=1)
                # vals = vals.reshape(num_ir_grid_points, num_band)
                # self._mode_kappa[i_sigma, i_temp, :, :, i] = vals
                w = diagonalize_collision_matrix(mat,
                                                 pinv_solver=self._pinv_solver,
                                                 log_level=self._log_level)
                if solver in [1, 2, 4, 5]:
                    mat = mat.T
                spectra = np.dot(mat.T, X) ** 2 * w
                for s, eigvec in zip(spectra, mat.T):
                    vals = s * (eigvec ** 2).reshape(-1, 3).sum(axis=1)
                    vals = vals.reshape(num_ir_grid_points, num_band)
                    self._mode_kappa[i_sigma, i_temp, :, :, i] += vals

        factor = self._conversion_factor * Kb * t ** 2
        self._mode_kappa[i_sigma, i_temp] *= factor
820180
private static void randomGaussian( Random rand , double sigma , int radius , Point2D_I32 pt ) {

		int x,y;

		while( true ) {
			x = (int)(rand.nextGaussian()*sigma);
			y = (int)(rand.nextGaussian()*sigma);
			if( Math.sqrt(x*x + y*y) < radius )
				break;
		}

		pt.set(x,y);
	}
891083
private void checkRange(long index) {
    Assert.indexNot(isEmpty(), "segment is empty");
    Assert.indexNot(index < firstIndex(), index + " is less than the first index in the segment");
    Assert.indexNot(index > lastIndex(), index + " is greater than the last index in the segment");
  }