Compute the gradients.

@param pij Desired affinity matrix
@param qij Projected affinity matrix
@param qij_isum Normalization factor
@param sol Current solution coordinates
@param meta Point metadata
507642
protected void computeGradient(AffinityMatrix pij, double[][] qij, double qij_isum, double[][] sol, double[] meta) {
    final int dim3 = dim * 3;
    int size = pij.size();
    for(int i = 0, off = 0; i < size; i++, off += dim3) {
      final double[] sol_i = sol[i], qij_i = qij[i];
      Arrays.fill(meta, off, off + dim, 0.); // Clear gradient only
      for(int j = 0; j < size; j++) {
        if(i == j) {
          continue;
        }
        final double[] sol_j = sol[j];
        final double qij_ij = qij_i[j];
        // Qij after scaling!
        final double q = MathUtil.max(qij_ij * qij_isum, MIN_QIJ);
        double a = 4 * (pij.get(i, j) - q); // SNE gradient
        for(int k = 0; k < dim; k++) {
          meta[off + k] += a * (sol_i[k] - sol_j[k]);
        }
      }
    }
  }
1049357
def compute_laplacian_matrix(affinity_matrix, method='auto', **kwargs):
    """"""
    if method == 'auto':
        method = 'geometric'
    return Laplacian.init(method, **kwargs).laplacian_matrix(affinity_matrix)
973071
public static int cudnnDivisiveNormalizationForward(
        cudnnHandle handle, 
        cudnnLRNDescriptor normDesc, 
        int mode, 
        Pointer alpha, 
        cudnnTensorDescriptor xDesc, /** same desc for means, temp, temp2 */
        Pointer x, 
        Pointer means, /** if NULL, means are assumed to be zero */
        Pointer temp, 
        Pointer temp2, 
        Pointer beta, 
        cudnnTensorDescriptor yDesc, 
        Pointer y)
    {
        return checkResult(cudnnDivisiveNormalizationForwardNative(handle, normDesc, mode, alpha, xDesc, x, means, temp, temp2, beta, yDesc, y));
    }
1129454
public boolean refine(List<CameraPinhole> calibration , DMatrix4x4 Q ) {
		if( calibration.size() != cameras.size )
			throw new RuntimeException("Calibration and cameras do not match");

		computeNumberOfCalibrationParameters();
		func = new ResidualK();

		if( func.getNumOfInputsN() > 6*calibration.size() )
			throw new IllegalArgumentException("Need more views to refine. eqs="+(3*calibration.size()+" unknowns="+func.getNumOfInputsN()));

		// Declared new each time to ensure all variables are properly zeroed


		// plane at infinity to the null space of Q
		ConvertDMatrixStruct.convert(Q,_Q);
		nullspace.process(_Q,1,p);
		CommonOps_DDRM.divide(p,p.get(3));
		p.numRows=3;

		// Convert input objects into an array which can be understood by optimization
		encode(calibration,p,param);

		// Configure optimization
//		optimizer.setVerbose(System.out,0);
		optimizer.setFunction(func,null); // Compute using a numerical Jacobian
		optimizer.initialize(param.data,converge.ftol,converge.gtol);

		// Tell it to run for at most 100 iterations
		if( !UtilOptimize.process(optimizer,converge.maxIterations) )
			return false;

		// extract solution
		decode(optimizer.getParameters(),calibration,p);
		recomputeQ(p,Q);

		return true;
	}
1084252
@Override
    public void setMetaCue(IMetaCue[] cuePoints) {
        Map<String, Object> cues = new HashMap<String, Object>();
        this.cuePoints = cuePoints;

        int j = 0;
        for (j = 0; j < this.cuePoints.length; j++) {
            cues.put(String.valueOf(j), this.cuePoints[j]);
        }

        //		"CuePoints", cuePointData
        //					'0',	MetaCue
        //							name, "test"
        //							type, "event"
        //							time, "0.1"
        //					'1',	MetaCue
        //							name, "test1"
        //							type, "event1"
        //							time, "0.5"

        this.put("cuePoints", cues);
    }
430886
private void computePixelFeatures() {
		for (int y = 0; y < derivX.height; y++) {
			int pixelIndex = y*derivX.width;
			int endIndex = pixelIndex+derivX.width;
			for (; pixelIndex < endIndex; pixelIndex++ ) {
				float dx = derivX.data[pixelIndex];
				float dy = derivY.data[pixelIndex];

				// angle from 0 to pi radians
				orientation.data[pixelIndex] = UtilAngle.atanSafe(dy,dx) + GrlConstants.F_PId2;
				// gradient magnitude
				magnitude.data[pixelIndex] = Math.sqrt(dx*dx + dy*dy);
			}
		}
	}
427828
def re_normalize_flux(self, kwargs_ps, norm_factor):
        
        for i, model in enumerate(self.point_source_type_list):
            if model == 'UNLENSED':
                kwargs_ps[i]['point_amp'] *= norm_factor
            elif model in ['LENSED_POSITION', 'SOURCE_POSITION']:
                if self._fixed_magnification_list[i] is True:
                    kwargs_ps[i]['source_amp'] *= norm_factor
                else:
                    kwargs_ps[i]['point_amp'] *= norm_factor
        return kwargs_ps
975961
private void parsePointMetadata(JsonParser jp) throws IOException, SQLException {
        jp.nextToken(); // FIELD_NAME coordinates        
        String coordinatesField = jp.getText();
        if (coordinatesField.equalsIgnoreCase(GeoJsonField.COORDINATES)) {
            jp.nextToken(); // START_ARRAY [ to parse the coordinate
            parseCoordinateMetadata(jp);            
        } else {
            throw new SQLException("Malformed GeoJSON file. Expected 'coordinates', found '" + coordinatesField + "'");
        }
    }
215339
@Override
    public boolean apply(final int[] mapping) {
        for (final int u : queryStereoIndices) {
            switch (queryTypes[u]) {
                case Tetrahedral:
                    if (!checkTetrahedral(u, mapping)) return false;
                    break;
                case Geometric:
                    if (!checkGeometric(u, otherIndex(u), mapping)) return false;
                    break;
            }
        }
        return true;
    }
378264
private float iterationSor(GrayF32 image1, GrayF32 deriv1X, GrayF32 deriv1Y,
							   int i, int ipx, int imx, int ipy, int imy) {
		float w = SOR_RELAXATION;

		// these variables could be precomputed once.  See equation 11
		float psid = psiData.data[i];
		float psig = gamma*psiGradient.data[i];

		float di = warpImage2.data[i] - image1.data[i];
		float dx2 = warpDeriv2X.data[i];
		float dy2 = warpDeriv2Y.data[i];
		float dxx2 = warpDeriv2XX.data[i];
		float dyy2 = warpDeriv2YY.data[i];
		float dxy2 = warpDeriv2XY.data[i];

		float Au = -psid*di*dx2 + alpha*divU.data[i]
				- psig*((dx2 - deriv1X.data[i])*warpDeriv2XX.data[i] + (dy2 - deriv1Y.data[i])*warpDeriv2XY.data[i]);

		float Av = -psid*di*dy2 + alpha*divV.data[i]
				- psig*((dx2 - deriv1X.data[i])*warpDeriv2XY.data[i] + (dy2 - deriv1Y.data[i])*warpDeriv2YY.data[i]);

		float Du = psid*dx2*dx2 + psig*(dxx2*dxx2 + dxy2*dxy2) + alpha*divD.data[i];
		float Dv = psid*dy2*dy2 + psig*(dyy2*dyy2 + dxy2*dxy2) + alpha*divD.data[i];
		float D = psid*dx2*dy2 + psig*(dxx2 + dyy2)*dxy2;

		// update the change in flow
		float psi_index = psiSmooth.data[i];
		float coef0 = 0.5f*(psiSmooth.data[ipx] + psi_index);
		float coef1 = 0.5f*(psiSmooth.data[imx] + psi_index);
		float coef2 = 0.5f*(psiSmooth.data[ipy] + psi_index);
		float coef3 = 0.5f*(psiSmooth.data[imy] + psi_index);

		float div_du = coef0 * du.data[ipx] + coef1 * du.data[imx] + coef2 * du.data[ipy] + coef3 * du.data[imy];
		float div_dv = coef0 * dv.data[ipx] + coef1 * dv.data[imx] + coef2 * dv.data[ipy] + coef3 * dv.data[imy] ;

		final float dui = du.data[i];
		final float dvi = dv.data[i];

		du.data[i] = (1f-w)*dui + w*(Au - D*dvi + alpha*div_du)/Du;
		dv.data[i] = (1f-w)*dvi + w*(Av - D*du.data[i] + alpha*div_dv)/Dv;

		return (du.data[i] - dui) * (du.data[i] - dui) + (dv.data[i] - dvi) * (dv.data[i] - dvi);
	}
1144997
private void parseParentGeometryMetadata(JsonParser jp) throws IOException, SQLException {
        if(jp.nextToken()!=JsonToken.VALUE_NULL){//START_OBJECT { in case of null geometry
        jp.nextToken(); // FIELD_NAME type     
        jp.nextToken(); //VALUE_STRING Point
        String geometryType = jp.getText();        
        parseGeometryMetadata(jp, geometryType);
        } 
        
    }
675501
private void computeGradientAndHessian(DMatrixRMaj param  )
    {
        // residuals = f(x) - y
        function.compute(param, residuals);

        computeNumericalJacobian(param,jacobian);

        CommonOps_DDRM.multTransA(jacobian, residuals, g);
        CommonOps_DDRM.multTransA(jacobian, jacobian,  H);

        CommonOps_DDRM.extractDiag(H,Hdiag);
    }
471876
protected static double[][] computePij(double[][] dist, double sigma, double initialScale) {
    final int size = dist.length;
    final double msigmasq = -.5 / (sigma * sigma);
    double[][] pij = new double[size][size];
    FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Computing affinities", size, LOG) : null;
    Duration timer = LOG.isStatistics() ? LOG.newDuration(GaussianAffinityMatrixBuilder.class.getName() + ".runtime.pijmatrix").begin() : null;
    MeanVariance mv = LOG.isStatistics() ? new MeanVariance() : null;
    for(int i = 0; i < size; i++) {
      double logP = computeH(i, dist[i], pij[i], msigmasq);
      if(mv != null) {
        mv.put(FastMath.exp(logP));
      }
      LOG.incrementProcessed(prog);
    }
    LOG.ensureCompleted(prog);
    if(LOG.isStatistics()) { // timer != null, mv != null
      LOG.statistics(timer.end());
      LOG.statistics(new DoubleStatistic(GaussianAffinityMatrixBuilder.class.getName() + ".perplexity.average", mv.getMean()));
      LOG.statistics(new DoubleStatistic(GaussianAffinityMatrixBuilder.class.getName() + ".perplexity.stddev", mv.getSampleStddev()));
    }
    // Scale pij to have the desired sum EARLY_EXAGGERATION
    double sum = 0.;
    for(int i = 1; i < size; i++) {
      final double[] pij_i = pij[i];
      for(int j = 0; j < i; j++) { // Nur über halbe Matrix!
        sum += (pij_i[j] += pij[j][i]); // Symmetrie herstellen
      }
    }
    // Scaling taken from original tSNE code:
    final double scale = initialScale / (2. * sum);
    for(int i = 1; i < size; i++) {
      final double[] pij_i = pij[i];
      for(int j = 0; j < i; j++) {
        pij_i[j] = pij[j][i] = MathUtil.max(pij_i[j] * scale, MIN_PIJ);
      }
    }
    return pij;
  }
1241877
def update(self, Pi, Tij):
         """
        from msmtools import analysis as msmana

        # update transition matrix by copy
        self._Tij = np.array(Tij)
        assert msmana.is_transition_matrix(self._Tij), 'Given transition matrix is not a stochastic matrix'
        assert self._Tij.shape[0] == self._nstates, 'Given transition matrix has unexpected number of states '
        # reset spectral decomposition
        self._spectral_decomp_available = False

        # check initial distribution
        assert np.all(Pi >= 0), 'Given initial distribution contains negative elements.'
        assert np.any(Pi > 0), 'Given initial distribution is zero'
        self._Pi = np.array(Pi) / np.sum(Pi)
860181
@Override
    public ZMatrixRMaj getQ(ZMatrixRMaj Q , boolean transposed ) {
        Q = UtilDecompositons_ZDRM.checkIdentity(Q,N,N);

        Arrays.fill(w,0,N*2,0);

        if( transposed ) {
            for( int j = N-2; j >= 0; j-- ) {
                QrHelperFunctions_ZDRM.extractHouseholderRow(QT,j,j+1,N,w,0);
                QrHelperFunctions_ZDRM.rank1UpdateMultL(Q, w, 0, gammas[j], j+1, j+1 , N);
            }
        } else {
            for( int j = N-2; j >= 0; j-- ) {
                QrHelperFunctions_ZDRM.extractHouseholderRow(QT,j,j+1,N,w,0);
                QrHelperFunctions_ZDRM.rank1UpdateMultR(Q, w, 0, gammas[j], j+1, j+1 , N, b);
            }
        }

        return Q;
    }
522472
boolean computeCalibrationPoints(SquareGrid grid) {
		calibrationPoints.reset();

		for (int row = 0; row < grid.rows-1; row++) {
			int offset = row%2;
			for (int col = offset; col < grid.columns; col += 2) {
				SquareNode a = grid.get(row,col);

				if( col > 0 ) {
					SquareNode b = grid.get(row+1,col-1);
					if( !setIntersection(a,b,calibrationPoints.grow()))
						return false;
				}

				if( col < grid.columns-1) {
					SquareNode b = grid.get(row+1,col+1);
					if( !setIntersection(a,b,calibrationPoints.grow()))
						return false;
				}
			}
		}

		return true;
	}
1283103
def get_pij_matrix(t, diag, A, A_inv):
    
    return A.dot(np.diag(np.exp(diag * t))).dot(A_inv)