# Copyright (C) 2003 by Intevation GmbH
# Authors:
# Frank Koormann <frank.koormann@intevation.de>
#
# This program is free software under the GPL (>=v2)
# Read the file COPYING coming with the software for details.

"""SciParam: Stochastic Helper Functions

   source: GREAT-ER Model Worker

    Implements:
    - LogNorm_FivePercentileToStddev()
    - LogNormToNorm()
    - NormToLogNorm()
"""

__version__ = "$Revision: 1.1 $"
# $Source: /greaterrepository/sciparam/SciParam/stochastic.py,v $
# $Id: stochastic.py,v 1.1 2003/02/05 17:48:17 tkoester Exp $

from math import sqrt, log, exp

def LogNorm_FivePercentileToStddev(descriptives):
    """Derive stddev for log-normal distribution from 5 percentile and mean."""

    # P( Z >= z ) = 0.05, standard-normal distribution
    z = -1.6449

    if len(descriptives) == 2:
        mean    = float(descriptives[0])
        perfive = float(descriptives[1])
        if ( mean == 0.0 or perfive == 0.0 ):
            return 0.0
        else:
            # Stddev of corresponding normal distribution
            sy = sqrt(z*z - 2 * log( perfive / mean )) + z

            # 
            sx = mean * sqrt( exp( sy*sy ) - 1 )
            
            return sx

    else:
        return 0.0


def LogNormToNorm(descriptives):
    """Transform log-normal descriptives to corresponding normal dist."""

    if len(descriptives) == 2:
        log_mean = float(descriptives[0])
        log_sdev = float(descriptives[1])

        if log_sdev == 0.0 or log_mean == 0.0:
            norm_sdev = 0.0
            if log_mean == 0.0:
                norm_mean = 0.0
            else:
                norm_mean = log( log_mean )
        else:  
            norm_var  = log( ((log_sdev*log_sdev) / (log_mean*log_mean)) + 1 )
            norm_sdev = sqrt ( norm_var )

            norm_mean = log( log_mean ) - 0.5 * norm_var

        return ( norm_mean, norm_sdev )

    else:
        return ( 0.0, 0.0 )


def NormToLogNorm(descriptives):
    """Transform normal descriptives to corresponding log-normal dist."""

    if len(descriptives) == 2:
        norm_mean = float(descriptives[0])
        norm_sdev = float(descriptives[1])

        norm_var = norm_sdev * norm_sdev

        log_mean = exp( norm_mean + 0.5 * norm_var )

        temp = exp ( norm_var ) - 1
  
        if temp > 0.0:
            log_sdev = log_mean * sqrt ( temp )
        else:
            log_sdev = 0.0;

        return (log_mean, log_sdev)

    else:
        return ( 0.0, 0.0 )