r"""For some classes (in particular $\Delta F=1$), only a subset of Wilson
    coefficients exist in WET-3 and WET-4. Therefore, depending on the number
    of flavours `f`, the dimensionality of the ADM has to be reduced.
322772
def get_permissible_wcs(classname, f):
    """
    # these are the problematic sectors
    classes = ['cu', 'db', 'sb', 'sd', 'mue', 'mutau', 'taue', 'dF0']
    if classname not in classes or f == 5:
        # for 5-flavour WET, nothing to do.
        # Neither for other classes (I, II, ...) because they exist either all
        # or not at all in WET-3 and WET-4 (they have specific flavours).
        return 'all'
    if f not in [3, 4]:
        raise ValueError("f must be 3, 4, or 5.")
    if classname == 'dF0':
        sector = 'dF=0'
    else:
        sector = classname
    perm_keys = wcxf.Basis['WET-{}'.format(f), 'JMS'].sectors[sector].keys()
    all_keys = coeffs[sector]
    return [i for i, c in enumerate(all_keys) if c in perm_keys]
557750
def sph_harmonic_tp(nrows, ncols, n, m):
    
    nuvec = pysphi.ynunm(n, m, n + 1)
    num = nuvec[1:] * (-1) ** m
    num = num[::-1]
    yvec = np.concatenate([num, nuvec])
    
    theta = np.linspace(0.0, np.pi, nrows)
    phi = np.linspace(0.0, 2.0 * np.pi - 2.0 * np.pi / ncols)
    nu = np.array(range(-n, n + 1), dtype=np.complex128)
    
    out = np.zeros([nrows, ncols], dtype=np.complex128)

    for nn in xrange(0, nrows):
        exp_vec = np.exp(1j * nu * theta[nn])
        for mm in xrange(0, ncols):          
            vl = np.sum(yvec * exp_vec)
            out[nn, mm] = 1j ** m * np.exp(1j * m * phi[mm]) * vl

    return out
201410
public static void multAddOuter( double alpha , DMatrix4x4 A , double beta , DMatrix4 u , DMatrix4 v , DMatrix4x4 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
    }
1277142
@Override
    public DMatrixRMaj getU(DMatrixRMaj U , boolean transpose , boolean compact ) {
        U = handleU(U, transpose, compact,m,n,min);
        CommonOps_DDRM.setIdentity(U);

        for( int i = 0; i < m; i++ ) u[i] = 0;

        for( int j = min-1; j >= 0; j-- ) {
            u[j] = 1;
            for( int i = j+1; i < m; i++ ) {
                u[i] = UBV.get(i,j);
            }
            if( transpose )
                QrHelperFunctions_DDRM.rank1UpdateMultL(U, u, gammasU[j], j, j, m);
            else
                QrHelperFunctions_DDRM.rank1UpdateMultR(U, u, gammasU[j], j, j, m, this.b);
        }

        return U;
    }
772988
def getUe(classname, eta_s, f, alpha_s, alpha_e, m_u, m_d, m_s, m_c, m_b, m_e, m_mu, m_tau):
    """"""
    args = f, m_u, m_d, m_s, m_c, m_b, m_e, m_mu, m_tau
    A = getattr(adm, 'adm_e_' + classname)(*args)
    perm_keys = get_permissible_wcs(classname, f)
    if perm_keys != 'all':
        # remove disallowed rows & columns if necessary
        A = A[perm_keys][:, perm_keys]
    w, v = admeig(classname, *args)
    b0s = 11 - 2 * f / 3
    a = w / (2 * b0s)
    K = np.linalg.inv(v) @ A.T @ v
    for i in range(K.shape[0]):
        for j in range(K.shape[1]):
            if a[i] - a[j] != 1:
                K[i, j] *= (eta_s**(a[j] + 1) - eta_s**a[i]) / (a[i] - a[j] - 1)
            else:
                K[i, j] *= eta_s**a[i] * log(1 / eta_s)
    return -alpha_e / (2 * b0s * alpha_s) * v @ K @ np.linalg.inv(v)
534290
def lm_deltaf(damping_times, modes):
    

    df = {}
    for lmn in modes:
        l, m, nmodes = int(lmn[0]), int(lmn[1]), int(lmn[2])
        for n in range(nmodes):
            df['%d%d%d' %(l,m,n)] = \
                1. / qnm_time_decay(damping_times['%d%d%d' %(l,m,n)], 1./1000)

    return min(df.values())
277094
@Override
  public MarginalSet computeMarginals(FactorGraph factorGraph) {
    VariableNumMap variables = factorGraph.getVariables();
    Preconditions.checkArgument(variables.getDiscreteVariables().size() == variables.size());

    // Initialize the mean field distribution to the uniform distribution over
    // all variables.
    long curTime = System.currentTimeMillis();
    System.out.println("init: " + (System.currentTimeMillis() - curTime)); 
    int numVars = variables.size();
    IndexedList<Integer> variableNums = new IndexedList<Integer>(variables.getVariableNums());
    List<DiscreteVariable> variableTypes = variables.getDiscreteVariables();
    List<Tensor> variableMarginals = Lists.newArrayListWithCapacity(numVars);
    for (int i = 0; i < numVars; i++) {
      int[] dimensions = new int[] { variableNums.get(i) };
      int[] sizes = new int[] { variableTypes.get(i).numValues() };
      variableMarginals.add(DenseTensor.constant(dimensions, sizes, 1.0 / sizes[0]));
    }
    System.out.println("logweights: "  + (System.currentTimeMillis() - curTime));
    // Get the log weights for each factor in the original factor graph.
    List<Tensor> logWeights = Lists.newArrayList();
    for (Factor factor : factorGraph.getFactors()) {
      logWeights.add(factor.coerceToDiscrete().getWeights().elementwiseLog());
    }
    System.out.println("done: " +  + (System.currentTimeMillis() - curTime)); 

    double updateL2 = Double.POSITIVE_INFINITY;
    int numIterations = 0;
    while (updateL2 > CONVERGENCE_DELTA) {
      updateL2 = 0.0;
      for (int i = 0; i < numVars; i++) {
        int curVarNum = variableNums.get(i);
        // Accumulate the messages from each factor containing this variable.
        DenseTensorBuilder messageAccumulator = new DenseTensorBuilder(
            variableMarginals.get(i).getDimensionNumbers(), variableMarginals.get(i).getDimensionSizes());
        Set<Integer> factorIndexes = factorGraph.getFactorsWithVariable(curVarNum);
        // System.out.println(i + " " + curVarNum + ": " + factorIndexes);
        for (int factorIndex : factorIndexes) {
          Tensor factorMessage = getFactorMessage(curVarNum, logWeights.get(factorIndex),
              variableMarginals, variableNums);
          // System.out.println("factorMessage: " +
          // Arrays.toString(factorMessage.getDimensionNumbers()));
          // System.out.println("messageAccumulator: " +
          // Arrays.toString(messageAccumulator.getDimensionNumbers()));
          messageAccumulator.increment(factorMessage);
        }

        // Update the marginal based on the inbound messages, setting
        // marginal equal to the logistic function of the accumulated messages.
        messageAccumulator.exp();
        double normalizingConstant = messageAccumulator.getTrace();
        messageAccumulator.multiply(1.0 / normalizingConstant);
        Tensor newMarginal = messageAccumulator.buildNoCopy();

        // Compute the convergence criteria
        Tensor delta = newMarginal.elementwiseAddition(
            variableMarginals.get(i).elementwiseProduct(SparseTensor.getScalarConstant(-1.0)));
        updateL2 += delta.getL2Norm();

        variableMarginals.set(i, newMarginal);
      }
      numIterations++;
      System.out.println("update L2: " + updateL2);
    }
    
    System.out.println("iterations:" + numIterations);
    // Format output as factors.
    List<Factor> marginals = Lists.newArrayList();
    for (int i = 0; i < numVars; i++) {
      marginals.add(new TableFactor(variables.intersection(variableNums.get(i)),
          variableMarginals.get(i)));
    }

    System.out.println("really done: " +  + (System.currentTimeMillis() - curTime));
    return new FactorMarginalSet(marginals, 1.0, factorGraph.getConditionedVariables(),
        factorGraph.getConditionedValues());
  }
651900
def check_semidefinite_positiveness(A):
    
    B = empty_like(A)
    B[:] = A
    B[diag_indices_from(B)] += sqrt(finfo(float).eps)
    try:
        cholesky(B)
    except LinAlgError:
        return False
    return True
231217
private static void gaussian(float a[][],
	                             int index[]) {
		int n = index.length;
		float c[] = new float[n];
		// Initialize the index
		for (int i = 0; i < n; ++i) {
			index[i] = i;
		}
		// Find the rescaling factors, one from each row
		for (int i = 0; i < n; ++i){
			float c1 = 0;
			for (int j = 0; j < n; ++j){
				float c0 = Math.abs(a[i][j]);
				if (c0 > c1) {
					c1 = c0;
				}
			}
			c[i] = c1;
		}
		// Search the pivoting element from each column
		int k = 0;
		for (int j = 0; j < n - 1; ++j){
			float pi1 = 0;
			for (int i = j; i < n; ++i){
				float pi0 = Math.abs(a[index[i]][j]);
				pi0 /= c[index[i]];
				if (pi0 > pi1) {

					pi1 = pi0;
					k = i;
				}
			}
			// Interchange rows according to the pivoting order
			int itmp = index[j];
			index[j] = index[k];
			index[k] = itmp;
			for (int i = j + 1; i < n; ++i){
				float pj = a[index[i]][j] / a[index[j]][j];
				// Record pivoting ratios below the diagonal
				a[index[i]][j] = pj;

				// Modify other elements accordingly
				for (int l = j + 1; l < n; ++l){
					a[index[i]][l] -= pj * a[index[j]][l];
				}
			}
		}
	}
776373
@Override
	public void computeHessian(DMatrixRMaj jacLeft , DMatrixRMaj jacRight) {
		A.reshape(jacLeft.numCols,jacLeft.numCols);
		B.reshape(jacLeft.numCols,jacRight.numCols);
		D.reshape(jacRight.numCols,jacRight.numCols);

		// take advantage of the inner product's symmetry when possible to reduce
		// the number of calculations
		MatrixMultProduct_DDRM.inner_reorder_lower(jacLeft,A);
		CommonOps_DDRM.symmLowerToFull(A);
		CommonOps_DDRM.multTransA(jacLeft,jacRight,B);
		MatrixMultProduct_DDRM.inner_reorder_lower(jacRight,D);
		CommonOps_DDRM.symmLowerToFull(D);
	}
1186035
public static void multAddOuter( double alpha , DMatrix6x6 A , double beta , DMatrix6 u , DMatrix6 v , DMatrix6x6 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a15 = alpha*A.a15 + beta*u.a1*v.a5;
        C.a16 = alpha*A.a16 + beta*u.a1*v.a6;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a25 = alpha*A.a25 + beta*u.a2*v.a5;
        C.a26 = alpha*A.a26 + beta*u.a2*v.a6;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a35 = alpha*A.a35 + beta*u.a3*v.a5;
        C.a36 = alpha*A.a36 + beta*u.a3*v.a6;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
        C.a45 = alpha*A.a45 + beta*u.a4*v.a5;
        C.a46 = alpha*A.a46 + beta*u.a4*v.a6;
        C.a51 = alpha*A.a51 + beta*u.a5*v.a1;
        C.a52 = alpha*A.a52 + beta*u.a5*v.a2;
        C.a53 = alpha*A.a53 + beta*u.a5*v.a3;
        C.a54 = alpha*A.a54 + beta*u.a5*v.a4;
        C.a55 = alpha*A.a55 + beta*u.a5*v.a5;
        C.a56 = alpha*A.a56 + beta*u.a5*v.a6;
        C.a61 = alpha*A.a61 + beta*u.a6*v.a1;
        C.a62 = alpha*A.a62 + beta*u.a6*v.a2;
        C.a63 = alpha*A.a63 + beta*u.a6*v.a3;
        C.a64 = alpha*A.a64 + beta*u.a6*v.a4;
        C.a65 = alpha*A.a65 + beta*u.a6*v.a5;
        C.a66 = alpha*A.a66 + beta*u.a6*v.a6;
    }
395334
def _vmomentdensity(self,R,z,n,m,o,nsigma=None,mc=False,nmc=10000,
                       _returnmc=False,_vrs=None,_vts=None,_vzs=None,
                       _rawgausssamples=False,
                       gl=False,ngl=_DEFAULTNGL,_returngl=False,_glqeval=None,
                       _return_actions=False,_jr=None,_lz=None,_jz=None,
                       _return_freqs=False,
                       _rg=None,_kappa=None,_nu=None,_Omega=None,
                       _sigmaR1=None,_sigmaz1=None,
                       **kwargs):
        """"""
        if isinstance(R,numpy.ndarray):
            return numpy.array([self._vmomentdensity(r,zz,n,m,o,nsigma=nsigma,
                                                    mc=mc,nmc=nmc,
                                                    gl=gl,ngl=ngl,**kwargs) for r,zz in zip(R,z)])
        if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                actionAngle.actionAngleAdiabaticGrid)):
            if n % 2 == 1. or o % 2 == 1.:
                return 0. #we know this must be the case
        if nsigma == None:
            nsigma= _NSIGMA
        if _sigmaR1 is None:
            sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr)
        else:
            sigmaR1= _sigmaR1
        if _sigmaz1 is None:
            sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz)
        else:
            sigmaz1= _sigmaz1
        thisvc= potential.vcirc(self._pot,R,use_physical=False)
        #Use the asymmetric drift equation to estimate va
        gamma= numpy.sqrt(0.5)
        va= sigmaR1**2./2./thisvc\
            *(gamma**2.-1. #Assume close to flat rotation curve, sigphi2/sigR2 =~ 0.5
               +R*(1./self._hr+2./self._hsr))
        if math.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center
        if gl:
            if ngl % 2 == 1:
                raise ValueError("ngl must be even")
            if not _glqeval is None and ngl != _glqeval.shape[0]:
                _glqeval= None
            #Use Gauss-Legendre integration for all
            if ngl == _DEFAULTNGL:
                glx, glw= self._glxdef, self._glwdef
                glx12, glw12= self._glxdef12, self._glwdef12
            elif ngl == _DEFAULTNGL2:
                glx, glw= self._glxdef2, self._glwdef2
                glx12, glw12= self._glxdef, self._glwdef
            else:
                glx, glw= numpy.polynomial.legendre.leggauss(ngl)
                glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2)
            #Evaluate everywhere
            if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                    actionAngle.actionAngleAdiabaticGrid)):
                vRgl= nsigma*sigmaR1/2.*(glx+1.)
                vzgl= nsigma*sigmaz1/2.*(glx+1.)
                vRglw= glw
                vzglw= glw
            else:
                vRgl= nsigma*sigmaR1/2.*(glx12+1.)
                #vRgl= 1.5/2.*(glx12+1.)
                vRgl= list(vRgl)
                vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.))
                #vRgl.extend(-1.5/2.*(glx12+1.))
                vRgl= numpy.array(vRgl)
                vzgl= nsigma*sigmaz1/2.*(glx12+1.)
                #vzgl= 1.5/2.*(glx12+1.)
                vzgl= list(vzgl)
                vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.))
                #vzgl.extend(-1.5/2.*(glx12+1.))
                vzgl= numpy.array(vzgl)
                vRglw= glw12
                vRglw= list(vRglw)
                vRglw.extend(glw12)
                vRglw= numpy.array(vRglw)
                vzglw= glw12
                vzglw= list(vzglw)
                vzglw.extend(glw12)
                vzglw= numpy.array(vzglw)
            if 'vTmax' in kwargs: vTmax = kwargs['vTmax']
            else:                 vTmax = 1.5
            vTgl= vTmax/2.*(glx+1.)
            #Tile everything
            vTgl= numpy.tile(vTgl,(ngl,ngl,1)).T
            vRgl= numpy.tile(numpy.reshape(vRgl,(1,ngl)).T,(ngl,1,ngl))
            vzgl= numpy.tile(vzgl,(ngl,ngl,1))
            vTglw= numpy.tile(glw,(ngl,ngl,1)).T #also tile weights
            vRglw= numpy.tile(numpy.reshape(vRglw,(1,ngl)).T,(ngl,1,ngl))
            vzglw= numpy.tile(vzglw,(ngl,ngl,1))
            #evaluate
            if _glqeval is None and _jr is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R+numpy.zeros(ngl*ngl*ngl),
                                           vRgl.flatten(),
                                           vTgl.flatten(),
                                           z+numpy.zeros(ngl*ngl*ngl),
                                           vzgl.flatten(),
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            elif not _jr is None and _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                                                 log=True,
                                                                 _return_actions=True,
                                                                 _return_freqs=True,
                                                                 use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            elif not _jr is None and not _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                           rg=_rg,kappa=_kappa,nu=_nu,
                                           Omega=_Omega,
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            else:
                logqeval= _glqeval
            if _returngl:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        logqeval)
            elif _return_actions and _return_freqs:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        jr,lz,jz,
                        rg,kappa,nu,Omega)
            elif _return_actions:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        jr,lz,jz)
            else:
                return numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                 *vTglw*vRglw*vzglw*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2)
        elif mc:
            mvT= (thisvc-va)/gamma/sigmaR1
            if _vrs is None:
                vrs= numpy.random.normal(size=nmc)
            else:
                vrs= _vrs
            if _vts is None:
                vts= numpy.random.normal(size=nmc)+mvT
            else:
                if _rawgausssamples:
                    vts= _vts+mvT
                else:
                    vts= _vts
            if _vzs is None:
                vzs= numpy.random.normal(size=nmc)
            else:
                vzs= _vzs
            Is= _vmomentsurfaceMCIntegrand(vzs,vrs,vts,numpy.ones(nmc)*R,
                                           numpy.ones(nmc)*z,
                                           self,sigmaR1,gamma,sigmaz1,mvT,
                                           n,m,o)
            if _returnmc:
                if _rawgausssamples:
                    return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o),
                        vrs,vts-mvT,vzs)
                else:
                    return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o),
                            vrs,vts,vzs)
            else:
                return numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
        else: #pragma: no cover because this is too slow; a warning is shown
            warnings.warn("Calculations using direct numerical integration using tplquad is not recommended and extremely slow; it has also not been carefully tested",galpyWarning)
            return integrate.tplquad(_vmomentsurfaceIntegrand,
                                     1./gamma*(thisvc-va)/sigmaR1-nsigma,
                                     1./gamma*(thisvc-va)/sigmaR1+nsigma,
                                     lambda x: 0., lambda x: nsigma,
                                     lambda x,y: 0., lambda x,y: nsigma,
                                     (R,z,self,sigmaR1,gamma,sigmaz1,n,m,o),
                                     **kwargs)[0]*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
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];
					}
				}
			}
		}
	}
63050
def unbound_dimensions(streams, kdims, no_duplicates=True):
    
    params = stream_parameters(streams, no_duplicates)
    return [d for d in kdims if d not in params]
866529
def dBinaryRochedy (r, D, q, F):
    
    return -r[1]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*r[1]*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5 + F*F*(1+q)*r[1]
37133
private int getAtomTypeXCount(IAtomContainer ac, IAtom atom) {
        List<IAtom> neighbours = ac.getConnectedAtomsList(atom);
        int nocounter = 0;
        IBond bond;
        for (IAtom neighbour : neighbours) {
            if ((neighbour.getSymbol().equals("N") || neighbour.getSymbol().equals("O"))
                    && !(Boolean) neighbour.getProperty("IS_IN_AROMATIC_RING")) {
                //if (ac.getMaximumBondOrder(neighbours[i]) == 1.0) {
                bond = ac.getBond(neighbour, atom);
                if (bond.getOrder() != IBond.Order.DOUBLE) {
                    nocounter += 1;
                }
            }
        }
        return nocounter;
    }
459067
def _Aterm(N, h, m, k, dW):
    """"""
    sqrt2h = np.sqrt(2.0/h)
    Xk = np.random.normal(0.0, 1.0, (N, m, 1))
    Yk = np.random.normal(0.0, 1.0, (N, m, 1))
    term1 = _dot(Xk, _t(Yk + sqrt2h*dW))
    term2 = _dot(Yk + sqrt2h*dW, _t(Xk))
    return (term1 - term2)/k