# Funciones de Julia (http://julialang.org/) para operar con 
# curvas de Lorenz discretas
# Funciones definidas:
#  lorenz_disc(n,y) 
#  lorenz(n,y,t)
#  lorenz_domin(na,ya,nb,yb)
#  grafica_lorenz(n,y)
#  grafica_dos_lorenz(na,ya,nb,yb)
# Ricard Torres <rtorres@itam.mx>

module lorenzgraficar
export lorenz_disc
export lorenz
export lorenz_domin
export grafica_lorenz
export grafica_dos_lorenz
using PyPlot 

"""
function lorenz_disc(n,y)\n
Curva de Lorenz discreta, es decir, frecuencias poblacionales\n
acumuladas y proporciones de ingreso acumuladas\n
Argumentos:\n
n: vector de numero de individuos por categoria de ingreso\n
y: vector de categorias de ingreso\n
Output: [F S]\n
F: frecuencias pobl acumuladas\n
S: proporciones ingreso acumuladas
"""
function lorenz_disc(n,y)
   ln = length(n)
   ly = length(y)
   if ln != ly 
      error("vectores 'n' e 'y' deben tener la misma longitud")
   end
   f = n / sum(n)
   # Ingreso agregado
   Y =  dot(n , y)
   # Proporciones del ingreso agregado por categoria
   s = (1/Y) * ( n .* y )
   # Vectores cumulativos (a~nadiendo un cero al inicio)
   F = zeros(1+ln)
   S = zeros(1+ly)
   F[2:end] = cumsum(f)
   S[2:end] = cumsum(s)
   return [F S]
end

"""
function lorenz(n,y,t)\n
Curva de Lorenz para proporciones arbitrarias (interpolacion lineal)\n
de una distribucion discreta de riqueza\n
Argumentos:\n 
n: vector de numero de individuos por categoria de ingreso\n
y: vector de categorias de ingreso\n
t: proporcion (o vector de proporciones)
"""
function lorenz(n,y,t)
  lt = length(t)
  ln = length(n)
  ly = length(y)
  if ln != ly 
      error("vectores 'n' e 'y' deben tener la misma longitud")
  end
  Z = lorenz_disc(n,y)
  lz = length(Z[:,1])
  lor = zeros(lt)
  for k = 1 : lt
      i1 = lz 
      for i = 2 : lz-1
          if t[k] < Z[i,1]
             i1 = i
             break
          end
      end 
      i2 = i1 - 1
      w = ( Z[i1,1] - t[k] ) / ( Z[i1,1] - Z[i2,1] )
      lor[k] =  w * Z[i2,2] + (1-w) * Z[i1,2]
  end
  if lt == 1
     return sum(lor) # convertir a escalar
  else
     return lor
  end
end

"""
function lorenz_domin(na,ya,nb,yb)\n
Comprobar si hay dominacion en el sentido de Lorenz entre\n
distribuciones discretas de riqueza (es suficiente comprobarlo\n
para los puntos que definen las curvas)\n
Retorno:\n 
  0: curvas coinciden\n 
  1: (na,ya) domina a (nb,yb)\n
 -1: (nb,yb) domina a (na,ya)\n 
  2: no dominacion
"""
function lorenz_domin(na,ya,nb,yb)
   la = length(na)
   lb = length(nb)
   if la != length(ya) || lb != length(yb)
      error("Vectores de poblacion e ingreso deben tener la misma longitud")
   end
   Za = lorenz_disc(na,ya)
   Zb = lorenz_disc(nb,yb)
   F = vcat( Za[2:end-1,1], Zb[2:end-1,1] ) # Union de rangos de las F
   F = sort(unique(F)) # Eliminar repeticiones, ordenar
   lF = length(F)
   lora = lorenz(na,ya,F)
   lorb = lorenz(nb,yb,F)
   tol = 1e-15 # Tolerancia para comparaciones
   a_encima = true 
   b_encima = true
   for i = 1 : lF
       if ( lora[i] > lorb[i] && lora[i]-lorb[i] > tol )
          b_encima = false
       elseif ( lorb[i] > lora[i] && lorb[i]-lora[i] > tol  )
          a_encima = false
       end
   end
   if a_encima && b_encima 
      return 0   # dominacion mutua -> igualdad
   elseif a_encima
      return 1   # (xa,pa) domina a (xb,pb)
   elseif b_encima
      return -1  # (xb,pb) domina a (xa,pa)
   else
      return 2   # ninguna domina a la otra
   end
end

"""
function grafica_lorenz(n, x)\n
Graficar la curva de Lorenz de una distribucion de ingresos finita (n, x)\n
n: son las frecuencias absolutas de cada categoria de ingreso\n
x: son los categorias de ingreso (en estricto orden ascendiente)
"""
function grafica_lorenz(n, x)
   lx = length(x)
   if lx != length(n)
      error("Categorias y frecuencias deben tener la misma longitud")
   end
   if any(x != sort(unique(x))) 
      error("Categorias deben estar en orden ascendiente estricto")
   end
   Z = lorenz_disc(n,x)
   figure()
   plot([0,1],[0,1],linestyle="--",linewidth=1,color="green")
   plot(Z[:,1],Z[:,2],linestyle="-",linewidth=2,color="blue")
   plot(Z[:,1],Z[:,2],linestyle="None",linewidth=2,marker="o") 
   axis([-0.05,1.05,-0.05,1.05])
   xlabel("Frecuencias acumuladas poblacion") 
   ylabel("Proporciones acumuladas ingreso")
   title("Curva de Lorenz")
end  # function

"""
function grafica_dos_lorenz(pa, xa, pb, xb)\n
Graficar las curvas de Lorenz de dos distribuciones de ingresos finitas\n
ni: son las frecuencias absolutas (numeros) de cada categoria\n
xi: son las categorias de ingreso (en estricto orden ascendiente)
"""
function grafica_dos_lorenz(na, xa, nb, xb)
   la = length(xa)
   lb = length(xb)
   if la != length(na) || lb != length(nb)
      error("Categorias y frecuencias deben tener la misma longitud")
   end
   if any(xa != sort(unique(xa))) || any(xb != sort(unique(xb))) 
      error("Categorias deben estar en orden ascendiente estricto")
   end
   Za = lorenz_disc(na,xa)
   Zb = lorenz_disc(nb,xb)
   figure()
   plot([0,1],[0,1],linestyle="--",linewidth=1,color="green")
   ha = plot(Za[:,1],Za[:,2],linestyle="-",linewidth=2,color="blue")
   plot(Za[:,1],Za[:,2],linestyle="None",linewidth=2,marker="o",color="blue") 
   hb = plot(Zb[:,1],Zb[:,2],linestyle="-",linewidth=2,color="red")
   plot(Zb[:,1],Zb[:,2],linestyle="None",linewidth=2,marker="o",color="red") 
   axis([-0.05,1.05,-0.05,1.05])
   legend([ha;hb],["Serie A","Serie B"],loc="upper left")
   xlabel("Frecuencias acumuladas poblacion") 
   ylabel("Proporciones acumuladas ingreso")
   title("Curvas de Lorenz")
end  # function

end  # module

using lorenzgraficar


# Ejemplo 1:
# Para generar datos, cargamos la libreria de distribuciones
# (eso toma su tiempo, pues de otro modo la ejecucion seria instantanea)
# using Distributions
# du = DiscreteUniform # Usar dist uniforme discreta para generar datos
# ly = 5 # Numero de categorias de ingreso
# y = sort(unique(rand(du(100,800),ly))) # ingreso en orden ascendiente
# ly = length(y) # por si habia repeticiones
# n = rand(du(50,400),ly) # Poblacion no tiene que estar ordenada
#  Graficos (usaremos PyPlot)
# using PyPlot
# Z = lorenz_disc(n,y)
# plot([0,1],[0,1],linestyle="--",color="green")
# plot(Z[:,1],Z[:,2],linestyle="-",color="blue",linewidth=2.0,marker="o")
# xlabel("Frecuencias acumuladas poblacion")
# ylabel("Proporciones acumuladas ingreso")
# axis([0,1,0,1])
# title("Curva de Lorenz")
#
# Ejemplo 2:
# na = [79,  62,  111,  62,  62]
# ya = [376,  429,  529,  618,  742]
# nb = [211,  56,  162,  134,  167,  235,  174,  167]
# yb = [109,  146,  233,  348,  428,  464,  513,  745]
# Za = lorenz_disc(na,ya)
# Zb = lorenz_disc(nb,yb)
# # using PyPlot
# # plot([0,1],[0,1],linestyle="--",color="green")
# # plot(Za[:,1],Za[:,2],linestyle="-",color="blue",linewidth=2.0,marker="o")
# # plot(Zb[:,1],Zb[:,2],linestyle="-",color="red",linewidth=2.0,marker="o")
# # legend(["Diagonal","Serie A","Serie B"],loc="upper left")
# # xlabel("Frecuencias acumuladas poblacion")
# # ylabel("Proporciones acumuladas ingreso")
# # title("Curvas de Lorenz")