Evaluate the Hamiltonian bny operator averaging

    :param theta:
    :param hamiltonian:
    :return:
548313
def objective_fun(theta, hamiltonian=None,
                  quantum_resource=QVMConnection(sync_endpoint='http://localhost:5000')):
    
    if hamiltonian is None:
        # Hamiltonian is Identity
        return 1.0

    if isinstance(hamiltonian, PauliSum):
        result = estimate_locally_commuting_operator(ucc_circuit(theta), hamiltonian,
                                                     1.0E-6, quantum_resource=quantum_resource)
        result = result[0][0].real  # first output is expected value, second is variance, third is shots
    elif isinstance(hamiltonian, np.ndarray) and isinstance(quantum_resource, QVMConnection):
        wf = quantum_resource.wavefunction(ucc_circuit(theta))
        wf = wf.amplitudes.reshape((-1, 1))
        result = np.conj(wf).T.dot(hamiltonian).dot(wf)[0, 0].real
        print(result)
    else:
        raise TypeError("type of hamiltonian or qvm is unrecognized")

    return result
1295444
def lnprior(self,theta):
        
        params,priors = self.params,self.priors
        kwargs = dict(list(zip(params,theta)))
        err = np.seterr(invalid='raise')
        try:
            lnprior = np.sum(np.log([priors[k](v) for k,v in list(kwargs.items())]))
        except (FloatingPointError,ValueError):
            lnprior = -np.inf
        np.seterr(**err)
        return lnprior
1008769
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]);
			}
		}
795457
def rdirichlet(theta, size=1):
    
    gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)])
    if size > 1 and np.size(theta) > 1:
        return (gammas.T / gammas.sum(1))[:-1].T
    elif np.size(theta) > 1:
        return (gammas[0] / gammas[0].sum())[:-1]
    else:
        return 1.
526865
def theta_limit(theta):
    
    while theta <= -1 * np.pi / 2:
        theta += np.pi
    while theta > np.pi / 2:
        theta -= np.pi
    return theta
909770
def orthonormality(V, ip_B=None):
    
    return norm(numpy.eye(V.shape[1]) - inner(V, V, ip_B=ip_B))
543612
def N_to_Ntriangles(N):
    

    theta = np.array([np.pi/2*(k-0.5)/N for k in range(1, N+1)])
    phi = np.array([[np.pi*(l-0.5)/Mk for l in range(1, Mk+1)] for Mk in np.array(1 + 1.3*N*np.sin(theta), dtype=int)])
    Ntri = 2*np.array([len(p) for p in phi]).sum()

    return Ntri
676261
def inter_spin_hamiltonian(self, u_int, J_coup):
        """"""
        J_coup *= u_int
        h_int  = (u_int - 2*J_coup)/2.*self.oper['sumSz2']
        h_int += J_coup*self.oper['sumSz-sp2']
        h_int -= J_coup/2.*self.oper['sumSz-or2']
        h_int -= J_coup*self.oper['Sfliphop']

        return h_int
897271
def get_averaged_interaction(self):
        

        # v[triplet, band0, band, band]
        v = self._interaction_strength
        w = self._weights_at_q
        v_sum = np.dot(w, v.sum(axis=2).sum(axis=2))
        return v_sum / np.prod(v.shape[2:])
210774
def _run_theta(cnv_info, data, work_dir, run_n3=True):
    
    out = {"caller": "theta"}
    max_normal = "0.9"
    opts = ["-m", max_normal]
    n2_result = _safe_run_theta(cnv_info["theta_input"], os.path.join(work_dir, "n2"), ".n2.results",
                                ["-n", "2"] + opts, data)
    if n2_result:
        out["estimate"] = n2_result
        if run_n3:
            n2_bounds = "%s.withBounds" % os.path.splitext(n2_result)[0]
            n3_result = _safe_run_theta(n2_bounds, os.path.join(work_dir, "n3"), ".n3.results",
                                        ["-n", "3", "--RESULTS", n2_result] + opts,
                                        data)
            if n3_result:
                best_result = _select_model(n2_bounds, n2_result, n3_result,
                                            os.path.join(work_dir, "n3"), data)
                out["estimate"] = best_result
                out["cnvs"] = _merge_theta_calls(n2_bounds, best_result, cnv_info["vrn_file"], data)
    return out
600535
def Ry_matrix(theta):
    """"""
    return np.array([
        [np.cos(theta), 0, np.sin(theta)],
        [0, 1, 0],
        [-np.sin(theta), 0, np.cos(theta)]
    ])
735329
private void hqr2() {
    // FIXME: does this fail on NaN/inf values?

    // This is derived from the Algol procedure hqr2,
    // by Martin and Wilkinson, Handbook for Auto. Comp.,
    // Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.

    // Initialize
    final int nn = this.n, low = 0, high = nn - 1;

    // Store roots isolated by balanc and compute matrix norm
    double norm = 0.;
    for(int i = 0; i < nn; i++) {
      // If eventually allowing low,high to be set, use this:
      // if(i < low || i > high) { d[i] = H[i][i]; e[i] = 0.; }
      for(int j = (i > 0 ? i - 1 : 0); j < nn; j++) {
        norm += abs(H[i][j]);
      }
    }

    // Outer loop over eigenvalue index
    {
      double exshift = 0.0;
      // double p = 0, q = 0, r = 0; // , s = 0, z = 0;
      int iter = 0;
      for(int n = nn - 1; n >= low;) {
        // Look for single small sub-diagonal element
        int l = n;
        for(; l > low; --l) {
          double s = abs(H[l - 1][l - 1]) + abs(H[l][l]);
          if(abs(H[l][l - 1]) < 0x1P-52 * (s == 0. ? norm : s)) {
            break;
          }
        }

        // Check for convergence
        if(l == n) {
          // One root found
          d[n] = (H[n][n] += exshift);
          e[n] = 0.;
          n--;
          iter = 0;
          continue;
        }
        final double[] Hn = H[n], Hnm1 = H[n - 1];
        if(l == n - 1) {
          // Two roots found
          double w = Hn[n - 1] * Hnm1[n];
          double p = (Hnm1[n - 1] - Hn[n]) * 0.5, q = p * p + w;
          double x = (Hn[n] += exshift), z = FastMath.sqrt(abs(q));
          Hnm1[n - 1] += exshift;

          if(q >= 0) {
            // Real pair
            z = (p >= 0) ? p + z : p - z;
            d[n - 1] = x + z;
            d[n] = (z != 0.) ? x - w / z : d[n - 1];
            e[n] = e[n - 1] = 0.;
            x = Hn[n - 1];
            double s = abs(x) + abs(z);
            p = x / s;
            q = z / s;
            double r = FastMath.hypot(p, q);
            p /= r;
            q /= r;

            // Row modification
            for(int j = n - 1; j < nn; j++) {
              double tmp = Hnm1[j];
              Hnm1[j] = q * tmp + p * Hn[j];
              Hn[j] = q * Hn[j] - p * tmp;
            }

            // Column modification
            for(int i = 0; i <= n; i++) {
              modifyQP(H[i], n, p, q);
            }

            // Accumulate transformations
            for(int i = low; i <= high; i++) {
              modifyQP(V[i], n, p, q);
            }
          }
          else {
            // Complex pair
            d[n] = d[n - 1] = x + p;
            e[n] = -(e[n - 1] = z);
          }
          n -= 2;
          iter = 0;
          continue;
        }
        // No convergence yet

        // Form shift
        double x = Hn[n], y = 0., w = 0.;
        y = Hnm1[n - 1];
        w = Hn[n - 1] * Hnm1[n];

        // Wilkinson's original ad hoc shift
        if(iter == 10) {
          exshift += x;
          for(int i = low; i <= n; i++) {
            H[i][i] -= x;
          }
          double s = abs(Hn[n - 1]) + abs(Hnm1[n - 2]);
          x = y = 0.75 * s;
          w = -0.4375 * s * s;
        }

        // MATLAB's new ad hoc shift
        if(iter == 30) {
          double s = (y - x) * 0.5;
          s = s * s + w;
          if(s > 0) {
            s = FastMath.sqrt(s);
            s = (y < x) ? -s : s;
            s = x - w / ((y - x) * 0.5 + s);
            for(int i = low; i <= n; i++) {
              H[i][i] -= s;
            }
            exshift += s;
            x = y = w = 0.964;
          }
        }

        ++iter; // (Could check iteration count here.)

        double p = 0., q = 0., r = 0.;
        // Look for two consecutive small sub-diagonal elements
        int m = n - 2;
        for(; m >= l; m--) {
          final double[] Hm = H[m], Hmp1 = H[m + 1];
          final double z = Hm[m], xz = x - z;
          p = (xz * (y - z) - w) / Hmp1[m] + Hm[m + 1];
          q = Hmp1[m + 1] - xz - y;
          r = H[m + 2][m + 1];
          double tmp = abs(p) + abs(q) + abs(r);
          p /= tmp;
          q /= tmp;
          r /= tmp;
          if(m == l || abs(Hm[m - 1]) * (abs(q) + abs(r)) < 0x1P-52 * (abs(p) * (abs(H[m - 1][m - 1]) + abs(z) + abs(Hmp1[m + 1])))) {
            break;
          }
        }

        for(int i = m + 2; i <= n; i++) {
          final double[] Hi = H[i];
          Hi[i - 2] = 0.;
          if(i > m + 2) {
            Hi[i - 3] = 0.;
          }
        }

        // Double QR step involving rows l:n and columns m:n
        for(int k = m, last = n - 1; k <= last; k++) {
          boolean notlast = (k != last);
          final double[] Hk = H[k], Hkp1 = H[k + 1],
              Hkp2 = notlast ? H[k + 2] : null;
          if(k != m) {
            p = Hk[k - 1];
            q = Hkp1[k - 1];
            r = notlast ? Hkp2[k - 1] : 0.;
            x = abs(p) + abs(q) + abs(r);
            if(x == 0.0) {
              continue; // Jama 1.0.3 fix
            }
            p /= x;
            q /= x;
            r /= x;
          }
          double s = FastMath.hypot(p, q, r);
          s = (p < 0) ? -s : s;
          if(s != 0) {
            if(k != m) {
              Hk[k - 1] = -s * x;
            }
            else if(l != m) {
              Hk[k - 1] = -Hk[k - 1];
            }
            p += s;
            x = p / s;
            y = q / s;
            double z = r / s;
            q /= p;
            r /= p;

            // Row modification
            for(int j = k; j < nn; j++) {
              double tmp = Hk[j] + q * Hkp1[j];
              if(notlast) {
                tmp += r * Hkp2[j];
                Hkp2[j] -= tmp * z;
              }
              Hk[j] -= tmp * x;
              Hkp1[j] -= tmp * y;
            }

            // Column modification
            for(int i = 0; i <= Math.min(n, k + 3); i++) {
              modifyQR(H[i], k, notlast, q, r, x, y, z);
            }

            // Accumulate transformations
            for(int i = low; i <= high; i++) {
              modifyQR(V[i], k, notlast, q, r, x, y, z);
            }
          } // (s != 0)
        } // k loop
        // check convergence
      } // while (n >= low)
    }

    if(norm == 0.0) {
      return;
    }

    // Backsubstitute to find vectors of upper triangular form
    for(int n = nn - 1; n >= 0; n--) {
      double p = d[n], q = e[n];
      if(q == 0) {
        hqr2BacksubstituteReal(n, p, norm);
      }
      else if(q < 0) {
        hqr2BacksubstituteComplex(n, p, q, norm);
      }
    }

    // Vectors of isolated roots
    // Only matters if the user can modify low, high:
    // for(int i = 0; i < nn; i++) {
    // if(i < low || i > high) { System.arraycopy(H[i], i, V[i], i, nn - i); }
    // }

    // Back transformation to get eigenvectors of original matrix
    hqr2BackTransformation(nn, low, high);
  }
1373221
def generate_random_ising_problem(solver, h_range=None, j_range=None):
    

    if h_range is None:
        h_range = solver.properties.get('h_range', [-1, 1])
    if j_range is None:
        j_range = solver.properties.get('j_range', [-1, 1])

    lin = {qubit: random.uniform(*h_range) for qubit in solver.nodes}
    quad = {edge: random.uniform(*j_range) for edge in solver.undirected_edges}

    return lin, quad
1003681
def v_theta(self):
        
        # compute theta
        norm = np.sqrt(np.sum(self.wxyz**2))
        theta = 2 * np.arccos(self.wxyz[0] / norm)

        # compute the unit vector
        v = np.array(self.wxyz[1:])
        v = v / np.sqrt(np.sum(v**2))

        return v, theta
641132
def hessian(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y
        dphi_dr = self._dphi_dr(x_shift, y_shift, theta_E, r_trunc)
        d2phi_dr2 = self._d2phi_dr2(x_shift, y_shift, theta_E, r_trunc)
        dr_dx, dr_dy = self._dr_dx(x, y)
        d2r_dx2, d2r_dy2, d2r_dxy = self._d2r_dx2(x_shift, y_shift)
        f_xx = d2r_dx2*dphi_dr + dr_dx**2*d2phi_dr2
        f_yy = d2r_dy2*dphi_dr + dr_dy**2*d2phi_dr2
        f_xy = d2r_dxy*dphi_dr + dr_dx*dr_dy*d2phi_dr2
        return f_xx, f_yy, f_xy
301605
function(theta) {
        _TEMP.identity();
        _TEMP.matrix[1][1] = _TEMP.matrix[2][2] = Math.cos(theta);
        _TEMP.matrix[2][1] = Math.sin(theta);
        _TEMP.matrix[1][2] = -_TEMP.matrix[2][1];
        return this.multiplySelf(_TEMP);
    }
546509
def make_pendulum(theta_max, l):
    
    g_l = 9.8 / l
    def deriv(t, y):
        """ Calculate [d/d theta(t), d/dt d/dt theta(t)]. """
        theta, dtheta_dt = y
        return np.array([dtheta_dt, - g_l * gv.sin(theta)])
    y0 = np.array([theta_max, 0.0])
    return gv.ode.Integrator(deriv=deriv).solution(0.0, y0)