"""
Photon scattering in quantum optical systems

This module includes a collection of functions for numerically computing photon
scattering in driven arbitrary systems coupled to some configuration of output
waveguides. The implementation of these functions closely follows the
mathematical treatment given in K.A. Fischer, et. al., Scattering of Coherent
Pulses from Quantum Optical Systems (2017, arXiv:1710.02875).
"""
# Author:  Ben Bartlett
# Contact: benbartlett@stanford.edu

import numpy as np
from itertools import product, combinations_with_replacement
from ..core import basis, tensor, zero_ket, Qobj, QobjEvo
from .propagator import propagator, Propagator

__all__ = ['temporal_basis_vector',
           'temporal_scattered_state',
           'scattering_probability']


def set_partition(collection, num_sets):
    """
    Enumerate all ways of partitioning collection into num_sets different
    lists, e.g. :
    list(set_partition([1,2], 2))
    >>> [[[1, 2], []], [[1], [2]], [[2], [1]], [[], [1, 2]]]

    Parameters
    ----------
    collection : iterable
        Collection to generate a set partition of.
    num_sets : int
        Number of sets to partition collection into.

    Returns
    -------
    partition : iterable
        The partitioning of collection into num_sets sets.
    """
    for partitioning in product(range(num_sets), repeat=len(collection)):
        partition = [[] for _ in range(num_sets)]
        for i, set_index in enumerate(partitioning):
            partition[set_index].append(collection[i])
        yield tuple(tuple(indices) for indices in partition)


def photon_scattering_amplitude(propagator, c_ops, tlist, taus, psi, psit):
    """
    Compute the scattering amplitude for a system emitting into multiple
    waveguides.

    Parameters
    ----------
    propagator : :class:Propagator
        Propagator
    c_ops : list
        list of collapse operators for each waveguide; these are assumed to
        include spontaneous decay rates, e.g.
        :math:`\\sigma = \\sqrt \\gamma \\cdot a`
    tlist : array_like
        List of times starting at the beginning and ending at the end of the
        evolution.
    taus : list-like
        List of (list of emission times) for each waveguide.
    psi : Qobj
        State at the start of the evolution
    psit : Qobj
        State at the end of the evolution.
    """
    # Extract the full list of taus
    tau_collapse = []
    for i, tau_wg in enumerate(taus):
        for t in tau_wg:
            tau_collapse.append((t, i))
    tau_collapse.sort(key=lambda tup: tup[0])  # sort tau_collapse by time

    tq = tlist[0]
    # Compute Prod Ueff(tq, tq-1)
    for tau in tau_collapse:
        tprev = tq
        tq, q = tau
        psi = c_ops[q] * propagator(tq, tprev) * psi

    psi = propagator(tlist[-1], tq) * psi
    return psit.overlap(psi)


def _temporal_basis_idx(waveguide_emission_indices, n_time_bins):
    """
    Generate a the global index for the excitation.
    """
    idx = []
    for i, wg_indices in enumerate(waveguide_emission_indices):
        for index in wg_indices:
            idx += [i * n_time_bins + index]
    idx = idx or [0]
    return tuple(idx)


def _temporal_basis_dims(waveguide_emission_indices, n_time_bins,
                         n_emissions=None):
    """
    Return the dims of the ``temporal_basis_vector``.
    """
    num_col = len(waveguide_emission_indices)
    if n_emissions is None:
        n_emissions = sum(
            [len(waveguide_indices)
             for waveguide_indices in waveguide_emission_indices])
    n_emissions = n_emissions or 1
    return [num_col * n_time_bins] * n_emissions


def temporal_basis_vector(waveguide_emission_indices, n_time_bins,
                          n_emissions=None):
    """
    Generate a temporal basis vector for emissions at specified time bins into
    specified waveguides.

    Parameters
    ----------
    waveguide_emission_indices : list or tuple
        List of indices where photon emission occurs for each waveguide,
        e.g. [[t1_wg1], [t1_wg2, t2_wg2], [], [t1_wg4, t2_wg4, t3_wg4]].
    n_time_bins : int
        Number of time bins; the range over which each index can vary.

    Returns
    -------
    temporal_basis_vector : :class: qutip.Qobj
        A basis vector representing photon scattering at the specified indices.
        If there are W waveguides, T times, and N photon emissions, then the
        basis vector has dimensionality (W*T)^N.
    """
    idx = _temporal_basis_idx(waveguide_emission_indices, n_time_bins)
    dims = _temporal_basis_dims(waveguide_emission_indices,
                                n_time_bins, n_emissions)
    return basis(dims, list(idx))


def _temporal_scattered_matrix(H, psi0, n_emissions, c_ops, tlist,
                               system_zero_state=None,
                               construct_effective_hamiltonian=True):
    """
    Compute the scattered n-photon state as an ndarray.
    """
    T = len(tlist)
    W = len(c_ops)
    em_dims = max(n_emissions, 1)
    phi_n = np.zeros([W * T] * em_dims, dtype=complex)

    if construct_effective_hamiltonian:
        # Construct an effective Hamiltonian from system hamiltonian and c_ops
        Heff = QobjEvo(H) - 1j / 2 * sum([op.dag() * op for op in c_ops])
    else:
        Heff = H

    evolver = Propagator(Heff, memoize=len(tlist))

    all_emission_indices = combinations_with_replacement(range(T), n_emissions)

    if system_zero_state is None:
        system_zero_state = psi0

    # Compute <omega_tau> for all combinations of tau
    for emission_indices in all_emission_indices:
        # Consider unique partitionings of emission times into waveguides
        partition = tuple(set(set_partition(emission_indices, W)))
        # Consider all possible partitionings of time bins by waveguide
        for indices in partition:
            taus = [[tlist[i] for i in wg_indices] for wg_indices in indices]
            phi_n_amp = photon_scattering_amplitude(
                evolver, c_ops, tlist,
                taus, psi0, system_zero_state
            )
            # Add scatter amplitude times temporal basis to overall state
            idx = _temporal_basis_idx(indices, T)
            phi_n[idx] = phi_n_amp
    return phi_n


def temporal_scattered_state(H, psi0, n_emissions, c_ops, tlist,
                             system_zero_state=None,
                             construct_effective_hamiltonian=True):
    """
    Compute the scattered n-photon state projected onto the temporal basis.

    Parameters
    ----------
    H : :class: qutip.Qobj or list
        System-waveguide(s) Hamiltonian or effective Hamiltonian in Qobj or
        list-callback format. If construct_effective_hamiltonian is not
        specified, an effective Hamiltonian is constructed from `H` and
        `c_ops`.
    psi0 : :class: qutip.Qobj
        Initial state density matrix :math:`\\rho(t_0)` or state vector
        :math:`\\psi(t_0)`.
    n_emissions : int
        Number of photon emissions to calculate.
    c_ops : list
        List of collapse operators for each waveguide; these are assumed to
        include spontaneous decay rates, e.g.
        :math:`\\sigma = \\sqrt \\gamma \\cdot a`
    tlist : array_like
        List of times for :math:`\\tau_i`. tlist should contain 0 and exceed
        the pulse duration / temporal region of interest.
    system_zero_state : :class: qutip.Qobj
        State representing zero excitations in the system. Defaults to
        :math:`\\psi(t_0)`
    construct_effective_hamiltonian : bool
        Whether an effective Hamiltonian should be constructed from H and
        c_ops:
        :math:`H_{eff} = H - \\frac{i}{2} \\sum_n \\sigma_n^\\dagger \\sigma_n`
        Default: True.

    Returns
    -------
    phi_n : :class: qutip.Qobj
        The scattered bath state projected onto the temporal basis given by
        tlist. If there are W waveguides, T times, and N photon emissions, then
        the state is a tensor product state with dimensionality T^(W*N).
    """
    T = len(tlist)
    W = len(c_ops)
    em_dims = max(n_emissions, 1)
    phi_n = _temporal_scattered_matrix(
        H, psi0, n_emissions, c_ops, tlist,
        system_zero_state, construct_effective_hamiltonian
    )
    return Qobj(phi_n.ravel(), dims=[[W * T] * em_dims, [1] * em_dims])


def scattering_probability(H, psi0, n_emissions, c_ops, tlist,
                           system_zero_state=None,
                           construct_effective_hamiltonian=True):
    """
    Compute the integrated probability of scattering n photons in an arbitrary
    system. This function accepts a nonlinearly spaced array of times.

    Parameters
    ----------
    H : :class: qutip.Qobj or list
        System-waveguide(s) Hamiltonian or effective Hamiltonian in Qobj or
        list-callback format. If construct_effective_hamiltonian is not
        specified, an effective Hamiltonian is constructed from H and
        `c_ops`.
    psi0 : :class: qutip.Qobj
        Initial state density matrix :math:`\\rho(t_0)` or state vector
        :math:`\\psi(t_0)`.
    n_emissions : int
        Number of photons emitted by the system (into any combination of
        waveguides).
    c_ops : list
        List of collapse operators for each waveguide; these are assumed to
        include spontaneous decay rates, e.g.
        :math:`\\sigma = \\sqrt \\gamma \\cdot a`.
    tlist : array_like
        List of times for :math:`\\tau_i`. tlist should contain 0 and exceed
        the pulse duration / temporal region of interest; tlist need not be
        linearly spaced.
    system_zero_state : :class: qutip.Qobj
        State representing zero excitations in the system. Defaults to
        `basis(systemDims, 0)`.
    construct_effective_hamiltonian : bool
        Whether an effective Hamiltonian should be constructed from H and
        c_ops:
        :math:`H_{eff} = H - \\frac{i}{2} \\sum_n \\sigma_n^\\dagger \\sigma_n`
        Default: True.

    Returns
    -------
    scattering_prob : float
        The probability of scattering n photons from the system over the time
        range specified.
    """
    phi_n = _temporal_scattered_matrix(H, psi0, n_emissions, c_ops, tlist,
                                       system_zero_state,
                                       construct_effective_hamiltonian)
    T = len(tlist)
    W = len(c_ops)

    # Compute <omega_tau> for all combinations of tau
    all_emission_indices = combinations_with_replacement(range(T), n_emissions)
    probs = np.zeros([T] * n_emissions)

    # Project scattered state onto temporal basis
    for emit_indices in all_emission_indices:
        # Consider unique emission time partitionings
        partition = tuple(set(set_partition(emit_indices, W)))
        # wg_indices_list = list(set_partition(indices, W))
        for wg_indices in partition:
            idx = _temporal_basis_idx(wg_indices, T)
            amplitude = phi_n[idx]
            probs[emit_indices] += np.real(amplitude.conjugate() * amplitude)

    # Iteratively integrate to obtain single value
    while probs.shape != ():
        probs = np.trapz(probs, x=tlist)
    return np.abs(probs)
