# -*- coding: utf-8 -*-
#
#  This file is part of Sequana software
#
#  Copyright (c) 2016-2020 - Sequana Development Team
#
#  File author(s):
#      <NAME> <<EMAIL>>
#
#  Distributed under the terms of the 3-clause BSD license.
#  The full license is in the LICENSE file, distributed with this software.
#
#  website: https://github.com/sequana/sequana
#  documentation: http://sequana.readthedocs.io
#
##############################################################################
"""Simple VST transformation"""
from sequana.lazy import numpy as np
import colorlog
logger = colorlog.getLogger(__name__)


__all__ = ['VST']

class VST():
    """

    ::
    
        v = VST()
        v.anscombe(X)

    Reference: <NAME> and <NAME>, "Optimal inversion of the
        generalized Anscombe transformation for Poisson-Gaussian noise",
        IEEE Trans. Image Process., doi:10.1109/TIP.2012.2202675

    """
    def __init__(self):
        pass

    @staticmethod
    def anscombe(x):
        r"""Compute the anscombe variance stabilizing transform.

        :param x: noisy Poisson-distributed data
        :return: data with variance approximately equal to 1.

        Reference: Anscombe, <NAME>. (1948), "The transformation of Poisson,
            binomial and negative-binomial data", Biometrika 35 (3-4): 246-254

        For Poisson distribution, the mean and variance are not independent. The
        anscombe transform aims at transforming the data so that the variance
        is about 1 for large enough mean; For mean zero, the varaince is still
        zero. So, it transform Poisson data to approximately Gaussian data with
        mean :math:`\sqrt{x+3/8} - 1/(4m^{1/2})`
        """

        #if np.mean(x) <4:
        #    logger.warning("Mean of input data below 4")
        try:
            # If a dataframe, we do not want to change it
            return 2.0 * np.sqrt(x + 3.0/8.0)
        except:
            return 2.0 * np.sqrt(np.array(x) + 3.0/8.0)

    @staticmethod
    def inverse_anscombe(x):
        """Compute the inverse transform

        This uses an approximation of the exact  unbiased inverse.
        Reference: <NAME>., & <NAME>. (2011). A closed-form
        approximation of the exact unbiased inverse of the Anscombe
        variance-stabilizing transformation. Image Processing.
        """
        # should be positive
        # should wotk for dataframe and numpy arrays 
        return (1.0/4.0 * np.power(x, 2) +
                1.0/4.0 * np.sqrt(3.0/2.0) * np.power(x, -1.0) -
                11.0/8.0 * np.power(x, -2.0) + 
                5.0/8.0 * np.sqrt(3.0/2.0) * np.power(x, -3.0) - 1.0 / 8.0)

    @staticmethod
    def generalized_anscombe(x, mu, sigma, gain=1.0):
        """Compute the generalized anscombe variance stabilizing transform

        Data should be a mixture of poisson and gaussian noise.

        The input signal  z  is assumed to follow the Poisson-Gaussian noise
        model::
    
            x = gain * p + n

        where gain is the camera gain and mu and sigma are the read noise
        mean and standard deviation. X should contain only positive values.  
        Negative values are ignored. Biased for low counts
        """
        try:
            # If a dataframe, we do not want to change it
            y = gain*x + (gain**2) * 3.0/8.0 + sigma**2 - gain*mu
            return (2.0 / gain) * np.sqrt(np.maximum(y, 0.0))
        except:
            x = np.array(x)
            y = gain*x + (gain**2) * 3.0/8.0 + sigma**2 - gain*mu
            return (2.0 / gain) * np.sqrt(np.maximum(y, 0.0))


# -*- coding: utf-8 -*-
#
#  This file is part of Sequana software
#
#  Copyright (c) 2016-2020 - Sequana Development Team
#
#  File author(s):
#      <NAME> <<EMAIL>>
#
#  Distributed under the terms of the 3-clause BSD license.
#  The full license is in the LICENSE file, distributed with this software.
#
#  website: https://github.com/sequana/sequana
#  documentation: http://sequana.readthedocs.io
#
##############################################################################
"""Simple VST transformation"""
from sequana.lazy import numpy as np
import colorlog
logger = colorlog.getLogger(__name__)


__all__ = ['VST']

class VST():
    """

    ::
    
        v = VST()
        v.anscombe(X)

    Reference: <NAME> and <NAME>, "Optimal inversion of the
        generalized Anscombe transformation for Poisson-Gaussian noise",
        IEEE Trans. Image Process., doi:10.1109/TIP.2012.2202675

    """
    def __init__(self):
        pass

    @staticmethod
    def anscombe(x):
        r"""Compute the anscombe variance stabilizing transform.

        :param x: noisy Poisson-distributed data
        :return: data with variance approximately equal to 1.

        Reference: Anscombe, <NAME>. (1948), "The transformation of Poisson,
            binomial and negative-binomial data", Biometrika 35 (3-4): 246-254

        For Poisson distribution, the mean and variance are not independent. The
        anscombe transform aims at transforming the data so that the variance
        is about 1 for large enough mean; For mean zero, the varaince is still
        zero. So, it transform Poisson data to approximately Gaussian data with
        mean :math:`\sqrt{x+3/8} - 1/(4m^{1/2})`
        """

        #if np.mean(x) <4:
        #    logger.warning("Mean of input data below 4")
        try:
            # If a dataframe, we do not want to change it
            return 2.0 * np.sqrt(x + 3.0/8.0)
        except:
            return 2.0 * np.sqrt(np.array(x) + 3.0/8.0)

    @staticmethod
    def inverse_anscombe(x):
        """Compute the inverse transform

        This uses an approximation of the exact  unbiased inverse.
        Reference: <NAME>., & <NAME>. (2011). A closed-form
        approximation of the exact unbiased inverse of the Anscombe
        variance-stabilizing transformation. Image Processing.
        """
        # should be positive
        # should wotk for dataframe and numpy arrays 
        return (1.0/4.0 * np.power(x, 2) +
                1.0/4.0 * np.sqrt(3.0/2.0) * np.power(x, -1.0) -
                11.0/8.0 * np.power(x, -2.0) + 
                5.0/8.0 * np.sqrt(3.0/2.0) * np.power(x, -3.0) - 1.0 / 8.0)

    @staticmethod
    def generalized_anscombe(x, mu, sigma, gain=1.0):
        """Compute the generalized anscombe variance stabilizing transform

        Data should be a mixture of poisson and gaussian noise.

        The input signal  z  is assumed to follow the Poisson-Gaussian noise
        model::
    
            x = gain * p + n

        where gain is the camera gain and mu and sigma are the read noise
        mean and standard deviation. X should contain only positive values.  
        Negative values are ignored. Biased for low counts
        """
        try:
            # If a dataframe, we do not want to change it
            y = gain*x + (gain**2) * 3.0/8.0 + sigma**2 - gain*mu
            return (2.0 / gain) * np.sqrt(np.maximum(y, 0.0))
        except:
            x = np.array(x)
            y = gain*x + (gain**2) * 3.0/8.0 + sigma**2 - gain*mu
            return (2.0 / gain) * np.sqrt(np.maximum(y, 0.0))


# -*- coding: utf-8 -*- # #  This file is part of Sequana software # #  Copyright (c) 2016-2020 - Sequana Development Team # #  File author(s): #      <NAME> <<EMAIL>> # #  Distributed under the terms of the 3-clause BSD license. #  The full license is in the LICENSE file, distributed with this software. # #  website: https://github.com/sequana/sequana #  documentation: http://sequana.readthedocs.io # ############################################################################## Simple VST transformation ::              v = VST()         v.anscombe(X)      Reference: <NAME> and <NAME>, "Optimal inversion of the         generalized Anscombe transformation for Poisson-Gaussian noise",         IEEE Trans. Image Process., doi:10.1109/TIP.2012.2202675 Compute the anscombe variance stabilizing transform.          :param x: noisy Poisson-distributed data         :return: data with variance approximately equal to 1.          Reference: Anscombe, <NAME>. (1948), "The transformation of Poisson,             binomial and negative-binomial data", Biometrika 35 (3-4): 246-254          For Poisson distribution, the mean and variance are not independent. The         anscombe transform aims at transforming the data so that the variance         is about 1 for large enough mean; For mean zero, the varaince is still         zero. So, it transform Poisson data to approximately Gaussian data with         mean :math:`\sqrt{x+3/8} - 1/(4m^{1/2})` #if np.mean(x) <4: #    logger.warning("Mean of input data below 4") # If a dataframe, we do not want to change it Compute the inverse transform          This uses an approximation of the exact  unbiased inverse.         Reference: <NAME>., & <NAME>. (2011). A closed-form         approximation of the exact unbiased inverse of the Anscombe         variance-stabilizing transformation. Image Processing. # should be positive # should wotk for dataframe and numpy arrays Compute the generalized anscombe variance stabilizing transform          Data should be a mixture of poisson and gaussian noise.          The input signal  z  is assumed to follow the Poisson-Gaussian noise         model::                  x = gain * p + n          where gain is the camera gain and mu and sigma are the read noise         mean and standard deviation. X should contain only positive values.           Negative values are ignored. Biased for low counts # If a dataframe, we do not want to change it