#!/usr/bin/python
# coding=utf-8

# Python program to produce the in PlanetMath encyclopedia entry:
# http://planetmath.org/encyclopedia/ChangeOfVariablesInIntegralOnMathbbRn.html
#
# Both the PyX and SciPy packages, along with Python 2.x, 
# are required to run this program.
#
# Author: Steve Cheng <steve@gold-saucer.org>

from pyx import canvas, path, unit, color, style, text, trafo, deco

# Returns the affine transform Mq + c in two-dimensional space
def affine_transform(M, q, c):
  return (M[0][0]*q[0] + M[0][1]*q[1] + c[0],
          M[1][0]*q[0] + M[1][1]*q[1] + c[1])

# Use Catmull-Rom interpolation to draw a smooth Bézier curve
# through the given points q.  The shape parameter controls the degree
# of curving.  The points q are assumed to form a loop (a closed curve)
# if loop is set to True.
def interp_curve(q, shape=1.0, loop=False):
  shape /= 6.0

  if loop:
    I = range(0, len(q)-2) + [-1]
  else:
    I = range(1, len(q)-2)
  
  segments = [ (  (q[i][0] + (q[i+1][0]-q[i-1][0])*shape,
                   q[i][1] + (q[i+1][1]-q[i-1][1])*shape),
                  (q[i+1][0] - (q[i+2][0]-q[i][0])*shape,
                   q[i+1][1] - (q[i+2][1]-q[i][1])*shape),
                  q[i+1] )
               for i in I ]

  if loop:
    pa = path.path(path.moveto(q[0][0], q[0][1]))
  else:
    pa = path.path(path.moveto(q[1][0], q[1][1]))

  pa.append(path.multicurveto_pt(
                     [ (unit.topt(p[0][0]), unit.topt(p[0][1]),
                        unit.topt(p[1][0]), unit.topt(p[1][1]),
                        unit.topt(p[2][0]), unit.topt(p[2][1]) )
                       for p in segments ]))
  
  if loop:
    pa.append(path.closepath())

  return pa

# The function to plot and its derivative
def g(x, y):
  return (x+y*y/2, 2*y/(x*x))
def Dg(x, y):
  return ( (1           , y        ) , 
           (-4*y/(x*x*x), 2.0/(x*x)) )


unit.set(uscale=1.25)
cvsImage = canvas.canvas()
cvs = cvsImage

# Function values on a grid mesh
x_range = [ 1.0 + 0.1*i for i in xrange(0, 20+3) ]
y_range = [ 1.0 + 0.1*j for j in xrange(0, 20+3) ]

# Draw grid lines for every K points
K = 4
x_grid  = [ x_range[i] for i in xrange(1, len(x_range)-1, K) ]
y_grid  = [ y_range[j] for j in xrange(1, len(y_range)-1, K) ]
for x in x_grid:
  cvs.stroke(interp_curve([ g(x,y) for y in y_range ]),
             attrs=[ style.linestyle.dotted ])
for y in y_grid:
  cvs.stroke(interp_curve([ g(x,y) for x in x_range ]),
             attrs=[ style.linestyle.dotted ])

# Which sub-rectangle in to highlight
m = 1*K+1
n = 2*K+1

# Shade the image sub-rectangle
#cvs.fill(interp_curve( 
#           [ g(x_range[i]  , y_range[n]  ) for i in xrange(m,m+K+1) ] +
#           [ g(x_range[m+K], y_range[j]  ) for j in xrange(n,n+K+1) ] +
#           [ g(x_range[i]  , y_range[n+K]) for i in xrange(m+K,m-1,-1) ] +
#           [ g(x_range[m]  , y_range[j]  ) for j in xrange(n+K,n-1,-1) ],
#           loop = True),
#         attrs=[ color.gray(0.8) ])

# Join the corner points of the image sub-rectangle
#cvs.stroke(path.path( 
#            path.moveto(*g(x_range[m]  , y_range[n]  )),
#            path.lineto(*g(x_range[m+K], y_range[n]  )),
#            path.lineto(*g(x_range[m+K], y_range[n+K])),
#            path.lineto(*g(x_range[m]  , y_range[n+K])),
#            path.closepath()))

# Show the Jacobian approximation and label it
cx = (x_range[m]+x_range[m+K])/2
cy = (y_range[n]+y_range[n+K])/2
dx = (x_range[m+K]-x_range[m])/2
dy = (y_range[n+K]-y_range[n])/2
L = Dg(cx, cy)
q = g(cx, cy)
cvs.fill(path.path(
             path.moveto(*affine_transform(L, ( dx, dy), q)),
             path.lineto(*affine_transform(L, (-dx, dy), q)),
             path.lineto(*affine_transform(L, (-dx,-dy), q)),
             path.lineto(*affine_transform(L, ( dx,-dy), q)),
             path.closepath()),
          attrs=[ color.cmyk.Orange ])
p = affine_transform(L, ( -dx*1.1, dy*1.1), q)
#cvs.text(p[0], p[1], r"linear approx. to image $g(Q)$ " +
#                     r"by $y+{\rm D} g(x)(Q-x)$",
#                     [ text.parbox(3.5), text.halign.flushright ])
cvs.text(p[0], p[1], r"approximation to image $g(Q)$")

# Draw and label the image centre point
cvs.fill(path.circle(q[0],q[1], 0.025),
          attrs=[ color.grey.black ])
cvs.text(q[0]+dx/2, q[1]-dy/2, r"$y$")


# Draw the grid in domain
cvsDomain = canvas.canvas()
cvs = cvsDomain

for x in x_grid:
  cvs.stroke(path.line(x, y_range[1], x, y_range[-2]),
             attrs=[ style.linestyle.dotted ])
for y in y_grid:
  cvs.stroke(path.line(x_range[1], y, x_range[-2], y),
             attrs=[ style.linestyle.dotted ])

# Highlight the domain sub-rectangle and label it
cvs.fill(path.path(
             path.moveto(cx+dx, cy+dy),
             path.lineto(cx-dx, cy+dy),
             path.lineto(cx-dx, cy-dy),
             path.lineto(cx+dx, cy-dy),
             path.closepath()),
          attrs=[ color.cmyk.Orange ])
cvs.text(cx+dx*1.1, cy+dy*1.1, r"small rectangle $Q$")

# Draw and label the domain centre point
cvs.fill(path.circle(cx,cy, 0.025),
          attrs=[ color.grey.black ])
cvs.text(cx+dx/2, cy-dy/2, r"$x$")


cvs = canvas.canvas()

# Put the domain and image depictions 
cvs.insert(cvsDomain, 
  [ trafo.translate(-cvsDomain.bbox().right()-0.30, 0) ] )
cvs.insert(cvsImage, 
  [ trafo.translate(-cvsImage.bbox().left()+0.25, 0) ] )

# Draw curved mapping arrow
cvs.stroke(path.curve(-0.75, 3.0, 
                      -0.25, 3.25, 
                       0.25, 3.25, 
                       0.75, 3.0), 
           [ deco.earrow ])
cvs.text( 0, 3.5, r"differentiable mapping $g$", [ text.halign.center ])
 
# Write encapsulated PostScript file
cvs.writeEPSfile('jacobian.eps')