Commit bb3335cc authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

repair and rename basic GISAS fit example

parent 94328883
#!/usr/bin/env python3
"""
Fake GISAS data according to sample model from gisas_fit2d, with added noise.
"""
import bornagain as ba
from bornagain import deg, angstrom, nm
import numpy as np
import gisas_fit2d as model
def fake_data():
"""
Generate fake "experimental" data, and save them as numpy array.
"""
params = {
'cylinder_height': 5*nm,
'cylinder_radius': 5*nm,
'prism_height': 5*nm,
'prism_base_edge': 5*nm
}
# Compute model distribution
simulation = model.get_simulation(params)
simulation.runSimulation()
data = simulation.result().array()
# Draw noisy data
np.random.seed(0)
noise_factor = 0.1
noisy = np.random.normal(data, noise_factor*np.sqrt(data))
noisy[noisy < 0.1] = 0.1
# Save to numpy
np.savetxt("gisas-model1.txt.gz", data)
if __name__ == '__main__':
fake_data()
......@@ -5,7 +5,7 @@ of substrate.
"""
import bornagain as ba
from bornagain import deg, angstrom, nm
from bornagain import deg, nm
import numpy as np
from matplotlib import pyplot as plt
......@@ -47,47 +47,19 @@ def get_simulation(params):
"""
Returns a GISAXS simulation with beam and detector defined
"""
simulation = ba.GISASSimulation()
simulation.setDetectorParameters(100, -1*deg, 1*deg, 100, 0, 2*deg)
simulation.setBeamParameters(1*angstrom, 0.2*deg, 0)
simulation.beam().setIntensity(1e+08)
simulation.setSample(get_sample(params))
return simulation
intensity = 1e8
beam = ba.Beam(intensity, 0.1*nm, ba.Direction(0.2*deg, 0))
det = ba.SphericalDetector(100, -1*deg, 1*deg, 100, 0, 2*deg)
sample = get_sample(params)
def create_real_data():
"""
Generating "experimental" data by running simulation with certain parameters.
The data is saved on disk in the form of numpy array.
"""
# default sample parameters
params = {
'cylinder_height': 5*nm,
'cylinder_radius': 5*nm,
'prism_height': 5*nm,
'prism_base_edge': 5*nm
}
# retrieving simulated data in the form of numpy array
simulation = get_simulation(params)
simulation.runSimulation()
real_data = simulation.result().array()
# spoiling simulated data with noise to produce "real" data
np.random.seed(0)
noise_factor = 0.1
noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
noisy[noisy < 0.1] = 0.1
np.savetxt("basic_fitting_tutorial_data.txt.gz", real_data)
return ba.GISASSimulation(beam, sample, det)
def load_real_data():
def load_data():
"""
Loads experimental data from file
"""
return np.loadtxt("basic_fitting_tutorial_data.txt.gz", dtype=float)
return np.loadtxt("gisas-model1.txt.gz", dtype=float)
def run_fitting():
......@@ -95,7 +67,7 @@ def run_fitting():
Setup simulation and fit
"""
real_data = load_real_data()
real_data = load_data()
fit_objective = ba.FitObjective()
fit_objective.addSimulationAndData(get_simulation, real_data)
......@@ -122,8 +94,8 @@ def run_fitting():
for fitPar in result.parameters():
print(fitPar.name(), fitPar.value, fitPar.error)
# saving simulation image corresponding to the best fit parameters
# np.savetxt("data.txt", fit_objective.simulationResult().array())
# Save simulation image corresponding to the best fit parameters
np.savetxt("fit.txt", fit_objective.simulationResult().array())
if __name__ == '__main__':
......
......@@ -13,7 +13,7 @@
# ************************************************************************** #
import bornagain as ba
import ba_plot
import ba_plot as bp
try: # workaround for build servers
import numpy as np
from matplotlib import pyplot as plt
......@@ -85,10 +85,10 @@ class PlotterGISAS(Plotter):
zmax = np.amax(arr) if self._zmax is None else self._zmax
zmin = zmax*1e-6 if self._zmin is None else self._zmin
ba.plot_colormap(real_data,
bp.plot_colormap(real_data,
title="Experimental data",
zmin=zmin,
zmax=zmax,
intensity_min=zmin,
intensity_max=zmax,
units=self._units,
xlabel=self._xlabel,
ylabel=self._ylabel,
......@@ -96,10 +96,10 @@ class PlotterGISAS(Plotter):
aspect=self._aspect)
self.make_subplot(2)
ba.plot_colormap(sim_data,
bp.plot_colormap(sim_data,
title="Simulated data",
zmin=zmin,
zmax=zmax,
intensity_min=zmin,
intensity_max=zmax,
units=self._units,
xlabel=self._xlabel,
ylabel=self._ylabel,
......@@ -107,10 +107,10 @@ class PlotterGISAS(Plotter):
aspect=self._aspect)
self.make_subplot(3)
ba.plot_colormap(diff,
bp.plot_colormap(diff,
title="Difference",
zmin=zmin,
zmax=zmax,
intensity_min=zmin,
intensity_max=zmax,
units=self._units,
xlabel=self._xlabel,
ylabel=self._ylabel,
......
......@@ -242,6 +242,11 @@ class ObserverCallbackWrapper(PyObserverCallback):
def finalize(self, minimizer_result):
return self.finalize_cpp(self.convert_result(minimizer_result))
def create_default_plotter(self):
import ba_fitmonitor as bafim
self.m_plotter = bafim.PlotterGISAS()
return self.m_plotter.plot
def initPlot(self, every_nth, callback = None):
if not callback:
callback = self.create_default_plotter()
......
.\" Automatically generated by Pod::Man 4.07 (Pod::Simple 3.32)
.\" Automatically generated by Pod::Man 4.14 (Pod::Simple 3.40)
.\"
.\" Standard preamble:
.\" ========================================================================
......@@ -54,16 +54,20 @@
.\" Avoid warning from groff about undefined register 'F'.
.de IX
..
.if !\nF .nr F 0
.if \nF>0 \{\
. de IX
. tm Index:\\$1\t\\n%\t"\\$2"
.nr rF 0
.if \n(.g .if rF .nr rF 1
.if (\n(rF:(\n(.g==0)) \{\
. if \nF \{\
. de IX
. tm Index:\\$1\t\\n%\t"\\$2"
..
. if !\nF==2 \{\
. nr % 0
. nr F 2
. if !\nF==2 \{\
. nr % 0
. nr F 2
. \}
. \}
.\}
.rr rF
.\"
.\" Accent mark definitions (@(#)ms.acc 1.5 88/02/08 SMI; from UCB 4.2).
.\" Fear. Run. Save yourself. No user-serviceable parts.
......@@ -129,7 +133,7 @@
.\" ========================================================================
.\"
.IX Title "BORNAGAIN 1"
.TH BORNAGAIN 1 "2017-09-04" "perl v5.24.1" "BornAgain manual"
.TH BORNAGAIN 1 "2021-03-09" "perl v5.32.1" "BornAgain manual"
.\" For nroff, turn off justification. Always turn off hyphenation; it makes
.\" way too many mistakes in technical documents.
.if n .ad l
......
......@@ -416,7 +416,7 @@ Controlled by the multi-threading machinery in ISimulation::runSingleSimulation(
C++ includes: DWBAComputation.h
";
%feature("docstring") DWBAComputation::DWBAComputation "DWBAComputation::DWBAComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, std::vector< SimulationElement >::iterator begin_it, std::vector< SimulationElement >::iterator end_it)
%feature("docstring") DWBAComputation::DWBAComputation "DWBAComputation::DWBAComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, std::vector< SimulationElement >::iterator begin_it, std::vector< SimulationElement >::iterator end_it, bool forcePolarized=false)
";
%feature("docstring") DWBAComputation::~DWBAComputation "DWBAComputation::~DWBAComputation() override
......@@ -896,7 +896,7 @@ Controlled by the multi-threading machinery in ISimulation::runSingleSimulation(
C++ includes: IComputation.h
";
%feature("docstring") IComputation::IComputation "IComputation::IComputation(const MultiLayer &sample, const SimulationOptions &options, ProgressHandler &progress)
%feature("docstring") IComputation::IComputation "IComputation::IComputation(const MultiLayer &sample, const SimulationOptions &options, ProgressHandler &progress, bool forcePolarized=false)
";
%feature("docstring") IComputation::~IComputation "IComputation::~IComputation()
......@@ -1034,11 +1034,6 @@ Run a simulation in a MPI environment.
%feature("docstring") ISimulation::setDetectorResolutionFunction "void ISimulation::setDetectorResolutionFunction(const IResolutionFunction2D &resolution_function)
";
%feature("docstring") ISimulation::setAnalyzerProperties "void ISimulation::setAnalyzerProperties(const kvector_t direction, double efficiency, double total_transmission)
Sets the polarization analyzer characteristics of the detector.
";
%feature("docstring") ISimulation::setSample "void ISimulation::setSample(const MultiLayer &sample)
The MultiLayer object will not be owned by the ISimulation object.
......@@ -2007,7 +2002,7 @@ Controlled by the multi-threading machinery in ISimulation::runSingleSimulation(
C++ includes: SpecularComputation.h
";
%feature("docstring") SpecularComputation::SpecularComputation "SpecularComputation::SpecularComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, SpecularElementIter begin_it, SpecularElementIter end_it)
%feature("docstring") SpecularComputation::SpecularComputation "SpecularComputation::SpecularComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, SpecularElementIter begin_it, SpecularElementIter end_it, bool forcePolarized=false)
";
%feature("docstring") SpecularComputation::~SpecularComputation "SpecularComputation::~SpecularComputation() override
......
......@@ -2561,6 +2561,9 @@ C++ includes: SpecularDetector1D.h
%feature("docstring") SpecularDetector1D::SpecularDetector1D "SpecularDetector1D::SpecularDetector1D(const IAxis &axis)
";
%feature("docstring") SpecularDetector1D::SpecularDetector1D "SpecularDetector1D::SpecularDetector1D()
";
%feature("docstring") SpecularDetector1D::~SpecularDetector1D "SpecularDetector1D::~SpecularDetector1D()
";
......@@ -2585,6 +2588,9 @@ Returns region of interest if exists.
Resets region of interest making whole detector plane available for the simulation.
";
%feature("docstring") SpecularDetector1D::setAxis "void SpecularDetector1D::setAxis(const IAxis &axis)
";
%feature("docstring") SpecularDetector1D::defaultAxesUnits "Axes::Units SpecularDetector1D::defaultAxesUnits() const override
Return default axes units.
......
......@@ -4685,6 +4685,9 @@ C++ includes: ISpecularStrategy.h
%feature("docstring") ISpecularStrategy::Execute "virtual coeffs_t ISpecularStrategy::Execute(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const =0
";
%feature("docstring") ISpecularStrategy::computeTopLayerR "virtual std::variant<complex_t, Eigen::Matrix2cd> ISpecularStrategy::computeTopLayerR(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const =0
";
// File: classLargeCylindersInDWBABuilder.xml
%feature("docstring") LargeCylindersInDWBABuilder "
......@@ -6403,7 +6406,7 @@ If the usage of average materials is requested, layers and particles are sliced
C++ includes: ProcessedSample.h
";
%feature("docstring") ProcessedSample::ProcessedSample "ProcessedSample::ProcessedSample(const MultiLayer &sample, const SimulationOptions &options)
%feature("docstring") ProcessedSample::ProcessedSample "ProcessedSample::ProcessedSample(const MultiLayer &sample, const SimulationOptions &options, bool forcePolarized=false)
";
%feature("docstring") ProcessedSample::~ProcessedSample "ProcessedSample::~ProcessedSample()
......@@ -7169,7 +7172,7 @@ C++ includes: Slice.h
%feature("docstring") Slice::setMaterial "void Slice::setMaterial(const Material &material)
";
%feature("docstring") Slice::material "Material Slice::material() const
%feature("docstring") Slice::material "const Material & Slice::material() const
";
%feature("docstring") Slice::thickness "double Slice::thickness() const
......@@ -7297,6 +7300,11 @@ Computes refraction angle reflection/transmission coefficients for given sliced
Computes refraction angle reflection/transmission coefficients for given sliced multilayer and a set of kz projections corresponding to each slice
";
%feature("docstring") SpecularMagneticStrategy::computeTopLayerR "std::variant< complex_t, Eigen::Matrix2cd > SpecularMagneticStrategy::computeTopLayerR(const std::vector< Slice > &slices, const std::vector< complex_t > &kzs) const override
Computes the Fresnel R coefficient for the top layer only Introduced in order to speed up pure reflectivity computations
";
// File: classSpecularMagneticStrategy__v1.xml
%feature("docstring") SpecularMagneticStrategy_v1 "
......@@ -7314,6 +7322,9 @@ Computes refraction angle reflection/transmission coefficients for given sliced
%feature("docstring") SpecularMagneticStrategy_v1::Execute "ISpecularStrategy::coeffs_t SpecularMagneticStrategy_v1::Execute(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const
";
%feature("docstring") SpecularMagneticStrategy_v1::computeTopLayerR "std::variant< complex_t, Eigen::Matrix2cd > SpecularMagneticStrategy_v1::computeTopLayerR(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const override
";
// File: classSpecularMagneticStrategy__v2.xml
%feature("docstring") SpecularMagneticStrategy_v2 "
......@@ -7335,6 +7346,9 @@ Computes refraction angle reflection/transmission coefficients for given sliced
Computes refraction angle reflection/transmission coefficients for given sliced multilayer and a set of kz projections corresponding to each slice
";
%feature("docstring") SpecularMagneticStrategy_v2::computeTopLayerR "std::variant< complex_t, Eigen::Matrix2cd > SpecularMagneticStrategy_v2::computeTopLayerR(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const override
";
// File: classSpecularMagneticTanhStrategy.xml
%feature("docstring") SpecularMagneticTanhStrategy "
......@@ -7378,6 +7392,11 @@ Computes refraction angles and transmission/reflection coefficients for given co
%feature("docstring") SpecularScalarStrategy::Execute "ISpecularStrategy::coeffs_t SpecularScalarStrategy::Execute(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const override
";
%feature("docstring") SpecularScalarStrategy::computeTopLayerR "std::variant< complex_t, Eigen::Matrix2cd > SpecularScalarStrategy::computeTopLayerR(const std::vector< Slice > &slices, const std::vector< complex_t > &kz) const override
Computes the Fresnel R coefficient for the top layer only Introduced in order to speed up pure reflectivity computations
";
// File: classSpecularScalarTanhStrategy.xml
%feature("docstring") SpecularScalarTanhStrategy "
......
......@@ -3045,6 +3045,11 @@ class FitObjective(object):
def finalize(self, minimizer_result):
return self.finalize_cpp(self.convert_result(minimizer_result))
def create_default_plotter(self):
import ba_fitmonitor as bafim
self.m_plotter = bafim.PlotterGISAS()
return self.m_plotter.plot
def initPlot(self, every_nth, callback = None):
if not callback:
callback = self.create_default_plotter()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment