Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix
512530
def get_delta_plus_epsilon(nat, i_now, dr_now, dr_old):
    
    XIJ = get_XIJ(nat, i_now, dr_now, dr_old)
    YIJ = get_YIJ(nat, i_now, dr_old)

    YIJ_invert = array_inverse(YIJ)

    # Perform sum_k X_ik Y_jk^-1
    epsilon = np.sum(XIJ.reshape(-1,3,1,3)*YIJ_invert.reshape(-1,1,3,3), axis=3)

    return epsilon
822350
def update_gradients_full(self, dL_dK, X, X2=None):
    #def dK_dtheta(self, dL_dK, X, X2, target):
        """"""
        X,slices = X[:,:-1],index_to_slices(X[:,-1])
        if X2 is None:
            X2,slices2 = X,slices
            K = np.zeros((X.shape[0], X.shape[0]))
        else:
            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
        
        vyt = self.variance_Yt
        vyx = self.variance_Yx

        lyt = 1./(2*self.lengthscale_Yt)
        lyx = 1./(2*self.lengthscale_Yx)

        a = self.a
        b = self.b
        c = self.c

        tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
        xdist = (X[:,1][:,None] - X2[:,1][None,:])**2
        #rdist = [tdist,xdist]
        ttdist = (X[:,0][:,None] - X2[:,0][None,:])
        
        rd=tdist.shape[0]

        dka = np.zeros([rd,rd])
        dkb = np.zeros([rd,rd])
        dkc = np.zeros([rd,rd])
        dkYdvart = np.zeros([rd,rd])
        dkYdvarx = np.zeros([rd,rd])
        dkYdlent = np.zeros([rd,rd])
        dkYdlenx = np.zeros([rd,rd])


        kyy = lambda tdist,xdist: np.exp(-lyt*(tdist) -lyx*(xdist))
        #k1 = lambda tdist: (lyt - lyt**2 * (tdist) )
        #k2 = lambda xdist: ( lyx**2 * (xdist)  - lyx )
        #k3 = lambda xdist: ( 3*lyx**2 - 6*xdist*lyx**3 + xdist**2*lyx**4 )
        #k4 = lambda tdist: -lyt*np.sqrt(tdist)

        k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )

        k2 = lambda xdist: ( 4*lyx**2 * (xdist)  - 2*lyx )

        k3 = lambda xdist: ( 3*4*lyx**2 - 6*8*xdist*lyx**3 + 16*xdist**2*lyx**4 )

        k4 = lambda ttdist: 2*lyt*(ttdist)

        dkyydlyx = lambda tdist,xdist: kyy(tdist,xdist)*(-xdist)
        dkyydlyt = lambda tdist,xdist: kyy(tdist,xdist)*(-tdist)

        dk1dlyt = lambda tdist: 2. - 4*2.*lyt*tdist
        dk2dlyx = lambda xdist: (4.*2.*lyx*xdist -2.)
        dk3dlyx = lambda xdist: (6.*4.*lyx - 18.*8*xdist*lyx**2 + 4*16*xdist**2*lyx**3)

        dk4dlyt = lambda ttdist: 2*(ttdist)

        for i, s1 in enumerate(slices):
            for j, s2 in enumerate(slices2):
                for ss1 in s1:
                    for ss2 in s2:
                        if i==0 and j==0:
                            dka[ss1,ss2] = 0
                            dkb[ss1,ss2] = 0
                            dkc[ss1,ss2] = 0
                            dkYdvart[ss1,ss2] = vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdvarx[ss1,ss2] = vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])
                        elif i==0 and j==1:
                            dka[ss1,ss2] = -k2(xdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkb[ss1,ss2] = k4(ttdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkc[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #dkYdvart[ss1,ss2] = 0
                            #dkYdvarx[ss1,ss2] = 0
                            #dkYdlent[ss1,ss2] = 0
                            #dkYdlenx[ss1,ss2] = 0
                            dkYdvart[ss1,ss2] = (-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdvarx[ss1,ss2] = (-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])* (-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)+\
                            vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*b*dk4dlyt(ttdist[ss1,ss2])
                            dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)+\
                            vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*dk2dlyx(xdist[ss1,ss2]))
                        elif i==1 and j==1:
                            dka[ss1,ss2] = (2*a*k3(xdist[ss1,ss2]) - 2*c*k2(xdist[ss1,ss2]))*vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkb[ss1,ss2] = 2*b*k1(tdist[ss1,ss2])*vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkc[ss1,ss2] = (-2*a*k2(xdist[ss1,ss2]) + 2*c )*vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdvart[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdvarx[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )*vyt* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])*( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 ) +\
                            vyx*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*b**2*dk1dlyt(tdist[ss1,ss2])
                            dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])*( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 ) +\
                            vyx*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])* (-2*a*c*dk2dlyx(xdist[ss1,ss2]) + a**2*dk3dlyx(xdist[ss1,ss2]) )
                        else:
                            dka[ss1,ss2] = -k2(xdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkb[ss1,ss2] = -k4(ttdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkc[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #dkYdvart[ss1,ss2] = 0
                            #dkYdvarx[ss1,ss2] = 0
                            #dkYdlent[ss1,ss2] = 0
                            #dkYdlenx[ss1,ss2] = 0
                            dkYdvart[ss1,ss2] = (-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdvarx[ss1,ss2] = (-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])* (-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)+\
                            vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*(-1)*b*dk4dlyt(ttdist[ss1,ss2])
                            dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)+\
                            vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*dk2dlyx(xdist[ss1,ss2])) 

        self.a.gradient = np.sum(dka * dL_dK)  

        self.b.gradient = np.sum(dkb * dL_dK) 

        self.c.gradient = np.sum(dkc * dL_dK)


        self.variance_Yt.gradient = np.sum(dkYdvart * dL_dK)  # Vy

        self.variance_Yx.gradient = np.sum(dkYdvarx * dL_dK)

        self.lengthscale_Yt.gradient = np.sum(dkYdlent*(-0.5*self.lengthscale_Yt**(-2)) * dL_dK)    #ly np.sum(dktheta2*(-self.lengthscale_Y**(-2)) * dL_dK) 

        self.lengthscale_Yx.gradient =  np.sum(dkYdlenx*(-0.5*self.lengthscale_Yx**(-2)) * dL_dK)
346347
def select_delta(dist_post_update, current_iteration,
                 clip_max, clip_min, d, theta, constraint):
  
  if current_iteration == 1:
    delta = 0.1 * (clip_max - clip_min)
  else:
    if constraint == 'l2':
      delta = np.sqrt(d) * theta * dist_post_update
    elif constraint == 'linf':
      delta = d * theta * dist_post_update

  return delta
686531
def update_gradients_full(self, dL_dK, X, X2=None):
        """"""
        if X2 is None: X2 = X
        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
        FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)

        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega))
        Lo = np.column_stack((self.basis_omega,self.basis_omega))
        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2))
        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
        Gint = self._int_computation( r,omega,phi, r,omega,phi)

        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]

        #dK_dvar
        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T)

        #dK_dlen
        da_dlen = [-1./self.lengthscale**2,0.]
        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega))
        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
        dGint_dlen = dGint_dlen + dGint_dlen.T
        dG_dlen = 1./2*Gint + self.lengthscale/2*dGint_dlen
        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T)

        #dK_dper
        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
        dFX2_dper = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X2,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X2)

        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period))
        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)

        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
        # SIMPLIFY!!!       IPPprim1 = (self.upper - self.lower)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
        IPPprim = np.where(np.logical_or(np.isnan(IPPprim1), np.isinf(IPPprim1)), IPPprim2, IPPprim1)


        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
        #IPPint2[0,0] = (self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)

        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period))
        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2))
        r2,omega2,phi2 = dLa_dper2.T,Lo[:,0:1],dLp_dper2.T

        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi)
        dGint_dper = dGint_dper + dGint_dper.T

        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]

        dG_dper = 1./self.variance*(self.lengthscale/2*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)))

        dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T)

        self.variance.gradient = np.sum(dK_dvar*dL_dK)
        self.lengthscale.gradient = np.sum(dK_dlen*dL_dK)
        self.period.gradient = np.sum(dK_dper*dL_dK)
1338767
def approximate_gradient(decision_function, sample, num_evals,
                         delta, constraint, shape, clip_min, clip_max):
  
  # Generate random vectors.
  noise_shape = [num_evals] + list(shape)
  if constraint == 'l2':
    rv = np.random.randn(*noise_shape)
  elif constraint == 'linf':
    rv = np.random.uniform(low=-1, high=1, size=noise_shape)

  axis = tuple(range(1, 1 + len(shape)))
  rv = rv / np.sqrt(np.sum(rv ** 2, axis=axis, keepdims=True))
  perturbed = sample + delta * rv
  perturbed = clip_image(perturbed, clip_min, clip_max)
  rv = (perturbed - sample) / delta

  # query the model.
  decisions = decision_function(perturbed)
  decision_shape = [len(decisions)] + [1] * len(shape)
  fval = 2 * decisions.astype(np_dtype).reshape(decision_shape) - 1.0

  # Baseline subtraction (when fval differs)
  if np.mean(fval) == 1.0:  # label changes.
    gradf = np.mean(rv, axis=0)
  elif np.mean(fval) == -1.0:  # label not change.
    gradf = - np.mean(rv, axis=0)
  else:
    fval = fval - np.mean(fval)
    gradf = np.mean(fval * rv, axis=0)

  # Get the gradient direction.
  gradf = gradf / np.linalg.norm(gradf)

  return gradf
1014780
def _initialize_stretching_matrix(self):
        

        self.S = np.zeros((self.nz, self.nz))

        if (self.nz==2) and (self.rd) and (self.delta):

            self.del1 = self.delta/(self.delta+1.)
            self.del2 = (self.delta+1.)**-1
            self.Us = self.Ubg[0]-self.Ubg[1]

            self.F1 = self.rd**-2 / (1.+self.delta)
            self.F2 = self.delta*self.F1
            self.S[0,0], self.S[0,1] = -self.F1,  self.F1
            self.S[1,0], self.S[1,1] =  self.F2, -self.F2

        else:

            for i in range(self.nz):

                if i == 0:
                    self.S[i,i]   = -self.f2/self.Hi[i]/self.gpi[i]
                    self.S[i,i+1] =  self.f2/self.Hi[i]/self.gpi[i]

                elif i == self.nz-1:
                    self.S[i,i]   = -self.f2/self.Hi[i]/self.gpi[i-1]
                    self.S[i,i-1] =  self.f2/self.Hi[i]/self.gpi[i-1]

                else:
                    self.S[i,i-1] = self.f2/self.Hi[i]/self.gpi[i-1]
                    self.S[i,i]   = -(self.f2/self.Hi[i]/self.gpi[i] +
                                        self.f2/self.Hi[i]/self.gpi[i-1])
                    self.S[i,i+1] = self.f2/self.Hi[i]/self.gpi[i]
981499
protected void computeNumericalJacobian( DenseMatrix64F param ,
                                             DenseMatrix64F pt ,
                                             DenseMatrix64F deriv )
    {
        double invDelta = 1.0/DELTA;

        func.compute(param,pt, temp0);

        // compute the jacobian by perturbing the parameters slightly
        // then seeing how it effects the results.
        for( int i = 0; i < param.numRows; i++ ) {
            param.data[i] += DELTA;
            func.compute(param,pt, temp1);
            // compute the difference between the two parameters and divide by the delta
            add(invDelta,temp1,-invDelta,temp0,temp1);
            // copy the results into the jacobian matrix
            System.arraycopy(temp1.data,0,deriv.data,i*pt.numRows,pt.numRows);

            param.data[i] -= DELTA;
        }
    }
1316092
def bound_perturbed_gmres(pseudo, p, epsilon, deltas):
    '''
    '''
    if not numpy.all(numpy.array(deltas) > epsilon):
        raise ArgumentError('all deltas have to be greater than epsilon')

    bound = []
    for delta in deltas:
        # get boundary paths
        paths = pseudo.contour_paths(delta)

        # get vertices on boundary
        vertices = paths.vertices()

        # evaluate polynomial
        supremum = numpy.max(numpy.abs(p(vertices)))

        # compute bound
        bound.append(epsilon/(delta-epsilon)
                     * paths.length()/(2*numpy.pi*delta)
                     * supremum)
    return bound
445314
public static void gradientInner(GrayF32 ii, double tl_x, double tl_y, double samplePeriod ,
									 int regionSize, double kernelSize,
									 float[] derivX, float derivY[])
	{
		// add 0.5 to c_x and c_y to have it round when converted to an integer pixel
		// this is faster than the straight forward method
		tl_x += 0.5;
		tl_y += 0.5;

		// round the kernel size
		int w = (int)(kernelSize+0.5);
		int r = w/2;
		if( r <= 0 ) r = 1;
		w = r*2+1;

		int i = 0;
		for( int y = 0; y < regionSize; y++ ) {
			int pixelsY = (int)(tl_y + y * samplePeriod);
			int indexRow1 = ii.startIndex + (pixelsY-r-1)*ii.stride - r - 1;
			int indexRow2 = indexRow1 + r*ii.stride;
			int indexRow3 = indexRow2 + ii.stride;
			int indexRow4 = indexRow3 + r*ii.stride;

			for( int x = 0; x < regionSize; x++ , i++) {
				int pixelsX = (int)(tl_x + x * samplePeriod);

				final int indexSrc1 = indexRow1 + pixelsX;
				final int indexSrc2 = indexRow2 + pixelsX;
				final int indexSrc3 = indexRow3 + pixelsX;
				final int indexSrc4 = indexRow4 + pixelsX;

				final float p0 = ii.data[indexSrc1];
				final float p1 = ii.data[indexSrc1+r];
				final float p2 = ii.data[indexSrc1+r+1];
				final float p3 = ii.data[indexSrc1+w];
				final float p11 = ii.data[indexSrc2];
				final float p4 = ii.data[indexSrc2+w];
				final float p10 = ii.data[indexSrc3];
				final float p5 = ii.data[indexSrc3+w];
				final float p9 = ii.data[indexSrc4];
				final float p8 = ii.data[indexSrc4+r];
				final float p7 = ii.data[indexSrc4+r+1];
				final float p6 = ii.data[indexSrc4+w];

				final float left = p8-p9-p1+p0;
				final float right = p6-p7-p3+p2;
				final float top = p4-p11-p3+p0;
				final float bottom = p6-p9-p5+p10;

				derivX[i] = right-left;
				derivY[i] = bottom-top;
//				System.out.printf("%2d %2d %2d %2d dx = %6.2f  dy = %6.2f\n",x,y,pixelsX,pixelsY,derivX[i],derivY[i]);
			}
		}
	}
495819
def get_D_square_min(atoms_now, atoms_old, i_now, j_now, delta_plus_epsilon=None):
    
    nat = len(atoms_now)
    assert len(atoms_now) == len(atoms_old)

    pos_now = atoms_now.positions
    pos_old = atoms_old.positions

    # Compute current and old distance vectors. Note that current distance
    # vectors cannot be taken from the neighbor calculation, because neighbors
    # are calculated from the sheared cell while these distance need to come
    # from the unsheared cell. Taking the distance from the unsheared cell
    # make periodic boundary conditions (and flipping of cell) a lot easier.
    dr_now = mic(pos_now[i_now] - pos_now[j_now], atoms_now.cell)
    dr_old = mic(pos_old[i_now] - pos_old[j_now], atoms_old.cell)

    # Sanity check: Shape needs to be identical!
    assert dr_now.shape == dr_old.shape

    if delta_plus_epsilon is None:
        # Get minimum strain tensor
        delta_plus_epsilon = get_delta_plus_epsilon(nat, i_now, dr_now, dr_old)

    # Spread epsilon out for each neighbor index
    delta_plus_epsilon_n = delta_plus_epsilon[i_now]

    # Compute D^2_min
    d_sq_n = np.sum(
        (
        dr_now-
        np.sum(delta_plus_epsilon_n.reshape(-1,3,3)*dr_old.reshape(-1,1,3),
               axis=2)
        )**2,
        axis=1)

    # For each atom, sum over all neighbors
    d_sq = np.bincount(i_now, weights=d_sq_n)

    return delta_plus_epsilon, d_sq
614211
def _D_constraint(self, neg_pairs, w):
    
    diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
    diff_sq = diff * diff
    dist = np.sqrt(diff_sq.dot(w))
    sum_deri1 = np.einsum('ij,i', diff_sq, 0.5 / np.maximum(dist, 1e-6))
    sum_deri2 = np.einsum(
        'ij,ik->jk',
        diff_sq,
        diff_sq / (-4 * np.maximum(1e-6, dist**3))[:,None]
    )
    sum_dist = dist.sum()
    return (
      np.log(sum_dist),
      sum_deri1 / sum_dist,
      sum_deri2 / sum_dist - np.outer(sum_deri1, sum_deri1) / (sum_dist * sum_dist)
    )
1035052
def _npdelta(self, a, delta):
        """"""
        for x in np.nditer(a, op_flags=["readwrite"]):
            delta += x
            x[...] = delta
        return a
789125
function find_mdw_colw(collw) {
	var delta = Math.abs(collw - cycle_width(collw)), _MDW = MDW;
	if(delta > 0.005) for(MDW=MIN_MDW; MDW<MAX_MDW; ++MDW) if(Math.abs(collw - cycle_width(collw)) <= delta) { delta = Math.abs(collw - cycle_width(collw)); _MDW = MDW; }
	MDW = _MDW;
}
1041459
def gradient(self):
        """"""
        result = np.zeros((self.numc, 3), float)
        for index1 in range(self.numc):
            result[index1] = self.gradient_component(index1)
        return result
463897
def gradF_t(X, Y, S, M_E, E, m0, rho):
    ''' 
    '''
    n, r = X.shape
    m, r = Y.shape

    XS = np.dot(X, S)
    YS = np.dot(Y, S.T)
    XSY = np.dot(XS, Y.T)

    Qx = np.dot(np.dot(X.T, ((M_E - XSY) * E)), YS) / n
    Qy = np.dot(np.dot(Y.T, ((M_E - XSY) * E).T), XS) / m

    W = np.dot((XSY - M_E) * E, YS) + np.dot(X, Qx) + rho * Gp(X, m0, r)
    Z = np.dot(((XSY - M_E) * E).T, XS) + np.dot(Y, Qy) + rho * Gp(Y, m0, r)

    return W, Z
877143
def from_deformation(cls, deformation):
        
        dfm = Deformation(deformation)
        return cls(0.5 * (np.dot(dfm.trans, dfm) - np.eye(3)))
452358
private void calcRadiation( int i, int j, WritableRaster demWR, WritableRaster sOmbraWR, WritableRaster insolationWR,
            double[] sunVector, WritableRaster gradientWR, double mr ) {
        double z = demWR.getSampleDouble(i, j, 0);
        double pressure = ATM * Math.exp(-0.0001184 * z);
        double ma = mr * pressure / ATM;
        double temp = 273 + pLapse * (z - 4000);
        double vap_psat = Math.exp(26.23 - 5416.0 / temp);
        double wPrec = 0.493 * pRH * vap_psat / temp;
        double taur = Math.exp((-.09030 * Math.pow(ma, 0.84)) * (1.0 + ma - Math.pow(ma, 1.01)));
        double d = pCmO3 * mr;
        double tauo = 1 - (0.1611 * d * Math.pow(1.0 + 139.48 * d, -0.3035) - 0.002715 * d)
                / (1.0 + 0.044 * d + 0.0003 * Math.pow(d, 2));
        double taug = Math.exp(-0.0127 * Math.pow(ma, 0.26));
        double tauw = 1 - 2.4959 * (wPrec * mr) / (1.0 + 79.034 * (wPrec * mr) * 0.6828 + 6.385 * (wPrec * mr));
        double taua = Math.pow((0.97 - 1.265 * Math.pow(pVisibility, (-0.66))), Math.pow(ma, 0.9));

        double In = 0.9751 * SOLARCTE * taur * tauo * taug * tauw * taua;

        double cosinc = scalarProduct(sunVector, gradientWR.getPixel(i, j, new double[3]));

        if (cosinc < 0) {
            cosinc = 0;
        }
        double tmp = insolationWR.getSampleDouble(i, j, 0);
        insolationWR.setSample(i, j, 0, In * cosinc * sOmbraWR.getSampleDouble(i, j, 0) / 1000 + tmp);
    }