diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index c91f2de2b3b3157b7baccd32953651baaabcfcdd..3d20364a3936d32cd98966eccd8605de20f47e4d 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -18,6 +18,7 @@
#include "Base/Axis/Scale.h"
#include "Base/Util/Assert.h"
#include <algorithm>
+#include <random>
#ifdef BORNAGAIN_PYTHON
@@ -26,6 +27,9 @@
namespace {
+auto const seed = 123;
+auto urbg = std::mt19937(seed);
+
PyObject* npExport(const Frame& frame, const std::vector<double>& flatData)
{
if (flatData.empty())
@@ -282,6 +286,11 @@ Datafield* Datafield::crop(double xmin, double xmax) const
#ifdef BORNAGAIN_PYTHON
+PyObject* Datafield::npXcenters() const
+{
+ return ::npExport(Frame(xAxis().clone()), xAxis().binCenters());
+}
+
PyObject* Datafield::npArray() const
{
return ::npExport(frame(), flatVector());
@@ -378,3 +387,17 @@ Datafield Datafield::flat() const
{
return {title(), frame().flat(), m_values, m_errSigmas};
}
+
+Datafield Datafield::noisy(double prefactor, double minimum) const
+{
+ std::vector<double> outval(size());
+ std::vector<double> errval(size());
+ for (size_t i = 0; i < size(); ++i) {
+ double mean = valAt(i);
+ double stdv = prefactor * sqrt(mean);
+ auto norm = std::normal_distribution<double>(mean, stdv);
+ outval[i] = std::max(norm(urbg), minimum);
+ errval[i] = stdv;
+ }
+ return {title(), frame().clone(), outval, errval};
+}
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index 1fbfb4675d347f2fb494b0fb9e5bf2cc6b512f25..14c1b766854743839475303613ba0a8c98107d1c 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -64,6 +64,7 @@ public:
double valAt(size_t i) const;
#ifdef BORNAGAIN_PYTHON
//! Returns data as Python numpy array.
+ PyObject* npXcenters() const;
PyObject* npArray() const;
PyObject* npErrors() const;
#endif
@@ -91,6 +92,7 @@ public:
Datafield plottableField(std::vector<std::string> labels = {}) const;
Datafield plottableField(std::string label) const;
Datafield flat() const;
+ Datafield noisy(double prefactor, double minimum) const;
//! Multiplies contents by constant factor, e.g. to shift curves in log plot.
void scale(double factor);
@@ -133,6 +135,8 @@ public:
//! Sets content of output data to specific value.
void setAllTo(const double& value); // Py: test only; rm when possible
+ const std::vector<double>& errorSigmas() const;
+
#ifndef SWIG
double& operator[](size_t i);
const double& operator[](size_t i) const;
@@ -141,7 +145,6 @@ public:
void setVector(const std::vector<double>& data_vector);
bool hasErrorSigmas() const;
- const std::vector<double>& errorSigmas() const;
std::vector<double>& errorSigmas();
private:
diff --git a/Wrap/Python/ba_fitmonitor.py b/Wrap/Python/ba_fitmonitor.py
index b199be1f11ca4f314e0b7db7bef9abd26b12d233..bd9c1c01de606fb9b7eaea875aa9e4b330dc0ebd 100644
--- a/Wrap/Python/ba_fitmonitor.py
+++ b/Wrap/Python/ba_fitmonitor.py
@@ -73,18 +73,18 @@ class PlotterGISAS(Plotter):
def plot(self, fit_objective):
Plotter.reset(self)
- real_data = fit_objective.experimentalData()
+ data = fit_objective.experimentalData()
sim_data = fit_objective.simulationResult()
diff = fit_objective.absoluteDifference()
self.make_subplot(1)
# same limits for both plots
- arr = real_data.npArray()
+ arr = data.npArray()
zmax = np.amax(arr) if self._zmax is None else self._zmax
zmin = zmax*1e-6 if self._zmin is None else self._zmin
- bp.plot_simres(real_data,
+ bp.plot_simres(data,
title="Experimental data",
intensity_min=zmin,
intensity_max=zmax,
diff --git a/auto/Examples/bayesian/likelihood_sampling.py b/auto/Examples/bayesian/likelihood_sampling.py
index e463127a6225d8a720744d9acba0e855758ca393..49ccddbce6b207bafbd3529d7dbd86b16532e1da 100755
--- a/auto/Examples/bayesian/likelihood_sampling.py
+++ b/auto/Examples/bayesian/likelihood_sampling.py
@@ -57,9 +57,9 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
- real_data = ba.readData2D(filepath).npArray()
- real_data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
- real_data[:, 2] = real_data[:, 1]*0.1 # arbitrary uncertainties of 10%
+ data = ba.readData2D(filepath).npArray()
+ data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
+ data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
def log_likelihood(P):
"""
@@ -70,9 +70,9 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
- q = real_data[:, 0]
- y = real_data[:, 1]
- dy = real_data[:, 2]
+ q = data[:, 0]
+ y = data[:, 1]
+ dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
@@ -104,14 +104,14 @@ if __name__ == '__main__':
plt.show()
# Plot and show MLE and data of reflectivity
- plt.errorbar(real_data[:, 0],
- real_data[:, 1],
- real_data[:, 2],
+ plt.errorbar(data[:, 0],
+ data[:, 1],
+ data[:, 2],
marker='.',
ls='')
plt.plot(
- real_data[:, 0],
- run_simulation(real_data[:, 0], *flat_samples.mean(axis=0)),
+ data[:, 0],
+ run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
diff --git a/auto/Examples/fit/scatter2d/consecutive_fitting.py b/auto/Examples/fit/scatter2d/consecutive_fitting.py
index 37b1dadb706e3dd018d5b05a6763ffe2f53fa8da..d53f842e29f67681b4f243d61672239c0d8022ab 100755
--- a/auto/Examples/fit/scatter2d/consecutive_fitting.py
+++ b/auto/Examples/fit/scatter2d/consecutive_fitting.py
@@ -50,7 +50,7 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, get_sample(P), detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -59,24 +59,17 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
def run_fitting():
"""
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
fit_objective.initPlot(10, observer)
diff --git a/auto/Examples/fit/scatter2d/custom_objective_function.py b/auto/Examples/fit/scatter2d/custom_objective_function.py
index bde5d2abe244199b63d66f25a9e8a5dd98c4730c..00d61fc9e3547e448e9c237556e7a8ba9df9feec 100755
--- a/auto/Examples/fit/scatter2d/custom_objective_function.py
+++ b/auto/Examples/fit/scatter2d/custom_objective_function.py
@@ -32,10 +32,10 @@ class MyObjective(ba.FitObjective):
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
objective = MyObjective()
- objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ objective.addFitPair(model.get_simulation, data, 1)
objective.initPrint(10)
P = ba.Parameters()
diff --git a/auto/Examples/fit/scatter2d/expfit_galaxi.py b/auto/Examples/fit/scatter2d/expfit_galaxi.py
index e5abf8669b3db63e4c27eec3c1449b912bde47a3..7de6f98915a85ef2e92d1dc9f32ac1fda5cd1773 100755
--- a/auto/Examples/fit/scatter2d/expfit_galaxi.py
+++ b/auto/Examples/fit/scatter2d/expfit_galaxi.py
@@ -122,19 +122,19 @@ def create_simulation(P):
return simulation
-def load_real_data(filename):
+def load_data(filename):
"""
Loads experimental data and returns numpy array.
"""
filepath = os.path.join(datadir, filename)
- return ba.readData2D(filepath).npArray()
+ return ba.readData2D(filepath)
def run_fitting():
- real_data = load_real_data("scatter2d/galaxi_data.tif.gz")
+ data = load_data("scatter2d/galaxi_data.tif.gz")
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(create_simulation, real_data, 1)
+ fit_objective.addFitPair(create_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
fit_objective.initPlot(10, observer)
diff --git a/auto/Examples/fit/scatter2d/find_background.py b/auto/Examples/fit/scatter2d/find_background.py
index 2edcbe882b60b5f99326b34ef7cc1e0c92efa69a..72565a617972f409f5a47aa5d92294e6018949cc 100755
--- a/auto/Examples/fit/scatter2d/find_background.py
+++ b/auto/Examples/fit/scatter2d/find_background.py
@@ -24,7 +24,7 @@ def get_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise, background and scale
to the simulated data.
@@ -42,18 +42,17 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- return result.npArray()
+ return result
if __name__ == '__main__':
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
- observer = ba_fitmonitor.PlotterGISAS()
+ observer = ba_fitmonitor.PlotterGISAS()
fit_objective.initPlot(10, observer)
P = ba.Parameters()
diff --git a/auto/Examples/fit/scatter2d/fit2d.py b/auto/Examples/fit/scatter2d/fit2d.py
index 72adfaf4706f775e604ed187d08bcf8fb104245d..00a7a3a4e5036df41ad8ac33d9de0e26d6876cd9 100755
--- a/auto/Examples/fit/scatter2d/fit2d.py
+++ b/auto/Examples/fit/scatter2d/fit2d.py
@@ -28,18 +28,19 @@ def get_simulation(P):
return simulation
-def real_data():
+def fake_data():
"""
Generating "experimental" data by running simulation with default parameters.
"""
simulation = get_simulation({'radius': 5 * nm})
result = simulation.simulate()
- return result.npArray()
+ return result
if __name__ == '__main__':
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data())
+ fit_objective.addFitPair(get_simulation, data)
fit_objective.initPrint(10)
P = ba.Parameters()
diff --git a/auto/Examples/fit/scatter2d/fit_along_slices.py b/auto/Examples/fit/scatter2d/fit_along_slices.py
index 6d90bb0b86f5286b2d425a238b3f0b051ec599fc..086302dcc8dc61ccc8824c0bde0e8916da4c3511 100755
--- a/auto/Examples/fit/scatter2d/fit_along_slices.py
+++ b/auto/Examples/fit/scatter2d/fit_along_slices.py
@@ -28,19 +28,18 @@ def get_masked_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
# initial values which we will have to find later during the fit
P = {'radius': 5*nm, 'height': 10*nm}
- # retrieving simulated data in the form of numpy array
simulation = model.get_simulation(P)
simulation.setBackground(ba.PoissonBackground())
result = simulation.simulate()
- return result.npArray()
+ return result.noisy(0.1, 0.1)
class PlotObserver:
@@ -57,7 +56,7 @@ class PlotObserver:
self.update(fit_objective)
@staticmethod
- def plot_real_data(data):
+ def plot_data(data):
"""
Plot experimental data as heatmap with horizontal/vertical lines
representing slices on top.
@@ -121,15 +120,15 @@ class PlotObserver:
"""
self.fig.clf()
- real_data = fit_objective.experimentalData()
+ data = fit_objective.experimentalData()
simul_data = fit_objective.simulationResult()
# plot real data
plt.subplot(2, 2, 1)
- self.plot_real_data(real_data)
+ self.plot_data(data)
# horizontal slices
- slices = [("real", real_data.xProjection(alpha_slice_value)),
+ slices = [("real", data.xProjection(alpha_slice_value)),
("simul", simul_data.xProjection(alpha_slice_value))]
title = ("Horizontal slice at alpha =" +
'{:3.1f}'.format(alpha_slice_value))
@@ -137,7 +136,7 @@ class PlotObserver:
self.plot_slices(slices, title)
# vertical slices
- slices = [("real", real_data.yProjection(phi_slice_value)),
+ slices = [("real", data.yProjection(phi_slice_value)),
("simul", simul_data.yProjection(phi_slice_value))]
title = "Vertical slice at phi =" + '{:3.1f}'.format(
phi_slice_value)
@@ -154,13 +153,13 @@ def run_fitting():
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_masked_simulation, real_data, 1)
+ fit_objective.addFitPair(get_masked_simulation, data, 1)
fit_objective.initPrint(10)
- # creating custom observer which will draw fit progress
+ # creating custom observer which will draw fit progress
plotter = PlotObserver()
fit_objective.initPlot(10, plotter)
diff --git a/auto/Examples/fit/scatter2d/fit_gisas.py b/auto/Examples/fit/scatter2d/fit_gisas.py
index 584c6338af4a0e51be9a7bdaf38fdea48443b815..253788863da5f4e72cdc0a39b0153d1607684548 100755
--- a/auto/Examples/fit/scatter2d/fit_gisas.py
+++ b/auto/Examples/fit/scatter2d/fit_gisas.py
@@ -13,13 +13,13 @@ from bornagain import ba_fitmonitor
from matplotlib import pyplot as plt
-def run_fitting():
+if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
- real_data = np.loadtxt(os.path.join(datadir, "scatter2d/faked_gisas1.dat.gz"),
- dtype=float)
+ fname = os.path.join(datadir, "scatter2d/faked_gisas1.dat.gz")
+ data = ba.readData2D(fname, ba.csv2D)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data)
+ fit_objective.addFitPair(model.get_simulation, data)
fit_objective.initPrint(10) # Print on every 10th iteration.
observer = ba_fitmonitor.PlotterGISAS()
@@ -39,8 +39,4 @@ def run_fitting():
# Save simulation image corresponding to the best fit parameters
np.savetxt("fit.txt", fit_objective.simulationResult().npArray())
-
-
-if __name__ == '__main__':
- run_fitting()
plt.show()
diff --git a/auto/Examples/fit/scatter2d/fit_with_masks.py b/auto/Examples/fit/scatter2d/fit_with_masks.py
index fdd070d5d9907926537aefe42868d70667c73e8d..95acd099683df91c9522b7f599ffa81379c0ef46 100755
--- a/auto/Examples/fit/scatter2d/fit_with_masks.py
+++ b/auto/Examples/fit/scatter2d/fit_with_masks.py
@@ -56,10 +56,10 @@ def add_mask_to_simulation(simulation):
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data().noisy(0.1, 0.1)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_masked_simulation, real_data, 1)
+ fit_objective.addFitPair(get_masked_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
fit_objective.initPlot(10, observer)
diff --git a/auto/Examples/fit/scatter2d/lmfit_basics.py b/auto/Examples/fit/scatter2d/lmfit_basics.py
index 4216709ddd27c811ebdea40cefaaecb5fcd0d4a4..ff4a7afe4978a3bfbb572452626fe2b78fa970ca 100755
--- a/auto/Examples/fit/scatter2d/lmfit_basics.py
+++ b/auto/Examples/fit/scatter2d/lmfit_basics.py
@@ -9,10 +9,10 @@ import model2_hexlattice as model
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ fit_objective.addFitPair(model.get_simulation, data, 1)
fit_objective.initPrint(10)
P = lmfit.Parameters()
diff --git a/auto/Examples/fit/scatter2d/lmfit_with_plotting.py b/auto/Examples/fit/scatter2d/lmfit_with_plotting.py
index 88e081b70a3810e45f9b9e06b60a8abc72dd81c8..4b2785ee77f4c8d231a8725b5f79fd5e4c9eba46 100755
--- a/auto/Examples/fit/scatter2d/lmfit_with_plotting.py
+++ b/auto/Examples/fit/scatter2d/lmfit_with_plotting.py
@@ -25,16 +25,16 @@ class LMFITPlotter:
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ fit_objective.addFitPair(model.get_simulation, data, 1)
fit_objective.initPrint(10)
P = lmfit.Parameters()
P.add('radius', value=7*nm, min=5*nm, max=8*nm)
P.add('length', value=10*nm, min=8*nm, max=14*nm)
-
+
plotter = LMFITPlotter(fit_objective)
result = lmfit.minimize(fit_objective.evaluate_residuals,
P, iter_cb=plotter)
diff --git a/auto/Examples/fit/scatter2d/minimizer_settings.py b/auto/Examples/fit/scatter2d/minimizer_settings.py
index a95e54328e9348d5c172e9effb52e6f6a19e478f..8d8e31f9d4f007c40c53ff3891f6e280e6014bfd 100755
--- a/auto/Examples/fit/scatter2d/minimizer_settings.py
+++ b/auto/Examples/fit/scatter2d/minimizer_settings.py
@@ -51,7 +51,7 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, get_sample(P), detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data from simulated image with default parameters.
"""
@@ -66,7 +66,7 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- return result.npArray()
+ return result
def run_fitting():
@@ -74,7 +74,7 @@ def run_fitting():
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
# prints info about available minimizers
print(ba.MinimizerFactory().catalogToString())
@@ -83,7 +83,7 @@ def run_fitting():
print(ba.MinimizerFactory().catalogDetailsToString())
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
P = ba.Parameters()
diff --git a/auto/Examples/fit/scatter2d/model1_cylinders.py b/auto/Examples/fit/scatter2d/model1_cylinders.py
index 789cc95a130fe664746def1931139b93480952e9..bcc2435f299c0e98ea231b61c495afb0bfa200e9 100755
--- a/auto/Examples/fit/scatter2d/model1_cylinders.py
+++ b/auto/Examples/fit/scatter2d/model1_cylinders.py
@@ -42,19 +42,13 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, sample, detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
P = {'radius': 5*nm, 'height': 10*nm}
- # retrieving simulated data in the form of numpy array
simulation = get_simulation(P)
result = simulation.simulate()
- real_data = result.npArray()
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
diff --git a/auto/Examples/fit/scatter2d/model2_hexlattice.py b/auto/Examples/fit/scatter2d/model2_hexlattice.py
index 46e5d9b8f260b31e6b5905c0bde6fc4acd410bf8..5cf760658db0df82d90cd1e57c24b6bfba6980e8 100755
--- a/auto/Examples/fit/scatter2d/model2_hexlattice.py
+++ b/auto/Examples/fit/scatter2d/model2_hexlattice.py
@@ -48,7 +48,7 @@ def get_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -56,12 +56,4 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # 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
- return noisy
+ return result.noisy(0.1, 0.1)
diff --git a/auto/Examples/fit/scatter2d/multiple_datasets.py b/auto/Examples/fit/scatter2d/multiple_datasets.py
index 24752b532f5dfd3294c01931765aa7bb510d2b3a..9def2912788efeb28804f15137e2bba9e2ecccdd 100755
--- a/auto/Examples/fit/scatter2d/multiple_datasets.py
+++ b/auto/Examples/fit/scatter2d/multiple_datasets.py
@@ -59,7 +59,7 @@ def simulation2(P):
return get_simulation(P)
-def create_real_data(incident_alpha):
+def fake_data(incident_alpha):
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -73,14 +73,7 @@ def create_real_data(incident_alpha):
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
class PlotObserver():
@@ -99,14 +92,14 @@ class PlotObserver():
@staticmethod
def plot_dataset(fit_objective, canvas):
for i_dataset in range(0, fit_objective.fitObjectCount()):
- real_data = fit_objective.experimentalData(i_dataset)
+ data = fit_objective.experimentalData(i_dataset)
simul_data = fit_objective.simulationResult(i_dataset)
chi2_map = fit_objective.relativeDifference(i_dataset)
- zmax = real_data.maxVal()
+ zmax = data.maxVal()
plt.subplot(canvas[i_dataset*3])
- bp.plot_simres(real_data,
+ bp.plot_simres(data,
title="\"Real\" data - #" +
str(i_dataset + 1),
intensity_min=1,
@@ -189,12 +182,12 @@ def run_fitting():
main function to run fitting
"""
- data1 = create_real_data(0.1*deg)
- data2 = create_real_data(0.4*deg)
+ data1 = fake_data(0.1*deg)
+ data2 = fake_data(0.4*deg)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(simulation1, data1, 1)
- fit_objective.addSimulationAndData(simulation2, data2, 1)
+ fit_objective.addFitPair(simulation1, data1, 1)
+ fit_objective.addFitPair(simulation2, data2, 1)
fit_objective.initPrint(10)
# creating custom observer which will draw fit progress
diff --git a/auto/Examples/fit/specular/Honeycomb_fit.py b/auto/Examples/fit/specular/Honeycomb_fit.py
index c382637f01ec51dcfce62a751af2b11e57ae3687..57cccad8437f825b970aa13f75349acf441bc84b 100755
--- a/auto/Examples/fit/specular/Honeycomb_fit.py
+++ b/auto/Examples/fit/specular/Honeycomb_fit.py
@@ -88,7 +88,7 @@ def run_simulation(qaxis, P, *, sign, T):
# Experimental data
####################################################################
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
# Load from disk.
filepath = os.path.join(datadir, filename)
with open(filepath, 'r') as f:
@@ -217,11 +217,11 @@ if __name__ == '__main__':
qmin = 0.08
qmax = 1.4
- expData = [
- get_experimental_data("specular/honeycomb300p.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb300m.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb150p.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb150m.dat", qmin, qmax)]
+ data = [
+ load_data("specular/honeycomb300p.dat", qmin, qmax),
+ load_data("specular/honeycomb300m.dat", qmin, qmax),
+ load_data("specular/honeycomb150p.dat", qmin, qmax),
+ load_data("specular/honeycomb150m.dat", qmin, qmax)]
simFunctions = [
lambda q, P: run_simulation(q, P, sign=+1, T=300),
@@ -234,14 +234,14 @@ if __name__ == '__main__':
# Plot data with initial model
simResults = [ f(qzs, P) for f in simFunctions ]
- plot(qzs, simResults, expData, [1, 1, 10, 10],
+ plot(qzs, simResults, data, [1, 1, 10, 10],
["300K $+$", "300K $-$", "150K $+$", "150K $-$"])
plot_sld_profile(P)
# Fit
- qaxes = [d[0] for d in expData]
- rdata = [d[1] for d in expData]
+ qaxes = [d[0] for d in data]
+ rdata = [d[1] for d in data]
def par_dict(*args):
return {name: value for name, value in zip(freeParNames, *args)} | fixedP
@@ -280,7 +280,7 @@ if __name__ == '__main__':
P = par_dict(result.x)
sim_results = [ f(qzs, P) for f in simFunctions ]
- plot(qzs, simResults, expData, [1, 1, 10, 10],
+ plot(qzs, simResults, data, [1, 1, 10, 10],
["300K $+$", "300K $-$", "150K $+$", "150K $-$"])
plot_sld_profile(P)
diff --git a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
index d2b77fe517a4130d7c86e4234c660cd9c7cd48dc..02b8a45b568af0cf16bc6f1dcd36d4f4525cc756 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
ax = fig.add_subplot(111)
for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
- ax.errorbar(exp[0],
- exp[1],
- xerr=exp[3],
- yerr=exp[2],
+ ax.errorbar(exp.xAxis().binCenters(),
+ exp.flatVector(),
+ # xerr=TODO i742,
+ yerr=exp.errorSigmas(),
fmt='.',
markersize=0.75,
linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
Plot the simulated spin asymmetry as well its
experimental counterpart with errorbars
"""
+ Yp = data_pp.npArray()
+ Ym = data_mm.npArray()
+ Ep = data_pp.npErrors()
+ Em = data_mm.npErrors()
# compute the errorbars of the spin asymmetry
- delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
- data_mm[1]**2 * data_pp[2]**2 ) /
- ( data_pp[1] + data_mm[1] )**4 )
+ delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
fig = plt.figure()
ax = fig.add_subplot(111)
- ax.errorbar(data_pp[0],
- (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
- xerr=data_pp[3],
+ ax.errorbar(data_pp.xAxis().binCenters(),
+ (Yp - Ym) / (Yp + Ym),
+ # xerr=TODO i742,
yerr=delta,
fmt='.',
markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
# Data Handling #
####################################################################
-def load_exp(fname):
- dat = numpy.loadtxt(fname)
- return numpy.transpose(dat)
+def load_data(fname):
+ flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+ return ba.readData1D(fname, ba.csv1D, flags)
####################################################################
# Main Function #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = load_exp(fname_stem + "pp.tab")
- expdata_mm = load_exp(fname_stem + "mm.tab")
+ expdata_pp = load_data(fname_stem + "pp.tab")
+ expdata_mm = load_data(fname_stem + "mm.tab")
fixedP = {
# parameters from our own fit run
diff --git a/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
index 3dc5a42df7dcb11db5a3a5cd8c5831fb15394cdc..5d7cff180c619151042766a1a0fe5233872108c0 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -32,36 +32,6 @@ def get_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
return simulation
-####################################################################
-# Fit Function #
-####################################################################
-
-def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
- startParams):
-
- fit_objective = ba.FitObjective()
- fit_objective.setObjectiveMetric("chi2")
-
- fit_objective.addSimulationAndData(
- lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
- r_uncertainty[0], 1)
- fit_objective.addSimulationAndData(
- lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
- r_uncertainty[1], 1)
-
- fit_objective.initPrint(10)
-
- P = ba.Parameters()
- for name, p in startParams.items():
- P.add(name, p[0], min=p[1], max=p[2])
-
- minimizer = ba.Minimizer()
-
- result = minimizer.minimize(fit_objective.evaluate, P)
- fit_objective.finalize(result)
-
- return {r.name(): r.value for r in result.parameters()}
-
####################################################################
# Main Function #
####################################################################
@@ -71,8 +41,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = psa.load_exp(fname_stem + "pp.tab")
- expdata_mm = psa.load_exp(fname_stem + "mm.tab")
+ data_pp = psa.load_data(fname_stem + "pp.tab")
+ data_mm = psa.load_data(fname_stem + "mm.tab")
fixedParams = {
# parameters can be moved here to keep them fixed
@@ -114,26 +84,43 @@ if __name__ == '__main__':
legend_pp = "$++$"
legend_mm = "$--$"
- psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+ psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
[legend_pp, legend_mm], [color_pp, color_mm])
- psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+ psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
+
+ fit_objective = ba.FitObjective()
+ fit_objective.setObjectiveMetric("chi2")
+ fit_objective.initPrint(10)
+
+ fit_objective.addFitPair(
+ lambda P: get_Simulation_pp(data_pp.npXcenters(), P),
+ data_pp, 1)
+ fit_objective.addFitPair(
+ lambda P: get_Simulation_mm(data_mm.npXcenters(), P),
+ data_mm, 1)
+
+ P = ba.Parameters()
+ for name, p in startParams.items():
+ P.add(name, p[0], min=p[1], max=p[2])
+
+ minimizer = ba.Minimizer()
+
+ result = minimizer.minimize(fit_objective.evaluate, P)
+ fit_objective.finalize(result)
+
+ fitResult = {r.name(): r.value for r in result.parameters()}
- fitResult = run_fit_ba([expdata_pp[0], expdata_mm[0]],
- [expdata_pp[1], expdata_mm[1]],
- [expdata_pp[2], expdata_mm[2]],
- [get_Simulation_pp, get_Simulation_mm],
- startParams)
print("Fit Result:")
print(fitResult)
q_pp, r_pp = psa.qr(get_Simulation_pp(qzs, fitResult).simulate())
q_mm, r_mm = psa.qr(get_Simulation_mm(qzs, fitResult).simulate())
- psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+ psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
[legend_pp, legend_mm], [color_pp, color_mm])
plt.draw()
- psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+ psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
plt.draw()
plt.show()
diff --git a/auto/Examples/fit/specular/Pt_layer_fit.py b/auto/Examples/fit/specular/Pt_layer_fit.py
index ccbcd4f5139ea49c1a60f8cdeb28941558ec974c..a9b82577a8b3665631eb793b3af5daecdedb1e41 100755
--- a/auto/Examples/fit/specular/Pt_layer_fit.py
+++ b/auto/Examples/fit/specular/Pt_layer_fit.py
@@ -60,7 +60,7 @@ def get_simulation(q_axis, P):
# Experimental data
####################################################################
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
filepath = os.path.join(datadir, filename)
data = np.genfromtxt(filepath, unpack=True)
@@ -144,7 +144,7 @@ if __name__ == '__main__':
qmin = 0.18
qmax = 2.4
- data = get_experimental_data("specular/RvsQ_36563_36662.dat.gz", qmin, qmax)
+ data = load_data("specular/RvsQ_36563_36662.dat.gz", qmin, qmax)
qzs = np.linspace(qmin, qmax, 1500)
diff --git a/auto/Examples/specular/PolarizedSpinAsymmetry.py b/auto/Examples/specular/PolarizedSpinAsymmetry.py
index d2b77fe517a4130d7c86e4234c660cd9c7cd48dc..02b8a45b568af0cf16bc6f1dcd36d4f4525cc756 100755
--- a/auto/Examples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
ax = fig.add_subplot(111)
for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
- ax.errorbar(exp[0],
- exp[1],
- xerr=exp[3],
- yerr=exp[2],
+ ax.errorbar(exp.xAxis().binCenters(),
+ exp.flatVector(),
+ # xerr=TODO i742,
+ yerr=exp.errorSigmas(),
fmt='.',
markersize=0.75,
linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
Plot the simulated spin asymmetry as well its
experimental counterpart with errorbars
"""
+ Yp = data_pp.npArray()
+ Ym = data_mm.npArray()
+ Ep = data_pp.npErrors()
+ Em = data_mm.npErrors()
# compute the errorbars of the spin asymmetry
- delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
- data_mm[1]**2 * data_pp[2]**2 ) /
- ( data_pp[1] + data_mm[1] )**4 )
+ delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
fig = plt.figure()
ax = fig.add_subplot(111)
- ax.errorbar(data_pp[0],
- (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
- xerr=data_pp[3],
+ ax.errorbar(data_pp.xAxis().binCenters(),
+ (Yp - Ym) / (Yp + Ym),
+ # xerr=TODO i742,
yerr=delta,
fmt='.',
markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
# Data Handling #
####################################################################
-def load_exp(fname):
- dat = numpy.loadtxt(fname)
- return numpy.transpose(dat)
+def load_data(fname):
+ flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+ return ba.readData1D(fname, ba.csv1D, flags)
####################################################################
# Main Function #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = load_exp(fname_stem + "pp.tab")
- expdata_mm = load_exp(fname_stem + "mm.tab")
+ expdata_pp = load_data(fname_stem + "pp.tab")
+ expdata_mm = load_data(fname_stem + "mm.tab")
fixedP = {
# parameters from our own fit run
diff --git a/auto/Examples/specular/VsGenx.py b/auto/Examples/specular/VsGenx.py
index 54486ebd1ecd512c0653fc35241f674515dce341..c44d0caebfcd3d9a3876b42468adc0bff67df124 100755
--- a/auto/Examples/specular/VsGenx.py
+++ b/auto/Examples/specular/VsGenx.py
@@ -16,7 +16,7 @@ def reference_data(filename):
"""
Loads and returns reference data from GenX simulation
"""
- ax_values, real_data = np.loadtxt(filename,
+ ax_values, data = np.loadtxt(filename,
usecols=(0, 1),
skiprows=3,
unpack=True)
@@ -24,7 +24,7 @@ def reference_data(filename):
# translate axis values from double incident angle to incident angle
ax_values *= 0.5
- return ax_values, real_data
+ return ax_values, data
def get_sample():
diff --git a/auto/MiniExamples/bayesian/likelihood_sampling.py b/auto/MiniExamples/bayesian/likelihood_sampling.py
index aeed815d4919212490203a26060a9699e485ad4e..9144431229378469ed89ad12dba6bc3a0a0165be 100755
--- a/auto/MiniExamples/bayesian/likelihood_sampling.py
+++ b/auto/MiniExamples/bayesian/likelihood_sampling.py
@@ -57,9 +57,9 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
- real_data = ba.readData2D(filepath).npArray()
- real_data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
- real_data[:, 2] = real_data[:, 1]*0.1 # arbitrary uncertainties of 10%
+ data = ba.readData2D(filepath).npArray()
+ data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
+ data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
def log_likelihood(P):
"""
@@ -70,9 +70,9 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
- q = real_data[:, 0]
- y = real_data[:, 1]
- dy = real_data[:, 2]
+ q = data[:, 0]
+ y = data[:, 1]
+ dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
@@ -104,14 +104,14 @@ if __name__ == '__main__':
# plt.show()
# Plot and show MLE and data of reflectivity
- plt.errorbar(real_data[:, 0],
- real_data[:, 1],
- real_data[:, 2],
+ plt.errorbar(data[:, 0],
+ data[:, 1],
+ data[:, 2],
marker='.',
ls='')
plt.plot(
- real_data[:, 0],
- run_simulation(real_data[:, 0], *flat_samples.mean(axis=0)),
+ data[:, 0],
+ run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
diff --git a/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py b/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py
index 939026e4e37538a98f0a7f32461ed01de6a41dbe..97c1adeb80325d4831912fde77a917ce39b2a13c 100755
--- a/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py
+++ b/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py
@@ -50,7 +50,7 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, get_sample(P), detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -59,24 +59,17 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
def run_fitting():
"""
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
plt.close() # (hide plot) fit_objective.initPlot(10, observer)
diff --git a/auto/MiniExamples/fit/scatter2d/custom_objective_function.py b/auto/MiniExamples/fit/scatter2d/custom_objective_function.py
index bde5d2abe244199b63d66f25a9e8a5dd98c4730c..00d61fc9e3547e448e9c237556e7a8ba9df9feec 100755
--- a/auto/MiniExamples/fit/scatter2d/custom_objective_function.py
+++ b/auto/MiniExamples/fit/scatter2d/custom_objective_function.py
@@ -32,10 +32,10 @@ class MyObjective(ba.FitObjective):
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
objective = MyObjective()
- objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ objective.addFitPair(model.get_simulation, data, 1)
objective.initPrint(10)
P = ba.Parameters()
diff --git a/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py b/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py
index 05a23a7696fab9198dfb6ba1736b626bde7e613f..b1f7d293c923b9b51a5835b56c081aa7d45af70c 100755
--- a/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py
+++ b/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py
@@ -122,19 +122,19 @@ def create_simulation(P):
return simulation
-def load_real_data(filename):
+def load_data(filename):
"""
Loads experimental data and returns numpy array.
"""
filepath = os.path.join(datadir, filename)
- return ba.readData2D(filepath).npArray()
+ return ba.readData2D(filepath)
def run_fitting():
- real_data = load_real_data("scatter2d/galaxi_data.tif.gz")
+ data = load_data("scatter2d/galaxi_data.tif.gz")
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(create_simulation, real_data, 1)
+ fit_objective.addFitPair(create_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
plt.close() # (hide plot) fit_objective.initPlot(10, observer)
diff --git a/auto/MiniExamples/fit/scatter2d/find_background.py b/auto/MiniExamples/fit/scatter2d/find_background.py
index 8d68baa690bb4a2037c8bce768ce3c3e6de5e6e7..adf22b4dee316aa098ba68960ba4c46e933f702a 100755
--- a/auto/MiniExamples/fit/scatter2d/find_background.py
+++ b/auto/MiniExamples/fit/scatter2d/find_background.py
@@ -24,7 +24,7 @@ def get_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise, background and scale
to the simulated data.
@@ -42,18 +42,17 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- return result.npArray()
+ return result
if __name__ == '__main__':
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
- observer = ba_fitmonitor.PlotterGISAS()
+ observer = ba_fitmonitor.PlotterGISAS()
plt.close() # (hide plot) fit_objective.initPlot(10, observer)
P = ba.Parameters()
diff --git a/auto/MiniExamples/fit/scatter2d/fit2d.py b/auto/MiniExamples/fit/scatter2d/fit2d.py
index e3c6ffed781197f9d7adf442807dd4ac993eaf9c..f937ac35dac9072cb63ee327297b537febbafb07 100755
--- a/auto/MiniExamples/fit/scatter2d/fit2d.py
+++ b/auto/MiniExamples/fit/scatter2d/fit2d.py
@@ -28,18 +28,19 @@ def get_simulation(P):
return simulation
-def real_data():
+def fake_data():
"""
Generating "experimental" data by running simulation with default parameters.
"""
simulation = get_simulation({'radius': 5 * nm})
result = simulation.simulate()
- return result.npArray()
+ return result
if __name__ == '__main__':
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data())
+ fit_objective.addFitPair(get_simulation, data)
fit_objective.initPrint(10)
P = ba.Parameters()
diff --git a/auto/MiniExamples/fit/scatter2d/fit_along_slices.py b/auto/MiniExamples/fit/scatter2d/fit_along_slices.py
index c4d3b3fd4573a6baade358e61182b87115e61bec..9383f16397d0ed7efda33f804e38dfd3de4778d0 100755
--- a/auto/MiniExamples/fit/scatter2d/fit_along_slices.py
+++ b/auto/MiniExamples/fit/scatter2d/fit_along_slices.py
@@ -28,19 +28,18 @@ def get_masked_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
# initial values which we will have to find later during the fit
P = {'radius': 5*nm, 'height': 10*nm}
- # retrieving simulated data in the form of numpy array
simulation = model.get_simulation(P)
simulation.setBackground(ba.PoissonBackground())
result = simulation.simulate()
- return result.npArray()
+ return result.noisy(0.1, 0.1)
class PlotObserver:
@@ -57,7 +56,7 @@ class PlotObserver:
self.update(fit_objective)
@staticmethod
- def plot_real_data(data):
+ def plot_data(data):
"""
Plot experimental data as heatmap with horizontal/vertical lines
representing slices on top.
@@ -121,15 +120,15 @@ class PlotObserver:
"""
self.fig.clf()
- real_data = fit_objective.experimentalData()
+ data = fit_objective.experimentalData()
simul_data = fit_objective.simulationResult()
# plot real data
plt.subplot(2, 2, 1)
- self.plot_real_data(real_data)
+ self.plot_data(data)
# horizontal slices
- slices = [("real", real_data.xProjection(alpha_slice_value)),
+ slices = [("real", data.xProjection(alpha_slice_value)),
("simul", simul_data.xProjection(alpha_slice_value))]
title = ("Horizontal slice at alpha =" +
'{:3.1f}'.format(alpha_slice_value))
@@ -137,7 +136,7 @@ class PlotObserver:
self.plot_slices(slices, title)
# vertical slices
- slices = [("real", real_data.yProjection(phi_slice_value)),
+ slices = [("real", data.yProjection(phi_slice_value)),
("simul", simul_data.yProjection(phi_slice_value))]
title = "Vertical slice at phi =" + '{:3.1f}'.format(
phi_slice_value)
@@ -154,13 +153,13 @@ def run_fitting():
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_masked_simulation, real_data, 1)
+ fit_objective.addFitPair(get_masked_simulation, data, 1)
fit_objective.initPrint(10)
- # creating custom observer which will draw fit progress
+ # creating custom observer which will draw fit progress
plotter = PlotObserver()
plt.close() # (hide plot) fit_objective.initPlot(10, plotter)
diff --git a/auto/MiniExamples/fit/scatter2d/fit_gisas.py b/auto/MiniExamples/fit/scatter2d/fit_gisas.py
index 3073a07e3aba00eb638ac5c107c029bfc4d5b807..d890ab2ed4265c2229f3d1c9ec7e5264353c6e0e 100755
--- a/auto/MiniExamples/fit/scatter2d/fit_gisas.py
+++ b/auto/MiniExamples/fit/scatter2d/fit_gisas.py
@@ -13,13 +13,13 @@ from bornagain import ba_fitmonitor
from matplotlib import pyplot as plt
-def run_fitting():
+if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
- real_data = np.loadtxt(os.path.join(datadir, "scatter2d/faked_gisas1.dat.gz"),
- dtype=float)
+ fname = os.path.join(datadir, "scatter2d/faked_gisas1.dat.gz")
+ data = ba.readData2D(fname, ba.csv2D)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data)
+ fit_objective.addFitPair(model.get_simulation, data)
fit_objective.initPrint(10) # Print on every 10th iteration.
observer = ba_fitmonitor.PlotterGISAS()
@@ -39,8 +39,4 @@ def run_fitting():
# Save simulation image corresponding to the best fit parameters
np.savetxt("fit.txt", fit_objective.simulationResult().npArray())
-
-
-if __name__ == '__main__':
- run_fitting()
plt.show()
diff --git a/auto/MiniExamples/fit/scatter2d/fit_with_masks.py b/auto/MiniExamples/fit/scatter2d/fit_with_masks.py
index 38fe262743b6284447cc442307638c9b1784932f..ca1130b852b4c877b4fafa3c355a0562a895cbd6 100755
--- a/auto/MiniExamples/fit/scatter2d/fit_with_masks.py
+++ b/auto/MiniExamples/fit/scatter2d/fit_with_masks.py
@@ -56,10 +56,10 @@ def add_mask_to_simulation(simulation):
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data().noisy(0.1, 0.1)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_masked_simulation, real_data, 1)
+ fit_objective.addFitPair(get_masked_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
plt.close() # (hide plot) fit_objective.initPlot(10, observer)
diff --git a/auto/MiniExamples/fit/scatter2d/lmfit_basics.py b/auto/MiniExamples/fit/scatter2d/lmfit_basics.py
index 4216709ddd27c811ebdea40cefaaecb5fcd0d4a4..ff4a7afe4978a3bfbb572452626fe2b78fa970ca 100755
--- a/auto/MiniExamples/fit/scatter2d/lmfit_basics.py
+++ b/auto/MiniExamples/fit/scatter2d/lmfit_basics.py
@@ -9,10 +9,10 @@ import model2_hexlattice as model
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ fit_objective.addFitPair(model.get_simulation, data, 1)
fit_objective.initPrint(10)
P = lmfit.Parameters()
diff --git a/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py b/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py
index 4c5903831bd633954fc55543d883026a21cd7b99..643fa35722e31a7fe6fd8dbe278b5b820c9aecfe 100755
--- a/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py
+++ b/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py
@@ -25,16 +25,16 @@ class LMFITPlotter:
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ fit_objective.addFitPair(model.get_simulation, data, 1)
fit_objective.initPrint(10)
P = lmfit.Parameters()
P.add('radius', value=7*nm, min=5*nm, max=8*nm)
P.add('length', value=10*nm, min=8*nm, max=14*nm)
-
+
result = lmfit.minimize(fit_objective.evaluate_residuals,
P)
diff --git a/auto/MiniExamples/fit/scatter2d/minimizer_settings.py b/auto/MiniExamples/fit/scatter2d/minimizer_settings.py
index e554661c94c00916bb2bd44f6c64be987d92542b..2093950d49840863a714e4cd1d05ff3b4b173f33 100755
--- a/auto/MiniExamples/fit/scatter2d/minimizer_settings.py
+++ b/auto/MiniExamples/fit/scatter2d/minimizer_settings.py
@@ -51,7 +51,7 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, get_sample(P), detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data from simulated image with default parameters.
"""
@@ -66,7 +66,7 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- return result.npArray()
+ return result
def run_fitting():
@@ -74,7 +74,7 @@ def run_fitting():
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
# prints info about available minimizers
print(ba.MinimizerFactory().catalogToString())
@@ -83,7 +83,7 @@ def run_fitting():
print(ba.MinimizerFactory().catalogDetailsToString())
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
P = ba.Parameters()
diff --git a/auto/MiniExamples/fit/scatter2d/model1_cylinders.py b/auto/MiniExamples/fit/scatter2d/model1_cylinders.py
index 4ebb3f6a961f374c4aa81c04e004cbe12d0d684c..b2bbf592a80e5762a1d962fc029bd13e82b3e19b 100755
--- a/auto/MiniExamples/fit/scatter2d/model1_cylinders.py
+++ b/auto/MiniExamples/fit/scatter2d/model1_cylinders.py
@@ -42,19 +42,13 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, sample, detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
P = {'radius': 5*nm, 'height': 10*nm}
- # retrieving simulated data in the form of numpy array
simulation = get_simulation(P)
result = simulation.simulate()
- real_data = result.npArray()
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
diff --git a/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py b/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
index eb30ce743834cf3ffc7e8396792bec6f31cc1bf8..e65c5f0bd9782637c6d6afca225d38915545e4f5 100755
--- a/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
+++ b/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
@@ -48,7 +48,7 @@ def get_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -56,12 +56,4 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # 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
- return noisy
+ return result.noisy(0.1, 0.1)
diff --git a/auto/MiniExamples/fit/scatter2d/multiple_datasets.py b/auto/MiniExamples/fit/scatter2d/multiple_datasets.py
index 896ce868dcc3a786eab5cd9d19e1b53c2687a592..aa1f4d20cc8d656edc88af848566790ed5e8bfed 100755
--- a/auto/MiniExamples/fit/scatter2d/multiple_datasets.py
+++ b/auto/MiniExamples/fit/scatter2d/multiple_datasets.py
@@ -59,7 +59,7 @@ def simulation2(P):
return get_simulation(P)
-def create_real_data(incident_alpha):
+def fake_data(incident_alpha):
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -73,14 +73,7 @@ def create_real_data(incident_alpha):
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
class PlotObserver():
@@ -99,14 +92,14 @@ class PlotObserver():
@staticmethod
def plot_dataset(fit_objective, canvas):
for i_dataset in range(0, fit_objective.fitObjectCount()):
- real_data = fit_objective.experimentalData(i_dataset)
+ data = fit_objective.experimentalData(i_dataset)
simul_data = fit_objective.simulationResult(i_dataset)
chi2_map = fit_objective.relativeDifference(i_dataset)
- zmax = real_data.maxVal()
+ zmax = data.maxVal()
plt.subplot(canvas[i_dataset*3])
- bp.plot_simres(real_data,
+ bp.plot_simres(data,
title="\"Real\" data - #" +
str(i_dataset + 1),
intensity_min=1,
@@ -189,12 +182,12 @@ def run_fitting():
main function to run fitting
"""
- data1 = create_real_data(0.1*deg)
- data2 = create_real_data(0.4*deg)
+ data1 = fake_data(0.1*deg)
+ data2 = fake_data(0.4*deg)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(simulation1, data1, 1)
- fit_objective.addSimulationAndData(simulation2, data2, 1)
+ fit_objective.addFitPair(simulation1, data1, 1)
+ fit_objective.addFitPair(simulation2, data2, 1)
fit_objective.initPrint(10)
# creating custom observer which will draw fit progress
diff --git a/auto/MiniExamples/fit/specular/Honeycomb_fit.py b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
index 11022ca19f68f7d2cdb56e6fa310d05a45567a46..e258f2b1599275ebafb344850783a0f0cd5373b6 100755
--- a/auto/MiniExamples/fit/specular/Honeycomb_fit.py
+++ b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
@@ -88,7 +88,7 @@ def run_simulation(qaxis, P, *, sign, T):
# Experimental data
####################################################################
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
# Load from disk.
filepath = os.path.join(datadir, filename)
with open(filepath, 'r') as f:
@@ -217,11 +217,11 @@ if __name__ == '__main__':
qmin = 0.08
qmax = 1.4
- expData = [
- get_experimental_data("specular/honeycomb300p.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb300m.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb150p.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb150m.dat", qmin, qmax)]
+ data = [
+ load_data("specular/honeycomb300p.dat", qmin, qmax),
+ load_data("specular/honeycomb300m.dat", qmin, qmax),
+ load_data("specular/honeycomb150p.dat", qmin, qmax),
+ load_data("specular/honeycomb150m.dat", qmin, qmax)]
simFunctions = [
lambda q, P: run_simulation(q, P, sign=+1, T=300),
@@ -234,14 +234,14 @@ if __name__ == '__main__':
# Plot data with initial model
simResults = [ f(qzs, P) for f in simFunctions ]
- plot(qzs, simResults, expData, [1, 1, 10, 10],
+ plot(qzs, simResults, data, [1, 1, 10, 10],
["300K $+$", "300K $-$", "150K $+$", "150K $-$"])
plot_sld_profile(P)
# Fit
- qaxes = [d[0] for d in expData]
- rdata = [d[1] for d in expData]
+ qaxes = [d[0] for d in data]
+ rdata = [d[1] for d in data]
def par_dict(*args):
return {name: value for name, value in zip(freeParNames, *args)} | fixedP
@@ -280,7 +280,7 @@ if __name__ == '__main__':
P = par_dict(result.x)
sim_results = [ f(qzs, P) for f in simFunctions ]
- plot(qzs, simResults, expData, [1, 1, 10, 10],
+ plot(qzs, simResults, data, [1, 1, 10, 10],
["300K $+$", "300K $-$", "150K $+$", "150K $-$"])
plot_sld_profile(P)
diff --git a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
index 37510dfc5bcc6471feee1f03dd63f50a1ab87967..d41fc7bd268e3cf251266026f335da70a08a64f2 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
ax = fig.add_subplot(111)
for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
- ax.errorbar(exp[0],
- exp[1],
- xerr=exp[3],
- yerr=exp[2],
+ ax.errorbar(exp.xAxis().binCenters(),
+ exp.flatVector(),
+ # xerr=TODO i742,
+ yerr=exp.errorSigmas(),
fmt='.',
markersize=0.75,
linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
Plot the simulated spin asymmetry as well its
experimental counterpart with errorbars
"""
+ Yp = data_pp.npArray()
+ Ym = data_mm.npArray()
+ Ep = data_pp.npErrors()
+ Em = data_mm.npErrors()
# compute the errorbars of the spin asymmetry
- delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
- data_mm[1]**2 * data_pp[2]**2 ) /
- ( data_pp[1] + data_mm[1] )**4 )
+ delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
fig = plt.figure()
ax = fig.add_subplot(111)
- ax.errorbar(data_pp[0],
- (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
- xerr=data_pp[3],
+ ax.errorbar(data_pp.xAxis().binCenters(),
+ (Yp - Ym) / (Yp + Ym),
+ # xerr=TODO i742,
yerr=delta,
fmt='.',
markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
# Data Handling #
####################################################################
-def load_exp(fname):
- dat = numpy.loadtxt(fname)
- return numpy.transpose(dat)
+def load_data(fname):
+ flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+ return ba.readData1D(fname, ba.csv1D, flags)
####################################################################
# Main Function #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = load_exp(fname_stem + "pp.tab")
- expdata_mm = load_exp(fname_stem + "mm.tab")
+ expdata_pp = load_data(fname_stem + "pp.tab")
+ expdata_mm = load_data(fname_stem + "mm.tab")
fixedP = {
# parameters from our own fit run
diff --git a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
index 0e75dd5f27032a0afaaab8969de02a27be6707fb..6c3a68be5210699ff3056b4610c3a7b12f31ad95 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -32,36 +32,6 @@ def get_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
return simulation
-####################################################################
-# Fit Function #
-####################################################################
-
-def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
- startParams):
-
- fit_objective = ba.FitObjective()
- fit_objective.setObjectiveMetric("chi2")
-
- fit_objective.addSimulationAndData(
- lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
- r_uncertainty[0], 1)
- fit_objective.addSimulationAndData(
- lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
- r_uncertainty[1], 1)
-
- fit_objective.initPrint(10)
-
- P = ba.Parameters()
- for name, p in startParams.items():
- P.add(name, p[0], min=p[1], max=p[2])
-
- minimizer = ba.Minimizer()
-
- result = minimizer.minimize(fit_objective.evaluate, P)
- fit_objective.finalize(result)
-
- return {r.name(): r.value for r in result.parameters()}
-
####################################################################
# Main Function #
####################################################################
@@ -71,8 +41,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = psa.load_exp(fname_stem + "pp.tab")
- expdata_mm = psa.load_exp(fname_stem + "mm.tab")
+ data_pp = psa.load_data(fname_stem + "pp.tab")
+ data_mm = psa.load_data(fname_stem + "mm.tab")
fixedParams = {
# parameters can be moved here to keep them fixed
@@ -114,26 +84,43 @@ if __name__ == '__main__':
legend_pp = "$++$"
legend_mm = "$--$"
- psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+ psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
[legend_pp, legend_mm], [color_pp, color_mm])
- psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+ psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
+
+ fit_objective = ba.FitObjective()
+ fit_objective.setObjectiveMetric("chi2")
+ fit_objective.initPrint(10)
+
+ fit_objective.addFitPair(
+ lambda P: get_Simulation_pp(data_pp.npXcenters(), P),
+ data_pp, 1)
+ fit_objective.addFitPair(
+ lambda P: get_Simulation_mm(data_mm.npXcenters(), P),
+ data_mm, 1)
+
+ P = ba.Parameters()
+ for name, p in startParams.items():
+ P.add(name, p[0], min=p[1], max=p[2])
+
+ minimizer = ba.Minimizer()
+
+ result = minimizer.minimize(fit_objective.evaluate, P)
+ fit_objective.finalize(result)
+
+ fitResult = {r.name(): r.value for r in result.parameters()}
- fitResult = run_fit_ba([expdata_pp[0], expdata_mm[0]],
- [expdata_pp[1], expdata_mm[1]],
- [expdata_pp[2], expdata_mm[2]],
- [get_Simulation_pp, get_Simulation_mm],
- startParams)
print("Fit Result:")
print(fitResult)
q_pp, r_pp = psa.qr(get_Simulation_pp(qzs, fitResult).simulate())
q_mm, r_mm = psa.qr(get_Simulation_mm(qzs, fitResult).simulate())
- psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+ psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
[legend_pp, legend_mm], [color_pp, color_mm])
#plt.draw()
- psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+ psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
#plt.draw()
#plt.show()
diff --git a/auto/MiniExamples/fit/specular/Pt_layer_fit.py b/auto/MiniExamples/fit/specular/Pt_layer_fit.py
index 8e0063aa1963cd7a2b8ddaaa58bc8d85b6452f23..b64d78e7e0813cd56a865fef7c1bfc5316ee46e1 100755
--- a/auto/MiniExamples/fit/specular/Pt_layer_fit.py
+++ b/auto/MiniExamples/fit/specular/Pt_layer_fit.py
@@ -60,7 +60,7 @@ def get_simulation(q_axis, P):
# Experimental data
####################################################################
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
filepath = os.path.join(datadir, filename)
data = np.genfromtxt(filepath, unpack=True)
@@ -144,7 +144,7 @@ if __name__ == '__main__':
qmin = 0.18
qmax = 2.4
- data = get_experimental_data("specular/RvsQ_36563_36662.dat.gz", qmin, qmax)
+ data = load_data("specular/RvsQ_36563_36662.dat.gz", qmin, qmax)
qzs = np.linspace(qmin, qmax, 1500)
diff --git a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
index 37510dfc5bcc6471feee1f03dd63f50a1ab87967..d41fc7bd268e3cf251266026f335da70a08a64f2 100755
--- a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
ax = fig.add_subplot(111)
for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
- ax.errorbar(exp[0],
- exp[1],
- xerr=exp[3],
- yerr=exp[2],
+ ax.errorbar(exp.xAxis().binCenters(),
+ exp.flatVector(),
+ # xerr=TODO i742,
+ yerr=exp.errorSigmas(),
fmt='.',
markersize=0.75,
linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
Plot the simulated spin asymmetry as well its
experimental counterpart with errorbars
"""
+ Yp = data_pp.npArray()
+ Ym = data_mm.npArray()
+ Ep = data_pp.npErrors()
+ Em = data_mm.npErrors()
# compute the errorbars of the spin asymmetry
- delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
- data_mm[1]**2 * data_pp[2]**2 ) /
- ( data_pp[1] + data_mm[1] )**4 )
+ delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
fig = plt.figure()
ax = fig.add_subplot(111)
- ax.errorbar(data_pp[0],
- (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
- xerr=data_pp[3],
+ ax.errorbar(data_pp.xAxis().binCenters(),
+ (Yp - Ym) / (Yp + Ym),
+ # xerr=TODO i742,
yerr=delta,
fmt='.',
markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
# Data Handling #
####################################################################
-def load_exp(fname):
- dat = numpy.loadtxt(fname)
- return numpy.transpose(dat)
+def load_data(fname):
+ flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+ return ba.readData1D(fname, ba.csv1D, flags)
####################################################################
# Main Function #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = load_exp(fname_stem + "pp.tab")
- expdata_mm = load_exp(fname_stem + "mm.tab")
+ expdata_pp = load_data(fname_stem + "pp.tab")
+ expdata_mm = load_data(fname_stem + "mm.tab")
fixedP = {
# parameters from our own fit run
diff --git a/auto/MiniExamples/specular/VsGenx.py b/auto/MiniExamples/specular/VsGenx.py
index e5b91dbb2c4f269e5eb6352a980d30651e308438..d0a990a43ab4e7b483208819f1fa9bed0a6f6cb8 100755
--- a/auto/MiniExamples/specular/VsGenx.py
+++ b/auto/MiniExamples/specular/VsGenx.py
@@ -16,7 +16,7 @@ def reference_data(filename):
"""
Loads and returns reference data from GenX simulation
"""
- ax_values, real_data = np.loadtxt(filename,
+ ax_values, data = np.loadtxt(filename,
usecols=(0, 1),
skiprows=3,
unpack=True)
@@ -24,7 +24,7 @@ def reference_data(filename):
# translate axis values from double incident angle to incident angle
ax_values *= 0.5
- return ax_values, real_data
+ return ax_values, data
def get_sample():
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index f4ee7229eda8294eab7be599dbb159d49a16e127..accb8d77c7a4857481ed01acc8264d5c5402293b 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2089,6 +2089,10 @@ class Datafield(object):
r"""valAt(Datafield self, size_t i) -> double"""
return _libBornAgainDevice.Datafield_valAt(self, i)
+ def npXcenters(self):
+ r"""npXcenters(Datafield self) -> PyObject *"""
+ return _libBornAgainDevice.Datafield_npXcenters(self)
+
def npArray(self):
r"""npArray(Datafield self) -> PyObject *"""
return _libBornAgainDevice.Datafield_npArray(self)
@@ -2148,6 +2152,10 @@ class Datafield(object):
r"""flat(Datafield self) -> Datafield"""
return _libBornAgainDevice.Datafield_flat(self)
+ def noisy(self, prefactor, minimum):
+ r"""noisy(Datafield self, double prefactor, double minimum) -> Datafield"""
+ return _libBornAgainDevice.Datafield_noisy(self, prefactor, minimum)
+
def scale(self, factor):
r"""scale(Datafield self, double factor)"""
return _libBornAgainDevice.Datafield_scale(self, factor)
@@ -2179,6 +2187,10 @@ class Datafield(object):
r"""setAllTo(Datafield self, double const & value)"""
return _libBornAgainDevice.Datafield_setAllTo(self, value)
+ def errorSigmas(self):
+ r"""errorSigmas(Datafield self) -> vdouble1d_t"""
+ return _libBornAgainDevice.Datafield_errorSigmas(self)
+
# Register Datafield in _libBornAgainDevice:
_libBornAgainDevice.Datafield_swigregister(Datafield)
class Beam(libBornAgainBase.ICloneable, libBornAgainParam.INode):
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index b9fe28cfe47e09f0d61fdfc3606d54ceded8ae4e..f02b47528cd895e8ce255178a16da14ec7e9f448 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -29680,6 +29680,39 @@ fail:
}
+SWIGINTERN PyObject *_wrap_Datafield_npXcenters(PyObject *self, PyObject *args) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = (Datafield *) 0 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ PyObject *swig_obj[1] ;
+ PyObject *result = 0 ;
+
+ if (!args) SWIG_fail;
+ swig_obj[0] = args;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 | 0 );
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_npXcenters" "', argument " "1"" of type '" "Datafield const *""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ {
+ try {
+ result = (PyObject *)((Datafield const *)arg1)->npXcenters();
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = result;
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
SWIGINTERN PyObject *_wrap_Datafield_npArray(PyObject *self, PyObject *args) {
PyObject *resultobj = 0;
Datafield *arg1 = (Datafield *) 0 ;
@@ -30283,6 +30316,54 @@ fail:
}
+SWIGINTERN PyObject *_wrap_Datafield_noisy(PyObject *self, PyObject *args) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = (Datafield *) 0 ;
+ double arg2 ;
+ double arg3 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ double val2 ;
+ int ecode2 = 0 ;
+ double val3 ;
+ int ecode3 = 0 ;
+ PyObject *swig_obj[3] ;
+ SwigValueWrapper< Datafield > result;
+
+ if (!SWIG_Python_UnpackTuple(args, "Datafield_noisy", 3, 3, swig_obj)) SWIG_fail;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 | 0 );
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_noisy" "', argument " "1"" of type '" "Datafield const *""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
+ if (!SWIG_IsOK(ecode2)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "Datafield_noisy" "', argument " "2"" of type '" "double""'");
+ }
+ arg2 = static_cast< double >(val2);
+ ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+ if (!SWIG_IsOK(ecode3)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "Datafield_noisy" "', argument " "3"" of type '" "double""'");
+ }
+ arg3 = static_cast< double >(val3);
+ {
+ try {
+ result = ((Datafield const *)arg1)->noisy(arg2,arg3);
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = SWIG_NewPointerObj((new Datafield(result)), SWIGTYPE_p_Datafield, SWIG_POINTER_OWN | 0 );
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
SWIGINTERN PyObject *_wrap_Datafield_scale(PyObject *self, PyObject *args) {
PyObject *resultobj = 0;
Datafield *arg1 = (Datafield *) 0 ;
@@ -30907,6 +30988,39 @@ fail:
}
+SWIGINTERN PyObject *_wrap_Datafield_errorSigmas(PyObject *self, PyObject *args) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = (Datafield *) 0 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ PyObject *swig_obj[1] ;
+ std::vector< double,std::allocator< double > > *result = 0 ;
+
+ if (!args) SWIG_fail;
+ swig_obj[0] = args;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 | 0 );
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_errorSigmas" "', argument " "1"" of type '" "Datafield const *""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ {
+ try {
+ result = (std::vector< double,std::allocator< double > > *) &((Datafield const *)arg1)->errorSigmas();
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(*result));
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
SWIGINTERN PyObject *Datafield_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
PyObject *obj;
if (!SWIG_Python_UnpackTuple(args, "swigregister", 1, 1, &obj)) return NULL;
@@ -42711,6 +42825,7 @@ static PyMethodDef SwigMethods[] = {
{ "Datafield_title", _wrap_Datafield_title, METH_O, "Datafield_title(Datafield self) -> std::string"},
{ "Datafield_setAt", _wrap_Datafield_setAt, METH_VARARGS, "Datafield_setAt(Datafield self, size_t i, double val)"},
{ "Datafield_valAt", _wrap_Datafield_valAt, METH_VARARGS, "Datafield_valAt(Datafield self, size_t i) -> double"},
+ { "Datafield_npXcenters", _wrap_Datafield_npXcenters, METH_O, "Datafield_npXcenters(Datafield self) -> PyObject *"},
{ "Datafield_npArray", _wrap_Datafield_npArray, METH_O, "Datafield_npArray(Datafield self) -> PyObject *"},
{ "Datafield_npErrors", _wrap_Datafield_npErrors, METH_O, "Datafield_npErrors(Datafield self) -> PyObject *"},
{ "Datafield_frame", _wrap_Datafield_frame, METH_O, "Datafield_frame(Datafield self) -> Frame"},
@@ -42728,6 +42843,7 @@ static PyMethodDef SwigMethods[] = {
"Datafield_plottableField(Datafield self, std::string label) -> Datafield\n"
""},
{ "Datafield_flat", _wrap_Datafield_flat, METH_O, "Datafield_flat(Datafield self) -> Datafield"},
+ { "Datafield_noisy", _wrap_Datafield_noisy, METH_VARARGS, "Datafield_noisy(Datafield self, double prefactor, double minimum) -> Datafield"},
{ "Datafield_scale", _wrap_Datafield_scale, METH_VARARGS, "Datafield_scale(Datafield self, double factor)"},
{ "Datafield_crop", _wrap_Datafield_crop, METH_VARARGS, "\n"
"Datafield_crop(Datafield self, double xmin, double ymin, double xmax, double ymax) -> Datafield\n"
@@ -42744,6 +42860,7 @@ static PyMethodDef SwigMethods[] = {
"Datafield_yProjection(Datafield self, double xlow, double xup) -> Datafield\n"
""},
{ "Datafield_setAllTo", _wrap_Datafield_setAllTo, METH_VARARGS, "Datafield_setAllTo(Datafield self, double const & value)"},
+ { "Datafield_errorSigmas", _wrap_Datafield_errorSigmas, METH_O, "Datafield_errorSigmas(Datafield self) -> vdouble1d_t"},
{ "Datafield_swigregister", Datafield_swigregister, METH_O, NULL},
{ "Datafield_swiginit", Datafield_swiginit, METH_VARARGS, NULL},
{ "new_Beam", _wrap_new_Beam, METH_VARARGS, "Beam(double intensity, double wavelength, double alpha, double phi=0)"},
diff --git a/rawEx/bayesian/likelihood_sampling.py b/rawEx/bayesian/likelihood_sampling.py
index 51f7e23937257f46adf7e764ce93212660f31b66..ce7f70578fd2690e35a7b13fc586ef6ab8af96c9 100755
--- a/rawEx/bayesian/likelihood_sampling.py
+++ b/rawEx/bayesian/likelihood_sampling.py
@@ -57,9 +57,9 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
- real_data = ba.readData2D(filepath).npArray()
- real_data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
- real_data[:, 2] = real_data[:, 1]*0.1 # arbitrary uncertainties of 10%
+ data = ba.readData2D(filepath).npArray()
+ data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
+ data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
def log_likelihood(P):
"""
@@ -70,9 +70,9 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
- q = real_data[:, 0]
- y = real_data[:, 1]
- dy = real_data[:, 2]
+ q = data[:, 0]
+ y = data[:, 1]
+ dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
@@ -104,14 +104,14 @@ if __name__ == '__main__':
<%= sm ? "# " : "" %>plt.show()
# Plot and show MLE and data of reflectivity
- plt.errorbar(real_data[:, 0],
- real_data[:, 1],
- real_data[:, 2],
+ plt.errorbar(data[:, 0],
+ data[:, 1],
+ data[:, 2],
marker='.',
ls='')
plt.plot(
- real_data[:, 0],
- run_simulation(real_data[:, 0], *flat_samples.mean(axis=0)),
+ data[:, 0],
+ run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
diff --git a/rawEx/fit/scatter2d/consecutive_fitting.py b/rawEx/fit/scatter2d/consecutive_fitting.py
index ea875fa9ae82b54b3bc62c6fd0347180396370a6..74ab218391ffbc0bd9c7d7db107b4c5e670999d1 100755
--- a/rawEx/fit/scatter2d/consecutive_fitting.py
+++ b/rawEx/fit/scatter2d/consecutive_fitting.py
@@ -50,7 +50,7 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, get_sample(P), detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -59,24 +59,17 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
def run_fitting():
"""
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
<%= sm ? "plt.close() # (hide plot) " :"" %>fit_objective.initPlot(10, observer)
diff --git a/rawEx/fit/scatter2d/custom_objective_function.py b/rawEx/fit/scatter2d/custom_objective_function.py
index bde5d2abe244199b63d66f25a9e8a5dd98c4730c..00d61fc9e3547e448e9c237556e7a8ba9df9feec 100755
--- a/rawEx/fit/scatter2d/custom_objective_function.py
+++ b/rawEx/fit/scatter2d/custom_objective_function.py
@@ -32,10 +32,10 @@ class MyObjective(ba.FitObjective):
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
objective = MyObjective()
- objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ objective.addFitPair(model.get_simulation, data, 1)
objective.initPrint(10)
P = ba.Parameters()
diff --git a/rawEx/fit/scatter2d/expfit_galaxi.py b/rawEx/fit/scatter2d/expfit_galaxi.py
index 41eb75ca250b06aa9b9e929488fb6755ca905c21..6e39dead1c3d2fa07546010d1aee49b31160e1c8 100755
--- a/rawEx/fit/scatter2d/expfit_galaxi.py
+++ b/rawEx/fit/scatter2d/expfit_galaxi.py
@@ -122,19 +122,19 @@ def create_simulation(P):
return simulation
-def load_real_data(filename):
+def load_data(filename):
"""
Loads experimental data and returns numpy array.
"""
filepath = os.path.join(datadir, filename)
- return ba.readData2D(filepath).npArray()
+ return ba.readData2D(filepath)
def run_fitting():
- real_data = load_real_data("scatter2d/galaxi_data.tif.gz")
+ data = load_data("scatter2d/galaxi_data.tif.gz")
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(create_simulation, real_data, 1)
+ fit_objective.addFitPair(create_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
<%= sm ? "plt.close() # (hide plot) " :"" %>fit_objective.initPlot(10, observer)
diff --git a/rawEx/fit/scatter2d/find_background.py b/rawEx/fit/scatter2d/find_background.py
index f864ec53d2ded5354201c18ee9bcc9b1e7745cb2..8720eb52d8d6f95aa9efa421f2fdd582a82129fb 100755
--- a/rawEx/fit/scatter2d/find_background.py
+++ b/rawEx/fit/scatter2d/find_background.py
@@ -24,7 +24,7 @@ def get_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise, background and scale
to the simulated data.
@@ -42,18 +42,17 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- return result.npArray()
+ return result
if __name__ == '__main__':
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
- observer = ba_fitmonitor.PlotterGISAS()
+ observer = ba_fitmonitor.PlotterGISAS()
<%= sm ? "plt.close() # (hide plot) " :"" %>fit_objective.initPlot(10, observer)
P = ba.Parameters()
diff --git a/rawEx/fit/scatter2d/fit2d.py b/rawEx/fit/scatter2d/fit2d.py
index 08734f6cab3fe559c71062cf3e914d290a0cbbb4..205bac5e0dc06059744647ca8aa9d25e2a237b61 100755
--- a/rawEx/fit/scatter2d/fit2d.py
+++ b/rawEx/fit/scatter2d/fit2d.py
@@ -28,18 +28,19 @@ def get_simulation(P):
return simulation
-def real_data():
+def fake_data():
"""
Generating "experimental" data by running simulation with default parameters.
"""
simulation = get_simulation({'radius': 5 * nm})
result = simulation.simulate()
- return result.npArray()
+ return result
if __name__ == '__main__':
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data())
+ fit_objective.addFitPair(get_simulation, data)
fit_objective.initPrint(10)
P = ba.Parameters()
diff --git a/rawEx/fit/scatter2d/fit_along_slices.py b/rawEx/fit/scatter2d/fit_along_slices.py
index 3de58250cdef41e1d8960be336ea11ba194b642a..da6989c7adb1df19204ddd275500dfc98288d484 100755
--- a/rawEx/fit/scatter2d/fit_along_slices.py
+++ b/rawEx/fit/scatter2d/fit_along_slices.py
@@ -28,19 +28,18 @@ def get_masked_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
# initial values which we will have to find later during the fit
P = {'radius': 5*nm, 'height': 10*nm}
- # retrieving simulated data in the form of numpy array
simulation = model.get_simulation(P)
simulation.setBackground(ba.PoissonBackground())
result = simulation.simulate()
- return result.npArray()
+ return result.noisy(0.1, 0.1)
class PlotObserver:
@@ -57,7 +56,7 @@ class PlotObserver:
self.update(fit_objective)
@staticmethod
- def plot_real_data(data):
+ def plot_data(data):
"""
Plot experimental data as heatmap with horizontal/vertical lines
representing slices on top.
@@ -121,15 +120,15 @@ class PlotObserver:
"""
self.fig.clf()
- real_data = fit_objective.experimentalData()
+ data = fit_objective.experimentalData()
simul_data = fit_objective.simulationResult()
# plot real data
plt.subplot(2, 2, 1)
- self.plot_real_data(real_data)
+ self.plot_data(data)
# horizontal slices
- slices = [("real", real_data.xProjection(alpha_slice_value)),
+ slices = [("real", data.xProjection(alpha_slice_value)),
("simul", simul_data.xProjection(alpha_slice_value))]
title = ("Horizontal slice at alpha =" +
'{:3.1f}'.format(alpha_slice_value))
@@ -137,7 +136,7 @@ class PlotObserver:
self.plot_slices(slices, title)
# vertical slices
- slices = [("real", real_data.yProjection(phi_slice_value)),
+ slices = [("real", data.yProjection(phi_slice_value)),
("simul", simul_data.yProjection(phi_slice_value))]
title = "Vertical slice at phi =" + '{:3.1f}'.format(
phi_slice_value)
@@ -154,13 +153,13 @@ def run_fitting():
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_masked_simulation, real_data, 1)
+ fit_objective.addFitPair(get_masked_simulation, data, 1)
fit_objective.initPrint(10)
- # creating custom observer which will draw fit progress
+ # creating custom observer which will draw fit progress
plotter = PlotObserver()
<%= sm ? "plt.close() # (hide plot) " :"" %>fit_objective.initPlot(10, plotter)
diff --git a/rawEx/fit/scatter2d/fit_gisas.py b/rawEx/fit/scatter2d/fit_gisas.py
index a255f1df1801bb4270efc726f1570f42d20504a9..bdba0108255aa4504e117a2830ae33ffe086fcfa 100755
--- a/rawEx/fit/scatter2d/fit_gisas.py
+++ b/rawEx/fit/scatter2d/fit_gisas.py
@@ -13,13 +13,13 @@ from bornagain import ba_fitmonitor
from matplotlib import pyplot as plt
-def run_fitting():
+if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
- real_data = np.loadtxt(os.path.join(datadir, "scatter2d/faked_gisas1.dat.gz"),
- dtype=float)
+ fname = os.path.join(datadir, "scatter2d/faked_gisas1.dat.gz")
+ data = ba.readData2D(fname, ba.csv2D)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data)
+ fit_objective.addFitPair(model.get_simulation, data)
fit_objective.initPrint(10) # Print on every 10th iteration.
observer = ba_fitmonitor.PlotterGISAS()
@@ -39,8 +39,4 @@ def run_fitting():
# Save simulation image corresponding to the best fit parameters
np.savetxt("fit.txt", fit_objective.simulationResult().npArray())
-
-
-if __name__ == '__main__':
- run_fitting()
plt.show()
diff --git a/rawEx/fit/scatter2d/fit_with_masks.py b/rawEx/fit/scatter2d/fit_with_masks.py
index 2eabe4707478dbdcf6dd57b88142fe1b8b4ad381..aaf0159d64602932b5ee0a1812ba5dff86f3bb0a 100755
--- a/rawEx/fit/scatter2d/fit_with_masks.py
+++ b/rawEx/fit/scatter2d/fit_with_masks.py
@@ -56,10 +56,10 @@ def add_mask_to_simulation(simulation):
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data().noisy(0.1, 0.1)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_masked_simulation, real_data, 1)
+ fit_objective.addFitPair(get_masked_simulation, data, 1)
fit_objective.initPrint(10)
observer = ba_fitmonitor.PlotterGISAS()
<%= sm ? "plt.close() # (hide plot) " :"" %>fit_objective.initPlot(10, observer)
diff --git a/rawEx/fit/scatter2d/lmfit_basics.py b/rawEx/fit/scatter2d/lmfit_basics.py
index 4216709ddd27c811ebdea40cefaaecb5fcd0d4a4..ff4a7afe4978a3bfbb572452626fe2b78fa970ca 100755
--- a/rawEx/fit/scatter2d/lmfit_basics.py
+++ b/rawEx/fit/scatter2d/lmfit_basics.py
@@ -9,10 +9,10 @@ import model2_hexlattice as model
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ fit_objective.addFitPair(model.get_simulation, data, 1)
fit_objective.initPrint(10)
P = lmfit.Parameters()
diff --git a/rawEx/fit/scatter2d/lmfit_with_plotting.py b/rawEx/fit/scatter2d/lmfit_with_plotting.py
index 2001282a54ce167b25179a3e1fc11d60dc7259d6..7e68bd8c35b0412ae8c948a426cc0d5e5b7cb0ab 100755
--- a/rawEx/fit/scatter2d/lmfit_with_plotting.py
+++ b/rawEx/fit/scatter2d/lmfit_with_plotting.py
@@ -25,16 +25,16 @@ class LMFITPlotter:
if __name__ == '__main__':
- real_data = model.create_real_data()
+ data = model.fake_data()
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
+ fit_objective.addFitPair(model.get_simulation, data, 1)
fit_objective.initPrint(10)
P = lmfit.Parameters()
P.add('radius', value=7*nm, min=5*nm, max=8*nm)
P.add('length', value=10*nm, min=8*nm, max=14*nm)
-
+
<%= sm ? "" :"plotter = LMFITPlotter(fit_objective)" %>
result = lmfit.minimize(fit_objective.evaluate_residuals,
P<%= sm ? "" :", iter_cb=plotter" %>)
diff --git a/rawEx/fit/scatter2d/minimizer_settings.py b/rawEx/fit/scatter2d/minimizer_settings.py
index c39b50349dea59a5692ec1cf8050bfb1a022b91f..48808a4725bb7535f7a249625808a6d3a065e84b 100755
--- a/rawEx/fit/scatter2d/minimizer_settings.py
+++ b/rawEx/fit/scatter2d/minimizer_settings.py
@@ -51,7 +51,7 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, get_sample(P), detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data from simulated image with default parameters.
"""
@@ -66,7 +66,7 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- return result.npArray()
+ return result
def run_fitting():
@@ -74,7 +74,7 @@ def run_fitting():
main function to run fitting
"""
- real_data = create_real_data()
+ data = fake_data()
# prints info about available minimizers
print(ba.MinimizerFactory().catalogToString())
@@ -83,7 +83,7 @@ def run_fitting():
print(ba.MinimizerFactory().catalogDetailsToString())
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(get_simulation, real_data, 1)
+ fit_objective.addFitPair(get_simulation, data, 1)
fit_objective.initPrint(10)
P = ba.Parameters()
diff --git a/rawEx/fit/scatter2d/model1_cylinders.py b/rawEx/fit/scatter2d/model1_cylinders.py
index a6d950027fa7b3d2ed21be576070d285c633f1ef..56d6d4e11bdc95d2d4cc7e0acb9ccb9a16ff37a8 100644
--- a/rawEx/fit/scatter2d/model1_cylinders.py
+++ b/rawEx/fit/scatter2d/model1_cylinders.py
@@ -42,19 +42,13 @@ def get_simulation(P):
return ba.ScatteringSimulation(beam, sample, detector)
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
P = {'radius': 5*nm, 'height': 10*nm}
- # retrieving simulated data in the form of numpy array
simulation = get_simulation(P)
result = simulation.simulate()
- real_data = result.npArray()
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
diff --git a/rawEx/fit/scatter2d/model2_hexlattice.py b/rawEx/fit/scatter2d/model2_hexlattice.py
index b29513ea87faa0e54b72b0aa1fb75ffa50df68f2..e6f82616e2a187b9a576e44b1f599de6230af5d2 100644
--- a/rawEx/fit/scatter2d/model2_hexlattice.py
+++ b/rawEx/fit/scatter2d/model2_hexlattice.py
@@ -48,7 +48,7 @@ def get_simulation(P):
return simulation
-def create_real_data():
+def fake_data():
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -56,12 +56,4 @@ def create_real_data():
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # 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
- return noisy
+ return result.noisy(0.1, 0.1)
diff --git a/rawEx/fit/scatter2d/multiple_datasets.py b/rawEx/fit/scatter2d/multiple_datasets.py
index 0c2d8df452c7da8235667ca3a683e37d1457dc55..5edec05979f6e0406da6049f0c9f81f5c4528860 100755
--- a/rawEx/fit/scatter2d/multiple_datasets.py
+++ b/rawEx/fit/scatter2d/multiple_datasets.py
@@ -59,7 +59,7 @@ def simulation2(P):
return get_simulation(P)
-def create_real_data(incident_alpha):
+def fake_data(incident_alpha):
"""
Generating "real" data by adding noise to the simulated data.
"""
@@ -73,14 +73,7 @@ def create_real_data(incident_alpha):
simulation = get_simulation(P)
result = simulation.simulate()
- # retrieving simulated data in the form of numpy array
- real_data = result.npArray()
-
- # spoiling simulated data with the noise to produce "real" data
- noise_factor = 0.1
- noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
- noisy[noisy < 0.1] = 0.1
- return noisy
+ return result.noisy(0.1, 0.1)
class PlotObserver():
@@ -99,14 +92,14 @@ class PlotObserver():
@staticmethod
def plot_dataset(fit_objective, canvas):
for i_dataset in range(0, fit_objective.fitObjectCount()):
- real_data = fit_objective.experimentalData(i_dataset)
+ data = fit_objective.experimentalData(i_dataset)
simul_data = fit_objective.simulationResult(i_dataset)
chi2_map = fit_objective.relativeDifference(i_dataset)
- zmax = real_data.maxVal()
+ zmax = data.maxVal()
plt.subplot(canvas[i_dataset*3])
- bp.plot_simres(real_data,
+ bp.plot_simres(data,
title="\"Real\" data - #" +
str(i_dataset + 1),
intensity_min=1,
@@ -189,12 +182,12 @@ def run_fitting():
main function to run fitting
"""
- data1 = create_real_data(0.1*deg)
- data2 = create_real_data(0.4*deg)
+ data1 = fake_data(0.1*deg)
+ data2 = fake_data(0.4*deg)
fit_objective = ba.FitObjective()
- fit_objective.addSimulationAndData(simulation1, data1, 1)
- fit_objective.addSimulationAndData(simulation2, data2, 1)
+ fit_objective.addFitPair(simulation1, data1, 1)
+ fit_objective.addFitPair(simulation2, data2, 1)
fit_objective.initPrint(10)
# creating custom observer which will draw fit progress
diff --git a/rawEx/fit/specular/Honeycomb_fit.py b/rawEx/fit/specular/Honeycomb_fit.py
index 39e0f863c6c9f170e387496c394f2f60963f8cff..bcb413ac4e6966ac1b85fb8d5a992da52d08157b 100755
--- a/rawEx/fit/specular/Honeycomb_fit.py
+++ b/rawEx/fit/specular/Honeycomb_fit.py
@@ -88,7 +88,7 @@ def run_simulation(qaxis, P, *, sign, T):
# Experimental data
####################################################################
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
# Load from disk.
filepath = os.path.join(datadir, filename)
with open(filepath, 'r') as f:
@@ -217,11 +217,11 @@ if __name__ == '__main__':
qmin = 0.08
qmax = 1.4
- expData = [
- get_experimental_data("specular/honeycomb300p.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb300m.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb150p.dat", qmin, qmax),
- get_experimental_data("specular/honeycomb150m.dat", qmin, qmax)]
+ data = [
+ load_data("specular/honeycomb300p.dat", qmin, qmax),
+ load_data("specular/honeycomb300m.dat", qmin, qmax),
+ load_data("specular/honeycomb150p.dat", qmin, qmax),
+ load_data("specular/honeycomb150m.dat", qmin, qmax)]
simFunctions = [
lambda q, P: run_simulation(q, P, sign=+1, T=300),
@@ -234,14 +234,14 @@ if __name__ == '__main__':
# Plot data with initial model
simResults = [ f(qzs, P) for f in simFunctions ]
- plot(qzs, simResults, expData, [1, 1, 10, 10],
+ plot(qzs, simResults, data, [1, 1, 10, 10],
["300K $+$", "300K $-$", "150K $+$", "150K $-$"])
plot_sld_profile(P)
# Fit
- qaxes = [d[0] for d in expData]
- rdata = [d[1] for d in expData]
+ qaxes = [d[0] for d in data]
+ rdata = [d[1] for d in data]
def par_dict(*args):
return {name: value for name, value in zip(freeParNames, *args)} | fixedP
@@ -280,7 +280,7 @@ if __name__ == '__main__':
P = par_dict(result.x)
sim_results = [ f(qzs, P) for f in simFunctions ]
- plot(qzs, simResults, expData, [1, 1, 10, 10],
+ plot(qzs, simResults, data, [1, 1, 10, 10],
["300K $+$", "300K $-$", "150K $+$", "150K $-$"])
plot_sld_profile(P)
diff --git a/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py b/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
index 782d35698dd97aa9ebb3a70e955f682baed42096..e10a8827c42083aedb8139073726a31919fb7593 100755
--- a/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -32,36 +32,6 @@ def get_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
return simulation
-####################################################################
-# Fit Function #
-####################################################################
-
-def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
- startParams):
-
- fit_objective = ba.FitObjective()
- fit_objective.setObjectiveMetric("chi2")
-
- fit_objective.addSimulationAndData(
- lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
- r_uncertainty[0], 1)
- fit_objective.addSimulationAndData(
- lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
- r_uncertainty[1], 1)
-
- fit_objective.initPrint(10)
-
- P = ba.Parameters()
- for name, p in startParams.items():
- P.add(name, p[0], min=p[1], max=p[2])
-
- minimizer = ba.Minimizer()
-
- result = minimizer.minimize(fit_objective.evaluate, P)
- fit_objective.finalize(result)
-
- return {r.name(): r.value for r in result.parameters()}
-
####################################################################
# Main Function #
####################################################################
@@ -71,8 +41,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = psa.load_exp(fname_stem + "pp.tab")
- expdata_mm = psa.load_exp(fname_stem + "mm.tab")
+ data_pp = psa.load_data(fname_stem + "pp.tab")
+ data_mm = psa.load_data(fname_stem + "mm.tab")
fixedParams = {
# parameters can be moved here to keep them fixed
@@ -114,26 +84,43 @@ if __name__ == '__main__':
legend_pp = "$++$"
legend_mm = "$--$"
- psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+ psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
[legend_pp, legend_mm], [color_pp, color_mm])
- psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+ psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
+
+ fit_objective = ba.FitObjective()
+ fit_objective.setObjectiveMetric("chi2")
+ fit_objective.initPrint(10)
+
+ fit_objective.addFitPair(
+ lambda P: get_Simulation_pp(data_pp.npXcenters(), P),
+ data_pp, 1)
+ fit_objective.addFitPair(
+ lambda P: get_Simulation_mm(data_mm.npXcenters(), P),
+ data_mm, 1)
+
+ P = ba.Parameters()
+ for name, p in startParams.items():
+ P.add(name, p[0], min=p[1], max=p[2])
+
+ minimizer = ba.Minimizer()
+
+ result = minimizer.minimize(fit_objective.evaluate, P)
+ fit_objective.finalize(result)
+
+ fitResult = {r.name(): r.value for r in result.parameters()}
- fitResult = run_fit_ba([expdata_pp[0], expdata_mm[0]],
- [expdata_pp[1], expdata_mm[1]],
- [expdata_pp[2], expdata_mm[2]],
- [get_Simulation_pp, get_Simulation_mm],
- startParams)
print("Fit Result:")
print(fitResult)
q_pp, r_pp = psa.qr(get_Simulation_pp(qzs, fitResult).simulate())
q_mm, r_mm = psa.qr(get_Simulation_mm(qzs, fitResult).simulate())
- psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+ psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
[legend_pp, legend_mm], [color_pp, color_mm])
<%= sm ? "#" : "" %>plt.draw()
- psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+ psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
<%= sm ? "#" : "" %>plt.draw()
<%= sm ? "#" : "" %>plt.show()
diff --git a/rawEx/fit/specular/Pt_layer_fit.py b/rawEx/fit/specular/Pt_layer_fit.py
index 17f91023a53e310fcf6ba83ca86323af5889388c..0d0b68d43839ac197d3c39e6b4e24ecd82905332 100755
--- a/rawEx/fit/specular/Pt_layer_fit.py
+++ b/rawEx/fit/specular/Pt_layer_fit.py
@@ -60,7 +60,7 @@ def get_simulation(q_axis, P):
# Experimental data
####################################################################
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
filepath = os.path.join(datadir, filename)
data = np.genfromtxt(filepath, unpack=True)
@@ -144,7 +144,7 @@ if __name__ == '__main__':
qmin = 0.18
qmax = 2.4
- data = get_experimental_data("specular/RvsQ_36563_36662.dat.gz", qmin, qmax)
+ data = load_data("specular/RvsQ_36563_36662.dat.gz", qmin, qmax)
qzs = np.linspace(qmin, qmax, 1500)
diff --git a/rawEx/specular/PolarizedSpinAsymmetry.py b/rawEx/specular/PolarizedSpinAsymmetry.py
index 15ae40211155a59284e76404ee97d1af7b7d289a..c3f6c0b2e36d74fc0a47df99622c1771efc57538 100755
--- a/rawEx/specular/PolarizedSpinAsymmetry.py
+++ b/rawEx/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
ax = fig.add_subplot(111)
for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
- ax.errorbar(exp[0],
- exp[1],
- xerr=exp[3],
- yerr=exp[2],
+ ax.errorbar(exp.xAxis().binCenters(),
+ exp.flatVector(),
+ # xerr=TODO i742,
+ yerr=exp.errorSigmas(),
fmt='.',
markersize=0.75,
linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
Plot the simulated spin asymmetry as well its
experimental counterpart with errorbars
"""
+ Yp = data_pp.npArray()
+ Ym = data_mm.npArray()
+ Ep = data_pp.npErrors()
+ Em = data_mm.npErrors()
# compute the errorbars of the spin asymmetry
- delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
- data_mm[1]**2 * data_pp[2]**2 ) /
- ( data_pp[1] + data_mm[1] )**4 )
+ delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
fig = plt.figure()
ax = fig.add_subplot(111)
- ax.errorbar(data_pp[0],
- (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
- xerr=data_pp[3],
+ ax.errorbar(data_pp.xAxis().binCenters(),
+ (Yp - Ym) / (Yp + Ym),
+ # xerr=TODO i742,
yerr=delta,
fmt='.',
markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
# Data Handling #
####################################################################
-def load_exp(fname):
- dat = numpy.loadtxt(fname)
- return numpy.transpose(dat)
+def load_data(fname):
+ flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+ return ba.readData1D(fname, ba.csv1D, flags)
####################################################################
# Main Function #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
datadir = os.getenv('BA_DATA_DIR', '')
fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
- expdata_pp = load_exp(fname_stem + "pp.tab")
- expdata_mm = load_exp(fname_stem + "mm.tab")
+ expdata_pp = load_data(fname_stem + "pp.tab")
+ expdata_mm = load_data(fname_stem + "mm.tab")
fixedP = {
# parameters from our own fit run
diff --git a/rawEx/specular/VsGenx.py b/rawEx/specular/VsGenx.py
index 8150c11d15f18bdbd357dc5f3db0a00c5e7bfcb4..dbd5a59f3604eda8f97083de288a0d2102e0e40b 100755
--- a/rawEx/specular/VsGenx.py
+++ b/rawEx/specular/VsGenx.py
@@ -16,7 +16,7 @@ def reference_data(filename):
"""
Loads and returns reference data from GenX simulation
"""
- ax_values, real_data = np.loadtxt(filename,
+ ax_values, data = np.loadtxt(filename,
usecols=(0, 1),
skiprows=3,
unpack=True)
@@ -24,7 +24,7 @@ def reference_data(filename):
# translate axis values from double incident angle to incident angle
ax_values *= 0.5
- return ax_values, real_data
+ return ax_values, data
def get_sample():