# Funciones de Julia (http://julialang.org/) para comprobar dominancia 
# estocastica (de 1er o 2o orden) entre distribuciones discretas finitas
# Definimos dos funciones: depo(pa,xa,pb,xb) y deso(pa,xa,pb,xb)
# Ricard Torres <rtorres@itam.mx>
"""
function depo(pa,xa,pb,xb)\n
Comprobar dominancia estocastica de primer orden entre las distribuciones\n
finitas discretas (pa,xa) y (pb,xb)\n
pi son las probabilidades asociadas a los distintos valores\n
xi son los valores que toman las v.a. (en orden ascendiente estricto)\n
Resultado:\n
  0 si hay dominacion mutua (las distribuciones coinciden)\n
  1 si (pa,xa) domina a (pb,xb)\n 
 -1 si (pa,xa) esta dominada por (pb,xb)\n
  2 si ninguna de las distribuciones domina a la otra
"""
function depo(pa,xa,pb,xb)
   la = length(xa)
   lb = length(xb)
   if la != length(pa) || lb != length(pb)
     error("Vectores de valores y probabilidades deben tener la misma longitud")
   end
   tol = 1e-10
   if max( abs(sum(pa)-1), abs(sum(pb)-1) ) > tol
      error("Vectores de probabilidad deben sumar 1")
   end
   if xa != sort(xa) || xb != sort(xb) || xa != unique(xa) || xb != unique(xb)
      error("Vectores de valores deben estar en orden ascendiente estricto")
   end
   xd = sort(unique(vcat(xa,xb))) # soporte comun de xa y xb
   ld = length(xd)
   fd = zeros(ld+1) # ponemos un cero al inicio
   ia = 1
   ib = 1
   for i = 2 : ld + 1
       fd[i] = fd[i-1] 
       if ia <= la && xd[i-1] == xa[ia]
          fd[i] += pa[ia]
          ia += 1
       end
       if ib <= lb && xd[i-1] == xb[ib]
          fd[i] -= pb[ib]
          ib += 1
       end
   end
   adomina = true
   bdomina = true
   for i = 2 : ld
       if fd[i] > tol
          adomina = false
       elseif fd[i] < - tol
          bdomina = false
       end
   end
   if adomina && bdomina 
      return 0   # dominacion mutua -> igualdad
   elseif adomina
      return 1   # (xa,pa) domina a (xb,pb)
   elseif bdomina
      return -1  # (xb,pb) domina a (xa,pa)
   else
      return 2   # ninguna domina a la otra
   end
end

"""
function deso(pa,xa,pb,xb)\n
Comprobar dominancia estocastica de segundo orden entre las distribuciones\n
finitas discretas (pa,xa) y (pb,xb)\n
pi son las probabilidades asociadas a los distintos valores\n
xi son los valores (en orden ascendiente estricto)\n
Resultado:\n
  0 si hay dominacion mutua (las distribuciones coinciden)\n
  1 si (pa,xa) domina a (pb,xb)\n 
 -1 si (pa,xa) esta dominada por (pb,xb)\n
  2 si ninguna distribucion domina a la otra
"""
function deso(pa,xa,pb,xb)
   la = length(xa)
   lb = length(xb)
   if la != length(pa) || lb != length(pb)
     error("Vectores de valores y probabilidades deben tener la misma longitud")
   end
   tol = 1e-10
   if max( abs(sum(pa)-1), abs(sum(pb)-1) ) > tol
      error("Vectores de probabilidad deben sumar 1")
   end
   if xa != sort(xa) || xb != sort(xb) || xa != unique(xa) || xb != unique(xb)
      error("Vectores de valores deben estar en orden ascendiente estricto")
   end
   xd = sort(unique(vcat(xa,xb))) #  soporte comun de xa y xb
   ld = length(xd)
   fd = zeros(ld+1) #  ponemos un cero al inicio
   ia = 1
   ib = 1
   for i = 2 : ld + 1
       fd[i] = fd[i-1] 
       if ia <= la && xd[i-1] == xa[ia]
          fd[i] += pa[ia]
          ia += 1
       end
       if ib <= lb && xd[i-1] == xb[ib]
          fd[i] -= pb[ib]
          ib += 1
       end
   end
   adomina = true
   bdomina = true
   integd = zeros(ld)
   for i = 2 : ld
       integd[i] = integd[i-1] + fd[i] * ( xd[i] - xd[i-1] )
       if integd[i] > tol
          adomina = false
       elseif integd[i] < - tol
          bdomina = false
       end
   end
   if adomina && bdomina 
      return 0   # dominacion mutua -> igualdad
   elseif adomina
      return 1   # (xa,pa) domina a (xb,pb)
   elseif bdomina
      return -1  # (xb,pb) domina a (xa,pa)
   else
      return 2   # ninguna domina a la otra
   end
end