| 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 |
|
Parent Directory
| 
Revision Log