Python torch.inverse() Examples

def compute_L_inverse(self,X,Y):
        N = X.size()[0] # num of points (along dim 0)
        # construct matrix K
        Xmat = X.expand(N,N)
        Ymat = Y.expand(N,N)
        P_dist_squared = torch.pow(Xmat-Xmat.transpose(0,1),2)+torch.pow(Ymat-Ymat.transpose(0,1),2)
        P_dist_squared[P_dist_squared==0]=1 # make diagonal 1 to avoid NaN in log computation
        K = torch.mul(P_dist_squared,torch.log(P_dist_squared))
        # construct matrix L
        O = torch.FloatTensor(N,1).fill_(1)
        Z = torch.FloatTensor(3,3).fill_(0)       
        P = torch.cat((O,X,Y),1)
        L = torch.cat((torch.cat((K,P),1),torch.cat((P.transpose(0,1),Z),1)),0)
        Li = torch.inverse(L)
        if self.use_cuda:
            Li = Li.cuda()
        return Li 

def LAFMagicFro(LAFs1, LAFs2, H1to2, xy_th  = 5.0, scale_log = 0.4):
    LHF2_in_1 = reprojectLAFs(LAFs2, torch.inverse(H1to2), True)
    LHF1 = LAFs_to_H_frames(LAFs1)
    idxs_in1, idxs_in_2 = get_closest_correspondences_idxs(LHF1, LHF2_in_1, xy_th, scale_log)
    if len(idxs_in1) == 0:
        print('Warning, no correspondences found')
        return None
    LHF1_good = LHF1[idxs_in1,:,:]
    LHF2_good = LHF2_in_1[idxs_in_2,:,:]
    scales1 = get_LHFScale(LHF1_good);
    scales2 = get_LHFScale(LHF2_good);
    max_scale = torch.max(scales1,scales2);
    min_scale = torch.min(scales1, scales2);
    mean_scale = 0.5 * (max_scale + min_scale)
    eps = 1e-12;
    dist_loss = (torch.sqrt((LHF1_good.view(-1,9) - LHF2_good.view(-1,9))**2 + eps) / V(mean_scale.data).view(-1,1).expand(LHF1_good.size(0),9)).mean(dim=1); 
    loss = dist_loss;
    #print dist_loss, scale_loss, shape_loss
    return loss, idxs_in1, idxs_in_2, LHF2_in_1[:,0:2,:] 

def _vertex_decimation(self, L):

        max_eigenvec = self._power_iteration(L)
        v_plus, v_minus = (max_eigenvec >= 0).squeeze(), (max_eigenvec < 0).squeeze()

        # print(v_plus, v_minus)

        # diagonal matrix, swap v_minus with v_plus not to incur in errors (does not change the matrix)
        if torch.sum(v_plus) == 0.:  # The matrix is diagonal, cannot reduce further
            if torch.sum(v_minus) == 0.:
                assert v_minus.shape[0] == L.shape[0], (v_minus.shape, L.shape)
                # I assumed v_minus should have ones, but this is not necessarily the case. So I added this if
                return torch.ones(v_minus.shape), L
            else:
                return v_minus, L

        L_plus_plus = L[v_plus][:, v_plus]
        L_plus_minus = L[v_plus][:, v_minus]
        L_minus_minus = L[v_minus][:, v_minus]
        L_minus_plus = L[v_minus][:, v_plus]

        L_new = L_plus_plus - torch.mm(torch.mm(L_plus_minus, torch.inverse(L_minus_minus)), L_minus_plus)

        return v_plus, L_new 

def affineAug(img, max_add = 0.5):
    img_s = img.squeeze()
    h,w = img_s.size()
    ### Generate A
    A = torch.eye(3)
    rand_add = max_add *(torch.rand(3,3) - 0.5) * 2.0
    ##No perspective change
    rand_add[2,0:2] = 0
    rand_add[2,2] = 0;
    A  = A + rand_add
    denormA = Grid2PxA(w,h)
    normA = Px2GridA(w, h)
    if img.is_cuda:
        A = A.cuda()
        denormA = denormA.cuda()
        normA = normA.cuda()
    grid = torch.nn.functional.affine_grid(A[0:2,:].unsqueeze(0), torch.Size((1,1,h,w)))
    H_Orig2New = torch.mm(torch.mm(denormA, torch.inverse(A)), normA)
    new_img = torch.nn.functional.grid_sample(img_s.float().unsqueeze(0).unsqueeze(0),  grid)  
    return new_img, H_Orig2New, 

def affineAug(img, max_add = 0.5):
    img_s = img.squeeze()
    h,w = img_s.size()
    ### Generate A
    A = torch.eye(3)
    rand_add = max_add *(torch.rand(3,3) - 0.5) * 2.0
    ##No perspective change
    rand_add[2,0:2] = 0
    rand_add[2,2] = 0;
    A  = A + rand_add
    denormA = Grid2PxA(w,h)
    normA = Px2GridA(w, h)
    if img.is_cuda:
        A = A.cuda()
        denormA = denormA.cuda()
        normA = normA.cuda()
    grid = torch.nn.functional.affine_grid(A[0:2,:].unsqueeze(0), torch.Size((1,1,h,w)))
    H_Orig2New = torch.mm(torch.mm(denormA, torch.inverse(A)), normA)
    new_img = torch.nn.functional.grid_sample(img_s.float().unsqueeze(0).unsqueeze(0),  grid)  
    return new_img, H_Orig2New, 

def LAFMagicFro(LAFs1, LAFs2, H1to2, xy_th  = 5.0, scale_log = 0.4):
    LHF2_in_1 = reprojectLAFs(LAFs2, torch.inverse(H1to2), True)
    LHF1 = LAFs_to_H_frames(LAFs1)
    idxs_in1, idxs_in_2 = get_closest_correspondences_idxs(LHF1, LHF2_in_1, xy_th, scale_log)
    if len(idxs_in1) == 0:
        print('Warning, no correspondences found')
        return None
    LHF1_good = LHF1[idxs_in1,:,:]
    LHF2_good = LHF2_in_1[idxs_in_2,:,:]
    scales1 = get_LHFScale(LHF1_good);
    scales2 = get_LHFScale(LHF2_good);
    max_scale = torch.max(scales1,scales2);
    min_scale = torch.min(scales1, scales2);
    mean_scale = 0.5 * (max_scale + min_scale)
    eps = 1e-12;
    dist_loss = (torch.sqrt((LHF1_good.view(-1,9) - LHF2_good.view(-1,9))**2 + eps) / V(mean_scale.data).view(-1,1).expand(LHF1_good.size(0),9)).mean(dim=1); 
    loss = dist_loss;
    #print dist_loss, scale_loss, shape_loss
    return loss, idxs_in1, idxs_in_2, LHF2_in_1[:,0:2,:] 

def compute_L_inverse(self,X,Y):
        N = X.size()[0] # num of points (along dim 0)
        # construct matrix K
        Xmat = X.expand(N,N)
        Ymat = Y.expand(N,N)
        P_dist_squared = torch.pow(Xmat-Xmat.transpose(0,1),2)+torch.pow(Ymat-Ymat.transpose(0,1),2)
        P_dist_squared[P_dist_squared==0]=1 # make diagonal 1 to avoid NaN in log computation
        K = torch.mul(P_dist_squared,torch.log(P_dist_squared))
        if self.reg_factor != 0:
            K+=torch.eye(K.size(0),K.size(1))*self.reg_factor
        # construct matrix L
        O = torch.FloatTensor(N,1).fill_(1)
        Z = torch.FloatTensor(3,3).fill_(0)       
        P = torch.cat((O,X,Y),1)
        L = torch.cat((torch.cat((K,P),1),torch.cat((P.transpose(0,1),Z),1)),0)
        Li = torch.inverse(L)
        if self.use_cuda:
            Li = Li.cuda()
        return Li 

def LAFMagicFro(LAFs1, LAFs2, H1to2, xy_th  = 5.0, scale_log = 0.4):
    LHF2_in_1 = reprojectLAFs(LAFs2, torch.inverse(H1to2), True)
    LHF1 = LAFs_to_H_frames(LAFs1)
    idxs_in1, idxs_in_2 = get_closest_correspondences_idxs(LHF1, LHF2_in_1, xy_th, scale_log)
    if len(idxs_in1) == 0:
        print('Warning, no correspondences found')
        return None
    LHF1_good = LHF1[idxs_in1,:,:]
    LHF2_good = LHF2_in_1[idxs_in_2,:,:]
    scales1 = get_LHFScale(LHF1_good);
    scales2 = get_LHFScale(LHF2_good);
    max_scale = torch.max(scales1,scales2);
    min_scale = torch.min(scales1, scales2);
    mean_scale = 0.5 * (max_scale + min_scale)
    eps = 1e-12;
    dist_loss = (torch.sqrt((LHF1_good.view(-1,9) - LHF2_good.view(-1,9))**2 + eps) / V(mean_scale.data).view(-1,1).expand(LHF1_good.size(0),9)).mean(dim=1); 
    loss = dist_loss;
    #print dist_loss, scale_loss, shape_loss
    return loss, idxs_in1, idxs_in_2, LHF2_in_1[:,0:2,:] 

def forward(self, inputs, cond_inputs=None, mode='direct'):
        if str(self.L_mask.device) != str(self.L.device):
            self.L_mask = self.L_mask.to(self.L.device)
            self.U_mask = self.U_mask.to(self.L.device)
            self.I = self.I.to(self.L.device)
            self.P = self.P.to(self.L.device)
            self.sign_S = self.sign_S.to(self.L.device)

        L = self.L * self.L_mask + self.I
        U = self.U * self.U_mask + torch.diag(
            self.sign_S * torch.exp(self.log_S))
        W = self.P @ L @ U

        if mode == 'direct':
            return inputs @ W, self.log_S.sum().unsqueeze(0).unsqueeze(
                0).repeat(inputs.size(0), 1)
        else:
            return inputs @ torch.inverse(
                W), -self.log_S.sum().unsqueeze(0).unsqueeze(0).repeat(
                    inputs.size(0), 1) 

def test_inverse_cache_is_used(self):
        self.transform.eval()
        self.transform.use_cache()

        self.transform.inverse(self.inputs)
        self.assertTrue(self.transform.weight_inverse.called)
        self.assertTrue(self.transform.logabsdet.called)
        self.transform.weight_inverse.reset_mock()
        self.transform.logabsdet.reset_mock()

        outputs, logabsdet = self.transform.inverse(self.inputs)
        # Cached values should be used.
        self.assertFalse(self.transform.weight_inverse.called)
        self.assertFalse(self.transform.logabsdet.called)
        self.assertEqual(outputs, self.outputs_inv)
        self.assertEqual(logabsdet, self.logabsdet_inv) 

def forward(self, inputs, cond_inputs=None, mode='direct', logdets=None):
        """ Performs a forward or backward pass for flow modules.
        Args:
            inputs: a tuple of inputs and logdets
            mode: to run direct computation or inverse
        """
        self.num_inputs = inputs.size(-1)

        if logdets is None:
            logdets = torch.zeros(inputs.size(0), 1, device=inputs.device)

        assert mode in ['direct', 'inverse']
        if mode == 'direct':
            for module in self._modules.values():
                inputs, logdet = module(inputs, cond_inputs, mode)
                logdets += logdet
        else:
            for module in reversed(self._modules.values()):
                inputs, logdet = module(inputs, cond_inputs, mode)
                logdets += logdet

        return inputs, logdets 

def setUp(self):
        features = 5
        batch_size = 10

        weight = torch.randn(features, features)
        inverse = torch.randn(features, features)
        logabsdet = torch.randn(1)
        self.transform = Linear(features)
        self.transform.bias.data = torch.randn(features)  # Just so bias isn't zero.

        self.inputs = torch.randn(batch_size, features)
        self.outputs_fwd = self.inputs @ weight.t() + self.transform.bias
        self.outputs_inv = (self.inputs - self.transform.bias) @ inverse.t()
        self.logabsdet_fwd = logabsdet * torch.ones(batch_size)
        self.logabsdet_inv = (-logabsdet) * torch.ones(batch_size)

        # Mocks for abstract methods.
        self.transform.forward_no_cache = MagicMock(
            return_value=(self.outputs_fwd, self.logabsdet_fwd))
        self.transform.inverse_no_cache = MagicMock(
            return_value=(self.outputs_inv, self.logabsdet_inv))
        self.transform.weight = MagicMock(return_value=weight)
        self.transform.weight_inverse = MagicMock(return_value=inverse)
        self.transform.logabsdet = MagicMock(return_value=logabsdet) 

def test_inverse(self):
        features = 100
        batch_size = 50

        for num_transforms in [1, 2, 11, 12]:
            with self.subTest(num_transforms=num_transforms):
                transform = orthogonal.HouseholderSequence(
                    features=features, num_transforms=num_transforms)
                matrix = transform.matrix()
                inputs = torch.randn(batch_size, features)
                outputs, logabsdet = transform.inverse(inputs)
                self.assert_tensor_is_good(outputs, [batch_size, features])
                self.assert_tensor_is_good(logabsdet, [batch_size])
                self.eps = 1e-5
                self.assertEqual(outputs, inputs @ matrix)
                self.assertEqual(logabsdet, utils.logabsdet(matrix) * torch.ones(batch_size)) 

def affineAug(img, max_add = 0.5):
    img_s = img.squeeze()
    h,w = img_s.size()
    ### Generate A
    A = torch.eye(3)
    rand_add = max_add *(torch.rand(3,3) - 0.5) * 2.0
    ##No perspective change
    rand_add[2,0:2] = 0
    rand_add[2,2] = 0;
    A  = A + rand_add
    denormA = Grid2PxA(w,h)
    normA = Px2GridA(w, h)
    if img.is_cuda:
        A = A.cuda()
        denormA = denormA.cuda()
        normA = normA.cuda()
    grid = torch.nn.functional.affine_grid(A[0:2,:].unsqueeze(0), torch.Size((1,1,h,w)))
    H_Orig2New = torch.mm(torch.mm(denormA, torch.inverse(A)), normA)
    new_img = torch.nn.functional.grid_sample(img_s.float().unsqueeze(0).unsqueeze(0),  grid)  
    return new_img, H_Orig2New, 

def _compute_w(self, XS, YS_inner):
        '''
            Use Newton's method to obtain w from support set XS, YS_inner
            https://github.com/bertinetto/r2d2/blob/master/fewshots/models/lrd2.py
        '''

        for i in range(self.iters):
            # use eta to store w_{i-1}^T X
            if i == 0:
                eta = torch.zeros_like(XS[:,0])  # support_size
            else:
                eta = (XS @ w).squeeze(1)

            mu = torch.sigmoid(eta)
            s = mu * (1 - mu)
            z = eta + (YS_inner - mu) / s
            Sinv = torch.diag(1.0/s)

            # Woodbury with regularization
            w = XS.t() @ torch.inverse(XS @ XS.t() + (10. ** self.lam) * Sinv) @ z.unsqueeze(1)

        return w 

def get_barycentric_coords(point, verts):
        if len(verts) == 2:
            diff = verts[1] - verts[0]
            diff_norm = torch.norm(diff)
            normalized_diff = diff / diff_norm
            u = torch.dot(verts[1] - point, normalized_diff) / diff_norm
            v = torch.dot(point - verts[0], normalized_diff) / diff_norm
            return u, v
        elif len(verts) == 3:
            # TODO Area method instead of LinAlg
            M = torch.cat([
                torch.cat([verts[0], verts[0].new_ones(1)]).unsqueeze(1),
                torch.cat([verts[1], verts[1].new_ones(1)]).unsqueeze(1),
                torch.cat([verts[2], verts[2].new_ones(1)]).unsqueeze(1),
            ], dim=1)
            invM = torch.inverse(M)
            uvw = torch.matmul(invM, torch.cat([point, point.new_ones(1)]).unsqueeze(1))
            return uvw
        else:
            raise ValueError('Barycentric coords only works for 2 or 3 points') 

def compute_L_inverse(self,X,Y):
        N = X.size()[0] # num of points (along dim 0)
        # construct matrix K
        Xmat = X.expand(N,N)
        Ymat = Y.expand(N,N)
        P_dist_squared = torch.pow(Xmat-Xmat.transpose(0,1),2)+torch.pow(Ymat-Ymat.transpose(0,1),2)
        P_dist_squared[P_dist_squared==0]=1 # make diagonal 1 to avoid NaN in log computation
        K = torch.mul(P_dist_squared,torch.log(P_dist_squared))
        # construct matrix L
        O = torch.FloatTensor(N,1).fill_(1)
        Z = torch.FloatTensor(3,3).fill_(0)       
        P = torch.cat((O,X,Y),1)
        L = torch.cat((torch.cat((K,P),1),torch.cat((P.transpose(0,1),Z),1)),0)
        Li = torch.inverse(L)
        if self.use_cuda:
            Li = Li.cuda()
        return Li 

def invH(H):
    """ Generate (H+damp)^{-1}, with predicted damping values
    :param approximate Hessian matrix JtWJ
    -----------
    :return the inverse of Hessian
    """
    # GPU is much slower for matrix inverse when the size is small (compare to CPU)
    # works (50x faster) than inversing the dense matrix in GPU
    if H.is_cuda:
        # invH = bpinv((H).cpu()).cuda()
        # invH = torch.inverse(H)
        invH = torch.inverse(H.cpu()).cuda()
    else:
        invH = torch.inverse(H)
    return invH 

def get_w_target_rr(data, vocab_size, ebd_model, w_target_lam):
    '''
        Compute the importance of every tokens in the support set
        Convert to Ridge Regression as it admits analytical solution.
        Using this explicit formula improve speed by 2x

        @return w: vocab_size * num_classes
    '''

    text_ebd = ebd_model(data)
    label = data['label'].clone()

    unique, inv_idx = torch.unique(label, sorted=True, return_inverse=True)
    new_label = torch.arange(len(unique), dtype=unique.dtype, device=unique.device)

    label = new_label[inv_idx]
    label_onehot = F.embedding(label,
            torch.eye(len(unique), dtype=torch.float, device=text_ebd.device))

    I = torch.eye(len(text_ebd), dtype=torch.float, device=text_ebd.device)

    w = text_ebd.t()\
            @ torch.inverse(text_ebd @ text_ebd.t() + w_target_lam * I)\
            @ label_onehot

    return w 

def _compute_w(self, XS, YS_onehot):
        '''
            Compute the W matrix of ridge regression
            @param XS: support_size x ebd_dim
            @param YS_onehot: support_size x way

            @return W: ebd_dim * way
        '''

        W = XS.t() @ torch.inverse(
                XS @ XS.t() + (10. ** self.lam) * self.I_support) @ YS_onehot

        return W 

def forward(self):

        constrainted_matrix = self.select_param()
        matrix_ = torch.squeeze(torch.squeeze(constrainted_matrix,dim=2),dim=2)
        matrix_t = torch.t(matrix_)
        matrixs = torch.mm(matrix_t,matrix_)
        trace_ = torch.trace(torch.mm(matrixs,torch.inverse(matrixs)))
        log_det = torch.logdet(matrixs)
        maha_loss = trace_ - log_det
        return maha_loss 

def transform(point, center, scale, resolution, invert=False):
    """Generate and affine transformation matrix.

    Given a set of points, a center, a scale and a targer resolution, the
    function generates and affine transformation matrix. If invert is ``True``
    it will produce the inverse transformation.

    Arguments:
        point {torch.tensor} -- the input 2D point
        center {torch.tensor or numpy.array} -- the center around which to perform the transformations
        scale {float} -- the scale of the face/object
        resolution {float} -- the output resolution

    Keyword Arguments:
        invert {bool} -- define wherever the function should produce the direct or the
        inverse transformation matrix (default: {False})
    """
    _pt = torch.ones(3)
    _pt[0] = point[0]
    _pt[1] = point[1]

    h = 200.0 * scale
    t = torch.eye(3)
    t[0, 0] = resolution / h
    t[1, 1] = resolution / h
    t[0, 2] = resolution * (-center[0] / h + 0.5)
    t[1, 2] = resolution * (-center[1] / h + 0.5)

    if invert:
        t = torch.inverse(t)

    new_point = (torch.matmul(t, _pt))[0:2]

    return new_point.int() 

def inverseLHFs(LHFs):
    LHF1_inv =torch.zeros(LHFs.size())
    if LHFs.is_cuda:
        LHF1_inv = LHF1_inv.cuda()
    for i in range(LHF1_inv.size(0)):
        LHF1_inv[i,:,:] = LHFs[i,:,:].inverse()
    return LHF1_inv 

def get_GT_correspondence_indexes_Fro(LAFs1,LAFs2, H1to2, dist_threshold = 4,
                                      skip_center_in_Fro = False):
    LHF2_in_1_pre = reprojectLAFs(LAFs2, torch.inverse(H1to2), True)
    LHF1_inv = inverseLHFs(LAFs_to_H_frames(LAFs1))
    frob_norm_dist = reproject_to_canonical_Frob_batched(LHF1_inv, LHF2_in_1_pre, batch_size = 2, skip_center = skip_center_in_Fro)
    min_dist, idxs_in_2 = torch.min(frob_norm_dist,1)
    plain_indxs_in1 = torch.arange(0, idxs_in_2.size(0))
    if LAFs1.is_cuda:
        plain_indxs_in1 = plain_indxs_in1.cuda()
    #print min_dist.min(), min_dist.max(), min_dist.mean()
    mask =  min_dist <= dist_threshold
    return min_dist[mask], plain_indxs_in1[mask], idxs_in_2[mask] 

def get_GT_correspondence_indexes(LAFs1, LAFs2, H1to2, dist_threshold = 4):    
    LHF2_in_1_pre = reprojectLAFs(LAFs2, torch.inverse(H1to2), True)
    just_centers1 = LAFs1[:,:,2];
    just_centers2_repr_to_1 = LHF2_in_1_pre[:,0:2,2];
    
    dist  = distance_matrix_vector(just_centers2_repr_to_1, just_centers1)
    min_dist, idxs_in_2 = torch.min(dist,1)
    plain_indxs_in1 = torch.arange(0, idxs_in_2.size(0))
    if LAFs1.is_cuda:
        plain_indxs_in1 = plain_indxs_in1.cuda()
    mask =  min_dist <= dist_threshold
    return min_dist[mask], plain_indxs_in1[mask], idxs_in_2[mask] 

def inverseLHFs(LHFs):
    LHF1_inv =torch.zeros(LHFs.size())
    if LHFs.is_cuda:
        LHF1_inv = LHF1_inv.cuda()
    for i in range(LHF1_inv.size(0)):
        LHF1_inv[i,:,:] = LHFs[i,:,:].inverse()
    return LHF1_inv 

def true_objective(theta1, theta2, ipd):
    p1 = torch.sigmoid(theta1)
    p2 = torch.sigmoid(theta2[[0,1,3,2,4]])
    p0 = (p1[0], p2[0])
    p = (p1[1:], p2[1:])
    # create initial laws, transition matrix and rewards:
    P0 = torch.stack(phi(*p0), dim=0).view(1,-1)
    P = torch.stack(phi(*p), dim=1)
    R = torch.from_numpy(ipd.payout_mat).view(-1,1).float()
    # the true value to optimize:
    objective = (P0.mm(torch.inverse(torch.eye(4) - hp.gamma*P))).mm(R)
    return -objective 

def true_objective(theta1, theta2, ipd):
    p1 = torch.sigmoid(theta1)
    p2 = torch.sigmoid(theta2[[0,1,3,2,4]])
    p0 = (p1[0], p2[0])
    p = (p1[1:], p2[1:])
    # create initial laws, transition matrix and rewards:
    P0 = torch.stack(phi(*p0), dim=0).view(1,-1)
    P = torch.stack(phi(*p), dim=1)
    R = torch.from_numpy(ipd.payout_mat).view(-1,1).float()
    # the true value to optimize:
    objective = (P0.mm(torch.inverse(torch.eye(4) - hp.gamma*P))).mm(R)
    return -objective 

def post_stabilization(self, world):
        v = world.get_v()
        M = world.M()
        Je = world.Je()
        Jc = None
        if world.contacts:
            Jc = world.Jc()
        ge = torch.matmul(Je, v)
        gc = None
        if Jc is not None:
            gc = torch.matmul(Jc, v) + torch.matmul(Jc, v) * -world.restitutions()

        u = torch.cat([Je.new_zeros(Je.size(1)), ge])
        if Jc is None:
            neq = Je.size(0) if Je.ndimension() > 0 else 0
            if neq > 0:
                P = torch.cat([torch.cat([M, -Je.t()], dim=1),
                               torch.cat([Je, Je.new_zeros(neq, neq)], dim=1)])
            else:
                P = M
            if self.cached_inverse is None:
                inv = torch.inverse(P)
            else:
                inv = self.cached_inverse
            x = torch.matmul(inv, u)
        else:
            v = gc
            Je = Je.unsqueeze(0)
            Jc = Jc.unsqueeze(0)
            h = u[:M.size(0)].unsqueeze(0)
            b = u[M.size(0):].unsqueeze(0)
            M = M.unsqueeze(0)
            v = v.unsqueeze(0)
            F = Jc.new_zeros(Jc.size(1), Jc.size(1)).unsqueeze(0)
            x = self.lcp_solver()(M, h, Jc, v, Je, b, F)
        dp = -x[:M.size(0)]
        return dp 

def to_norm_affine(affine, src_size, dst_size, align_corners: bool = False):
    """
    Given ``affine`` defined for coordinates in the pixel space, compute the corresponding affine
    for the normalized coordinates.

    Args:
        affine (torch Tensor): Nxdxd batched square matrix
        src_size (sequence of int): source image spatial shape
        dst_size (sequence of int): target image spatial shape
        align_corners: if True, consider -1 and 1 to refer to the centers of the
            corner pixels rather than the image corners.
            See also: https://pytorch.org/docs/stable/nn.functional.html#torch.nn.functional.grid_sample

    Raises:
        ValueError: affine must be a tensor
        ValueError: affine must be Nxdxd, got {tuple(affine.shape)}
        ValueError: affine suggests a {sr}-D transform, but the sizes are src_size={src_size}, dst_size={dst_size}

    """
    if not torch.is_tensor(affine):
        raise ValueError("affine must be a tensor")
    if affine.ndim != 3 or affine.shape[1] != affine.shape[2]:
        raise ValueError(f"affine must be Nxdxd, got {tuple(affine.shape)}")
    sr = affine.shape[1] - 1
    if sr != len(src_size) or sr != len(dst_size):
        raise ValueError(
            f"affine suggests a {sr}-D transform, but the sizes are src_size={src_size}, dst_size={dst_size}"
        )

    src_xform = normalize_transform(src_size, affine.device, affine.dtype, align_corners)
    dst_xform = normalize_transform(dst_size, affine.device, affine.dtype, align_corners)
    new_affine = src_xform @ affine @ torch.inverse(dst_xform)
    return new_affine