Contents of /trunk/ripley/test/python/cook1.py

Revision 3479 - (show annotations)
Wed Mar 23 04:03:22 2011 UTC (10 years ago) by gross
File MIME type: text/x-python
File size: 6382 byte(s)

 1 import numpy as np 2 import matplotlib.pyplot as plt 3 import random 4 5 from esys.escript import * 6 from esys.escript.linearPDEs import LinearPDE 7 from esys.escript.pdetools import Locator 8 from esys.finley import Rectangle,Brick 9 10 import scipy.optimize as so 11 12 NE=200 13 NE=100 14 H=1. 15 padding=0.0 16 numContours=25 17 numPlotsGrid=25 18 19 class PoissonTest(object): 20 def __init__(self, domain, a, b, xc): 21 self.domain=domain 22 self.a=a 23 self.b=b 24 self.xc=xc 25 self.pde=LinearPDE(domain) 26 self.pde.setSymmetryOn() 27 x=self.domain.getX() 28 DIM=domain.getDim() 29 self.pde.setValue(A=kronecker(DIM), q=whereZero(x[DIM-1]-inf(x[DIM-1]))) 30 self.u_ref=self.getU(a, b, xc) 31 self.omega=1. 32 33 def getU(self, a, b, xc): 34 x=Function(self.domain).getX() 35 self.pde.setValue(Y=exp(a-(b*length(x-xc))**2)) 36 u=self.pde.getSolution() 37 print "PoissonTest.getU: a,b,xc, u =",a,b,xc,inf(u),sup(u) 38 return u 39 40 def getCostFunction(self, a, b, xc): 41 u=interpolate(self.getU(a, b, xc), FunctionOnBoundary(self.domain)) 42 C=0.5*integrate( self.omega * (u-self.u_ref)**2) 43 return C 44 45 class PoissonTestTopData(PoissonTest): 46 def __init__(self, domain, a, b, xc): 47 PoissonTest.__init__(self, domain, a, b, xc) 48 x=FunctionOnBoundary(self.domain).getX() 49 DIM=self.domain.getDim() 50 self.omega=whereZero(x[DIM-1]-sup(x[DIM-1])) 51 class PoissonTestDiscreteData(PoissonTest): 52 def __init__(self, N, sigma, domain, a, b, xc): 53 PoissonTest.__init__(self, domain, a, b, xc) 54 self.N=N 55 self.sigma=sigma 56 x=domain.getX() 57 x_min=inf(x[0]) 58 x_max=sup(x[0]) 59 z_max=sup(x[self.domain.getDim()-1]) 60 dx=(x_max-x_min)/N 61 self.loc=Locator(self.domain, [ (dx/2+dx*i, z_max) for i in xrange(N) ]) 62 uu=self.loc(self.u_ref) 63 s=self.sigma*(max(uu)+min(uu))/2. 64 print self.loc.getX() 65 self.data = numpy.array([ random.gauss(d, s) for d in self.loc(self.u_ref) ]) 66 print self.data 67 def getCostFunction(self, a, b, xc): 68 u=numpy.array(self.loc(self.getU(a, b, xc))) 69 C=(0.5/self.N)*numpy.sum((u-self.data)**2) 70 return C 71 72 73 74 def plotContour(p, filename, title): 75 x0=np.linspace(-padding, 1+padding, numPlotsGrid) 76 X0, Y0 = np.meshgrid(x0, H*x0) 77 78 plt.figure() 79 C=X0.copy() 80 for i in range(C.shape[0]): 81 for j in range(C.shape[1]): 82 C[i,j]=p.getCostFunction(p.a,p.b, (X0[i,j], Y0[i,j])) 83 84 C_max=np.amax(C) 85 C_min=np.amin(C) 86 levels = np.linspace(C_min, C_max, numContours) 87 cset1 = plt.contourf(X0, Y0, C, levels, cmap=plt.cm.get_cmap('jet', len(levels)-1)) 88 cset2 = plt.contour(X0, Y0, C, cset1.levels, colors = 'k', hold='on') 89 for c in cset2.collections: c.set_linestyle('solid') 90 plt.title('Cost function (%s)'%title) 91 plt.xlabel("xc_0") 92 plt.ylabel("xc_1") 93 plt.colorbar(cset1) 94 plt.savefig(filename) 95 96 97 98 plotContour(PoissonTestDiscreteData(5, 0.10, Rectangle(NE, int(NE*H),l1=H), 0,15,(0.5,0.5)), "cost_passion_discrete_w_cc.png","b=10, xc=(0.5,0.2)") 99 100 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.2)), "cost_passion_top_w_cb.png","b=10, xc=(0.5,0.2)") 101 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.5)), "cost_passion_top_w_cc.png","b=10, xc=(0.5,0.5)") 102 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.8)), "cost_passion_top_w_ct.png","b=10, xc=(0.5,0.8)") 103 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.2)), "cost_passion_top_w_lb.png","b=10, xc=(0.2,0.2)") 104 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.5)), "cost_passion_top_w_lc.png","b=10, xc=(0.2,0.5)") 105 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.8)), "cost_passion_top_w_lt.png","b=10, xc=(0.2,0.8)") 106 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.2)), "cost_passion_all_w_cb.png","b=10, xc=(0.5,0.2)") 107 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.5)), "cost_passion_all_w_cc.png","b=10, xc=(0.5,0.5)") 108 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.8)), "cost_passion_all_w_ct.png","b=10, xc=(0.5,0.8)") 109 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.2)), "cost_passion_all_w_lb.png","b=10, xc=(0.2,0.2)") 110 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.5)), "cost_passion_all_w_lc.png","b=10, xc=(0.2,0.5)") 111 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.8)), "cost_passion_all_w_lt.png","b=10, xc=(0.2,0.8)") 112 113 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.2)), "cost_passion_top_n_cb.png","b=100, xc=(0.5,0.2)") 114 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.5)), "cost_passion_top_n_cc.png","b=100, xc=(0.5,0.5)") 115 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.8)), "cost_passion_top_n_ct.png","b=100, xc=(0.5,0.8)") 116 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.2)), "cost_passion_top_n_lb.png","b=100, xc=(0.2,0.2)") 117 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.5)), "cost_passion_top_n_lc.png","b=100, xc=(0.2,0.5)") 118 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.8)), "cost_passion_top_n_lt.png","b=100, xc=(0.2,0.8)") 119 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.2)), "cost_passion_all_n_cb.png","b=100, xc=(0.5,0.2)") 120 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.5)), "cost_passion_all_n_cc.png","b=100, xc=(0.5,0.5)") 121 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.8)), "cost_passion_all_n_ct.png","b=100, xc=(0.5,0.8)") 122 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.2)), "cost_passion_all_n_lb.png","b=100, xc=(0.2,0.2)") 123 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.5)), "cost_passion_all_n_lc.png","b=100, xc=(0.2,0.5)") 124 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.8)), "cost_passion_all_n_lt.png","b=100, xc=(0.2,0.8)") 125 126 127 128 129