Source Code for Module reflectometry.model1d.examples.demo_chisq_surfaces

 1  #!/usr/bin/env python 
 2  # This program is public domain. 
 3   
 4  """Pure numeric python implementation of Parratt's recursion for computing 
 5  reflectivity as a function of Q and lambda.""" 
 6   
 7  import pylab as Plot 
 8  import numpy as N 
 9  import parratt 
10   
11 -def theory(Q,d,rho): 
12      """Compute reflectometry theory from using parratt formalism""" 
13      r = parratt.refl(Q,d,rho) 
14      R = r*N.conj(r) 
15      return R.real 
16   
17 -def sampledata(Q,d,rho): 
18      """Generate sample data given the model Q,d,rho""" 
19      # Fresnel weighting 
20      dR = N.sqrt(theory(Q,[0,0],[rho[0],rho[-1]])) 
21      R = theory(Q,d,rho) + N.random.normal(scale=dR) 
22      return R,dR 
23   
24 -class Chisq: 
25      """Compute chisq(p) given f(p),y,dy""" 
26 -    def __init__(self,f,y,dy): 
27          self.y,self.dy,self.f = y,dy,f 
28 -    def __call__(self,p): 
29          return N.sum(((self.f(p)-self.y)/self.dy)**2) 
30   
31 -def build_landscape(fn,x,y,dim=(50,50)): 
32      """Construct a 2D landscape for the function f((x,y)) over 
33         x=[lo, hi], y=[lo,hi] with the number of steps dim=(nx,ny) 
34         defaulting to 50x50 
35      """ 
36      d1 = N.linspace(x[0],x[1],dim[0]) 
37      d2 = N.linspace(y[0],y[1],dim[1]) 
38      Z = N.zeros((len(d1),len(d2))) 
39      for i in range(len(d1)): 
40          for j in range(len(d2)): 
41              Z[i,j] = fn((d1[i],d2[j])) 
42      return d1,d2,Z 
43   
44 -def plot_mesh(x,y,Z): 
45      """Plot a regular mesh z with grid spacing x,y""" 
46      X,Y = N.meshgrid(x,y) 
47      #Plot.pcolor(X,Y,Z) 
48      Plot.contourf(X,Y,Z) 
49      Plot.show() 
50   
51 -class Demo_Model: 
52      """Example model.  Defines a reflectometry model object for a structure and a 
53         a dataset. 
54   
55            model(p) returns f(x;p) corresponding to dataset measured points. 
56            model.goodness(p) returns chisq for the given parameters. 
57            model.x indicates the measured points 
58            model.y, model.dy is the measurement and uncertainty at the points 
59            model.rho,model.d are the parameters to the reflectometry model 
60      """ 
61 -    def __init__(self): 
62          """Create a new instance of the demo model""" 
63          self.d = N.array([0,20,130,500,20,0],'d') 
64          self.rho = N.array([0,6e-4,5e-4,9e-4,10e-4,2e-4],'d') 
65          self.x = N.linspace(0.001,0.3,100) 
66          self.y,self.dy = sampledata(self.x,self.d,self.rho) 
67          self.goodness = Chisq(self,self.y,self.dy) 
68   
69 -    def __call__(self,p): 
70          """Evaluate the model at parameter vector p, returning f(x;p)""" 
71          # Protect the base values 
72          d = self.d.copy() 
73          rho = self.rho.copy() 
74   
75          # Substitute the fitting parameters 
76          d[2],d[3] = p 
77   
78          # Evaluate the function 
79          return theory(self.x,d,rho) 
80   
81  demo_model = Demo_Model() 
82   
83   
84 -def demo_landscape(): 
85      """Plot the chisq landscape of the demo model""" 
86      plot_mesh(*build_landscape(demo_model.goodness,[100,400],[300,900],dim=(100,100))) 
87   
88  if __name__ == "__main__": 
89      demo_landscape() 
90   
91  __id__ = "$Id: demo_parratt.py 37 2006-02-10 22:13:51Z pkienzle $" 
92