gisas model with cylinders only

Examples/fit/gisas2d/fake.py → Examples/fit/gisas2d/gisas_fake1.py

@@ -3,10 +3,8 @@
 Fake GISAS data according to sample model from gisas_fit2d, with added noise.
 """
 
-import bornagain as ba
-from bornagain import deg, angstrom, nm
+import gisas_model1 as model
 import numpy as np
-import gisas_fit2d as model
 
 
 def fake_data():
@@ -14,12 +12,7 @@ def fake_data():
     Generate fake "experimental" data, and save them as numpy array.
     """
 
-    params = {
-        'lg(intensity)': 6,
-        'lg(background)': -99, # background is not part of theory, but added below
-        'cylinder_height': 4*nm,
-        'cylinder_radius': 5*nm,
-    }
+    params, bg = model.model_parameters_and_bg()
 
     # Compute model distribution
     simulation = model.get_simulation(params)
@@ -28,7 +21,7 @@ def fake_data():
 
     # Draw noisy data
     np.random.seed(0)
-    data = np.random.normal(theory, np.sqrt(theory)) + np.random.poisson(10**0.5, theory.shape)
+    data = np.random.normal(theory, np.sqrt(theory)) + np.random.poisson(bg, theory.shape)
 
     # Save to numpy
     np.savetxt("gisas-model1.txt.gz", data)

Examples/fit/gisas2d/gisas_fit2d.py

 #!/usr/bin/env python3
 """
-Fitting example: 4 parameters fit for mixture of cylinders and prisms on top
-of substrate.
+Basic GISAS 2D fitting example.
+The model is in gisas_model1.
+Fake experimental data are generated by gisas_fake1.
 """
 
+import gisas_model1 as model
 import bornagain as ba
-from bornagain import deg, nm
 import numpy as np
 from matplotlib import pyplot as plt
 
 
-def get_sample(params):
-    """
-    Returns a sample with uncorrelated cylinders and prisms on a substrate.
-    """
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
-
-    # defining materials
-    m_vacuum = ba.HomogeneousMaterial("Vacuum", 0, 0)
-    m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
-    m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
-    # collection of particles
-    ff = ba.FormFactorCylinder(cylinder_radius, cylinder_height)
-    cylinder = ba.Particle(m_particle, ff)
-    layout = ba.ParticleLayout()
-    layout.addParticle(cylinder)
-
-    # vacuum layer with particles and substrate form multi layer
-    vacuum_layer = ba.Layer(m_vacuum)
-    vacuum_layer.addLayout(layout)
-    substrate_layer = ba.Layer(m_substrate, 0)
-    multi_layer = ba.MultiLayer()
-    multi_layer.addLayer(vacuum_layer)
-    multi_layer.addLayer(substrate_layer)
-    return multi_layer
-
-
-def get_simulation(params):
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    beam = ba.Beam(10**params['lg(intensity)'], 0.1*nm, ba.Direction(0.2*deg, 0))
-    det = ba.SphericalDetector(100, -1.5*deg, 1.5*deg, 100, 0, 3*deg)
-    sample = get_sample(params)
-
-    simulation = ba.GISASSimulation(beam, sample, det)
-    simulation.setBackground(ba.ConstantBackground(10**params['lg(background)']))
-
-    return simulation
-
-
 def run_fitting():
-    """
-    Setup simulation and fit
-    """
-
     real_data = np.loadtxt("gisas-model1.txt.gz", dtype=float)
 
     fit_objective = ba.FitObjective()
-    fit_objective.addSimulationAndData(get_simulation, real_data)
+    fit_objective.addSimulationAndData(model.get_simulation, real_data)
 
-    # Print fit progress on every n-th iteration.
-    fit_objective.initPrint(10)
-
-    # Plot fit progress on every n-th iteration. Will slow down fit.
-    fit_objective.initPlot(10)
-
-    params = ba.Parameters()
-    params.add("lg(intensity)", 5)
-    params.add("lg(background)", 1)
-    params.add("cylinder_height", 6.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
+    fit_objective.initPrint(10) # Print on every 10th iteration.
+    fit_objective.initPlot(10) # Plot on every 10th iteration. Will slow down the fit.
 
     minimizer = ba.Minimizer()
+    params = model.start_parameters_1()
     result = minimizer.minimize(fit_objective.evaluate, params)
 
     fit_objective.finalize(result)
 
-    print("Fitting completed.")
+    print("Fit completed.")
     print("chi2:", result.minValue())
     for fitPar in result.parameters():
         print(fitPar.name(), fitPar.value, fitPar.error)
@@ -89,8 +36,5 @@ def run_fitting():
 
 
 if __name__ == '__main__':
-    # uncomment line below to regenerate "experimental" data file
-    # create_real_data()
-
     run_fitting()
     plt.show()

Examples/fit/gisas2d/gisas_model1.py

+"""
+Parametric model of a GISAS simulation.
+Dilute cylinders on a substrate.
+"""
+
+import bornagain as ba
+from bornagain import deg, nm
+
+
+mat_vacuum = ba.HomogeneousMaterial("Vacuum", 0, 0)
+mat_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
+mat_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
+
+
+def get_sample(params):
+    cylinder_height = params["cylinder_height"]
+    cylinder_radius = params["cylinder_radius"]
+
+    # collection of particles
+    ff = ba.FormFactorCylinder(cylinder_radius, cylinder_height)
+    cylinder = ba.Particle(mat_particle, ff)
+    layout = ba.ParticleLayout()
+    layout.addParticle(cylinder)
+
+    # vacuum layer with particles and substrate form multi layer
+    vacuum_layer = ba.Layer(mat_vacuum)
+    vacuum_layer.addLayout(layout)
+    substrate_layer = ba.Layer(mat_substrate, 0)
+    multi_layer = ba.MultiLayer()
+    multi_layer.addLayer(vacuum_layer)
+    multi_layer.addLayer(substrate_layer)
+    return multi_layer
+
+
+def get_simulation(params):
+    beam = ba.Beam(10**params['lg(intensity)'], 0.1*nm, ba.Direction(0.2*deg, 0))
+    det = ba.SphericalDetector(100, -1.5*deg, 1.5*deg, 100, 0, 3*deg)
+    sample = get_sample(params)
+
+    simulation = ba.GISASSimulation(beam, sample, det)
+    if 'lg(background)' in params:
+        simulation.setBackground(ba.ConstantBackground(10**params['lg(background)']))
+
+    return simulation
+
+
+def model_parameters_and_bg():
+    return ({
+        'lg(intensity)': 6,
+        'cylinder_height': 4*nm,
+        'cylinder_radius': 5*nm,
+    }, 10**0.5)
+
+
+def start_parameters_1():
+    params = ba.Parameters()
+    params.add("lg(intensity)", 5)
+    params.add("lg(background)", 1)
+    params.add("cylinder_height", 6.*nm, min=0.01)
+    params.add("cylinder_radius", 6.*nm, min=0.01)
+    return params