<commit_before>#!/usr/bin/python
# -*- coding: utf-8 -*-
##
# SMC.py: Tomgraphic models module
##
# © 2012 <NAME> (<EMAIL>) and
#        <NAME> (<EMAIL>)
#     
# This file is a part of the Qinfer project.
# Licensed under the AGPL version 3.
##
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

## FEATURES ###################################################################

from __future__ import division

## IMPORTS ####################################################################

import numpy as np
import scipy.linalg as la

from qinfer.abstract_model import Model

## CLASSES ####################################################################

class QubitStatePauliModel(Model):
    """
    Represents an experimental system with unknown quantum state,
    and a limited visibility projective measurement.
    
    States are represented in the Pauli representation.
    """    
    
    @property
    def n_modelparams(self):
        return 3
        
    @property
    def expparams_dtype(self):
        # return 'float' <---- This implies a two-index array of scalars,
        #                      but we need a one-index array of records.
        return [('axis', '3f4'), ('vis', 'float')]
        #                 ^<commit_msg>Started refactoring tomography in terms of QuTiP.<commit_after>#!/usr/bin/python
# -*- coding: utf-8 -*-
##
# SMC.py: Tomgraphic models module
##
# © 2012 <NAME> (<EMAIL>) and
#        <NAME>. Granade (<EMAIL>)
#     
# This file is a part of the Qinfer project.
# Licensed under the AGPL version 3.
##
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

## TODO #######################################################################
# - Refactor modelparams <-> Qobj into common functions that memoize over
#   bases for each nq.
# - Write QST, QPT models using new distributions.
# - Deprecate old functionality.
# - Write unit tests for new mps <-> Qobj conversions.

## FEATURES ###################################################################

from __future__ import division

## IMPORTS ####################################################################

import numpy as np
import scipy.linalg as la

from qinfer.abstract_model import Model
from qinfer.distributions import Distribution, SingleSampleMixin

try:
    import qutip as qt
    from distutils.version import LooseVersion
    _qt_version = LooseVersion(qt.version.version)
    if _qt_version < LooseVersion('3.1'):
        qt = None
except ImportError:
    qt = None

## FUNCTIONS #################################################################

def require_qutip():
    if qt is None:
        raise ImportError("QuTiP version 3.1.0 or later required.")

def conjugate(B, X):
    return np.dot(B, np.dot(X, B.conj().transpose()))

def _qubit_superdims(nq):
    ds = [2] * nq
    return [[ds, ds]] * 2

## CLASSES ####################################################################

class GinibreQubitDistribution(SingleSampleMixin, Distribution):
    """
    Creates a new uniform prior over density operators on :math:`n` qubits,
    using the rank-:math:`k` Ginibre distribution. Sampled density operators
    are represented by vectors in the Pauli basis excluding the traceful
    element.
    
    This class is implemented using QuTiP (v3.1.0 or later), and thus will not
    work unless QuTiP is installed.


    :param int nq: Number of qubits.
    :param int rank: Rank of density operators to be sampled. If ``None``,
        full-rank density operators are sampled.
    """
    
    def __init__(self, nq=1, rank=None):
        require_qutip()
        self._nq = nq
        self._rank = rank
        self._dim = 2**nq

        # This arcane line takes the vectorized Pauli basis transformation
        # used by QuTiP and converts into a form that's useful for us.
        # Notably an array (not a matrix!) that excludes the traceful parts.
        self._paulis = qt.visualization._pauli_basis(nq).data.todense().H[1:, :].view(np.ndarray)

    @property
    def n_rvs(self):
        return self._dim ** 2 - 1
        
    def _sample(self):
        # Generate and flatten a density operator, so that we can multiply it
        # by the transformation defined above.
        rho = qt.rand_dm_ginibre(self._dim, self._rank).data.todense().view(np.ndarray).flatten(order='C')
        return np.real(np.dot(self._paulis, rho))

class BCSZQubitDistribution(SingleSampleMixin, Distribution):
    """
    Represents the BCSZ prior over CPTP maps of a given Choi (Kraus) rank.
    The returned model parameters are a vectorization of a supermatrix
    in the Pauli basis, with the first row omitted.

    :param int nq: Number of qubits on which random variate CPTP maps act.
    :param int rank: Choi/Kraus rank of the random variate maps. For
        ``rank = None``, full-rank will be assumed, such that the BCSZ
        prior is the projection of the Hilbert-Schmidt prior onto the
        CPTP constraint.
    :param bool enforce_tp: If ``False``, relaxes the constraint that
        sampled channels are trace-preserving.
    """
    def __init__(self, nq=1, rank=None, enforce_tp=True):
        require_qutip()
        self._nq = nq
        self._basis = (
            qt.visualization._pauli_basis(nq) / np.sqrt(2**nq)
        ).data.todense().H.view(np.ndarray)
        self._dims = 2**nq
        self._superdims = 4**nq
        self._rank = rank
        self._enforce_tp = enforce_tp

    @property
    def n_rvs(self):
        # The top "row" of the superoperator in the Pauli basis is completely
        # determined by the assumption of a CPTP map (S_{0i} = delta_{0i}),
        # and so we strip it off, removing d² params.
        return self._superdims ** 2 - self._superdims

    def _sample(self):
        # Since the dense representation of S is a matrix and not an array,
        # flattening in FORTRAN order gives MATLAB-like behavior (implicitly
        # two-index). Thus, we need to index away the first axis.
        return np.real(
            conjugate(self._basis,
                qt.rand_super_bcsz(
                    self._dims, enforce_tp=self._enforce_tp, rank=self._rank
                ).data.todense().view(np.ndarray)
            )
        )[1:, :].flatten(order='F')

    def to_qobj(self, modelparams):
        return [
            qt.Qobj(
                conjugate(
                    self._basis.conj().transpose(),
                     np.vstack((
                         np.hstack((
                             [1], np.zeros((self._superdims - 1, ))
                         )),
                        modelparams.reshape((self._superdims - 1, self._superdims), order='F')
                    ))
                ),
                dims=_qubit_superdims(self._nq)
            )
            for x in modelparams
        ]

class QubitStatePauliModel(Model):
    """
    Represents an experimental system with unknown quantum state,
    and a limited visibility projective measurement.
    
    States are represented in the Pauli representation.
    """    
    
    @property
    def n_modelparams(self):
        return 3
        
    @property
    def expparams_dtype(self):
        # return 'float' <---- This implies a two-index array of scalars,
        #                      but we need a one-index array of records.
        return [('axis', '3f4'), ('vis', 'float')]
        #                 ^
{'fraction_non_alphanumeric': 0.14752763234438626, 'fraction_numerical': 0.005351948807446189, 'mean_word_length': 2.6735042735042733, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3592809', 'n_tokens_mistral': 2458, 'n_tokens_neox': 2340, 'n_words': 869}