from __future__ import division

from scipy.optimize import *
from scipy.misc import derivative
from numpy import inf,linspace
from scipy.integrate import quad,quadrature
#from numpy.testing import *
from csv import DictWriter

#normalize transportation cost to t=1
t=.5
gu=.9
gb=1.1

q=[1.,1.]
p=[.5,.5]

def v(q,gamma,p):
  return max(q*gamma-q**2/2-p,0)

def d0(gamma,q,p):
  if v(q[0],gamma,p[0])+v(q[1],gamma,p[1])<t:
    return v(q[0],gamma,p[0])/t
  elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])>t:
    return 1
  elif abs(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))<=t:
    return 1/2+(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))/(2*t)
  elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])<-t:
    return 0

def d1(gamma,q,p):
  if v(q[0],gamma,p[0])+v(q[1],gamma,p[1])<t:
    return v(q[1],gamma,p[1])/t
  elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])>t:
    return 0
  elif abs(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))<=t:
    return 1/2+(v(q[1],gamma,p[1])-v(q[0],gamma,p[0]))/(2*t)
  elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])<-t:
    return 1


def br(p,q):
  def pdf(gamma):
    return 1/(gb-gu)
  def pi0(p0):
    def integrand(gamma):
      return max(p0-q[0]**2/2,0)*d0(gamma,q,[p0,p[1]])
    value, error = quad(integrand,gu,gb)
    return -value
  def pi1(p1):
    def integrand(gamma):
      return max(p1-q[1]**2/2,0)*d1(gamma,q,[p[0],p1])
    value, error = quad(integrand,gu,gb)
    return -value
  pmax = v(q[0],gb,0)
  x0=p[0]
  p0 = float(fmin_l_bfgs_b(pi0,x0,bounds=[(0,pmax)],approx_grad=1)[0])
  print fmin_l_bfgs_b(pi0,x0,approx_grad=1),pi0(p[0])
  x0=p[1]
  p1 = float(fmin_l_bfgs_b(pi1,x0,bounds=[(0,pmax)],approx_grad=1)[0])
  print fmin_l_bfgs_b(pi1,x0,approx_grad=1),pi1(p[1])
  return [p0,p1]

def eqP(p,q):
  print fixed_point(br,x0=p,args=(q,))
  return fixed_point(br,x0=p,args=(q,))

print eqP(p,q)