#!/usr/bin/python3
#exec(open('/path/to/file/ConsMaxM.py').read())
"""
Nel problema max u(x,y) su px+qy=m
trova il massimo m tale che x=0
Nel problema max u(x,y) su px+qy=m con preferenze quasilineari del tipo
u(x,y)=ax+bv(y) con v concava Inada
u_x/u_y = a/bv'(y) con x=0 diventa a/bv'(m/q) che per m che va da zero a infinito
passa da sotto a sopra p/q quindi per m piccolo x=0 e dopo il threshold m definito da
a/bv'(m/q)=p/q avremo x>0 nel punto di ottimo.
Nell'esempio di sotto v(y)=sqrt(y) che può anche essere ln(y) ecc.
"""
from sympy import *

x = Symbol('x')
y = Symbol('y')
m = Symbol('m')

#funzione di utilità 
#funzione
def u(x,y):
	a=1; b=10
	return a*x+b*sqrt(y)
#utilità marginali
ux=diff(u(x,y),x)
uy=diff(u(x,y),y)
print('il rapporto u_x/u_y è uguale a',ux/uy)
#prezzi
[p,q]=[1,2]

#soluzione problema
S=solve([ux/uy-p/q, p*x+q*y-m,q*y-m])
[x0,y0,M]=[float(S[0][x]),float(S[0][y]),float(S[0][m])]
print('la soluzione è [x,y,m]=',[x0,y0,M])

#plot
from numpy import linspace, meshgrid, arange
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.lines as lines
import numpy as np

fig = plt.figure()
#plot budget line
t=linspace(0,M/p,400)
b=M/q-(p/q)*t
plt.plot(t,b,'r')
#plot three indifference curves
xbar=2.5*M/p
ybar=1.6*M/q
X = linspace(0, xbar,400)
Y = linspace(0, ybar,400)
X,Y=meshgrid(X,Y)
uu=lambdify((x,y), u(x,y), modules=['numpy'])
z=uu(X,Y)
U=u(x0,y0)
levels=[0.8*U,U,1.2*U]
cs=plt.contour(X,Y,z,levels)
plt.clabel(cs,fontsize=12)
plt.xlabel('Bene x',fontsize=14)
plt.ylabel('Bene y',fontsize=14)
plt.xlim([0, xbar])
plt.ylim([0, ybar])
plt.title('Problema del consumatore',fontsize=18)
# save figure if you want
#fig.savefig('ConsMaxM.pdf')
plt.show()