Finds optimal angles with the quantum variational eigensolver method.

        Stored VQE result

        :returns: A tuple of the beta angles and the gamma angles for the optimal solution.
359614
def get_angles(self) -> Tuple[List[float], List[float]]:
        
        stacked_params = np.hstack((self.betas, self.gammas))
        vqe = VQE(self.minimizer, minimizer_args=self.minimizer_args,
                  minimizer_kwargs=self.minimizer_kwargs)
        cost_ham = reduce(lambda x, y: x + y, self.cost_ham)
        # maximizing the cost function!
        param_prog = self.get_parameterized_program()
        result = vqe.vqe_run(param_prog, cost_ham, stacked_params, qc=self.qc,
                             **self.vqe_options)
        self.result = result
        betas = result.x[:self.steps]
        gammas = result.x[self.steps:]
        return betas, gammas
961416
def vertical_modes(self):
        

        evals,evecs = np.linalg.eig(-self.S)

        asort = evals.argsort()

        # deformation wavenumbers and radii
        self.kdi2 = evals[asort]
        self.radii = np.zeros_like(self.kdi2)
        self.radii[0] = self.g*self.H/np.abs(self.f) # barotropic def. radius
        self.radii[1:] = 1./np.sqrt(self.kdi2[1:])

        # eigenstructure
        self.pmodes = evecs[:,asort]

        # normalize to have unit L2-norm
        Ai = (self.H / (self.Hi[:,np.newaxis]*(self.pmodes**2)).sum(axis=0))**0.5
        self.pmodes = Ai[np.newaxis,:]*self.pmodes
106408
def _get_v_angle_guess(self, case):
        
        v_angle = array([bus.v_angle * (pi / 180.0)
                         for bus in case.connected_buses])
        return v_angle
706640
def _optimal_orientation_from_detector(detector_name, tc):
    

    d = Detector(detector_name)
    ra, dec = d.optimal_orientation(tc)

    return ra, dec
499781
def optimise(args, **kwargs):
    

    expected_output_files = kwargs.pop("expected_output_files", None)
    if not expected_output_files:
        expected_output_files = ["optimised.pkl"]
        if args.plotting:
            expected_output_files.extend([
                "projection-estimate.{}".format(args.plot_format),
                "projection-optimised.{}".format(args.plot_format)
            ])
        
    # Estimate the model parameters, unless they are already specified.
    model = sick.models.Model(args.model)
    initial_theta = model._configuration.get("initial_theta", {})
    if len(set(model.parameters).difference(initial_theta)) == 0:
        model, data, metadata = _pre_solving(args, expected_output_files)

    else:
        model, data, metadata, initial_theta = estimate(args, 
            expected_output_files=expected_output_files, __return_result=True)

    try:
        theta, chisq, dof, model_fluxes = model.optimise(data, 
            initial_theta=initial_theta, full_output=True, debug=args.debug)

    except:
        logger.exception("Failed to optimise model parameters")
        raise

    metadata["optimised"] = {
        "theta": theta,
        "chi_sq": chisq,
        "dof": dof,
        "r_chi_sq": chisq/dof
    }

    logger.info("Optimised model parameters are:")
    _announce_theta(theta)
    logger.info("With a chi-sq value of {0:.1f} (reduced {1:.1f}; DOF {2:.1f})"\
        .format(chisq, chisq/dof, dof))

    if args.plotting:
        fig = sick.plot.spectrum(data, model_flux=model_fluxes)
        filename = _prefix(args, "projection-optimised.{}".format(
            args.plot_format))
        fig.savefig(filename)
        logger.info("Created figure {}".format(filename))

    if kwargs.pop("__return_result", False):
        return (model, data, metadata, theta)

    # Write the results to file.
    _write_output(_prefix(args, "optimised.pkl"), metadata)
    return None
505948
def get_optimization_gradients(self):
        """"""
        grad_array = self.fields.get("Opt point       1 Gradient at each geome")
        if grad_array is None:
            return []
        else:
            return np.reshape(grad_array, (-1, len(self.molecule.numbers), 3))
1066380
def solution_path(self, min_lambda, max_lambda, lambda_bins, verbose=0):
        ''''''
        lambda_grid = np.exp(np.linspace(np.log(max_lambda), np.log(min_lambda), lambda_bins))
        aic_trace = np.zeros(lambda_grid.shape) # The AIC score for each lambda value
        aicc_trace = np.zeros(lambda_grid.shape) # The AICc score for each lambda value (correcting for finite sample size)
        bic_trace = np.zeros(lambda_grid.shape) # The BIC score for each lambda value
        dof_trace = np.zeros(lambda_grid.shape) # The degrees of freedom of each final solution
        log_likelihood_trace = np.zeros(lambda_grid.shape)
        beta_trace = []
        best_idx = None
        best_plateaus = None

        # Solve the series of lambda values with warm starts at each point
        for i, lam in enumerate(lambda_grid):
            if verbose:
                print('#{0} Lambda = {1}'.format(i, lam))

            # Fit to the final values
            beta = self.solve(lam)

            if verbose:
                print('Calculating degrees of freedom')

            # Count the number of free parameters in the grid (dof)
            plateaus = calc_plateaus(beta, self.edges)
            dof_trace[i] = len(plateaus)

            if verbose:
                print('Calculating AIC')

            # Get the negative log-likelihood
            log_likelihood_trace[i] = self.log_likelihood(beta)

            # Calculate AIC = 2k - 2ln(L)
            aic_trace[i] = 2. * dof_trace[i] - 2. * log_likelihood_trace[i]
            
            # Calculate AICc = AIC + 2k * (k+1) / (n - k - 1)
            aicc_trace[i] = aic_trace[i] + 2 * dof_trace[i] * (dof_trace[i]+1) / (len(beta) - dof_trace[i] - 1.)

            # Calculate BIC = -2ln(L) + k * (ln(n) - ln(2pi))
            bic_trace[i] = -2 * log_likelihood_trace[i] + dof_trace[i] * (np.log(len(beta)) - np.log(2 * np.pi))

            # Track the best model thus far
            if best_idx is None or bic_trace[i] < bic_trace[best_idx]:
                best_idx = i
                best_plateaus = plateaus

            # Save the trace of all the resulting parameters
            beta_trace.append(np.array(beta))
            
            if verbose:
                print('DoF: {0} AIC: {1} AICc: {2} BIC: {3}'.format(dof_trace[i], aic_trace[i], aicc_trace[i], bic_trace[i]))

        if verbose:
            print('Best setting (by BIC): lambda={0} [DoF: {1}, AIC: {2}, AICc: {3} BIC: {4}]'.format(lambda_grid[best_idx], dof_trace[best_idx], aic_trace[best_idx], aicc_trace[best_idx], bic_trace[best_idx]))

        return {'aic': aic_trace,
                'aicc': aicc_trace,
                'bic': bic_trace,
                'dof': dof_trace,
                'loglikelihood': log_likelihood_trace,
                'beta': np.array(beta_trace),
                'lambda': lambda_grid,
                'best_idx': best_idx,
                'best': beta_trace[best_idx],
                'plateaus': best_plateaus}
470058
def get_matrix(theta, phi, psi, dx, dy, dz,
                   matrix=np.zeros((4, 4), dtype=DTYPE),
                   angles=np.zeros(3, dtype=DTYPE)):
        
        # NB!: matrix and angles by default are being overwritten on each call
        # thus, only created once at compile time.

        angles[0] = theta
        angles[1] = phi
        angles[2] = psi

        cx, cy, cz = np.cos(angles)
        sx, sy, sz = np.sin(angles)

        rotation = matrix[:3, :3]
        rotation.flat = (cx * cz - sx * cy * sz,
                         cx * sz + sx * cy * cz, sx * sy,
                         -sx * cz - cx * cy * sz,
                         -sx * sz + cx * cy * cz, cx * sy,
                         sy * sz,
                         -sy * cz, cy)

        # Translation component
        matrix[:3, 3] = dx, dy, dz
        matrix[3, 3] = 1.
        return matrix
717062
def axis_angle(self):
        
        qw, qx, qy, qz = self.quaternion
        theta = 2 * np.arccos(qw)
        omega = np.array([1,0,0])
        if theta > 0:
            rx = qx / np.sqrt(1.0 - qw**2)
            ry = qy / np.sqrt(1.0 - qw**2)
            rz = qz / np.sqrt(1.0 - qw**2)
            omega = np.array([rx, ry, rz])
        return theta * omega
23070
public DMatrixRMaj getQ() {
        Equation eq = new Equation();

        DMatrixRMaj Q = CommonOps_DDRM.identity(QR.numRows);
        DMatrixRMaj u = new DMatrixRMaj(QR.numRows,1);

        int N = Math.min(QR.numCols,QR.numRows);

        eq.alias(u,"u",Q,"Q",QR,"QR",QR.numRows,"r");

        // compute Q by first extracting the householder vectors from the columns of QR and then applying it to Q
        for( int j = N-1; j>= 0; j-- ) {
            eq.alias(j,"j",gammas[j],"gamma");

            eq.process("u(j:,0) = [1 ; QR((j+1):,j)]");
            eq.process("Q=(eye(r)-gamma*u*u')*Q");
        }

        return Q;
    }
708211
def channel_angle(chan0: Channel, chan1: Channel) -> bk.BKTensor:
    """"""
    return fubini_study_angle(chan0.vec, chan1.vec)
585271
def optimise(self, data, initial_theta=None, full_output=False, **kwargs):
        

        data = self._format_data(data)

        if initial_theta is None:
            initial_theta = self.estimate(data)

        # Which parameters will be optimised, and which will be fixed?
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)
        #parameters = set(self.parameters).difference(ignore_parameters)
        parameters = [p for p in self.parameters if p not in ignore_parameters]

        # What model wavelength ranges will be required?
        wavelengths_required = []
        for channel, spectrum in zip(matched_channels, data):
            if channel is None: continue
            z = initial_theta.get("z",
                initial_theta.get("z_{}".format(channel), 0))
            wavelengths_required.append(
                [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])

        # Create the spectrum approximator/interpolator.
        closest_point = [initial_theta[p] for p in self.grid_points.dtype.names]
        subset_bounds = self._initialise_approximator(
            closest_point=closest_point, 
            wavelengths_required=wavelengths_required, **kwargs)
        
        # Get the optimisation keyword arguments.
        op_kwargs = self._configuration.get("optimise", {}).copy()
        op_kwargs.update(kwargs)

        # Get fixed keywords.
        fixed = op_kwargs.pop("fixed", {})
        if fixed:
            # Remove non-parameters from the 'fixed' keywords.
            keys = set(fixed).intersection(parameters)
            # If the 'fixed' value is provided, use that. Otherwise if it is
            # None then use the initial_theta value.
            fixed = dict(zip(keys, 
                [(fixed[k], initial_theta.get(k, None))[fixed[k] is None] \
                    for k in keys]))

            logger.info("Fixing keyword arguments (these will not be optimised)"
                ": {}".format(fixed))

        # Remove fixed parameters from the parameters to be optimised
        #parameters = list(set(parameters).difference(fixed))
        parameters = [p for p in parameters if p not in fixed]

        # Translate input bounds.
        nbs = (None, None) # No boundaries.
        input_bounds = op_kwargs.pop("bounds", {})
        op_kwargs["bounds"] = [input_bounds.get(p, subset_bounds.get(p, nbs)) \
            for p in parameters]
        
        # Apply data masks now so we don't have to do it on the fly.
        masked_data, pixels_affected = self._apply_data_mask(data)

        # Prepare the convolution functions.
        self._create_convolution_functions(matched_channels, data, parameters,
            fixed_parameters=fixed)

        logger.info("Optimising parameters: {0}".format(", ".join(parameters)))
        logger.info("Optimisation keywords: {0}".format(op_kwargs))

        # Create the objective function.
        debug = kwargs.get("debug", False)
        def nlp(theta):
            # Apply fixed keywords
            t_ = theta.copy()
            p_ = [] + parameters
            for parameter, value in fixed.iteritems():
                p_.append(parameter)
                t_ = np.append(theta, value)

            return -inference.ln_probability(t_, p_, self, data, debug,
                matched_channels=matched_channels)

        # Do the optimisation.
        p0 = np.array([initial_theta[p] for p in parameters])

        x_opt = op.minimise(nlp, p0, **op_kwargs)

        # Put the result into a usable form.
        x_opt_theta = OrderedDict(zip(parameters, x_opt))
        x_opt_theta.update(fixed)

        if full_output:
            # Create model fluxes and calculate some metric.
            chi_sq, dof, model_fluxes = self._chi_sq(x_opt_theta, data)


            # Remove any prepared convolution functions.
            self._destroy_convolution_functions()
    
            return (x_opt_theta, chi_sq, dof, model_fluxes)

        # Remove any prepared convolution functions.
        self._destroy_convolution_functions()

        return x_opt_theta
1500
def optimizeAngle(angle):
    
    # First, we put the new angle in the range ]-360, 360[.
    # The modulo operator yields results with the sign of the
    # divisor, so for negative dividends, we preserve the sign
    # of the angle.
    if angle < 0:
        angle %= -360
    else:
        angle %= 360
    # 720 degrees is unnecessary, as 360 covers all angles.
    # As "-x" is shorter than "35x" and "-xxx" one character
    # longer than positive angles <= 260, we constrain angle
    # range to [-90, 270[ (or, equally valid: ]-100, 260]).
    if angle >= 270:
        angle -= 360
    elif angle < -90:
        angle += 360
    return angle
1087765
def solid_angles(self, permutation=None):
        
        if permutation is None:
            return self._solid_angles
        else:
            return [self._solid_angles[ii] for ii in permutation]
904651
def _determine_minimizer(self):
        
        if self.constraints:
            return SLSQP
        elif any([bound is not None for pair in self.model.bounds for bound in pair]):
            # If any bound is set
            return LBFGSB
        else:
            return BFGS
235257
def vecangle(v1, v2, deg=True):
    
    if np.array_equal(v1, v2):
        return 0.0
    dm = np.dot(v1, v2)
    cm = np.linalg.norm(v1) * np.linalg.norm(v2)
    angle = np.arccos(round(dm / cm, 10))  # Round here to prevent floating point errors
    return np.degrees([angle, ])[0] if deg else angle
1107735
def get_scipy_eig_vec(self):
    
    if not self.params['has_conv']:
      matrix_m = self.sess.run(self.dual_object.matrix_m)
      min_eig_vec_val, estimated_eigen_vector = eigs(matrix_m, k=1, which='SR',
                                                     tol=1E-4)
      min_eig_vec_val = np.reshape(np.real(min_eig_vec_val), [1, 1])
      return np.reshape(estimated_eigen_vector, [-1, 1]), min_eig_vec_val
    else:
      dim = self.dual_object.matrix_m_dimension
      input_vector = tf.placeholder(tf.float32, shape=(dim, 1))
      output_vector = self.dual_object.get_psd_product(input_vector)

      def np_vector_prod_fn(np_vector):
        np_vector = np.reshape(np_vector, [-1, 1])
        output_np_vector = self.sess.run(output_vector, feed_dict={input_vector:np_vector})
        return output_np_vector
      linear_operator = LinearOperator((dim, dim), matvec=np_vector_prod_fn)
      # Performing shift invert scipy operation when eig val estimate is available
      min_eig_vec_val, estimated_eigen_vector = eigs(linear_operator,
                                                     k=1, which='SR', tol=1E-4)
      min_eig_vec_val = np.reshape(np.real(min_eig_vec_val), [1, 1])
      return np.reshape(estimated_eigen_vector, [-1, 1]), min_eig_vec_val