A special bidiagonal _matrix in the QTT-format
    M = IPAS(D, A)
    Generates I+a*S_{-1} _matrix in the QTT-format:
    1 0 0 0
    a 1 0 0
    0 a 1 0
    0 0 a 1
    Convenient for Crank-Nicolson and time gradient matrices
501535
def IpaS(d, a, tt_instance=True):
    '''
    '''

    if d == 1:
        M = _np.array([[1, 0], [a, 1]]).reshape((1, 2, 2, 1), order='F')
    else:
        M = [None] * d
        M[0] = _np.zeros((1, 2, 2, 2))
        M[0][0, :, :, 0] = _np.array([[1, 0], [a, 1]])
        M[0][0, :, :, 1] = _np.array([[0, a], [0, 0]])
        for i in xrange(1, d - 1):
            M[i] = _np.zeros((2, 2, 2, 2))
            M[i][:, :, 0, 0] = _np.eye(2)
            M[i][:, :, 1, 0] = _np.array([[0, 0], [1, 0]])
            M[i][:, :, 1, 1] = _np.array([[0, 1], [0, 0]])
        M[d - 1] = _np.zeros((2, 2, 2, 1))
        M[d - 1][:, :, 0, 0] = _np.eye(2)
        M[d - 1][:, :, 1, 0] = _np.array([[0, 0], [1, 0]])
    if tt_instance:
        M = _matrix.matrix.from_list(M)
    return M
452328
def funTransQuadD(k, s):
    
    sqrtk = np.sqrt(complex(k))
    a = np.cosh(sqrtk * s)
    b = np.sinh(sqrtk * s) / sqrtk
    c = sqrtk * np.sinh(sqrtk * s)
    d = np.cosh(sqrtk * s)
    return np.matrix([[a.real, b.real], [c.real, d.real]], dtype=np.double)
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;
    }
  }
1228098
def matrix2sheet(self,float_row,float_col):
        
        xoffset = float_col*self.__xstep
        if isinstance(self.lbrt[0], datetime_types):
            xoffset = np.timedelta64(int(round(xoffset)), self._time_unit)
        x = self.lbrt[0] + xoffset
        yoffset = float_row*self.__ystep
        if isinstance(self.lbrt[3], datetime_types):
            yoffset = np.timedelta64(int(round(yoffset)), self._time_unit)
        y = self.lbrt[3] - yoffset
        return x, y
916924
def build_diagonals(self):
    

    ##########################################################
    # INCORPORATE BOUNDARY CONDITIONS INTO COEFFICIENT ARRAY #
    ##########################################################

    # Roll to keep the proper coefficients at the proper places in the
    # arrays: Python will naturally just do vertical shifts instead of 
    # diagonal shifts, so this takes into account the horizontal compoent 
    # to ensure that boundary values are at the right place.
    self.l2 = np.roll(self.l2, -2)
    self.l1 = np.roll(self.l1, -1)
    self.r1 = np.roll(self.r1, 1)
    self.r2 = np.roll(self.r2, 2)

    # Then assemble these rows: this is where the periodic boundary condition 
    # can matter.
    if self.coeff_matrix is not None:
      pass
    elif self.BC_E == 'Periodic' and self.BC_W == 'Periodic':
      # In this case, the boundary-condition-related stacking has already 
      # happened inside b.c.-handling function. This is because periodic
      # boundary conditions require extra diagonals to exist on the edges of 
      # the solution array
      pass
    else:
      self.diags = np.vstack((self.l2,self.l1,self.c0,self.r1,self.r2))
      self.offsets = np.array([-2,-1,0,1,2])

    # Everybody now (including periodic b.c. cases)
    self.coeff_matrix = spdiags(self.diags, self.offsets, self.nx, self.nx, format='csr')
129669
def symbolic_Rz_matrix(symbolic_theta):
    """"""
    return sympy.Matrix([
        [sympy.cos(symbolic_theta), -sympy.sin(symbolic_theta), 0],
        [sympy.sin(symbolic_theta), sympy.cos(symbolic_theta), 0],
        [0, 0, 1]
    ])
403996
def GaussSeidelStep(self, S, X, ATYU, rho, lcw, W2):
        """"""

        Xss = np.zeros_like(S, dtype=self.dtype)
        for ax in self.axes:
            Xss += sl.zpad(X[(slice(None),)*ax + (slice(0, -1),)],
                           (1, 0), ax)
            Xss += sl.zpad(X[(slice(None),)*ax + (slice(1, None),)],
                           (0, 1), ax)
        return (rho*(Xss + ATYU) + W2*S) / (W2 + rho*lcw)
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]
695968
def PCA(Y, components):
	

	N,D = Y.shape
	sv = linalg.svd(Y, full_matrices=0);
	[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
	v = s0.std(axis=0)
	s0 /= v;
	w0 *= v;
	return [s0, w0]

	if N>D:
		sv = linalg.svd(Y, full_matrices=0);
		[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
		v = s0.std(axis=0)
		s0 /= v;
		w0 *= v;
		return [s0, w0]
	else:
		K=np.cov(Y)
		sv = linalg.eigh(K)
		std_var = np.sqrt(sv[0])
		pc = sv[1]*std_var[np.newaxis(),0]
		#import ipdb
		#ipdb.set_trace()
		return [pc,std_var]
1100694
def cublasStpsv(handle, uplo, trans, diag, n, AP, x, incx):
    
    
    status = _libcublas.cublasStpsv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, int(AP), int(x), incx)
    cublasCheckStatus(status)
785859
private void tql2 (double[][] V) {

   //  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 < n; i++) {
         e[i-1] = e[i];
      }
      e[n-1] = 0.0;
   
      double f = 0.0;
      double tst1 = 0.0;
      double eps = Math.pow(2.0,-52.0);
      for (int l = 0; l < n; l++) {

         // Find small subdiagonal element
   
         tst1 = Math.max(tst1,Math.abs(d[l]) + Math.abs(e[l]));
         int m = l;
         while (m < n) {
            if (Math.abs(e[m]) <= eps*tst1) {
               break;
            }
            m++;
         }
   
         // If m == l, d[l] is an eigenvalue,
         // otherwise, iterate.
   
         if (m > l) {
            int iter = 0;
            do {
               iter = iter + 1;  // (Could check iteration count here.)
   
               // Compute implicit shift
   
               double g = d[l];
               double p = (d[l+1] - g) / (2.0 * e[l]);
               double r = Math.hypot(p,1.0);
               if (p < 0) {
                  r = -r;
               }
               d[l] = e[l] / (p + r);
               d[l+1] = e[l] * (p + r);
               double dl1 = d[l+1];
               double h = g - d[l];
               for (int i = l+2; i < n; i++) {
                  d[i] -= h;
               }
               f = f + h;
   
               // Implicit QL transformation.
   
               p = d[m];
               double c = 1.0;
               double c2 = c;
               double c3 = c;
               double el1 = e[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 * e[i];
                  h = c * p;
                  r = Math.hypot(p,e[i]);
                  e[i+1] = s * r;
                  s = e[i] / r;
                  c = p / r;
                  p = c * d[i] - s * g;
                  d[i+1] = h + s * (c * g + s * d[i]);
   
                  // Accumulate transformation.
   
                  for (int k = 0; k < n; k++) {
                     h = V[k][i+1];
                     V[k][i+1] = s * V[k][i] + c * h;
                     V[k][i] = c * V[k][i] - s * h;
                  }
               }
               p = -s * s2 * c3 * el1 * e[l] / dl1;
               e[l] = s * p;
               d[l] = c * p;
   
               // Check for convergence.
   
            } while (Math.abs(e[l]) > eps*tst1);
         }
         d[l] = d[l] + f;
         e[l] = 0.0;
      }
     
      // Sort eigenvalues and corresponding vectors.
   
      for (int i = 0; i < n-1; i++) {
         int k = i;
         double p = d[i];
         for (int j = i+1; j < n; j++) {
            if (d[j] < p) {
               k = j;
               p = d[j];
            }
         }
         if (k != i) {
            d[k] = d[i];
            d[i] = p;
            for (int j = 0; j < n; j++) {
               p = V[j][i];
               V[j][i] = V[j][k];
               V[j][k] = p;
            }
         }
      }
   }
1028992
def numbaGaussian2d(psf, sy, sx):
    '''
    
    '''
    ps0, ps1 = psf.shape
    c0,c1 = ps0//2, ps1//2
    ssx = 2*sx**2
    ssy = 2*sy**2
    for i in range(ps0):
        for j in range(ps1):
            psf[i,j]=exp( -( (i-c0)**2/ssy
                            +(j-c1)**2/ssx) )
    psf/=psf.sum()
693611
private static void tred(DenseMatrix V, double[] d, double[] e) {
        int n = V.nrows();
        for (int i = 0; i < n; i++) {
            d[i] = V.get(n - 1, i);
        }

        // Householder reduction to tridiagonal form.
        for (int i = n - 1; i > 0; i--) {

            // Scale to avoid under/overflow.
            double scale = 0.0;
            double h = 0.0;
            for (int k = 0; k < i; k++) {
                scale = scale + Math.abs(d[k]);
            }
            if (scale == 0.0) {
                e[i] = d[i - 1];
                for (int j = 0; j < i; j++) {
                    d[j] = V.get(i - 1, j);
                    V.set(i, j, 0.0);
                    V.set(j, i, 0.0);
                }
            } else {

                // Generate Householder vector.
                for (int k = 0; k < i; k++) {
                    d[k] /= scale;
                    h += d[k] * d[k];
                }
                double f = d[i - 1];
                double g = Math.sqrt(h);
                if (f > 0) {
                    g = -g;
                }
                e[i] = scale * g;
                h = h - f * g;
                d[i - 1] = f - g;
                for (int j = 0; j < i; j++) {
                    e[j] = 0.0;
                }

                // Apply similarity transformation to remaining columns.
                for (int j = 0; j < i; j++) {
                    f = d[j];
                    V.set(j, i, f);
                    g = e[j] + V.get(j, j) * f;
                    for (int k = j + 1; k <= i - 1; k++) {
                        g += V.get(k, j) * d[k];
                        e[k] += V.get(k, j) * f;
                    }
                    e[j] = g;
                }
                f = 0.0;
                for (int j = 0; j < i; j++) {
                    e[j] /= h;
                    f += e[j] * d[j];
                }
                double hh = f / (h + h);
                for (int j = 0; j < i; j++) {
                    e[j] -= hh * d[j];
                }
                for (int j = 0; j < i; j++) {
                    f = d[j];
                    g = e[j];
                    for (int k = j; k <= i - 1; k++) {
                        V.sub(k, j, (f * e[k] + g * d[k]));
                    }
                    d[j] = V.get(i - 1, j);
                    V.set(i, j, 0.0);
                }
            }
            d[i] = h;
        }

        for (int j = 0; j < n; j++) {
            d[j] = V.get(j, j);
        }
        e[0] = 0.0;
    }
1327131
private void _symmetricTred2 ()
  {
    // This is derived from the Algol procedures tred2 by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
    }

    // Householder reduction to tridiagonal form.

    for (int i = m_nDim - 1; i > 0; i--)
    {
      // Scale to avoid under/overflow.
      double scale = 0.0;
      double h = 0.0;
      for (int k = 0; k < i; k++)
        scale += MathHelper.abs (m_aEVd[k]);
      if (scale == 0.0)
      {
        m_aEVe[i] = m_aEVd[i - 1];
        for (int j = 0; j < i; j++)
        {
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
          m_aEigenVector[j][i] = 0.0;
        }
      }
      else
      {
        // Generate Householder vector.

        for (int k = 0; k < i; k++)
        {
          m_aEVd[k] /= scale;
          h += m_aEVd[k] * m_aEVd[k];
        }
        double f = m_aEVd[i - 1];
        double g = Math.sqrt (h);
        if (f > 0)
        {
          g = -g;
        }
        m_aEVe[i] = scale * g;
        h -= f * g;
        m_aEVd[i - 1] = f - g;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] = 0.0;
        }

        // Apply similarity transformation to remaining columns.
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          m_aEigenVector[j][i] = f;
          g = m_aEVe[j] + m_aEigenVector[j][j] * f;
          for (int k = j + 1; k <= i - 1; k++)
          {
            g += m_aEigenVector[k][j] * m_aEVd[k];
            m_aEVe[k] += m_aEigenVector[k][j] * f;
          }
          m_aEVe[j] = g;
        }
        f = 0.0;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] /= h;
          f += m_aEVe[j] * m_aEVd[j];
        }
        final double hh = f / (h + h);
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] -= hh * m_aEVd[j];
        }
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          g = m_aEVe[j];
          for (int k = j; k <= i - 1; k++)
          {
            m_aEigenVector[k][j] -= (f * m_aEVe[k] + g * m_aEVd[k]);
          }
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
        }
      }
      m_aEVd[i] = h;
    }

    // Accumulate transformations.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      m_aEigenVector[m_nDim - 1][i] = m_aEigenVector[i][i];
      m_aEigenVector[i][i] = 1.0;
      final double h = m_aEVd[i + 1];
      if (h != 0.0)
      {
        for (int k = 0; k <= i; k++)
        {
          m_aEVd[k] = m_aEigenVector[k][i + 1] / h;
        }
        for (int j = 0; j <= i; j++)
        {
          double g = 0.0;
          for (int k = 0; k <= i; k++)
          {
            g += m_aEigenVector[k][i + 1] * m_aEigenVector[k][j];
          }
          for (int k = 0; k <= i; k++)
          {
            m_aEigenVector[k][j] -= g * m_aEVd[k];
          }
        }
      }
      for (int k = 0; k <= i; k++)
      {
        m_aEigenVector[k][i + 1] = 0.0;
      }
    }
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
      m_aEigenVector[m_nDim - 1][j] = 0.0;
    }
    m_aEigenVector[m_nDim - 1][m_nDim - 1] = 1.0;
    m_aEVe[0] = 0.0;
  }
473464
def canonical(self):
        
        Ac1 = np.hstack((self.deltah, np.zeros((self.nh, self.nz))))
        Ac2 = np.hstack((np.zeros((self.nz, self.nh)), self.a22))
        Ac = np.vstack((Ac1, Ac2))
        Bc = np.vstack((self.thetah, np.zeros((self.nz, self.nc))))
        Cc = np.vstack((np.zeros((self.nh, self.nw)), self.c2))
        Rc1 = np.hstack((self.llambda.T.dot(self.llambda), -
                         self.llambda.T.dot(self.ub)))
        Rc2 = np.hstack((-self.ub.T.dot(self.llambda), self.ub.T.dot(self.ub)))
        Rc = np.vstack((Rc1, Rc2))
        Qc = self.pih.T.dot(self.pih)
        Nc = np.hstack(
            (self.pih.T.dot(self.llambda), -self.pih.T.dot(self.ub)))

        lq_aux = LQ(Qc, Rc, Ac, Bc, N=Nc, beta=self.beta)

        P1, F1, d1 = lq_aux.stationary_values()

        self.F_b = F1[:, 0:self.nh]
        self.F_f = F1[:, self.nh:]

        self.pihat = np.linalg.cholesky(self.pih.T.dot(
            self.pih) + self.beta.dot(self.thetah.T).dot(P1[0:self.nh, 0:self.nh]).dot(self.thetah)).T
        self.llambdahat = self.pihat.dot(self.F_b)
        self.ubhat = - self.pihat.dot(self.F_f)

        return
876339
def numpartition_qaoa(n_step, nums, minimizer=None, sampler=None):
    
    hamiltonian = pauli.Expr.zero()
    for i, x in enumerate(nums):
        hamiltonian += pauli.Z[i] * x
    hamiltonian = (hamiltonian ** 2).simplify()

    return vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, n_step), minimizer, sampler)
1325771
public static double[][] computeWeightMatrix(final int quanth, final int quants, final int quantb) {
    final int dim = quanth * quants * quantb;
    final DoubleWrapper tmp = new DoubleWrapper(); // To return cosine
    assert (dim > 0);
    final double[][] m = new double[dim][dim];
    for(int x = 0; x < dim; x++) {
      final int hx = x / (quantb * quants);
      final int sx = (x / quantb) % quants;
      final int bx = x % quantb;
      for(int y = x; y < dim; y++) {
        final int hy = y / (quantb * quants);
        final int sy = (y / quantb) % quants;
        final int by = y % quantb;

        final double shx = FastMath.sinAndCos((hx + .5) / quanth * MathUtil.TWOPI, tmp);
        final double chx = tmp.value;
        final double shy = FastMath.sinAndCos((hy + .5) / quanth * MathUtil.TWOPI, tmp);
        final double chy = tmp.value;
        final double cos = chx * (sx + .5) / quants - chy * (sy + .5) / quants;
        final double sin = shx * (sx + .5) / quants - shy * (sy + .5) / quants;
        final double db = (bx - by) / (double) quantb;
        final double val = 1. - FastMath.sqrt((db * db + sin * sin + cos * cos) / 5);
        m[x][y] = m[y][x] = val;
      }
    }
    return m;
  }