Adjoint of the linear operator.

        Note that this implementation uses an approximation that is only
        valid for small displacements.
497263
def adjoint(self):
        
        # TODO allow users to select what method to use here.
        div_op = Divergence(domain=self.displacement.space, method='forward',
                            pad_mode='symmetric')
        jacobian_det = self.domain.element(
            np.exp(-div_op(self.displacement)))

        return jacobian_det * self.inverse
149480
def interpolation_linear(x, x1, x2, y1, y2):
    
    m = (y2 - y1) / (x2 - x1)
    t = (x - x1)
    return m * t + y1
1333974
def matmul_adjoint_x(dz, x, y, transpose_a, transpose_b):
  """"""
  if not transpose_a and not transpose_b:
    return tf.matmul(dz, y, transpose_b=True)
  elif not transpose_a and transpose_b:
    return tf.matmul(dz, y)
  elif transpose_a and not transpose_b:
    return tf.matmul(y, dz, transpose_b=True)
  else:  # transpose_a and transpose_b
    return tf.matmul(y, dz, transpose_a=True, transpose_b=True)
393543
def inverse(self):
        """"""
        sign = '-' if self.sign == '+' else '+'
        return DiscreteFourierTransform(
            domain=self.range, range=self.domain, axes=self.axes,
            halfcomplex=self.halfcomplex, sign=sign)
582764
private static Operator substitute(Operator oper, List[] equatedIds) {
    Selector op1 = oper.getOperands()[0];
    Selector op2 = (oper.getOperands().length == 1) ? null : oper.getOperands()[1];
    if (op1 instanceof Identifier)
      op1 = substitute((Identifier) op1, equatedIds);
    else if (op1 instanceof Operator)
      op1 = substitute((Operator) op1, equatedIds);
    if (op1 == null)
      return null;
    if (op2 != null) {
      if (op2 instanceof Identifier)
        op2 = substitute((Identifier) op2, equatedIds);
      else if (op2 instanceof Operator)
        op2 = substitute((Operator) op2, equatedIds);
      if (op2 == null)
        return null;
    }
    if (op1 == oper.getOperands()[0] && (op2 == null || op2 == oper.getOperands()[1]))
      return oper;
    else if (oper instanceof LikeOperator) {
      LikeOperatorImpl loper = (LikeOperatorImpl) oper;
        return new LikeOperatorImpl(loper.getOp(), op1, loper.getInternalPattern(), loper.getPattern(), loper.isEscaped(), loper.getEscape());
    }
    else
      return (op2 == null) ? new OperatorImpl(oper.getOp(), op1) : new OperatorImpl(oper.getOp(), op1, op2);
  }
1172222
def adjoint(self):
        
        return PointwiseInner(vfspace=self.range, vecfield=self.vecfield,
                              weighting=self.weights)
169894
def relaxation_operators(p):
    
    k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - p)]])
    k1 = np.array([[0.0, np.sqrt(p)], [0.0, 0.0]])
    return k0, k1
1179028
def expression(self):
        
        return operatorPrecedence(self.relation, [
            (self.unary_op, 1, opAssoc.RIGHT),
            (self.binary_op_p1, 2, opAssoc.LEFT),
            (self.binary_op_p2, 2, opAssoc.LEFT)])
1007463
def matmul_adjoint_y(dz, x, y, transpose_a, transpose_b):
  """"""
  if not transpose_a and not transpose_b:
    return tf.matmul(x, dz, transpose_a=True)
  elif not transpose_a and transpose_b:
    return tf.matmul(dz, x, transpose_a=True)
  elif transpose_a and not transpose_b:
    return tf.matmul(x, dz)
  else:  # transpose_a and transpose_b
    return tf.matmul(dz, x, transpose_a=True, transpose_b=True)
1092191
private function hqr2()
    {
        //  Initialize
        $nn = $this->n;
        $n  = $nn - 1;
        $low = 0;
        $high = $nn - 1;
        $eps = pow(2.0, -52.0);
        $exshift = 0.0;
        $p = $q = $r = $s = $z = 0;
        // Store roots isolated by balanc and compute matrix norm
        $norm = 0.0;

        for ($i = 0; $i < $nn; ++$i) {
            if (($i < $low) or ($i > $high)) {
                $this->d[$i] = $this->H[$i][$i];
                $this->e[$i] = 0.0;
            }
            for ($j = max($i-1, 0); $j < $nn; ++$j) {
                $norm = $norm + abs($this->H[$i][$j]);
            }
        }

        // Outer loop over eigenvalue index
        $iter = 0;
        while ($n >= $low) {
            // Look for single small sub-diagonal element
            $l = $n;
            while ($l > $low) {
                $s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
                if ($s == 0.0) {
                    $s = $norm;
                }
                if (abs($this->H[$l][$l-1]) < $eps * $s) {
                    break;
                }
                --$l;
            }
            // Check for convergence
            // One root found
            if ($l == $n) {
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->d[$n] = $this->H[$n][$n];
                $this->e[$n] = 0.0;
                --$n;
                $iter = 0;
            // Two roots found
            } elseif ($l == $n-1) {
                $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
                $p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
                $q = $p * $p + $w;
                $z = sqrt(abs($q));
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
                $x = $this->H[$n][$n];
                // Real pair
                if ($q >= 0) {
                    if ($p >= 0) {
                        $z = $p + $z;
                    } else {
                        $z = $p - $z;
                    }
                    $this->d[$n-1] = $x + $z;
                    $this->d[$n] = $this->d[$n-1];
                    if ($z != 0.0) {
                        $this->d[$n] = $x - $w / $z;
                    }
                    $this->e[$n-1] = 0.0;
                    $this->e[$n] = 0.0;
                    $x = $this->H[$n][$n-1];
                    $s = abs($x) + abs($z);
                    $p = $x / $s;
                    $q = $z / $s;
                    $r = sqrt($p * $p + $q * $q);
                    $p = $p / $r;
                    $q = $q / $r;
                    // Row modification
                    for ($j = $n-1; $j < $nn; ++$j) {
                        $z = $this->H[$n-1][$j];
                        $this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
                        $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
                    }
                    // Column modification
                    for ($i = 0; $i <= $n; ++$i) {
                        $z = $this->H[$i][$n-1];
                        $this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
                        $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
                    }
                    // Accumulate transformations
                    for ($i = $low; $i <= $high; ++$i) {
                        $z = $this->V[$i][$n-1];
                        $this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
                        $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
                    }
                // Complex pair
                } else {
                    $this->d[$n-1] = $x + $p;
                    $this->d[$n] = $x + $p;
                    $this->e[$n-1] = $z;
                    $this->e[$n] = -$z;
                }
                $n = $n - 2;
                $iter = 0;
            // No convergence yet
            } else {
                // Form shift
                $x = $this->H[$n][$n];
                $y = 0.0;
                $w = 0.0;
                if ($l < $n) {
                    $y = $this->H[$n-1][$n-1];
                    $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
                }
                // Wilkinson's original ad hoc shift
                if ($iter == 10) {
                    $exshift += $x;
                    for ($i = $low; $i <= $n; ++$i) {
                        $this->H[$i][$i] -= $x;
                    }
                    $s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
                    $x = $y = 0.75 * $s;
                    $w = -0.4375 * $s * $s;
                }
                // MATLAB's new ad hoc shift
                if ($iter == 30) {
                    $s = ($y - $x) / 2.0;
                    $s = $s * $s + $w;
                    if ($s > 0) {
                        $s = sqrt($s);
                        if ($y < $x) {
                            $s = -$s;
                        }
                        $s = $x - $w / (($y - $x) / 2.0 + $s);
                        for ($i = $low; $i <= $n; ++$i) {
                            $this->H[$i][$i] -= $s;
                        }
                        $exshift += $s;
                        $x = $y = $w = 0.964;
                    }
                }
                // Could check iteration count here.
                $iter = $iter + 1;
                // Look for two consecutive small sub-diagonal elements
                $m = $n - 2;
                while ($m >= $l) {
                    $z = $this->H[$m][$m];
                    $r = $x - $z;
                    $s = $y - $z;
                    $p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
                    $q = $this->H[$m+1][$m+1] - $z - $r - $s;
                    $r = $this->H[$m+2][$m+1];
                    $s = abs($p) + abs($q) + abs($r);
                    $p = $p / $s;
                    $q = $q / $s;
                    $r = $r / $s;
                    if ($m == $l) {
                        break;
                    }
                    if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
                        $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
                        break;
                    }
                    --$m;
                }
                for ($i = $m + 2; $i <= $n; ++$i) {
                    $this->H[$i][$i-2] = 0.0;
                    if ($i > $m+2) {
                        $this->H[$i][$i-3] = 0.0;
                    }
                }
                // Double QR step involving rows l:n and columns m:n
                for ($k = $m; $k <= $n-1; ++$k) {
                    $notlast = ($k != $n-1);
                    if ($k != $m) {
                        $p = $this->H[$k][$k-1];
                        $q = $this->H[$k+1][$k-1];
                        $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
                        $x = abs($p) + abs($q) + abs($r);
                        if ($x != 0.0) {
                            $p = $p / $x;
                            $q = $q / $x;
                            $r = $r / $x;
                        }
                    }
                    if ($x == 0.0) {
                        break;
                    }
                    $s = sqrt($p * $p + $q * $q + $r * $r);
                    if ($p < 0) {
                        $s = -$s;
                    }
                    if ($s != 0) {
                        if ($k != $m) {
                            $this->H[$k][$k-1] = -$s * $x;
                        } elseif ($l != $m) {
                            $this->H[$k][$k-1] = -$this->H[$k][$k-1];
                        }
                        $p = $p + $s;
                        $x = $p / $s;
                        $y = $q / $s;
                        $z = $r / $s;
                        $q = $q / $p;
                        $r = $r / $p;
                        // Row modification
                        for ($j = $k; $j < $nn; ++$j) {
                            $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
                            if ($notlast) {
                                $p = $p + $r * $this->H[$k+2][$j];
                                $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
                            }
                            $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
                            $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
                        }
                        // Column modification
                        for ($i = 0; $i <= min($n, $k+3); ++$i) {
                            $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
                            if ($notlast) {
                                $p = $p + $z * $this->H[$i][$k+2];
                                $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
                            }
                            $this->H[$i][$k] = $this->H[$i][$k] - $p;
                            $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
                        }
                        // Accumulate transformations
                        for ($i = $low; $i <= $high; ++$i) {
                            $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
                            if ($notlast) {
                                $p = $p + $z * $this->V[$i][$k+2];
                                $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
                            }
                            $this->V[$i][$k] = $this->V[$i][$k] - $p;
                            $this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
                        }
                    }  // ($s != 0)
                }  // k loop
            }  // check convergence
        }  // while ($n >= $low)

        // Backsubstitute to find vectors of upper triangular form
        if ($norm == 0.0) {
            return;
        }

        for ($n = $nn-1; $n >= 0; --$n) {
            $p = $this->d[$n];
            $q = $this->e[$n];
            // Real vector
            if ($q == 0) {
                $l = $n;
                $this->H[$n][$n] = 1.0;
                for ($i = $n-1; $i >= 0; --$i) {
                    $w = $this->H[$i][$i] - $p;
                    $r = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $s = $r;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0.0) {
                            if ($w != 0.0) {
                                $this->H[$i][$n] = -$r / $w;
                            } else {
                                $this->H[$i][$n] = -$r / ($eps * $norm);
                            }
                        // Solve real equations
                        } else {
                            $x = $this->H[$i][$i+1];
                            $y = $this->H[$i+1][$i];
                            $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
                            $t = ($x * $s - $z * $r) / $q;
                            $this->H[$i][$n] = $t;
                            if (abs($x) > abs($z)) {
                                $this->H[$i+1][$n] = (-$r - $w * $t) / $x;
                            } else {
                                $this->H[$i+1][$n] = (-$s - $y * $t) / $z;
                            }
                        }
                        // Overflow control
                        $t = abs($this->H[$i][$n]);
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n] = $this->H[$j][$n] / $t;
                            }
                        }
                    }
                }
            // Complex vector
            } elseif ($q < 0) {
                $l = $n-1;
                // Last vector component imaginary so matrix is triangular
                if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
                    $this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
                    $this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
                } else {
                    $this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
                    $this->H[$n-1][$n-1] = $this->cdivr;
                    $this->H[$n-1][$n]   = $this->cdivi;
                }
                $this->H[$n][$n-1] = 0.0;
                $this->H[$n][$n] = 1.0;
                for ($i = $n-2; $i >= 0; --$i) {
                    // double ra,sa,vr,vi;
                    $ra = 0.0;
                    $sa = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
                        $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    $w = $this->H[$i][$i] - $p;
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $r = $ra;
                        $s = $sa;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0) {
                            $this->cdiv(-$ra, -$sa, $w, $q);
                            $this->H[$i][$n-1] = $this->cdivr;
                            $this->H[$i][$n]   = $this->cdivi;
                        } else {
                            // Solve complex equations
                            $x = $this->H[$i][$i+1];
                            $y = $this->H[$i+1][$i];
                            $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
                            $vi = ($this->d[$i] - $p) * 2.0 * $q;
                            if ($vr == 0.0 & $vi == 0.0) {
                                $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
                            }
                            $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
                            $this->H[$i][$n-1] = $this->cdivr;
                            $this->H[$i][$n]   = $this->cdivi;
                            if (abs($x) > (abs($z) + abs($q))) {
                                $this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
                                $this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
                            } else {
                                $this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
                                $this->H[$i+1][$n-1] = $this->cdivr;
                                $this->H[$i+1][$n]   = $this->cdivi;
                            }
                        }
                        // Overflow control
                        $t = max(abs($this->H[$i][$n-1]), abs($this->H[$i][$n]));
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
                                $this->H[$j][$n]   = $this->H[$j][$n] / $t;
                            }
                        }
                    } // end else
                } // end for
            } // end else for complex case
        } // end for

        // Vectors of isolated roots
        for ($i = 0; $i < $nn; ++$i) {
            if ($i < $low | $i > $high) {
                for ($j = $i; $j < $nn; ++$j) {
                    $this->V[$i][$j] = $this->H[$i][$j];
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix
        for ($j = $nn-1; $j >= $low; --$j) {
            for ($i = $low; $i <= $high; ++$i) {
                $z = 0.0;
                for ($k = $low; $k <= min($j, $high); ++$k) {
                    $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
                }
                $this->V[$i][$j] = $z;
            }
        }
    }
127032
def operator_is(u):
    """"""
    global _aux
    if np.ndim(u) == 2:
        P = _P2
    elif np.ndim(u) == 3:
        P = _P3
    else:
        raise ValueError("u has an invalid number of dimensions "
                         "(should be 2 or 3)")
    
    if u.shape != _aux.shape[1:]:
        _aux = np.zeros((len(P),) + u.shape)
    
    for _aux_i, P_i in zip(_aux, P):
        _aux_i[:] = binary_dilation(u, P_i)
    
    return _aux.min(0)
57043
private void orthes() {
    // FIXME: does this fail on NaN/inf values?

    // This is derived from the Algol procedures orthes and ortran,
    // by Martin and Wilkinson, Handbook for Auto. Comp.,
    // Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutines in EISPACK.

    final int low = 0, high = n - 1;
    for(int m = low + 1; m <= high - 1; m++) {
      // Scale column.
      double scale = 0.0;
      for(int i = m; i <= high; i++) {
        scale += abs(H[i][m - 1]);
      }
      if(scale > 0.0 || scale < 0.0) {
        // Compute Householder transformation.
        double h = 0.0;
        for(int i = high; i >= m; i--) {
          double oi = ort[i] = H[i][m - 1] / scale;
          h += oi * oi;
        }
        double g = FastMath.sqrt(h), om = ort[m];
        g = (om > 0) ? -g : g;
        h -= om * g;
        ort[m] = om - g;

        // Apply Householder similarity transformation
        // H = (I-u*u'/h)*H*(I-u*u')/h)

        for(int j = m; j < n; j++) {
          double f = 0.0;
          for(int i = high; i >= m; i--) {
            f += ort[i] * H[i][j];
          }
          f /= h;
          for(int i = m; i <= high; i++) {
            H[i][j] -= f * ort[i];
          }
        }

        for(int i = 0; i <= high; i++) {
          final double[] Hi = H[i];
          double f = 0.0;
          for(int j = high; j >= m; j--) {
            f += ort[j] * Hi[j];
          }
          f /= h;
          for(int j = m; j <= high; j++) {
            Hi[j] -= f * ort[j];
          }
        }
        ort[m] *= scale;
        H[m][m - 1] = scale * g;
      }
    }

    // Accumulate transformations (Algol's ortran).
    for(int i = 0; i < n; i++) {
      Arrays.fill(V[i], 0.);
      V[i][i] = 1.;
    }

    for(int m = high - 1; m >= low + 1; m--) {
      final double[] H_m = H[m];
      if(H_m[m - 1] != 0.0) {
        for(int i = m + 1; i <= high; i++) {
          ort[i] = H[i][m - 1];
        }
        final double ort_m = ort[m];
        for(int j = m; j <= high; j++) {
          double g = 0.0;
          for(int i = m; i <= high; i++) {
            g += ort[i] * V[i][j];
          }
          // Double division avoids possible underflow
          g = (g / ort_m) / H_m[m - 1];
          for(int i = m; i <= high; i++) {
            V[i][j] += g * ort[i];
          }
        }
      }
    }
  }
779400
function adjoint(out, a) {
  // Caching this value is nessecary if out == a
  var a0 =  a[0]
  out[0] =  a[3]
  out[1] = -a[1]
  out[2] = -a[2]
  out[3] =  a0

  return out
}
1018201
def inverse(self):
        """"""
        sign = '+' if self.sign == '-' else '-'
        return DiscreteFourierTransformInverse(
            domain=self.range, range=self.domain, axes=self.axes,
            halfcomplex=self.halfcomplex, sign=sign)
1167745
def is_hermitian(operator):
    
    if isinstance(operator, qt.Qobj):
        return (operator.dag() - operator).norm(FROBENIUS) / operator.norm(FROBENIUS) < EPS
    if isinstance(operator, np.ndarray):
        return np.linalg.norm(operator.T.conj() - operator) / np.linalg.norm(operator) < EPS
    return spnorm(operator.H - operator) / spnorm(operator) < EPS
844282
def operate(self, point):
        
        affine_point = np.array([point[0], point[1], point[2], 1])
        return np.dot(self.affine_matrix, affine_point)[0:3]
86176
public Vector4d hermite(Vector4dc t0, Vector4dc v1, Vector4dc t1, double t, Vector4d dest) {
        double t2 = t * t;
        double t3 = t2 * t;
        dest.x = (x + x - v1.x() - v1.x() + t1.x() + t0.x()) * t3 + (3.0 * v1.x() - 3.0 * x - t0.x() - t0.x() - t1.x()) * t2 + x * t + x;
        dest.y = (y + y - v1.y() - v1.y() + t1.y() + t0.y()) * t3 + (3.0 * v1.y() - 3.0 * y - t0.y() - t0.y() - t1.y()) * t2 + y * t + y;
        dest.z = (z + z - v1.z() - v1.z() + t1.z() + t0.z()) * t3 + (3.0 * v1.z() - 3.0 * z - t0.z() - t0.z() - t1.z()) * t2 + z * t + z;
        dest.w = (w + w - v1.w() - v1.w() + t1.w() + t0.w()) * t3 + (3.0 * v1.w() - 3.0 * w - t0.w() - t0.w() - t1.w()) * t2 + w * t + w;
        return dest;
    }