Sobolev type FoM enforcing both gradient and absolute similarity.
962750
def fom(reco, true_image):
    """"""
    gradient = odl.Gradient(reco.space)
    return (gradient(reco - true_image).norm() +
            reco.space.dist(reco, true_image))
645123
def smooth(data, fw):
    """"""
    if fw == 0:
        fdata = data
    else:
        fdata = lfilter(np.ones(fw)/fw, 1, data)
    return fdata
371479
function renormalize(wt) {

    // If a fixed tree size is specified, scale x and y based on the extent.
    // Compute the left-most, right-most, and depth-most nodes for extents.
    if (size != null) {
      var left = wt,
          right = wt,
          bottom = wt;
      var toVisit = [wt],
          node;
      while (node = toVisit.pop()) {
        var t = node.t;
        if (t.x < left.t.x) left = node;
        if (t.x > right.t.x) right = node;
        if (t.depth > bottom.t.depth) bottom = node;
        if (node.children) 
          toVisit = toVisit.concat(node.children);
      }

      var sep = separation == null ? 0.5 : separation(left.t, right.t)/2;
      var tx = sep - left.t.x;
      var kx = size[0] / (right.t.x + sep + tx);
      var ky = size[1] / (bottom.t.depth > 0 ? bottom.t.depth : 1);

      toVisit = [wt];
      while (node = toVisit.pop()) {
        var t = node.t;
        t.x = (t.x + tx) * kx;
        t.y = t.depth * ky;
        if (setNodeSizes) {
          t.x_size *= kx;
          t.y_size *= ky;
        }
        if (node.children) 
          toVisit = toVisit.concat(node.children);
      }
    }

    // Else either a fixed node size, or node size function was specified.
    // In this case, we translate such that the root node is at x = 0.
    else {
      var rootX = wt.t.x;
      moveRight(wt, -rootX);
    }
  }
375536
private void computePsiSmooth(GrayF32 ux , GrayF32 uy , GrayF32 vx , GrayF32 vy ,
								  GrayF32 psiSmooth ) {
		int N = derivFlowUX.width * derivFlowUX.height;

		for( int i = 0; i < N; i++ ) {
			float vux = ux.data[i];
			float vuy = uy.data[i];
			float vvx = vx.data[i];
			float vvy = vy.data[i];

			float mu = vux*vux + vuy*vuy;
			float mv = vvx*vvx + vvy*vvy;

			psiSmooth.data[i] = (float)(1.0/(2.0*Math.sqrt(mu + mv + EPSILON*EPSILON)));
		}
	}
329850
def update_gradients_full(self, dL_dK, X, X2=None):
                """"""
                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

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

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

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

                dkdubias = np.zeros([rd,rd])

                kyy = lambda tdist: np.exp(-lyt*(tdist))
                dkyydlyt = lambda tdist: kyy(tdist)*(-tdist)

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

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

                dk1dlyt = lambda tdist: 2. - 4*2.*lyt*tdist

                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:
                                    dkYdvart[ss1,ss2] = kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])
                                    dkdubias[ss1,ss2] = 0
                                elif i==0 and j==1:
                                    dkYdvart[ss1,ss2] = (k4(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    #dkYdvart[ss1,ss2] = ((2*lyt*ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (k4(ttdist[ss1,ss2])+1.)+\
                                    vyt*kyy(tdist[ss1,ss2])*(dk4dlyt(ttdist[ss1,ss2]))
                                    #dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (2*lyt*(ttdist[ss1,ss2])+1.)+\
                                    #vyt*kyy(tdist[ss1,ss2])*(2*ttdist[ss1,ss2])
                                    dkdubias[ss1,ss2] = 0
                                elif i==1 and j==1:
                                    dkYdvart[ss1,ss2] = (k1(tdist[ss1,ss2]) + 1. )* kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])*( k1(tdist[ss1,ss2]) + 1. ) +\
                          			vyt*kyy(tdist[ss1,ss2])*dk1dlyt(tdist[ss1,ss2])
                                    dkdubias[ss1,ss2] = 1
                                else:
                                    dkYdvart[ss1,ss2] = (-k4(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    #dkYdvart[ss1,ss2] = (-2*lyt*(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (-k4(ttdist[ss1,ss2])+1.)+\
                                    vyt*kyy(tdist[ss1,ss2])*(-dk4dlyt(ttdist[ss1,ss2]) )
                                    dkdubias[ss1,ss2] = 0
                                    #dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (-2*lyt*(ttdist[ss1,ss2])+1.)+\
                                    #vyt*kyy(tdist[ss1,ss2])*(-2)*(ttdist[ss1,ss2])
   

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

                self.lengthscale_Yt.gradient =  np.sum(dkYdlent*(-0.5*self.lengthscale_Yt**(-2)) * dL_dK)

                self.ubias.gradient = np.sum(dkdubias * dL_dK)
804933
public static void hessianNaive(GrayF32 integral, int skip , int size ,
									GrayF32 intensity)
	{
		final int w = intensity.width;
		final int h = intensity.height;

		// get convolution kernels for the second order derivatives
		IntegralKernel kerXX = DerivativeIntegralImage.kernelDerivXX(size,null);
		IntegralKernel kerYY = DerivativeIntegralImage.kernelDerivYY(size,null);
		IntegralKernel kerXY = DerivativeIntegralImage.kernelDerivXY(size,null);

		float norm = 1.0f/(size*size);

		for( int y = 0; y < h; y++ ) {
			for( int x = 0; x < w; x++ ) {

				int xx = x*skip;
				int yy = y*skip;

				computeHessian(integral, intensity, kerXX, kerYY, kerXY, norm, y, yy, x, xx);
			}
		}
	}
645002
def adapt(self, d, x):
        
        y = np.dot(self.w, x)
        e = d - y
        self.eps = self.eps - self.ro * self.mu * e * self.last_e * \
            np.dot(x, self.last_x) / \
            (np.dot(self.last_x, self.last_x) + self.eps)**2
        nu = self.mu / (self.eps + np.dot(x, x))
        self.w += nu * e * x        
        self.last_e = e
83464
def gradient(self):
        """"""
        functional = self

        class KLCrossEntCCGradient(Operator):

            """The gradient operator of this functional."""

            def __init__(self):
                """Initialize a new instance."""
                super(KLCrossEntCCGradient, self).__init__(
                    functional.domain, functional.domain, linear=False)

            def _call(self, x):
                """Apply the gradient operator to the given point."""
                if functional.prior is None:
                    return self.domain.element(np.exp(x))
                else:
                    return functional.prior * np.exp(x)

        return KLCrossEntCCGradient()
300078
def gradient_component(self, index1):
        """"""
        result = np.zeros(3, float)
        for index2 in range(self.numc):
            if self.scaling[index1, index2] > 0:
                for (se, ve), (sg, vg) in zip(self.yield_pair_energies(index1, index2), self.yield_pair_gradients(index1, index2)):
                    result += (sg*self.directions[index1, index2]*ve + se*vg)*self.scaling[index1, index2]
        return result
1327593
def _msge_with_gradient_overdetermined(data, delta, xvschema, skipstep, p):
    
    t, m, l = data.shape
    d = None
    l, k = 0, 0
    nt = np.ceil(t / skipstep)
    for trainset, testset in xvschema(t, skipstep):

        a, b = _construct_var_eqns(atleast_3d(data[trainset, :, :]), p)
        c, d = _construct_var_eqns(atleast_3d(data[testset, :, :]), p)

        e = sp.linalg.inv(np.eye(a.shape[1]) * delta ** 2 + a.T.dot(a))

        ba = b.transpose().dot(a)
        dc = d.transpose().dot(c)
        bae = ba.dot(e)
        baee = bae.dot(e)
        baecc = bae.dot(c.transpose().dot(c))

        l += np.sum(baecc * bae - 2 * bae * dc) + np.sum(d ** 2)
        k += np.sum(baee * dc - baecc * baee) * 4 * delta

    return l / (nt * d.size), k / (nt * d.size)
1366201
def _twofilter_smoothing_ON2(self, t, ti, info, phi, lwinfo):
        
        sp, sw = 0., 0.
        upb = lwinfo.max() + self.wgt[t].lw.max()
        if hasattr(self.model, 'upper_bound_trans'):
            upb += self.model.upper_bound_trans(t + 1)
        # Loop over n, to avoid having in memory a NxN matrix
        for n in range(self.N):
            omegan = np.exp(lwinfo + self.wgt[t].lw[n] - upb
                            + self.model.logpt(t + 1, self.X[t][n],
                                               info.hist.X[ti]))
            sp += np.sum(omegan * phi(self.X[t][n], info.hist.X[ti]))
            sw += np.sum(omegan)
        return sp / sw
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]);
			}
		}
	}
1048578
def dstep(self, W):
        

        # Compute residual X D - S in frequency domain
        Ryf = sl.inner(self.Zf, self.Df, axis=self.cri.axisM) - self.Sf
        # Transform to spatial domain, apply mask, and transform back to
        # frequency domain
        Ryf[:] = sl.rfftn(W * sl.irfftn(Ryf, self.cri.Nv, self.cri.axisN),
                          None, self.cri.axisN)
        # Compute gradient
        gradf = sl.inner(np.conj(self.Zf), Ryf, axis=self.cri.axisK)

        # If multiple channel signal, single channel dictionary
        if self.cri.C > 1 and self.cri.Cd == 1:
            gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True)

        # Update gradient step
        self.eta = self.eta_a / (self.j + self.eta_b)

        # Compute gradient descent
        self.Gf[:] = self.Df - self.eta * gradf
        self.G = sl.irfftn(self.Gf, self.cri.Nv, self.cri.axisN)

        # Eval proximal operator
        self.Dprv[:] = self.D
        self.D[:] = self.Pcn(self.G)
941733
def smeft_evolve_continuous(C_in, scale_in, scale_out, newphys=True, **kwargs):
    """"""
    sol = _smeft_evolve(C_in, scale_in, scale_out, newphys=newphys,
                        dense_output=True, **kwargs)
    @np.vectorize
    def _rge_solution(scale):
        t = log(scale)
        y = sol.sol(t).view(complex)
        yd = C_array2dict(y)
        yw = arrays2wcxf_nonred(yd)
        return yw
    def rge_solution(scale):
        # this is to return a scalar if the input is scalar
        return _rge_solution(scale)[()]
    return rge_solution
690710
def _smeft_evolve(C_in, scale_in, scale_out, newphys=True, **kwargs):
    """"""
    def fun(t0, y):
        return beta.beta_array(C=C_array2dict(y.view(complex)),
                               newphys=newphys).view(float) / (16 * pi**2)
    y0 = C_dict2array(C_in).view(float)
    sol = solve_ivp(fun=fun,
                    t_span=(log(scale_in), log(scale_out)),
                    y0=y0, **kwargs)
    return sol
975116
def hybrid_forward(self, F, x, sampled_values, label, weight, bias):
        """"""
        sampled_candidates, _, _ = sampled_values
        # (batch_size,)
        label = F.reshape(label, shape=(-1,))
        # (num_sampled+batch_size,)
        ids = F.concat(sampled_candidates, label, dim=0)
        # lookup weights and biases
        # (num_sampled+batch_size, dim)
        w_all = F.Embedding(data=ids, weight=weight,
                            input_dim=self._num_classes, output_dim=self._in_unit,
                            sparse_grad=self._sparse_grad)
        # (num_sampled+batch_size, 1)
        b_all = F.take(bias, indices=ids)
        return self._dense(x, sampled_values, label, w_all, b_all)
1214304
def grad_eigh(ans, x, UPLO='L'):
    """"""
    N = x.shape[-1]
    w, v = ans              # Eigenvalues, eigenvectors.
    def vjp(g):
        wg, vg = g          # Gradient w.r.t. eigenvalues, eigenvectors.
        w_repeated = anp.repeat(w[..., anp.newaxis], N, axis=-1)
        off_diag = anp.ones((N, N)) - anp.eye(N)
        F = off_diag / (T(w_repeated) - w_repeated + anp.eye(N))
        return _dot(v * wg[..., anp.newaxis, :] + _dot(v, F * _dot(T(v), vg)), T(v))
    return vjp