Single application of the adjoint projection.
966146
def _apply_adj(self, a):
        # is projection the zero operator?
        if self.V.shape[1] == 0:
            return numpy.zeros(a.shape)
        ''''''
        c = inner(self.V, a, ip_B=self.ip_B)
        if self.Q is not None and self.R is not None:
            c = self.Q.dot(scipy.linalg.solve_triangular(self.R.T.conj(), c,
                                                         lower=True))
        return self.W.dot(c)
1172222
def adjoint(self):
        
        return PointwiseInner(vfspace=self.range, vecfield=self.vecfield,
                              weighting=self.weights)
699930
def _multicomplex2(f, fx, x, h):
        """"""
        n = len(x)
        ee = np.diag(h)
        hess = np.outer(h, h)
        cmplx_wrap = Bicomplex.__array_wrap__
        for i in range(n):
            for j in range(i, n):
                zph = Bicomplex(x + 1j * ee[i, :], ee[j, :])
                hess[i, j] = cmplx_wrap(f(zph)).imag12 / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess
154216
def _one_iteration(self, F, Ybus, V, Vm, Va, pv, pq, pvpq):
        
        J = self._build_jacobian(Ybus, V, pv, pq, pvpq)

        # Update step.
        dx = -1 * spsolve(J, F)
#        dx = -1 * linalg.lstsq(J.todense(), F)[0]

        # Update voltage vector.
        npv = len(pv)
        npq = len(pq)
        if npv > 0:
            Va[pv] = Va[pv] + dx[range(npv)]
        if npq > 0:
            Va[pq] = Va[pq] + dx[range(npv, npv + npq)]
            Vm[pq] = Vm[pq] + dx[range(npv + npq, npv + npq + npq)]

        V = Vm * exp(1j * Va)
        Vm = abs(V) # Avoid wrapped round negative Vm.
        Va = angle(V)

        return V, Vm, Va
371693
private void _symmetricTql2 ()
  {
    // This is derived from the Algol procedures tql2, by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    for (int i = 1; i < m_nDim; i++)
      m_aEVe[i - 1] = m_aEVe[i];
    m_aEVe[m_nDim - 1] = 0.0;

    double f = 0.0;
    double tst1 = 0.0;
    for (int l = 0; l < m_nDim; l++)
    {
      // Find small subdiagonal element
      tst1 = Math.max (tst1, MathHelper.abs (m_aEVd[l]) + MathHelper.abs (m_aEVe[l]));
      int m = l;
      while (m < m_nDim)
      {
        if (MathHelper.abs (m_aEVe[m]) <= EPSILON * tst1)
          break;
        m++;
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.

      if (m > l)
      {
        do
        {
          // Compute implicit shift
          double g = m_aEVd[l];
          double p = (m_aEVd[l + 1] - g) / (2.0 * m_aEVe[l]);
          double r = MathHelper.hypot (p, 1.0);
          if (p < 0)
          {
            r = -r;
          }
          m_aEVd[l] = m_aEVe[l] / (p + r);
          m_aEVd[l + 1] = m_aEVe[l] * (p + r);
          final double dl1 = m_aEVd[l + 1];
          double h = g - m_aEVd[l];
          for (int i = l + 2; i < m_nDim; i++)
          {
            m_aEVd[i] -= h;
          }
          f += h;

          // Implicit QL transformation.
          p = m_aEVd[m];
          double c = 1.0;
          double c2 = c;
          double c3 = c;
          final double el1 = m_aEVe[l + 1];
          double s = 0.0;
          double s2 = 0.0;
          for (int i = m - 1; i >= l; i--)
          {
            c3 = c2;
            c2 = c;
            s2 = s;
            g = c * m_aEVe[i];
            h = c * p;
            r = MathHelper.hypot (p, m_aEVe[i]);
            m_aEVe[i + 1] = s * r;
            s = m_aEVe[i] / r;
            c = p / r;
            p = c * m_aEVd[i] - s * g;
            m_aEVd[i + 1] = h + s * (c * g + s * m_aEVd[i]);

            // Accumulate transformation.

            for (int k = 0; k < m_nDim; k++)
            {
              h = m_aEigenVector[k][i + 1];
              m_aEigenVector[k][i + 1] = s * m_aEigenVector[k][i] + c * h;
              m_aEigenVector[k][i] = c * m_aEigenVector[k][i] - s * h;
            }
          }
          p = -s * s2 * c3 * el1 * m_aEVe[l] / dl1;
          m_aEVe[l] = s * p;
          m_aEVd[l] = c * p;

          // Check for convergence.

        } while (MathHelper.abs (m_aEVe[l]) > EPSILON * tst1);
      }
      m_aEVd[l] = m_aEVd[l] + f;
      m_aEVe[l] = 0.0;
    }

    // Sort eigenvalues and corresponding vectors.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      int k = i;
      double p = m_aEVd[i];
      for (int j = i + 1; j < m_nDim; j++)
      {
        if (m_aEVd[j] < p)
        {
          k = j;
          p = m_aEVd[j];
        }
      }
      if (k != i)
      {
        m_aEVd[k] = m_aEVd[i];
        m_aEVd[i] = p;
        for (int j = 0; j < m_nDim; j++)
        {
          p = m_aEigenVector[j][i];
          m_aEigenVector[j][i] = m_aEigenVector[j][k];
          m_aEigenVector[j][k] = p;
        }
      }
    }
  }
20912
def adjoint(self):
        """"""
        if not self.is_linear:
            raise NotImplementedError('this operator is not linear and '
                                      'thus has no adjoint')

        forward_op = self

        class ResizingOperatorAdjoint(ResizingOperatorBase):

            """Adjoint of `ResizingOperator`.

            See `the online documentation
            <https://odlgroup.github.io/odl/math/resizing_ops.html>`_
            on resizing operators for mathematical details.
            """

            def _call(self, x, out):
                """Implement ``self(x, out)``."""
                with writable_array(out) as out_arr:
                    resize_array(x.asarray(), self.range.shape,
                                 offset=self.offset, pad_mode=self.pad_mode,
                                 pad_const=0, direction='adjoint',
                                 out=out_arr)

            @property
            def adjoint(self):
                """Adjoint of the adjoint, i.e. the original operator."""
                return forward_op

            @property
            def inverse(self):
                """(Pseudo-)Inverse of this operator.

                Note that in axes where ``self`` extends, the returned operator
                acts as a proper inverse, while in restriction axes, the
                operation is not invertible.
                """
                return ResizingOperatorAdjoint(
                    domain=self.range, range=self.domain,
                    pad_mode=self.pad_mode)

        return ResizingOperatorAdjoint(domain=self.range, range=self.domain,
                                       pad_mode=self.pad_mode)
1127369
def vector_projection(v1, v2):
    '''
    '''
    return scalar_projection(v1, v2) * v2 / np.linalg.norm(v2)
352964
def r(w, Xs):
    '''
    '''
    n = len(Xs)
    P = projection_matrix(w)
    c = C(w, Xs)

    return np.sqrt(sum(np.dot(c - X, np.dot(P, c - X)) for X in Xs) / n)
17299
def spm(x, model, y=None, n_samples=None, dx_min=-0.1,
        dx_max=0.1, n_dxs=5, dy_min=-0.1, dy_max=0.1, n_dys=5,
        angle_min=-30, angle_max=30, n_angles=31, black_border_size=0):
  
  if y is None:
    preds = model.get_probs(x)
    # Using model predictions as ground truth to avoid label leaking
    preds_max = reduce_max(preds, 1, keepdims=True)
    y = tf.to_float(tf.equal(preds, preds_max))
    y = tf.stop_gradient(y)
    del preds
  y = y / reduce_sum(y, 1, keepdims=True)

  # Define the range of transformations
  dxs = np.linspace(dx_min, dx_max, n_dxs)
  dys = np.linspace(dy_min, dy_max, n_dys)
  angles = np.linspace(angle_min, angle_max, n_angles)

  if n_samples is None:
    import itertools
    transforms = list(itertools.product(*[dxs, dys, angles]))
  else:
    sampled_dxs = np.random.choice(dxs, n_samples)
    sampled_dys = np.random.choice(dys, n_samples)
    sampled_angles = np.random.choice(angles, n_samples)
    transforms = zip(sampled_dxs, sampled_dys, sampled_angles)
  transformed_ims = parallel_apply_transformations(
      x, transforms, black_border_size)

  def _compute_xent(x):
    preds = model.get_logits(x)
    return tf.nn.softmax_cross_entropy_with_logits_v2(
        labels=y, logits=preds)

  all_xents = tf.map_fn(
      _compute_xent,
      transformed_ims,
      parallel_iterations=1)  # Must be 1 to avoid keras race conditions

  # Return the adv_x with worst accuracy

  # all_xents is n_total_samples x batch_size (SB)
  all_xents = tf.stack(all_xents)  # SB

  # We want the worst case sample, with the largest xent_loss
  worst_sample_idx = tf.argmax(all_xents, axis=0)  # B

  batch_size = tf.shape(x)[0]
  keys = tf.stack([
      tf.range(batch_size, dtype=tf.int32),
      tf.cast(worst_sample_idx, tf.int32)
  ], axis=1)
  transformed_ims_bshwc = tf.einsum('sbhwc->bshwc', transformed_ims)
  after_lookup = tf.gather_nd(transformed_ims_bshwc, keys)  # BHWC
  return after_lookup
412662
def _apply_projection(self, Av):
        '''
        PAv, UAp = self.projection.apply_complement(Av, return_Ya=True)
        self._UAps.append(UAp)
        c = UAp.copy()
        rhos = self.rhos
        if self.iter > 0:
            c -= (1 + rhos[-1]/rhos[-2])*self._UAps[-2]
        if self.iter > 1:
            c += rhos[-2]/rhos[-3]*self._UAps[-3]
        c *= ((-1)**self.iter) / numpy.sqrt(rhos[-1])
        if self.iter > 0:
            c -= numpy.sqrt(rhos[-2]/rhos[-1]) * self.C[:, [-1]]
        self.C = numpy.c_[self.C, c]
        return PAv
657239
def xstep(self):
        
        """

        self.YU[:] = self.Y - self.U
        self.block_sep0(self.YU)[:] += self.S
        Zf = sl.rfftn(self.YU, None, self.cri.axisN)
        Z0f = self.block_sep0(Zf)
        Z1f = self.block_sep1(Zf)

        DZ0f = np.conj(self.Df) * Z0f
        DZ0fBQ = sl.dot(self.B.dot(self.Q).T, DZ0f, axis=self.cri.axisC)
        Z1fQ = sl.dot(self.Q.T, Z1f, axis=self.cri.axisC)
        b = DZ0fBQ + Z1fQ

        Xh = sl.solvedbd_sm(self.gDf, (self.mu / self.rho) * self.GHGf + 1.0,
                            b, self.c, axis=self.cri.axisM)
        self.Xf[:] = sl.dot(self.Q, Xh, axis=self.cri.axisC)
        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            DDXf = np.conj(self.Df) *  sl.inner(self.Df, self.Xf,
                                                axis=self.cri.axisM)
            DDXfBB = sl.dot(self.B.T.dot(self.B), DDXf, axis=self.cri.axisC)
            ax = self.rho * (DDXfBB + self.Xf) + \
                 self.mu * self.GHGf * self.Xf
            b = self.rho * (sl.dot(self.B.T, DZ0f, axis=self.cri.axisC)
                            + Z1f)
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
423540
function hm(n, v, incv, offv, tau, b, offb) {
    // Apply householder transformation b = b - tau v v' b
    // LAPACK: DLARF

    if (tau === 0.0) {
      return;
    }

    // compute v' b
    var d = b[offb + 0];
    for (var i = 1; i < n; ++i) {
      d += b[offb + i] * v[offv + incv*i];
    }

    // compute b = b - tau v v' b
    b[offb + 0] -= tau * d;
    for (var i = 1; i < n; ++i) {
      b[offb + i] -= tau * v[offv + incv*i] * d;
    }
  }
643781
def fit_transform(
        self,
        X,
        projection="sum",
        scaler=preprocessing.MinMaxScaler(),
        distance_matrix=False,
    ):
        

        projections = projection
        scalers = scaler
        distance_matrices = distance_matrix

        # Turn single projection arguments into a pipeline
        if isinstance(projection, list) and isinstance(projection[0], int):
            projections = [projection]

        if not isinstance(projection, list):
            projections = [projection]

        # Turn single scaler arguments into a pipeline
        if not isinstance(scaler, list):
            scalers = [scaler]

        # Turn single distance matrix arguments into a pipeline
        if not isinstance(distance_matrix, list):
            distance_matrices = [distance_matrix]

        # set defaults to first list item, if not (correctly) set by the user
        if len(scalers) != len(projections):
            scalers = [scalers[0]] * len(projections)

        if len(distance_matrices) != len(projections):
            distance_matrices = [distance_matrices[0]] * len(projections)

        if self.verbose > 0:
            print("..Composing projection pipeline of length %s:" % (len(projections)))
            print("\tProjections: %s" % ("\n\t\t".join(map(str, projections))))
            print("\tDistance matrices: %s" % ("\n".join(map(str, distance_matrices))))
            print("\tScalers: %s" % ("\n".join(map(str, scalers))))

        # Pipeline Stack the projection functions
        lens = X
        for projection, scaler, distance_matrix in zip(
            projections, scalers, distance_matrices
        ):
            lens = self.project(
                lens,
                projection=projection,
                scaler=scaler,
                distance_matrix=distance_matrix,
            )

        return lens
602614
def laplacian_eigenmaps(self, num_dims=None, normed=True, val_thresh=1e-8):
    '''
    '''
    L = self.laplacian(normed=normed)
    return _null_space(L, num_dims, val_thresh, overwrite=True)
1304808
private SimpleMatrix computeAffineH( List<AssociatedPair> observations ,
										 DMatrixRMaj H , DMatrixRMaj Hzero ) {

		SimpleMatrix A = new SimpleMatrix(observations.size(),3);
		SimpleMatrix b = new SimpleMatrix(A.numRows(),1);

		Point2D_F64 c = new Point2D_F64();
		Point2D_F64 k = new Point2D_F64();

		for( int i = 0; i < observations.size(); i++ ) {
			AssociatedPair a = observations.get(i);

			GeometryMath_F64.mult(Hzero, a.p1, k);
			GeometryMath_F64.mult(H,a.p2,c);

			A.setRow(i,0,k.x,k.y,1);
			b.set(i,0,c.x);
		}

		SimpleMatrix x = A.solve(b);

		SimpleMatrix Ha = SimpleMatrix.identity(3);
		Ha.setRow(0,0,x.getDDRM().data);

		return Ha;
	}
738950
static void computePseudo(DMatrixRMaj A , DMatrixRMaj output ) {
		final int N = A.numRows;

		DMatrix2x2 m = new DMatrix2x2();
		DMatrix2x2 m_inv = new DMatrix2x2();

		for (int i = 0, index=0; i < N; i++) {
			double a_i1 = A.data[index++];
			double a_i2 = A.data[index++];
			m.a11 += a_i1*a_i1;
			m.a12 += a_i1*a_i2;
			m.a22 += a_i2*a_i2;
		}
		m.a21 = m.a12;
		CommonOps_DDF2.invert(m,m_inv);

		output.reshape(2,N);
		for (int i = 0,index=0; i < N; i++) {
			output.data[i] = A.data[index++]*m_inv.a11 + A.data[index++]*m_inv.a12;
		}
		int end = 2*N;
		for (int i = N,A_index=0; i < end; i++) {
			output.data[i] = A.data[A_index++]*m_inv.a21 + A.data[A_index++]*m_inv.a22;
		}
	}
1365855
def convex_conj(self):
        """"""
        if isinstance(self.domain, ProductSpace):
            norm = GroupL1Norm(self.domain, 2)
        else:
            norm = L1Norm(self.domain)

        return FunctionalQuadraticPerturb(norm.convex_conj,
                                          quadratic_coeff=self.gamma / 2)