module Distribution
  #
  # Ruby version implements three methods on this module:
  # * Genz:: Used by default, with improvement to calculate p on rho > 0.95
  # * Hull:: Port from a C++ code
  # * Jantaravareerat:: Iterative (slow and buggy)
  #

  module BivariateNormal
    module Ruby_
      class << self
        SIDE = 0.1 # :nodoc:
        LIMIT = 5 # :nodoc:
        # Return the partial derivative of cdf over x, with y and rho constant
        # Reference:
        # * Tallis, 1962, p.346, cited by Olsson, 1979
        def partial_derivative_cdf_x(x, y, rho)
          Distribution::Normal.pdf(x) * Distribution::Normal.cdf((y - rho * x).quo(Math.sqrt(1 - rho**2)))
        end
        alias_method :pd_cdf_x,  :partial_derivative_cdf_x
        # Probability density function for a given x, y and rho value.
        #
        # Source: http://en.wikipedia.org/wiki/Multivariate_normal_distribution
        def pdf(x, y, rho, s1 = 1.0, s2 = 1.0)
          1.quo(2 * Math::PI * s1 * s2 * Math.sqrt(1 - rho**2)) * (Math.exp(-(1.quo(2 * (1 - rho**2))) *
            ((x**2.quo(s1)) + (y**2.quo(s2)) - (2 * rho * x * y).quo(s1 * s2))))
        end

        def f(x, y, aprime, bprime, rho)
          r = aprime * (2 * x - aprime) + bprime * (2 * y - bprime) + 2 * rho * (x - aprime) * (y - bprime)
          Math.exp(r)
        end

        # CDF for a given x, y and rho value.
        # Uses Genz algorithm (cdf_genz method).
        #
        def cdf(a, b, rho)
          cdf_genz(a, b, rho)
        end

        def sgn(x)
          if (x >= 0)
            1
          else
            -1
          end
        end

        # Normal cumulative distribution function (cdf) for a given x, y and rho.
        # Based on Hull (1993, cited by Arne, 2003)
        #
        # References:
        # * Arne, B.(2003). Financial Numerical Recipes in C ++. Available on  http://finance.bi.no/~bernt/gcc_prog/recipes/recipes/node23.html
        def cdf_hull(a, b, rho)
          # puts "a:#{a} - b:#{b} - rho:#{rho}"
          if a <= 0 && b <= 0 && rho <= 0
            # puts "ruta 1"
            aprime = a.quo(Math.sqrt(2.0 * (1.0 - rho**2)))
            bprime = b.quo(Math.sqrt(2.0 * (1.0 - rho**2)))
            aa = [0.3253030, 0.4211071, 0.1334425, 0.006374323]
            bb = [0.1337764, 0.6243247, 1.3425378, 2.2626645]
            sum = 0
            4.times do |i|
              4.times do |j|
                sum += aa[i] * aa[j] * f(bb[i], bb[j], aprime, bprime, rho)
              end
            end
            sum *= (Math.sqrt(1.0 - rho**2).quo(Math::PI))
            return sum
          elsif (a * b * rho <= 0.0)

            # puts "ruta 2"
            if a <= 0 && b >= 0 && rho >= 0
              return Distribution::Normal.cdf(a) - cdf(a, -b, -rho)
            elsif a >= 0.0 && b <= 0.0 && rho >= 0
              return Distribution::Normal.cdf(b) - cdf(-a, b, -rho)
            elsif a >= 0.0 && b >= 0.0 && rho <= 0
              return Distribution::Normal.cdf(a) + Distribution::Normal.cdf(b) - 1.0 + cdf(-a, -b, rho)
            end
          elsif (a * b * rho >= 0.0)
            # puts "ruta 3"
            denum = Math.sqrt(a**2 - 2 * rho * a * b + b**2)
            rho1 = ((rho * a - b) * sgn(a)).quo(denum)
            rho2 = ((rho * b - a) * sgn(b)).quo(denum)
            delta = (1.0 - sgn(a) * sgn(b)).quo(4)
            # puts "#{rho1} - #{rho2}"
            return cdf(a, 0.0, rho1) + cdf(b, 0.0, rho2) - delta
          end
          fail "Should'nt be here! #{a} - #{b} #{rho}"
        end

        # CDF. Iterative method by Jantaravareerat (n/d)
        #
        # Reference:
        # * Jantaravareerat, M. & Thomopoulos, N. (n/d). Tables for standard bivariate normal distribution

        def cdf_jantaravareerat(x, y, rho, s1 = 1, s2 = 1)
          # Special cases
          return 1 if x > LIMIT && y > LIMIT
          return 0 if x < -LIMIT || y < -LIMIT
          return Distribution::Normal.cdf(y) if  x > LIMIT
          return Distribution::Normal.cdf(x) if  y > LIMIT

          # puts "x:#{x} - y:#{y}"
          x = -LIMIT if x < -LIMIT
          x = LIMIT if x > LIMIT
          y = -LIMIT if y < -LIMIT
          y = LIMIT if y > LIMIT

          x_squares = ((LIMIT + x) / SIDE).to_i
          y_squares = ((LIMIT + y) / SIDE).to_i
          sum = 0
          x_squares.times do |i|
            y_squares.times do |j|
              z1 = -LIMIT + (i + 1) * SIDE
              z2 = -LIMIT + (j + 1) * SIDE
              # puts " #{z1}-#{z2}"
              h = (pdf(z1, z2, rho, s1, s2) + pdf(z1 - SIDE, z2, rho, s1, s2) + pdf(z1, z2 - SIDE, rho, s1, s2) + pdf(z1 - SIDE, z2 - SIDE, rho, s1, s2)).quo(4)
              sum += (SIDE**2) * h # area
            end
          end
          sum
        end
        # Normal cumulative distribution function (cdf) for a given x, y and rho.
        # Ported from Fortran code by Alan Genz
        #
        # Original documentation
        #    DOUBLE PRECISION FUNCTION BVND( DH, DK, R )
        #    A function for computing bivariate normal probabilities.
        #
        #       Alan Genz
        #       Department of Mathematics
        #       Washington State University
        #       Pullman, WA 99164-3113
        #       Email : alangenz_AT_wsu.edu
        #
        #    This function is based on the method described by
        #        Drezner, Z and G.O. Wesolowsky, (1989),
        #        On the computation of the bivariate normal integral,
        #        Journal of Statist. Comput. Simul. 35, pp. 101-107,
        #    with major modifications for double precision, and for |R| close to 1.
        #
        # Original location:
        # * http://www.math.wsu.edu/faculty/genz/software/fort77/tvpack.f
        def cdf_genz(x, y, rho)
          dh = -x
          dk = -y
          r = rho
          twopi = 6.283185307179586

          w = 11.times.collect { [nil] * 4 }
          x = 11.times.collect { [nil] * 4 }

          data = [
            0.1713244923791705E+00, -0.9324695142031522E+00,
            0.3607615730481384E+00, -0.6612093864662647E+00,
            0.4679139345726904E+00, -0.2386191860831970E+00]

          (1..3).each do|i|
            w[i][1] = data[(i - 1) * 2]
            x[i][1] = data[(i - 1) * 2 + 1]
          end
          data = [
            0.4717533638651177E-01, -0.9815606342467191E+00,
            0.1069393259953183E+00, -0.9041172563704750E+00,
            0.1600783285433464E+00, -0.7699026741943050E+00,
            0.2031674267230659E+00, -0.5873179542866171E+00,
            0.2334925365383547E+00, -0.3678314989981802E+00,
            0.2491470458134029E+00, -0.1252334085114692E+00]
          (1..6).each do|i|
            w[i][2] = data[(i - 1) * 2]
            x[i][2] = data[(i - 1) * 2 + 1]
          end
          data = [
            0.1761400713915212E-01, -0.9931285991850949E+00,
            0.4060142980038694E-01, -0.9639719272779138E+00,
            0.6267204833410906E-01, -0.9122344282513259E+00,
            0.8327674157670475E-01, -0.8391169718222188E+00,
            0.1019301198172404E+00, -0.7463319064601508E+00,
            0.1181945319615184E+00, -0.6360536807265150E+00,
            0.1316886384491766E+00, -0.5108670019508271E+00,
            0.1420961093183821E+00, -0.3737060887154196E+00,
            0.1491729864726037E+00, -0.2277858511416451E+00,
            0.1527533871307259E+00, -0.7652652113349733E-01]

          (1..10).each do|i|
            w[i][3] = data[(i - 1) * 2]
            x[i][3] = data[(i - 1) * 2 + 1]
          end

          if  r.abs < 0.3
            ng = 1
            lg = 3
          elsif  r.abs < 0.75
            ng = 2
            lg = 6
          else
            ng = 3
            lg = 10
          end

          h = dh
          k = dk
          hk = h * k
          bvn = 0
          if  r.abs < 0.925
            if  r.abs > 0
              hs = (h * h + k * k).quo(2)
              asr = Math.asin(r)
              (1..lg).each do |i|
                [-1, 1].each do |is|
                  sn = Math.sin(asr * (is * x[i][ng] + 1).quo(2))
                  bvn += w[i][ng] * Math.exp((sn * hk - hs).quo(1 - sn * sn))
                end # do
              end # do
              bvn *= asr.quo(2 * twopi)
            end # if
            bvn += Distribution::Normal.cdf(-h) * Distribution::Normal.cdf(-k)

          else # r.abs
            if  r < 0
              k = -k
              hk = -hk
            end

            if  r.abs < 1
              as = (1 - r) * (1 + r)
              a = Math.sqrt(as)
              bs = (h - k)**2
              c = (4 - hk).quo(8)
              d = (12 - hk).quo(16)
              asr = -(bs.quo(as) + hk).quo(2)
              if  asr > -100
                bvn = a * Math.exp(asr) * (1 - c * (bs - as) * (1 - d * bs.quo(5)).quo(3) + c * d * as * as.quo(5))
              end
              if  -hk < 100
                b = Math.sqrt(bs)
                bvn -= Math.exp(-hk.quo(2)) * Math.sqrt(twopi) * Distribution::Normal.cdf(-b.quo(a)) * b *
                       (1 - c * bs * (1 - d * bs.quo(5)).quo(3))
              end

              a = a.quo(2)
              (1..lg).each do |i|
                [-1, 1].each do |is|
                  xs = (a * (is * x[i][ng] + 1))**2
                  rs = Math.sqrt(1 - xs)
                  asr = -(bs / xs + hk).quo(2)
                  if  asr > -100
                    bvn += a * w[i][ng] * Math.exp(asr) *
                           (Math.exp(-hk * (1 - rs).quo(2 * (1 + rs))) .quo(rs) - (1 + c * xs * (1 + d * xs)))
                  end
                end
              end
              bvn = -bvn / twopi
            end

            if  r > 0
              bvn += Distribution::Normal.cdf(-[h, k].max)
            else
              bvn = -bvn
              if  k > h
                bvn = bvn + Distribution::Normal.cdf(k) - Distribution::Normal.cdf(h)
              end
            end
          end
          bvn
        end
        private :f, :sgn
      end
    end
  end
end
