from pylab import *

x = linspace(0., 1., 2001)

def truth(L1, L2, do_plot=False):
  p = exp(-L1*(x - 0.5)**2 - L2*sin(2.*pi*x/0.5)**2)
  Z = trapz(p, x=x)
  H = trapz(p/Z*log(p/Z), x=x)

  logZ = 200*log(Z)
  H *= 200

  if do_plot:
    plot(x, p)
    title([logZ, H])
    show()

  return [logZ, H]
  

def grid():
  # The two lagrange multipliers
  L1 = linspace(0., 20., 101)
  L2 = linspace(0., 20., 101)
  [L1, L2] = meshgrid(L1, L2)
  L2 = L2[::-1, :]

  logZ = L1.copy()
  info = L1.copy()

  for i in xrange(0, 101):
    for j in xrange(0, 101):
      [logZ[i, j], info[i, j]] = truth(L1[i, j], L2[i, j])
    print(i+1)

  subplot(1, 2, 1)
  imshow(logZ)

  subplot(1, 2, 2)
  imshow(info)
  show()


if __name__ == '__main__':
  truth(10., 1., do_plot=True)
  grid()



@@ -318,44 +318,75 @@
   # 
-The two lagrange multiplier
+Calculate log(Z) and H for some canonical distribution
 s%0A  
-L
+T
 1 = 
+exp(
 lins
@@ -390,29 +390,30 @@
 inspace(
-0., 2
+-2., 1
 0., 101)
 %0A  L2 = 
@@ -404,24 +404,29 @@
 0., 101)
+)
 %0A  
-L
+T
 2 = 
+exp(
 linspace
@@ -430,21 +430,22 @@
 ace(
-0., 2
+-2., 1
 0., 101)
 %0A  %5B
@@ -444,17 +444,18 @@
 101)
+)
 %0A  %5B
-L1, L
+T1, T
 2%5D =
@@ -468,24 +468,24 @@
 rid(
-L1, L
+T1, T
 2)%0A  
-L
+T
 2 = 
-L
+T
 2%5B::
@@ -491,17 +491,16 @@
 :-1, :%5D%0A
-%0A
   logZ =
@@ -500,17 +500,17 @@
  logZ = 
-L
+T
 1.copy()
@@ -516,16 +516,13 @@
 )%0A  
-info = L
+H = T
 1.co
@@ -602,20 +602,17 @@
 %5Bi, j%5D, 
-info
+H
 %5Bi, j%5D%5D 
@@ -619,17 +619,20 @@
 = truth(
-L
+1./T
 1%5Bi, j%5D,
@@ -632,17 +632,20 @@
 %5Bi, j%5D, 
-L
+1./T
 2%5Bi, j%5D)
@@ -661,16 +661,28 @@
 t(i+1)%0A%0A
+  figure(1)%0A
   subplo
@@ -709,44 +709,326 @@
 logZ
-)%0A%0A  subplot(1, 2, 2)%0A  imshow(info)
+, extent=(log(T1.min()), log(T1.max()), log(T2.min()), log(T2.max())))%0A  xlabel('ln($T_1$)')%0A  ylabel('ln($T_2$)')%0A  title('log(Z)')%0A%0A  subplot(1, 2, 2)%0A  imshow(H, extent=(log(T1.min()), log(T1.max()), log(T2.min()), log(T2.max())))%0A  xlabel('ln($T_1$)')%0A  ylabel('ln($T_2$)')%0A  title('H, max=%7Bm%7D'.format(m=H.max()))%0A
 %0A  s
@@ -1034,17 +1034,16 @@
 show()%0A%0A
-%0A
 if __nam

Bring into line with showresults plots