# Source code for qutip.utilities

    """
    This module contains utility functions that are commonly needed in other
    qutip modules.
    """

    __all__ = ['n_thermal', 'clebsch', 'convert_unit']

    import numpy as np

    def n_thermal(w, w_th):
        """
        Return the number of photons in thermal equilibrium for an harmonic
        oscillator mode with frequency 'w', at the temperature described by
        'w_th' where :math:`\\omega_{\\rm th} = k_BT/\\hbar`.

        Parameters
        ----------

        w : *float* or *array*
            Frequency of the oscillator.

        w_th : *float*
            The temperature in units of frequency (or the same units as `w`).

        Returns
        -------

        n_avg : *float* or *array*

            Return the number of average photons in thermal equilibrium for a
            an oscillator with the given frequency and temperature.

        """

        if isinstance(w, np.ndarray):
            return 1.0 / (np.exp(w / w_th) - 1.0)

        else:
            if (w_th > 0) and np.exp(w / w_th) != 1.0:
                return 1.0 / (np.exp(w / w_th) - 1.0)
            else:
                return 0.0

    def _factorial_prod(N, arr):
        arr[:int(N)] += 1

    def _factorial_div(N, arr):
        arr[:int(N)] -= 1

    def _to_long(arr):
        prod = 1
        for i, v in enumerate(arr):
            prod *= (i+1)**int(v)
        return prod

    def clebsch(j1, j2, j3, m1, m2, m3):
        """Calculates the Clebsch-Gordon coefficient
        for coupling (j1,m1) and (j2,m2) to give (j3,m3).

        Parameters
        ----------
        j1 : float
            Total angular momentum 1.

        j2 : float
            Total angular momentum 2.

        j3 : float
            Total angular momentum 3.

        m1 : float
            z-component of angular momentum 1.

        m2 : float
            z-component of angular momentum 2.

        m3 : float
            z-component of angular momentum 3.

        Returns
        -------
        cg_coeff : float
            Requested Clebsch-Gordan coefficient.

        """
        if m3 != m1 + m2:
            return 0
        vmin = int(np.max([-j1 + j2 + m3, -j1 + m1, 0]))
        vmax = int(np.min([j2 + j3 + m1, j3 - j1 + j2, j3 + m3]))

        c_factor = np.zeros((int(j1 + j2 + j3 + 1)), np.int32)
        _factorial_prod(j3 + j1 - j2, c_factor)
        _factorial_prod(j3 - j1 + j2, c_factor)
        _factorial_prod(j1 + j2 - j3, c_factor)
        _factorial_prod(j3 + m3, c_factor)
        _factorial_prod(j3 - m3, c_factor)
        _factorial_div(j1 + j2 + j3 + 1, c_factor)
        _factorial_div(j1 - m1, c_factor)
        _factorial_div(j1 + m1, c_factor)
        _factorial_div(j2 - m2, c_factor)
        _factorial_div(j2 + m2, c_factor)
        C = np.sqrt((2.0 * j3 + 1.0)*_to_long(c_factor))

        s_factors = np.zeros(((vmax + 1 - vmin), (int(j1 + j2 + j3))), np.int32)
        sign = (-1) ** (vmin + j2 + m2)
        for i,v in enumerate(range(vmin, vmax + 1)):
            factor = s_factors[i,:]
            _factorial_prod(j2 + j3 + m1 - v, factor)
            _factorial_prod(j1 - m1 + v, factor)
            _factorial_div(j3 - j1 + j2 - v, factor)
            _factorial_div(j3 + m3 - v, factor)
            _factorial_div(v + j1 - j2 - m3, factor)
            _factorial_div(v, factor)
        common_denominator = -np.min(s_factors, axis=0)
        numerators = s_factors + common_denominator
        S = sum([(-1)**i * _to_long(vec) for i,vec in enumerate(numerators)]) * \
            sign / _to_long(common_denominator)
        return C * S

    # -----------------------------------------------------------------------------
    # Functions for unit conversions
    #
    _e = 1.602176565e-19  # C
    _kB = 1.3806488e-23   # J/K
    _h = 6.62606957e-34   # Js

    _unit_factor_tbl = {
        #   "unit": "factor that convert argument from unit 'unit' to Joule"
        "J": 1.0,
        "eV": _e,
        "meV": 1.0e-3 * _e,
        "GHz": 1.0e9 * _h,
        "mK": 1.0e-3 * _kB,
    }

    def convert_unit(value, orig="meV", to="GHz"):
        """
        Convert an energy from unit `orig` to unit `to`.

        Parameters
        ----------
        value : float / array
            The energy in the old unit.

        orig : string
            The name of the original unit ("J", "eV", "meV", "GHz", "mK")

        to : string
            The name of the new unit ("J", "eV", "meV", "GHz", "mK")

        Returns
        -------
        value_new_unit : float / array
            The energy in the new unit.
        """
        if orig not in _unit_factor_tbl:
            raise TypeError("Unsupported unit %s" % orig)

        if to not in _unit_factor_tbl:
            raise TypeError("Unsupported unit %s" % to)

        return value * (_unit_factor_tbl[orig] / _unit_factor_tbl[to])

    def convert_GHz_to_meV(w):
        """
        Convert an energy from unit GHz to unit meV.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # 1 GHz = 4.1357e-6 eV = 4.1357e-3 meV
        w_meV = w * 4.1357e-3
        return w_meV

    def convert_meV_to_GHz(w):
        """
        Convert an energy from unit meV to unit GHz.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # 1 meV = 1.0/4.1357e-3 GHz
        w_GHz = w / 4.1357e-3
        return w_GHz

    def convert_J_to_meV(w):
        """
        Convert an energy from unit J to unit meV.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # 1 eV = 1.602e-19 J
        w_meV = 1000.0 * w / _e
        return w_meV

    def convert_meV_to_J(w):
        """
        Convert an energy from unit meV to unit J.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # 1 eV = 1.602e-19 J
        w_J = 0.001 * w * _e
        return w_J

    def convert_meV_to_mK(w):
        """
        Convert an energy from unit meV to unit mK.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # 1 mK = 0.0000861740 meV
        w_mK = w / 0.0000861740
        return w_mK

    def convert_mK_to_meV(w):
        """
        Convert an energy from unit mK to unit meV.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # 1 mK = 0.0000861740 meV
        w_meV = w * 0.0000861740
        return w_meV

    def convert_GHz_to_mK(w):
        """
        Convert an energy from unit GHz to unit mK.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.
        """
        # h v [Hz] = kB T [K]
        # h 1e9 v [GHz] = kB 1e-3 T [mK]
        # T [mK] = 1e12 * (h/kB) * v [GHz]
        w_mK = w * 1.0e12 * (_h / _kB)
        return w_mK

    def convert_mK_to_GHz(w):
        """
        Convert an energy from unit mK to unit GHz.

        Parameters
        ----------
        w : float / array
            The energy in the old unit.

        Returns
        -------
        w_new_unit : float / array
            The energy in the new unit.

        """
        w_GHz = w * 1.0e-12 * (_kB / _h)
        return w_GHz

    def _version2int(version_string):
        str_list = version_string.split(
            "-dev")[0].split("rc")[0].split("a")[0].split("b")[0].split(
            "post")[0].split('.')
        return sum([int(d if len(d) > 0 else 0) * (100 ** (3 - n))
                    for n, d in enumerate(str_list[:3])])

    def _blas_info():
        config = np.__config__
        if hasattr(config, 'blas_ilp64_opt_info'):
            blas_info = config.blas_ilp64_opt_info
        elif hasattr(config, 'blas_opt_info'):
            blas_info = config.blas_opt_info
        else:
            blas_info = {}
        _has_lib_key = 'libraries' in blas_info.keys()
        blas = None
        if hasattr(config,'mkl_info') or \
                (_has_lib_key and any('mkl' in lib for lib in blas_info['libraries'])):
            blas = 'INTEL MKL'
        elif hasattr(config,'openblas_info') or \
                (_has_lib_key and any('openblas' in lib for lib in blas_info['libraries'])):
            blas = 'OPENBLAS'
        elif 'extra_link_args' in blas_info.keys() and ('-Wl,Accelerate' in blas_info['extra_link_args']):
            blas = 'Accelerate'
        else:
            blas = 'Generic'
        return blas

    def available_cpu_count():
        """
        Get the number of cpus.
        It tries to only get the number available to qutip.
        """
        import os
        import multiprocessing
        try:
            import psutil
        except ImportError:
            psutil = None
        num_cpu = 0

        if 'QUTIP_NUM_PROCESSES' in os.environ:
            # We consider QUTIP_NUM_PROCESSES=0 as unset.
            num_cpu = int(os.environ['QUTIP_NUM_PROCESSES'])

        if num_cpu == 0 and 'SLURM_CPUS_PER_TASK' in os.environ:
            num_cpu = int(os.environ['SLURM_CPUS_PER_TASK'])

        if num_cpu == 0 and hasattr(os, 'sched_getaffinity'):
            num_cpu = len(os.sched_getaffinity(0))

        if (
            num_cpu == 0
            and psutil is not None
            and hasattr(psutil.Process(), "cpu_affinity")
        ):
            num_cpu = len(psutil.Process().cpu_affinity())

        if num_cpu == 0:
            try:
                num_cpu = multiprocessing.cpu_count()
            except NotImplementedError:
                pass

        return num_cpu or 1