mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2026-10-28 07:51:28 +01:00
CI/CD: Integrate pre-commit hooks and GitHub Actions workflow (#303)
To improve codebase quality and reduce human error, this PR introduces the pre-commit framework. This ensures that all code adheres to project standards before it is even committed, maintaining a consistent style and catching common mistakes early. Key Changes: - Code Formatting: Automated C++ formatting using clang-format (based on the project's .clang-format file). - Syntax Validation: Basic checks for file integrity and syntax. - Spell Check: Automated scanning for typos in source code and comments. - CMake Formatting: Standardization of CMakeLists.txt and .cmake configuration files. - GitHub Workflow: Added a CI action that validates every Pull Request against the pre-commit configuration to ensure compliance. The configuration includes a [ci] block to handle automated fixes within the PR. Currently, this is disabled. If we want the CI to automatically commit formatting fixes back to the PR branch, this can be toggled to true in .pre-commit-config.yaml. ```yaml ci: autofix_commit_msg: [pre-commit] auto fixes from pre-commit hooks autofix_prs: false autoupdate_schedule: monthly ``` The last large commit with the fit functions, for example, was not formatted according to the clang-format rules. This PR would allow to avoid similar mistakes in the future. Python fomat with `ruff` for tests and sanitiser for `.ipynb` notebooks can be added as well.
This commit is contained in:
lunin_l authored and 
GitHub committed
1 parent
85098d2008
commit
8f8173feb6
36 files changed
+4104
-3853
No files matched your search
+35
-45
@@ -1,37 +1,38 @@
|
||||
# SPDX-License-Identifier: MPL-2.0
|
||||
|
||||
find_package (Python 3.10 COMPONENTS Interpreter Development.Module REQUIRED)
|
||||
find_package(
|
||||
Python 3.10
|
||||
COMPONENTS Interpreter Development.Module
|
||||
REQUIRED)
|
||||
set(PYBIND11_FINDPYTHON ON) # Needed for RH8
|
||||
|
||||
# Download or find pybind11 depending on configuration
|
||||
if(AARE_FETCH_PYBIND11)
|
||||
FetchContent_Declare(
|
||||
pybind11
|
||||
GIT_REPOSITORY https://github.com/pybind/pybind11
|
||||
GIT_TAG v2.13.6
|
||||
)
|
||||
FetchContent_MakeAvailable(pybind11)
|
||||
FetchContent_Declare(
|
||||
pybind11
|
||||
GIT_REPOSITORY https://github.com/pybind/pybind11
|
||||
GIT_TAG v2.13.6)
|
||||
FetchContent_MakeAvailable(pybind11)
|
||||
else()
|
||||
find_package(pybind11 2.13 REQUIRED)
|
||||
find_package(pybind11 2.13 REQUIRED)
|
||||
endif()
|
||||
|
||||
# Add the compiled python extension
|
||||
pybind11_add_module(
|
||||
_aare # name of the module
|
||||
src/module.cpp # source file
|
||||
)
|
||||
|
||||
set_target_properties(_aare PROPERTIES
|
||||
LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
|
||||
_aare # name of the module
|
||||
src/module.cpp # source file
|
||||
)
|
||||
|
||||
set_target_properties(_aare PROPERTIES LIBRARY_OUTPUT_DIRECTORY
|
||||
${CMAKE_BINARY_DIR})
|
||||
target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
|
||||
|
||||
target_include_directories(_aare SYSTEM PRIVATE
|
||||
$<TARGET_PROPERTY:Minuit2::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>
|
||||
)
|
||||
target_include_directories(
|
||||
_aare SYSTEM
|
||||
PRIVATE $<TARGET_PROPERTY:Minuit2::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
|
||||
|
||||
# List of python files to be copied to the build directory
|
||||
set( PYTHON_FILES
|
||||
set(PYTHON_FILES
|
||||
aare/__init__.py
|
||||
aare/CtbRawFile.py
|
||||
aare/ClusterFinder.py
|
||||
@@ -42,43 +43,32 @@ set( PYTHON_FILES
|
||||
aare/RawFile.py
|
||||
aare/transform.py
|
||||
aare/ScanParameters.py
|
||||
aare/utils.py
|
||||
)
|
||||
|
||||
aare/utils.py)
|
||||
|
||||
# Copy the python files to the build directory
|
||||
foreach(FILE ${PYTHON_FILES})
|
||||
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} )
|
||||
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE})
|
||||
endforeach(FILE ${PYTHON_FILES})
|
||||
|
||||
set_target_properties(_aare PROPERTIES
|
||||
LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/aare
|
||||
)
|
||||
|
||||
set(PYTHON_EXAMPLES
|
||||
examples/play.py
|
||||
examples/fits.py
|
||||
)
|
||||
set_target_properties(_aare PROPERTIES LIBRARY_OUTPUT_DIRECTORY
|
||||
${CMAKE_BINARY_DIR}/aare)
|
||||
|
||||
set(PYTHON_EXAMPLES examples/play.py examples/fits.py)
|
||||
|
||||
# Copy the python examples to the build directory
|
||||
foreach(FILE ${PYTHON_EXAMPLES})
|
||||
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} )
|
||||
message(STATUS "Copying ${FILE} to ${CMAKE_BINARY_DIR}/${FILE}")
|
||||
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE})
|
||||
message(STATUS "Copying ${FILE} to ${CMAKE_BINARY_DIR}/${FILE}")
|
||||
endforeach(FILE ${PYTHON_EXAMPLES})
|
||||
|
||||
|
||||
if(AARE_INSTALL_PYTHONEXT)
|
||||
install(
|
||||
TARGETS _aare
|
||||
EXPORT "${TARGETS_EXPORT_NAME}"
|
||||
LIBRARY DESTINATION aare
|
||||
COMPONENT python
|
||||
)
|
||||
install(
|
||||
TARGETS _aare
|
||||
EXPORT "${TARGETS_EXPORT_NAME}"
|
||||
LIBRARY DESTINATION aare COMPONENT python)
|
||||
|
||||
install(
|
||||
FILES ${PYTHON_FILES}
|
||||
DESTINATION aare
|
||||
COMPONENT python
|
||||
)
|
||||
endif()
|
||||
install(
|
||||
FILES ${PYTHON_FILES}
|
||||
DESTINATION aare
|
||||
COMPONENT python)
|
||||
endif()
|
||||
+132
-112
@@ -5,72 +5,88 @@
|
||||
#include <pybind11/stl.h>
|
||||
#include <pybind11/stl_bind.h>
|
||||
|
||||
#include "aare/Chi2.hpp"
|
||||
#include "aare/Fit.hpp"
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/Models.hpp"
|
||||
#include "aare/Chi2.hpp"
|
||||
|
||||
namespace py = pybind11;
|
||||
using namespace pybind11::literals;
|
||||
|
||||
template <typename Model, typename FCN>
|
||||
py::object fit_dispatch(
|
||||
const aare::FitModel<Model>& model,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj,
|
||||
int n_threads);
|
||||
py::object
|
||||
fit_dispatch(const aare::FitModel<Model> &model,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj, int n_threads);
|
||||
|
||||
template <typename Model>
|
||||
void bind_fit_model(py::module& m, const char* name) {
|
||||
template <typename Model> void bind_fit_model(py::module &m, const char *name) {
|
||||
using FM = aare::FitModel<Model>;
|
||||
using FCN = aare::func::Chi2Model1DGrad<Model>;
|
||||
py::class_<FM>(m, name)
|
||||
.def(py::init<unsigned int, unsigned int, double, bool>(),
|
||||
py::arg("strategy") = 0,
|
||||
py::arg("max_calls") = 100,
|
||||
py::arg("tolerance") = 0.5,
|
||||
py::arg("compute_errors") = false)
|
||||
.def("SetParLimits", py::overload_cast<unsigned int, double, double>(&FM::SetParLimits), py::arg("idx"), py::arg("lo"), py::arg("hi"))
|
||||
.def("SetParLimits", py::overload_cast<const std::string&, double, double>(&FM::SetParLimits), py::arg("idx"), py::arg("lo"), py::arg("hi"))
|
||||
.def("FixParameter", py::overload_cast<unsigned int, double>(&FM::FixParameter), py::arg("idx"), py::arg("val"))
|
||||
.def("FixParameter", py::overload_cast<const std::string&, double>(&FM::FixParameter), py::arg("idx"), py::arg("val"))
|
||||
.def("ReleaseParameter", py::overload_cast<unsigned int>(&FM::ReleaseParameter), py::arg("idx"))
|
||||
.def("ReleaseParameter", py::overload_cast<const std::string&>(&FM::ReleaseParameter), py::arg("idx"))
|
||||
.def("SetParameter", py::overload_cast<unsigned int, double>(&FM::SetParameter), py::arg("idx"), py::arg("val"))
|
||||
.def("SetParameter", py::overload_cast<const std::string&, double>(&FM::SetParameter), py::arg("idx"), py::arg("val"))
|
||||
py::arg("strategy") = 0, py::arg("max_calls") = 100,
|
||||
py::arg("tolerance") = 0.5, py::arg("compute_errors") = false)
|
||||
.def("SetParLimits",
|
||||
py::overload_cast<unsigned int, double, double>(&FM::SetParLimits),
|
||||
py::arg("idx"), py::arg("lo"), py::arg("hi"))
|
||||
.def("SetParLimits",
|
||||
py::overload_cast<const std::string &, double, double>(
|
||||
&FM::SetParLimits),
|
||||
py::arg("idx"), py::arg("lo"), py::arg("hi"))
|
||||
.def("FixParameter",
|
||||
py::overload_cast<unsigned int, double>(&FM::FixParameter),
|
||||
py::arg("idx"), py::arg("val"))
|
||||
.def("FixParameter",
|
||||
py::overload_cast<const std::string &, double>(&FM::FixParameter),
|
||||
py::arg("idx"), py::arg("val"))
|
||||
.def("ReleaseParameter",
|
||||
py::overload_cast<unsigned int>(&FM::ReleaseParameter),
|
||||
py::arg("idx"))
|
||||
.def("ReleaseParameter",
|
||||
py::overload_cast<const std::string &>(&FM::ReleaseParameter),
|
||||
py::arg("idx"))
|
||||
.def("SetParameter",
|
||||
py::overload_cast<unsigned int, double>(&FM::SetParameter),
|
||||
py::arg("idx"), py::arg("val"))
|
||||
.def("SetParameter",
|
||||
py::overload_cast<const std::string &, double>(&FM::SetParameter),
|
||||
py::arg("idx"), py::arg("val"))
|
||||
.def("GetParName", &FM::GetParName, py::arg("idx"))
|
||||
.def("GetParNames", &FM::GetParNames)
|
||||
.def_property_readonly("par_names", &FM::GetParNames)
|
||||
.def_property_readonly("n_par", [](py::object /*cls*/) { return Model::npar; })
|
||||
.def_property_readonly("n_par",
|
||||
[](py::object /*cls*/) { return Model::npar; })
|
||||
.def_property("max_calls", &FM::max_calls, &FM::SetMaxCalls)
|
||||
.def_property("tolerance", &FM::tolerance, &FM::SetTolerance)
|
||||
.def_property("compute_errors", &FM::compute_errors, &FM::SetComputeErrors)
|
||||
.def("__call__",
|
||||
[](const FM& /*self*/,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> par)
|
||||
{
|
||||
.def_property("compute_errors", &FM::compute_errors,
|
||||
&FM::SetComputeErrors)
|
||||
.def(
|
||||
"__call__",
|
||||
[](const FM & /*self*/,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast>
|
||||
par) {
|
||||
auto x_view = make_view_1d(x);
|
||||
auto p_view = make_view_1d(par);
|
||||
|
||||
std::vector<double> pvec(p_view.begin(), p_view.end());
|
||||
|
||||
auto* result = new aare::NDArray<double, 1>({x_view.size()});
|
||||
auto *result = new aare::NDArray<double, 1>({x_view.size()});
|
||||
for (ssize_t i = 0; i < x_view.size(); ++i)
|
||||
(*result)(i) = Model::eval(x_view[i], pvec);
|
||||
|
||||
return return_image_data(result);
|
||||
},
|
||||
py::arg("x"), py::arg("par"))
|
||||
.def("fit",
|
||||
[](const FM& self,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj,
|
||||
int n_threads) -> py::object
|
||||
{
|
||||
return fit_dispatch<Model, FCN>(self, x, y, y_err_obj, n_threads);
|
||||
.def(
|
||||
"fit",
|
||||
[](const FM &self,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj, int n_threads) -> py::object {
|
||||
return fit_dispatch<Model, FCN>(self, x, y, y_err_obj,
|
||||
n_threads);
|
||||
},
|
||||
R"doc(
|
||||
Fit this model to 1D or 3D data using Minuit2.
|
||||
@@ -86,24 +102,21 @@ void bind_fit_model(py::module& m, const char* name) {
|
||||
n_threads : int
|
||||
Number of threads for 3D parallel loop.
|
||||
)doc",
|
||||
py::arg("x"),
|
||||
py::arg("y"),
|
||||
py::arg("y_err") = py::none(),
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err") = py::none(),
|
||||
py::arg("n_threads") = 4);
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
py::dict pack_1d_result_dict(const aare::NDArray<double, 1>& result,
|
||||
bool compute_errors)
|
||||
{
|
||||
py::dict pack_1d_result_dict(const aare::NDArray<double, 1> &result,
|
||||
bool compute_errors) {
|
||||
constexpr std::size_t npar = Model::npar;
|
||||
|
||||
auto res = result.view();
|
||||
|
||||
auto par_out = new NDArray<double, 1>({npar}, 0.0);
|
||||
auto chi2_out = new NDArray<double, 1>({1}, 0.0);
|
||||
auto par_out = new NDArray<double, 1>({npar}, 0.0);
|
||||
auto chi2_out = new NDArray<double, 1>({1}, 0.0);
|
||||
|
||||
auto par_view = par_out->view();
|
||||
auto par_view = par_out->view();
|
||||
auto chi2_view = chi2_out->view();
|
||||
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
@@ -111,7 +124,7 @@ py::dict pack_1d_result_dict(const aare::NDArray<double, 1>& result,
|
||||
}
|
||||
|
||||
if (compute_errors) {
|
||||
auto err_out = new NDArray<double, 1>({npar}, 0.0);
|
||||
auto err_out = new NDArray<double, 1>({npar}, 0.0);
|
||||
auto err_view = err_out->view();
|
||||
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
@@ -120,58 +133,59 @@ py::dict pack_1d_result_dict(const aare::NDArray<double, 1>& result,
|
||||
|
||||
chi2_view(0) = res(2 * npar);
|
||||
|
||||
return py::dict(
|
||||
"par"_a = return_image_data(par_out),
|
||||
"par_err"_a = return_image_data(err_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
"par_err"_a = return_image_data(err_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
} else {
|
||||
chi2_view(0) = res(npar);
|
||||
|
||||
return py::dict(
|
||||
"par"_a = return_image_data(par_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
}
|
||||
}
|
||||
|
||||
// Helper: typed dispatch for one Model, handles 1D/3D + y_err logic
|
||||
template <typename Model, typename FCN>
|
||||
py::object fit_dispatch(
|
||||
const aare::FitModel<Model>& model,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj,
|
||||
int n_threads)
|
||||
{
|
||||
py::object
|
||||
fit_dispatch(const aare::FitModel<Model> &model,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj, int n_threads) {
|
||||
constexpr std::size_t npar = Model::npar;
|
||||
|
||||
if (y.ndim() == 3) {
|
||||
auto par_out = new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
|
||||
auto chi2_out= new NDArray<double, 2>({y.shape(0), y.shape(1)}, 0.0);
|
||||
auto par_out =
|
||||
new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
|
||||
auto chi2_out = new NDArray<double, 2>({y.shape(0), y.shape(1)}, 0.0);
|
||||
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_3d(y);
|
||||
|
||||
if (!y_err_obj.is_none()) {
|
||||
auto y_err = py::cast<py::array_t<double,
|
||||
py::array::c_style | py::array::forcecast>>(y_err_obj);
|
||||
auto y_err = py::cast<
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast>>(
|
||||
y_err_obj);
|
||||
|
||||
if (y_err.ndim() != 3) {
|
||||
throw std::runtime_error("For 3D input y, y_err must also be 3D.");
|
||||
throw std::runtime_error(
|
||||
"For 3D input y, y_err must also be 3D.");
|
||||
}
|
||||
|
||||
auto err_out = new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
|
||||
|
||||
auto err_out =
|
||||
new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
|
||||
auto y_view_err = make_view_3d(y_err);
|
||||
|
||||
aare::fit_3d<Model, FCN>(model, x_view, y_view, y_view_err,
|
||||
par_out->view(), err_out->view(), chi2_out->view(), n_threads);
|
||||
|
||||
aare::fit_3d<Model, FCN>(model, x_view, y_view, y_view_err,
|
||||
par_out->view(), err_out->view(),
|
||||
chi2_out->view(), n_threads);
|
||||
|
||||
if (model.compute_errors()) {
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
"par_err"_a = return_image_data(err_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
} else {
|
||||
delete err_out;
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
}
|
||||
} else {
|
||||
@@ -179,9 +193,10 @@ py::object fit_dispatch(
|
||||
NDView<double, 3> dummy_err{};
|
||||
NDView<double, 3> dummy_err_out{};
|
||||
|
||||
aare::fit_3d<Model, FCN>(model, x_view, y_view, dummy_err,
|
||||
par_out->view(), dummy_err_out, chi2_out->view(), n_threads);
|
||||
|
||||
aare::fit_3d<Model, FCN>(model, x_view, y_view, dummy_err,
|
||||
par_out->view(), dummy_err_out,
|
||||
chi2_out->view(), n_threads);
|
||||
|
||||
return py::dict("par"_a = return_image_data(par_out),
|
||||
"chi2"_a = return_image_data(chi2_out));
|
||||
}
|
||||
@@ -192,15 +207,18 @@ py::object fit_dispatch(
|
||||
auto y_view = make_view_1d(y);
|
||||
|
||||
if (!y_err_obj.is_none()) {
|
||||
auto y_err = py::cast<py::array_t<double,
|
||||
py::array::c_style | py::array::forcecast>>(y_err_obj);
|
||||
auto y_err = py::cast<
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast>>(
|
||||
y_err_obj);
|
||||
|
||||
if (y_err.ndim() != 1) {
|
||||
throw std::runtime_error("For 1D input y, y_err must also be 1D.");
|
||||
throw std::runtime_error(
|
||||
"For 1D input y, y_err must also be 1D.");
|
||||
}
|
||||
|
||||
auto y_view_err = make_view_1d(y_err);
|
||||
result = aare::fit_pixel<Model, FCN>(model, x_view, y_view, y_view_err);
|
||||
result =
|
||||
aare::fit_pixel<Model, FCN>(model, x_view, y_view, y_view_err);
|
||||
} else {
|
||||
result = aare::fit_pixel<Model, FCN>(model, x_view, y_view);
|
||||
}
|
||||
@@ -403,7 +421,6 @@ void define_fit_bindings(py::module &m) {
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
|
||||
m.def(
|
||||
"fit_pol1",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
@@ -575,7 +592,6 @@ void define_fit_bindings(py::module &m) {
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
|
||||
m.def(
|
||||
"fit_scurve2",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
@@ -661,7 +677,6 @@ void define_fit_bindings(py::module &m) {
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
|
||||
// ── Bind model classes ──────────────────────────────────────────
|
||||
bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
|
||||
bind_fit_model<aare::model::RisingScurve>(m, "RisingScurve");
|
||||
@@ -669,48 +684,57 @@ void define_fit_bindings(py::module &m) {
|
||||
bind_fit_model<aare::model::Pol1>(m, "Pol1");
|
||||
bind_fit_model<aare::model::Pol2>(m, "Pol2");
|
||||
|
||||
m.def("fit",
|
||||
m.def(
|
||||
"fit",
|
||||
[](py::object model_obj,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj,
|
||||
int n_threads) -> py::object
|
||||
{
|
||||
py::object y_err_obj, int n_threads) -> py::object {
|
||||
using namespace aare::model;
|
||||
using namespace aare::func;
|
||||
|
||||
|
||||
// ── Polynomial of degree 1 ───────
|
||||
if(py::isinstance< aare::FitModel<Pol1> >(model_obj)) {
|
||||
const auto& mdl = model_obj.cast< const aare::FitModel<Pol1>& >();
|
||||
return fit_dispatch<Pol1, Chi2Pol1>(mdl, x, y, y_err_obj, n_threads);
|
||||
if (py::isinstance<aare::FitModel<Pol1>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<Pol1> &>();
|
||||
return fit_dispatch<Pol1, Chi2Pol1>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
}
|
||||
|
||||
// ── Polynomial of degree 2 ───────
|
||||
if(py::isinstance< aare::FitModel<Pol2> >(model_obj)) {
|
||||
const auto& mdl = model_obj.cast< const aare::FitModel<Pol2>& >();
|
||||
return fit_dispatch<Pol2, Chi2Pol2>(mdl, x, y, y_err_obj, n_threads);
|
||||
if (py::isinstance<aare::FitModel<Pol2>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<Pol2> &>();
|
||||
return fit_dispatch<Pol2, Chi2Pol2>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
}
|
||||
// ── Gaussian ───────
|
||||
if(py::isinstance< aare::FitModel<Gaussian> >(model_obj)) {
|
||||
const auto& mdl = model_obj.cast< const aare::FitModel<Gaussian>& >();
|
||||
return fit_dispatch<Gaussian, Chi2Gaussian>(mdl, x, y, y_err_obj, n_threads);
|
||||
if (py::isinstance<aare::FitModel<Gaussian>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<Gaussian> &>();
|
||||
return fit_dispatch<Gaussian, Chi2Gaussian>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
}
|
||||
|
||||
// ── Rising Scurve ───────
|
||||
if(py::isinstance< aare::FitModel<RisingScurve> >(model_obj)) {
|
||||
const auto& mdl = model_obj.cast< const aare::FitModel<RisingScurve>& >();
|
||||
return fit_dispatch<RisingScurve, Chi2RisingScurve>(mdl, x, y, y_err_obj, n_threads);
|
||||
if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<RisingScurve> &>();
|
||||
return fit_dispatch<RisingScurve, Chi2RisingScurve>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
}
|
||||
|
||||
// ── Falling Scurve ───────
|
||||
if(py::isinstance< aare::FitModel<FallingScurve> >(model_obj)) {
|
||||
const auto& mdl = model_obj.cast< const aare::FitModel<FallingScurve>& >();
|
||||
return fit_dispatch<FallingScurve, Chi2FallingScurve>(mdl, x, y, y_err_obj, n_threads);
|
||||
if (py::isinstance<aare::FitModel<FallingScurve>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<FallingScurve> &>();
|
||||
return fit_dispatch<FallingScurve, Chi2FallingScurve>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
}
|
||||
|
||||
throw std::runtime_error(
|
||||
"Unknown model type. Expected Pol1, Pol2, Gaussian, RisingScurve or FallingScurve."
|
||||
);
|
||||
"Unknown model type. Expected Pol1, Pol2, Gaussian, "
|
||||
"RisingScurve or FallingScurve.");
|
||||
},
|
||||
R"(
|
||||
Fit a model to 1D or 3D data using Minuit2.
|
||||
@@ -741,10 +765,6 @@ void define_fit_bindings(py::module &m) {
|
||||
"par_err" : (rows, cols, npar) parameter errors (if compute_errors).
|
||||
"chi2" : (rows, cols) chi-squared per pixel.
|
||||
)",
|
||||
py::arg("model"),
|
||||
py::arg("x"),
|
||||
py::arg("y"),
|
||||
py::arg("y_err") = py::none(),
|
||||
py::arg("n_threads") = 4
|
||||
);
|
||||
py::arg("model"), py::arg("x"), py::arg("y"),
|
||||
py::arg("y_err") = py::none(), py::arg("n_threads") = 4);
|
||||
}
|
||||
+341
-336
@@ -1,339 +1,344 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "a9463197",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from aare import Interpolator\n",
|
||||
"import numpy as np \n",
|
||||
"import matplotlib.pyplot as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "7fdc9ff6",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def calculate_rosenblatt_lookup_table(etacube, etabins_x, etabins_y, etabins_e):\n",
|
||||
"\n",
|
||||
" \"\"\" \n",
|
||||
" Rosenblatt transform: \n",
|
||||
" transforms joint probability P_X,Y to uniform probability between [0,1] U\n",
|
||||
" U_0 = F_X(x) = P(x' < x) marginal CDF of X\n",
|
||||
" U_1 = F_Y|X(y|x) = P(y' < y | x) conditional CDF\n",
|
||||
" \"\"\"\n",
|
||||
"\n",
|
||||
" #marginal CDF for eta_x\n",
|
||||
" pEtaX = etacube.sum(axis=1) # marginal distribution eta_x\n",
|
||||
" CDF_EtaX = np.concatenate([np.zeros(shape=[1, len(etabins_e)]), np.cumsum(pEtaX, axis=0)], axis=0)\n",
|
||||
" CDF_EtaX = CDF_EtaX / CDF_EtaX[-1] #normalize \n",
|
||||
"\n",
|
||||
" # conditional CDF F_Y|X (a CDF in y direction for each etaX bin)\n",
|
||||
" CDF_EtaY_cond = np.cumsum(etacube, axis=1)\n",
|
||||
" np.divide(CDF_EtaY_cond, np.repeat(pEtaX[:, None, :], len(etabins_y), axis=1), out=np.zeros_like(CDF_EtaY_cond, dtype=np.float64), where=np.repeat(pEtaX[:, None, :], len(etabins_y), axis=1)!=0)\n",
|
||||
"\n",
|
||||
" CDF_EtaY_cond = np.concatenate([np.zeros(shape=[len(etabins_x), 1, len(etabins_e)]), CDF_EtaY_cond], axis=1)\n",
|
||||
" #CDF_EtaY_cond = CDF_EtaY_cond / np.where(CDF_EtaY_cond[:,-1:,:] == 0, 1, CDF_EtaY_cond[:,-1:,:]) #normalize\n",
|
||||
"\n",
|
||||
" # calculate (u,v) at the center of each original bin\n",
|
||||
" \"\"\" \n",
|
||||
" Xc_1d = 0.5 * (F_EtaX_edges[:-1] + F_EtaX_edges[1:]) # why the average between two values?\n",
|
||||
" Yc_2d = 0.5 * (F_EtaY_cond_edges[:, :-1] + F_EtaY_cond_edges[:, 1:])\n",
|
||||
"\n",
|
||||
" X_tab = np.repeat(Xc_1d[:, None], nBinY, axis=1) # (nBinX, nBinY)\n",
|
||||
" Y_tab = Yc_2d # (nBinX, nBinY)\n",
|
||||
" \"\"\" \n",
|
||||
" X_tab = np.repeat(CDF_EtaX[:, None, :], etabins_y.shape[0], axis=1) #TODO: actually does not depend on eta_y at all \n",
|
||||
" Y_tab = CDF_EtaY_cond \n",
|
||||
"\n",
|
||||
" print(\"shape X_tab: \", X_tab.shape)\n",
|
||||
" print(\"shape Y_tab: \", Y_tab.shape)\n",
|
||||
"\n",
|
||||
" return X_tab, Y_tab\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "616d04d7",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Text(0, 0.5, 'y')"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
"cells":
|
||||
[
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "a9463197",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from aare import Interpolator\n",
|
||||
"import numpy as np \n",
|
||||
"import matplotlib.pyplot as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "7fdc9ff6",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def calculate_rosenblatt_lookup_table(etacube, etabins_x, etabins_y, etabins_e):\n",
|
||||
"\n",
|
||||
" \"\"\" \n",
|
||||
" Rosenblatt transform: \n",
|
||||
" transforms joint probability P_X,Y to uniform probability between [0,1] U\n",
|
||||
" U_0 = F_X(x) = P(x' < x) marginal CDF of X\n",
|
||||
" U_1 = F_Y|X(y|x) = P(y' < y | x) conditional CDF\n",
|
||||
" \"\"\"\n",
|
||||
"\n",
|
||||
" #marginal CDF for eta_x\n",
|
||||
" pEtaX = etacube.sum(axis=1) # marginal distribution eta_x\n",
|
||||
" CDF_EtaX = np.concatenate([np.zeros(shape=[1, len(etabins_e)]), np.cumsum(pEtaX, axis=0)], axis=0)\n",
|
||||
" CDF_EtaX = CDF_EtaX / CDF_EtaX[-1] #normalize \n",
|
||||
"\n",
|
||||
" # conditional CDF F_Y|X (a CDF in y direction for each etaX bin)\n",
|
||||
" CDF_EtaY_cond = np.cumsum(etacube, axis=1)\n",
|
||||
" np.divide(CDF_EtaY_cond, np.repeat(pEtaX[:, None, :], len(etabins_y), axis=1), out=np.zeros_like(CDF_EtaY_cond, dtype=np.float64), where=np.repeat(pEtaX[:, None, :], len(etabins_y), axis=1)!=0)\n",
|
||||
"\n",
|
||||
" CDF_EtaY_cond = np.concatenate([np.zeros(shape=[len(etabins_x), 1, len(etabins_e)]), CDF_EtaY_cond], axis=1)\n",
|
||||
" #CDF_EtaY_cond = CDF_EtaY_cond / np.where(CDF_EtaY_cond[:,-1:,:] == 0, 1, CDF_EtaY_cond[:,-1:,:]) #normalize\n",
|
||||
"\n",
|
||||
" # calculate (u,v) at the center of each original bin\n",
|
||||
" \"\"\" \n",
|
||||
" Xc_1d = 0.5 * (F_EtaX_edges[:-1] + F_EtaX_edges[1:]) # why the average between two values?\n",
|
||||
" Yc_2d = 0.5 * (F_EtaY_cond_edges[:, :-1] + F_EtaY_cond_edges[:, 1:])\n",
|
||||
"\n",
|
||||
" X_tab = np.repeat(Xc_1d[:, None], nBinY, axis=1) # (nBinX, nBinY)\n",
|
||||
" Y_tab = Yc_2d # (nBinX, nBinY)\n",
|
||||
" \"\"\" \n",
|
||||
" X_tab = np.repeat(CDF_EtaX[:, None, :], etabins_y.shape[0], axis=1) #TODO: actually does not depend on eta_y at all \n",
|
||||
" Y_tab = CDF_EtaY_cond \n",
|
||||
"\n",
|
||||
" print(\"shape X_tab: \", X_tab.shape)\n",
|
||||
" print(\"shape Y_tab: \", Y_tab.shape)\n",
|
||||
"\n",
|
||||
" return X_tab, Y_tab\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "616d04d7",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Text(0, 0.5, 'y')"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAAAbwAAAHFCAYAAACJnUzWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAA6OdJREFUeJzs/Xl4HNWdLo6/rd60tVuyZK1ty4BsGSxkBzExi4kdSJwhJCgEE25IyALkhvGQDOEmM2HyZU2exzMJYZw7iUmYhBACNyE2mYhfFkDEmInCkiDjRXZsLBvb6la35G5JpW5Jvap+f7RO6dTpc6qqW63Nqvd5eISrq+qcOnXqsy8WWZZlmDBhwoQJE+c4CuZ6AiZMmDBhwsRswGR4JkyYMGFiUcBkeCZMmDBhYlHAZHgmTJgwYWJRwGR4JkyYMGFiUcBkeCZMmDBhYlHAZHgmTJgwYWJRwGR4JkyYMGFiUcBkeCZMmDBhYlHAZHgmcsYbb7yBm266CbW1tXA4HKipqcHWrVvx+uuvZ3WfBx98EBaLJac57N27FxaLBXv37s3peqPYvHkzNm/ebOg8i8UCi8WCgoICuFwuNDY24qabbsLu3bsxMTGRcc3KlSvxuc99Lqv5vPbaa3jwwQcxPDyc1XXsWGT9du/endV9tDA2NoYHH3yQ+06efPJJWCwWnDp1Km/jZYNEIoE1a9bg3/7t3+ZkfBpDQ0MoKyvDb37zm7meyuKBbMJEDvi///f/ygUFBfJll10mP/XUU/Krr74q//znP5cvu+wyuaCgQP7P//xPw/fq7e2VX3/99ZzmIUmS/Prrr8uSJOV0vVFs2rRJ3rRpk6Hzzj//fPn111+XX3/9dfnll1+W/+u//ku+7rrrZADyVVddJQ8PD6uu2bdvn9zT05PVfL7zne/IAOR33303q+vYsV555RUZgLxr166s7qOFs2fPygDkBx54IOO3gYEB+fXXX5ej0WjexssGO3bskKuqquRIJDIn47N48MEH5cbGRjkWi831VBYFTIZnImt0dnbKBQUF8kc+8hE5kUiofkskEvJHPvIRuaCgQO7s7NS8z+jo6ExOM6/IhuGtXbuW+9sTTzwhA5A/8YlPTHs+2TK8sbEx7vHZZnhziUQiIdfX18tf//rX53oqCgKBgGyz2eRnnnlmrqeyKGAyPBNZ47rrrpOtVqvc29vL/f3MmTOy1WqVP/KRjyjHHnjgARmA3NXVJd94441yWVmZXFNTo/qNRjQale+55x65urpaLioqkq+66ir5rbfekhsaGuTPfvazynmEYL/yyivKsc9+9rNySUmJfPz4cfnaa6+VS0pKZI/HI99zzz0ZmsWDDz4ov/e975XLy8tll8slv+c975F//OMfyxMTE6rz8sHwZFmWP/zhD8sWi0U+deqUcox9plQqJX/zm9+UV69eLRcWFsput1u++OKL5R07dqjWi/2PrEFDQ4N83XXXyc8995y8fv162el0yv/yL//CHYus389//nP5K1/5ilxdXS0XFhbK73vf++R9+/YZWoPPfvazckNDgyzLsvzuu+9y50bG/OlPf8pl1D/5yU/klpYW2el0yuXl5fLHPvYx+ciRIxnjGH2vPDz33HMyAPnw4cPKsfHxcXn9+vXyBRdcoNK8/X6/XF1dLW/atElOJpO696Zx9dVXy01NTRl7aGJiQr7gggvkD3/4w6rj1157rXzVVVdlNYaJ3GD68ExkhVQqhVdeeQWXXnopPB4P95zly5ejtbUVe/bsQSqVUv328Y9/HI2Njdi1axd++MMfCsf5/Oc/jx07duDzn/882tvbceONN+KGG24w7LNKJBK4/vrrcc0116C9vR233XYb/uM//gP//u//rjrv1KlT+OIXv4hf/epX+PWvf42Pf/zj+NKXvoRvfvObhsbJFtdffz1kWcaf/vQn4Tnf/va38eCDD+KTn/wkfve73+HZZ5/F7bffrjz7HXfcgS996UsAgF//+td4/fXX8frrr+OSSy5R7rFv3z587Wtfw5e//GW88MILuPHGGzXn9a//+q84efIkfvzjH+PHP/4x+vr6sHnzZpw8eTKr56utrcULL7wAALj99tuVud13333Ca7Zv347bb78da9euxa9//Wt873vfw8GDB3H55Zfj+PHjqnONvlcefve736GqqgoXXXSRcqywsBC/+tWvMDAwgNtuuw0AMDExgU996lOQZRm/+MUvYLVas1qDf/qnf8KxY8fwxz/+UXX8D3/4A06cOIF//Md/VB3fvHkz/vznP2ftjzWRA+aa45pYWAgEAjIA+X/9r/+led7NN98sA5D7+/tlWZ7SSu6///6Mc1kN7/DhwzIARSsh+MUvfqHSFmRZrOEBkH/1q1+prv/whz8sNzU1CeecSqXkRCIhP/zww3JFRYVKQs+XhveHP/xBBiD/+7//u3KM1bo+8pGPyOvXr9ccR8uk2dDQIFutVvnYsWPc33jrd8kll6ie99SpU7LdbpfvuOMO1bPpaXiyrG3SZDW8oaEhuaioKEPrOXPmjOx0OuVbbrlFNU4u75XgwgsvlP/+7/+e+9uzzz4rA5B37Ngh33///XJBQYH80ksv6d6Th1QqJZ9//vlyW1ub6vi1114rX3DBBRmaX0dHhwxA/sMf/pDTeCaMw9TwTMwI5Mk2i2z0pZ6mAQCvvvoqAOATn/iE6vjWrVths9kMjW+xWPDRj35UdaylpQWnT59WHduzZw8+8IEPwO12w2q1wm634/7770coFMLAwIChsbKBbKD95Hvf+14cOHAA27Ztw4svvoiRkZGsx2lpacHq1asNn3/LLbeo3lVDQwOuuOIKvPLKK1mPnQ1ef/11jI+PZ0SpLl++HFdffXWGlmT0vfLQ19eHqqoq7m+f+MQn8A//8A/42te+hm9961v413/9V3zwgx/M7mEmUVBQgLvuugu//e1vcebMGQDAiRMn8MILL2Dbtm0Z3wSZk8/ny2k8E8ZhMjwTWaGyshLFxcV49913Nc87deoUiouLsXTpUtXx2tpa3TFCoRAAoLq6WnXcZrOhoqLC0DyLi4tRWFioOuZ0OhGNRpV//+Uvf8GWLVsAAP/1X/+FP//5z/jrX/+Kb3zjGwCA8fFxQ2NlA0KY6+rqhOfce++9eOSRR/DGG2/g2muvRUVFBa655hq89dZbhscxss40ampquMfIu5gpkPvz5ltXV5cxvpH3KsL4+HjGtTRuu+02JBIJ2Gw2fPnLXzYyfc17FRUVKWb7H/zgBygqKlLMpjTInGZiv5lQw2R4JrKC1WrF+9//frz11lvwer3cc7xeL7q6unD11Vdn+D+M5NsRptbf3686nkwm80qAf/nLX8Jut+O3v/0tPvGJT+CKK67ApZdemrf78/D888/DYrHgfe97n/Acm82Ge+65B/v27cPg4CB+8YtfoLe3Fx/60IcwNjZmaJxs8xoDgQD3GC1gFBYWIhaLZZwXDAazGosGub/f78/4ra+vD5WVlTnfm0VlZSUGBwe5v42OjuLWW2/F6tWrUVRUhDvuuGNaY7ndbnz2s5/Fj3/8YwwODuKnP/0pbrnlFpSVlWWcS+aUz2c1wYfJ8ExkjXvvvReyLGPbtm0ZQSmpVAr/8A//AFmWce+99+Z0f8IMnn32WdXx3bt3I5lM5jZpDiwWC2w2m4opj4+P4+c//3nexqDx05/+FH/4wx/wyU9+EitWrDB0TVlZGbZu3Yp//Md/xODgoJKw7XQ6lfnmA7/4xS9U5tbTp0/jtddeUyXbr1y5Eu+8846K6YVCIbz22muqe2Uzt8svvxxFRUV4+umnVce9Xi/27NmDa665JpfH4WLNmjU4ceIE97c777wTZ86cwa9//Wv85Cc/wfPPP4//+I//mNZ4X/7ylxEMBrF161YMDw/jrrvu4p5HAoPoYBoTMwNjDhETJihceeWV2LFjB+6++25s3LgRd911F1asWIEzZ87gBz/4Ad58803s2LEDV1xxRU73X7t2LT75yU/iu9/9LqxWK66++mocPnwY3/3ud+F2u1FQkB857brrrsOjjz6KW265Bf/7f/9vhEIhPPLIIwrBzhXj4+N44403lP8/efIkfvOb3+C3v/0tNm3apBmdCgAf/ehH0dzcjEsvvRTLli3D6dOnsWPHDjQ0NGDVqlUAgIsvvhgA8L3vfQ+f/exnYbfb0dTUBJfLldOcBwYGcMMNN+ALX/gCJEnCAw88gMLCQpXQcuutt+JHP/oRPv3pT+MLX/gCQqEQvv3tb2PJkiWqe7lcLjQ0NKC9vR3XXHMNli5disrKSqxcuTJj3LKyMtx3333413/9V3zmM5/BJz/5SYRCITz00EMoLCzEAw88kNPz8LB582Y8/PDDGBsbQ3FxsXL8xz/+MZ5++mn89Kc/xdq1a7F27Vrcdddd+Jd/+RdceeWVeO973wsgbaY/77zz8NnPfhZPPvmk7nirV6/G3//93+MPf/gDNm7ciHXr1nHPe+ONN1BRUaG8UxMziDkNmTGxoPH666/LW7dulaurq2WbzSZXVVXJH//4x+XXXnst41wSiXn27FnhbzRIHl5VVZVcWFgoX3bZZfLrr78uu91u+Stf+YpynlYenpFxnnjiCbmpqUl2Op3y+eefL2/fvl3+yU9+khEBmU2UJqj8s5KSEvn888+Xt27dKu/atUtOpVIZ17CRk9/97nflK664Qq6srJQdDoe8YsUK+fbbb1fl7smyLN97771yXV2dXFBQwM3D40ErD+/LX/6yvGzZMtnpdCp5jyx+9rOfyRdeeKFcWFgoX3TRRfKzzz6bEaUpy7L88ssvy+95z3tkp9NpKA/vxz/+sdzS0iI7HA7Z7XbLbW1tqnw5Wc7uvfLQ09MjWywWVZTnwYMH5aKiItWayHJ6/7W2tsorV66Uh4aGZFmW5UOHDskAskpcf/LJJ2UA8i9/+Uvu7xMTE3JDQ4P8pS99yfA9TeQOiywbCBLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## joint conditional distribution \n",
|
||||
"rng = np.random.default_rng(42)\n",
|
||||
"n = 10_000_000 \n",
|
||||
" \n",
|
||||
"## circular distribution x,y dependant \n",
|
||||
"theta = rng.uniform(0, 4*np.pi, size=n)\n",
|
||||
"r = rng.uniform(0.3, 0.4) \n",
|
||||
"#r = 0.3 + 0.4 * (theta / (4*np.pi)) # 螺旋半径从 0.3 到 0.7\n",
|
||||
"x_raw = r * np.cos(theta) + 0.1 * rng.normal(size=n)\n",
|
||||
"y_raw = r * np.sin(theta) + 0.1 * rng.normal(size=n)\n",
|
||||
" \n",
|
||||
"# scale to [0, 1]\n",
|
||||
"x = (x_raw - x_raw.min()) / (x_raw.max() - x_raw.min())\n",
|
||||
"y = (y_raw - y_raw.min()) / (y_raw.max() - y_raw.min())\n",
|
||||
" \n",
|
||||
"plt.scatter(x[:5000], y[:5000], s=1, alpha=0.5)\n",
|
||||
"plt.title(\"Original Distribution (x, y)\")\n",
|
||||
"plt.gca().set_aspect('equal')\n",
|
||||
"plt.xlabel(\"x\")\n",
|
||||
"plt.ylabel(\"y\")\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "a713be97",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def plot(x, y, CDF_X, CDF_Y, u, v): \n",
|
||||
" #plot\n",
|
||||
" fig, axs = plt.subplots(1, 4, figsize=(15, 4))\n",
|
||||
" \n",
|
||||
" axs[0].scatter(x[:5000], y[:5000], s=1, alpha=0.5)\n",
|
||||
" axs[0].set_title(\"Original (x, y)\")\n",
|
||||
" axs[0].set_aspect('equal')\n",
|
||||
" axs[0].set_xlabel(\"x\")\n",
|
||||
" axs[0].set_ylabel(\"y\")\n",
|
||||
" \n",
|
||||
" axs[1].imshow(CDF_X[:,:,0], origin='lower', extent=[0,1,0,1], cmap='viridis')\n",
|
||||
" axs[1].set_title(\"normalized CDF_X (u = F_X(x))\")\n",
|
||||
" axs[1].set_xlabel(\"y\")\n",
|
||||
" axs[1].set_ylabel(\"P_X(x < x')\")\n",
|
||||
"\n",
|
||||
" axs[2].imshow(CDF_Y[:,:,0], origin='lower', extent=[0,1,0,1], cmap='viridis')\n",
|
||||
" axs[2].set_title(\"normalized conditional CDF_Y (u = F_Y|X(y|x))\")\n",
|
||||
" axs[2].set_xlabel(\"P_Y|X(y < y'|x)\")\n",
|
||||
" axs[2].set_ylabel(\"x\")\n",
|
||||
"\n",
|
||||
" axs[3].scatter(u[:5000], v[:5000], s=1, alpha=0.5, color='blue')\n",
|
||||
" axs[3].set_title(\"Mapped (u, v)\")\n",
|
||||
" axs[3].set_aspect('equal')\n",
|
||||
" axs[3].set_xlim(0,1); axs[2].set_ylim(0,1)\n",
|
||||
" axs[3].set_xlabel(\"u\")\n",
|
||||
" axs[3].set_ylabel(\"v\")\n",
|
||||
" \n",
|
||||
" plt.tight_layout()\n",
|
||||
" plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "0e87346a",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"(300, 300, 2)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"#calculate joint probability density function\n",
|
||||
"nBins = 300\n",
|
||||
"range = np.array([[x.min(),x.max()], [y.min(),y.max()]])\n",
|
||||
"bin_width = (range[:,0] - range[:,1])/nBins\n",
|
||||
"H, xedges, yedges = np.histogram2d(x, y, bins=nBins, range=range)\n",
|
||||
"H = H.astype(int) \n",
|
||||
"H = np.repeat(H[:,:,None], 2, axis=2) #add energy axis\n",
|
||||
"\n",
|
||||
"print(H.shape)\n",
|
||||
"ebins = [0,1,2] #dummy energy bins"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "d63e9295",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAABdIAAAGWCAYAAAB8cHHOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXl8FEX6xp/qniMBQlCOEAynCIjiQXBVUAFxdRFZdWXlhxeKuiAKq6yK6CIi66LiwXqAFwqsirpeqHhhBF0UlEMFFkUBOQIJkSskQJKZ6fr9UV093T3dMz2TyQXv188409XV1dXVPT3kqbefl3HOOQiCIAiCIAiCIAiCIAiCIAiCcESp6w4QBEEQBEEQBEEQBEEQBEEQRH2GhHSCIAiCIAiCIAiCIAiCIAiCiAMJ6QRBEARBEARBEARBEARBEAQRBxLSCYIgCIIgCIIgCIIgCIIgCCIOJKQTBEEQBEEQBEEQBEEQBEEQRBxISCcIgiAIgiAIgiAIgiAIgiCIOJCQThAEQRAEQRAEQRAEQRAEQRBxICGdIAiCIAiCIAiCIAiCIAiCIOJAQjpBEARBEARBEARBEARBEARBxIGEdIIgCIIgCIIgCIIgCIKo58yePRuMMTDGsHjx4pj1nHN07twZjDH069ev1vtXU8jj3rx5s6f6I0aMwB/+8Iea7ZRHfv75ZwQCAaxataquu0KkARLSCYIgCIIgCIIgCIIgCKKBkJWVhVmzZsWUf/HFF9i4cSOysrLqoFf1g++++w5z5szBP/7xj7ruCgCgS5cuuPLKK3HbbbfVdVeINEBCOkEQBEEQBEEQBEEQBEE0EIYOHYq33noL+/fvt5TPmjULZ555Jtq1a1dHPat7HnzwQfzud79Dr1696rorBrfccgu+/PJLfP3113XdFaKakJBOEARBEARBEARBEARBEA2EYcOGAQDmzZtnlJWWluKtt97CiBEjHLeZPHkyTj/9dBx99NFo2rQpevbsiVmzZoFzbqnXoUMHXHTRRXjnnXdw0kknISMjA506dcITTzxhqbd48WIwxvDyyy9j3LhxaN26NTIzM9G3b1989913MftfsWIF/vjHP+Loo49GRkYGTj31VLzxxhsx9ZYtW4Y+ffogIyMDbdq0wYQJExAKhTyNy86dO/HOO+/g6quvtpS7WcPIY3CyyYlHKBRCq1atYvYDAPv27UNmZibGjRtnlOXn5+P444/HM888k9R+iPoHCekEQRAEQRAEQRAEQRAE0UBo2rQphgwZghdffNEomzdvHhRFwdChQx232bx5M0aOHIk33ngDb7/9Nv70pz9hzJgxmDJlSkzd77//Hrfeeituu+02vPPOO+jduzf++te/4pFHHompe/fdd2PTpk144YUX8MILL2DHjh3o168fNm3aZNRZtGgR+vTpg3379uGZZ57B/Pnzccopp2Do0KGYPXu2UW/dunUYMGAA9u3bh9mzZ+OZZ57Bd99959mm5dNPP0UoFEL//v091U8Vv9+Pq666yvGpgHnz5qGiogLXXXedpbxfv3746KOPYiYuiIaFr647QBAEQRAEQRAEQRAEQRCEd0aMGIH+/fvjf//7H0444QS8+OKL+POf/+zqj/7SSy8ZnzVNQ79+/cA5x7/+9S9MnDgRjDFj/Y4dO/Ddd9/h5JNPBgAMHDgQJSUlmDJlCkaPHo1GjRoZdVu2bIl33nnH2P6ss87Ccccdh6lTp+L5558HAIwePRonnHACPv/8c/h8Qoq84IILsGvXLtx999245pproCgK7r//fnDO8fnnnyMnJwcAMGjQIJx44omexmTp0qXIzMxEt27dvA5jylx33XV4/PHH8frrr+PGG280ymfPno38/Hz06NHDUr9nz56YOXMm1q9fXyv9I2oGikgnCIIgCIIgCIIgCIIgiAZE3759ceyxx+LFF1/EmjVrsHz5cldbFwD4/PPPcd555yE7OxuqqsLv9+Pee+/F7t27UVJSYql7wgknGCK65IorrsD+/fuxatWqmHKzCN++fXv07t0bixYtAgBs2LABP/30E6688koAQDgcNl4XXnghioqKsH79egAicn3AgAGGiA4Aqqq6Rtnb2bFjB1q2bGnpT03Ro0cP5OfnWyYofvzxR3z77beO56FVq1YAgO3bt9d434iag4R0giAIgiAIgiAIgiAIgmhAMMZw3XXX4eWXX8YzzzyDLl264Oyzz3as++233+L8888HADz//PP46quvsHz5ctxzzz0AgEOHDlnqt27dOqYNWbZ7925PdWW9nTt3AgBuv/12+P1+y2v06NEAgF27dhltx9t3Ig4dOoSMjAxPddPBiBEjsHTpUvz0008ARNR/MBg0POzNyH7Zx5poWJCQThAEQRAEQRAEQRAEQRANjGuvvRa7du3CM888E+PJbea1116D3+/HBx98gMsvvxy9e/dGr169XOsXFxe7ljVv3txTXVmvRYsWAIAJEyZg+fLljq9TTjnFaDvevhPRokUL7NmzJ6ZcitiVlZWWcingp8qwYcMQDAYxe/ZsRCIR/Pvf/8Yll1yCo446Kqau7JccD6JhQkI6QRAEQRAEQRAEQRAEQTQwjjnmGNxxxx0YPHgwhg8f7lqPMQafzwdVVY2yQ4cO4d///rdj/f/973/44YcfLGWvvvoqsrKy0LNnT0v5vHnzLAk0t2zZgq+//hr9+vUDAHTt2hXHHXccfvjhB/Tq1cvxJX3d+/fvj4KCAiOKHQAikQhef/11T+PRrVs37N69G6WlpZbyDh06AABWr15tKX/vvfc8tevGUUcdhUsuuQRz587FBx98gOLiYld7nU2bNkFRFHTt2rVa+yTqFhLSCYIgCIIgCIIgCIIgCKIB8uCDD+Ldd99Fbm6ua51BgwahvLwcV1xxBRYuXIjXXnsNZ599NoLBoGP9Nm3a4I9//CNeeuklfPzxx7jqqquwcOFC/P3vf7ckGgWAkpISXHrppViwYAFeffVVnHfeecjIyMCECROMOs8++ywKCgpwwQUXYN68efjyyy/x7rvvYurUqfjzn/9s1Pv73/8OADj33HPx+uuv4/3338egQYNw4MABT2MhE6h+8803lvLTTjsNXbt2xe2334558+bh448/xsiRI7FkyZKYNhYvXgzGGO677z5P+xwxYgSKiopwyy23IC8vD+edd55jvWXLluGUU05xjFYnGg4kpBMEQRAEQRAEQRAEQRDEYcq5555rJCUdPHgw7rnnHgwZMgR33XWXY/1TTjkFjz32GB599FFcfPHF+Oqrr/DYY4/hzjvvjKn7z3/+E+3bt8d1112HESNGIDc3F4sWLcKxxx5r1Onfvz++/fZbNGvWDLfeeivOO+883HTTTfjss88swvOJJ56Izz77DE2bNsXw4cPxl7/8BSeddBImTpzo6Tj79OmDDh06YP78+ZZyVVXx/vvvo1u3bhg1ahSuueYaBINBPPXUUzFtlJeXA0DciQkz5513Htq2bYvCwkIMHz4cihIrtZaXl6OgoMBIuEo0XBg3P39BEARBEARBEARBEARBEMQRSYcOHXDiiSfigw8+iFtv8eLF6N+/P/7zn/9gyJAhtdS7xDz66KN44IEHsH37dmRmZia9/Z133ol58+bhl19+SVvi0lmzZuGvf/0rtm3bRhHpDRyKSCcIgiAIgiAIgiAIgiAIosFz8803Izs7G08//XRK2y9atAgTJ05Mm4geDofx0EMPYcKECSSiHwb46roDBEEQBEEQBEEQBEEQBEEQ1SUjIwP//ve/8d1336W0/fLly9Pan23btuGqq67C3/72t7S2S9QNFJFOEARBEIcx1157rZGlXtKhQwdce+21tdqPzZs3gzGG2bNne6q/adMm3HLLLejSpQsyMzPRqFEjnHDCCfj73/+O7du3G/WuvfZaMMaMV+PGjdGhQwcjOVJlZWVM2/369bNsY36tXbvWU//Wr1+PRo0a4YorrohZt3fvXhxzzDE4/fTTEYlEErY1d+5ctGzZEmVlZZ72Xdu4jVWLFi08tzF37lwwxvDcc8/FrPv666+hqipuv/12S3koFEK3bt3w4IMPJt3niRMnomfPntA0zVK+d+9eNGvWDO+++27SbRJ1T0O9n9UXnMbPnkxt3bp1uO+++7B582ZP29clySSC279/Px544AH06tULTZs2RTAYRIcOHTBixAisWrXKqDd79mzLfS4jIwOtW7dG//79MXXqVJSUlMS0fd9997neJ528d+3s27cPeXl5rr8ZS5YsgaqqlqR9bhzuvyfTpk0DYwzvvfee4/oLLrgARx99NHbs2AFAnBv7NTtixAgEg0GsWbMmZvsHH3wQjDG8//77lvKNGzciGAxi6dKlnvppRiYtdPpOJaKgoABNmjSx/LtHcs455+DWW29Nuk2C8MLmzZsT2roA0cSe9cnWRXLWWWdhzJgxdd0NAEDHjh1x7733pi3CnahjOEHUAkuXLuVDhgzhrVu35n6/n+fk5PDLLruMf/3110m1M2nSJJ7qZbto0SIOgC9atCil7b3St29f3rdvX091q6qqeNeuXfnUqVNrtE9e2LNnD8/OzubvvPNOXXeFIIg0Mnz4cN6+fXtL2apVq/iGDRtqtR+//vorB8BfeumlhHXff/993rhxY96+fXs+bdo0/tlnn/GCggI+ffp0ftJJJ/FTTjnFqDt8+HCemZnJly5dypcuXco///xzPmfOHP5///d/XFVVfsIJJ/Bt27ZZ2u/bty/v1KmTsY35dfDgQc/H9Pjjj3MA/M0337SU/9///R/PyMjgP/74Y8I2Dhw4wI855hg+bdo0z/utbQDwIUOGxIzVihUrkmrn4osv5k2aNOG//vqrUVZeXs47d+7Mu3Xrxg8dOmSpP336dN6qVSteXl6edJ/37dvHmzVrxl988cWYdffddx/v3Lkzr6ysTLpdom5piPez+oTT+C1dutRyj/zPf/7j+u/lDRs28FLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1500x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## check that interpolator transforms to uniform distribution \n",
|
||||
"interpolator = Interpolator(H, xedges[:-1], yedges[:-1], ebins[:-1]) #dummy energy bins\n",
|
||||
"CDF_X = interpolator.get_ietax() #lookup table x \n",
|
||||
"CDF_Y = interpolator.get_ietay() #lookup table y \n",
|
||||
"\n",
|
||||
"x_indices = ((x - range[0,0])/bin_width[0]).astype(int)\n",
|
||||
"y_indices = ((y - range[1,0])/bin_width[1]).astype(int)\n",
|
||||
"\n",
|
||||
"u = CDF_X[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"v = CDF_Y[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"\n",
|
||||
"plot(x,y,CDF_X,CDF_Y, u,v)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "85823394",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"shape X_tab: (301, 300, 2)\n",
|
||||
"shape Y_tab: (300, 301, 2)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAABdIAAAGWCAYAAAB8cHHOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXlg1GT+xp83maOFliJXKZZLERDxAlwRVEBcFZFVV5QfqKCoC+KxLh6ILoKyLt6yHuCFAh54rqLihQi6KCgUFRFFAREKLZVCSwu0nZm8vz/evJkkk8xkplNa4PvZHWcmefNeSVP6vN88X8Y55yAIgiAIgiAIgiAIgiAIgiAIwhGlvjtAEARBEARBEARBEARBEARBEA0ZEtIJgiAIgiAIgiAIgiAIgiAIIg4kpBMEQRAEQRAEQRAEQRAEQRBEHEhIJwiCIAiCIAiCIAiCIAiCIIg4kJBOEARBEARBEARBEARBEARBEHEgIZ0gCIIgCIIgCIIgCIIgCIIg4kBCOkEQBEEQBEEQBEEQBEEQBEHEgYR0giAIgiAIgiAIgiAIgiAIgogDCekEQRAEQRAEQRAEQRAEQRAEEQcS0gmCIAiCIAiCIAiCIAiigTN79mwwxsAYw5IlS2L2c87RqVMnMMbQv3///d6/ukKOe9OmTZ7Kjx49Guecc07ddsojv/zyCwKBAFatWlXfXSHSAAnpBEEQBEEQBEEQBEEQBHGAkJ2djVmzZsVs//zzz7FhwwZkZ2fXQ68aBt9++y3mzJmDf/3rX/XdFQBA586dcemll+If//hHfXeFSAMkpBMEQRAEQRAEQRAEQRDEAcKwYcPw1ltvYffu3Zbts2bNwimnnIJ27drVU8/qn/vuuw9/+tOf0KtXr/ruisH111+PL774Al999VV9d4WoJSSkEwRBEARBEARBEARBEMQBwvDhwwEA8+bNM7aVl5fjrbfewujRox2Pufvuu3HyySejWbNmaNKkCXr06IFZs2aBc24p16FDB5x33nl4++23cdxxxyEjIwNHHHEEHnvsMUu5JUuWgDGGl156CePHj0fr1q2RmZmJfv364dtvv41pf+XKlfjLX/6CZs2aISMjAyeeeCJef/31mHLLly9H3759kZGRgTZt2mDixIkIhUKe5mX79u14++23cfnll1u2u1nDyDE42eTEIxQKoVWrVjHtAEBZWRkyMzMxfvx4Y1vPnj1x9NFH46mnnkqqHaLhQUI6QRAEQRAEQRAEQRAEQRwgNGnSBEOHDsXzzz9vbJs3bx4URcGwYcMcj9m0aRPGjBmD119/Hf/973/x17/+FTfccAOmTp0aU/a7777DTTfdhH/84x94++230adPH/z973/HQw89FFP2jjvuwMaNG/Hcc8/hueeew7Zt29C/f39s3LjRKLN48WL07dsXZWVleOqppzB//nyccMIJGDZsGGbPnm2UW7t2LQYOHIiysjLMnj0bTz31FL799lvPNi2ffPIJQqEQBgwY4Kl8qvj9flx22WWOTwXMmzcPVVVVuPLKKy3b+/fvjw8//DBm4YI4sPDVdwcIgiAIgiAIgiAIgiAIgvDO6NGjMWDAAPz444845phj8Pzzz+Piiy929Ud/4YUXjM+apqF///7gnOM///kPJk2aBMaYsX/btm349ttvcfzxxwMABg0ahJKSEkydOhXjxo1Do0aNjLItW7bE22+/bRx/6qmn4qijjsK0adPw7LPPAgDGjRuHY445Bp999hl8PiFFnn322dixYwfuuOMOjBw5Eoqi4J577gHnHJ999hlyc3MBAIMHD0b37t09zcmyZcuQmZmJrl27ep3GlLnyyivx6KOP4rXXXsM111xjbJ89ezZ69uyJY4891lK+R48emDlzJtatW7df+kfUDRSRThAEQRAEQRAEQRAEQRAHEP369cORRx6J559/Hj/88ANWrFjhausCAJ999hnOPPNM5OTkQFVV+P1+3HXXXSgtLUVJSYml7DHHHGOI6JIRI0Zg9+7dWLVqVcx2swjfvn179OnTB4sXLwYArF+/Hj///DMuvfRSAEA4HDZe5557LoqKirBu3ToAInJ94MCBhogOAKqqukbZ29m2bRtatmxp6U9dceyxx6Jnz56WBYqffvoJ33zzjeN5aNWqFQBg69atdd43ou4gIZ0gCIIgCIIgCIIgCIIgDiAYY7jyyivx0ksv4amnnkLnzp1x2mmnOZb95ptvcNZZZwEAnn32WXz55ZdYsWIF7rzzTgDAvn37LOVbt24dU4fcVlpa6qmsLLd9+3YAwC233AK/3295jRs3DgCwY8cOo+54bSdi3759yMjI8FQ2HYwePRrLli3Dzz//DEBE/QeDQcPD3ozsl32uiQMLEtIJgiAIgiAIgiAIgiAI4gDjiiuuwI4dO/DUU0/FeHKbefXVV+H3+/H+++/jkksuQZ8+fdCrVy/X8sXFxa7bmjdv7qmsLNeiRQsAwMSJE7FixQrH1wknnGDUHa/tRLRo0QI7d+6M2S5F7Orqast2KeCnyvDhwxEMBjF79mxEIhG8+OKLuOCCC3DYYYfFlJX9kvNBHJiQkE4QBEEQBEEQBEEQBEEQBxiHH344br31VgwZMgSjRo1yLccYg8/ng6qqxrZ9+/bhxRdfdCz/448/4vvvv7dse+WVV5CdnY0ePXpYts+bN8+SQPP333/HV199hf79+wMAunTpgqOOOgrff/89evXq5fiSvu4DBgzAokWLjCh2AIhEInjttdc8zUfXrl1RWlqK8vJyy/YOHToAAFavXm3Z/u6773qq143DDjsMF1xwAebOnYv3338fxcXFrvY6GzduhKIo6NKlS63aJOoXEtIJgiAIgiAIgiAIgiAI4gDkvvvuwzvvvIO8vDzXMoMHD0ZlZSVGjBiBhQsX4tVXX8Vpp52GYDDoWL5Nmzb4y1/+ghdeeAEfffQRLrvsMixcuBD//Oc/LYlGAaCkpAQXXnghFixYgFdeeQVnnnkmMjIyMHHiRKPM008/jUWLFuHss8/GvHnz8MUXX+Cdd97BtGnTcPHFFxvl/vnPfwIAzjjjDLz22mt47733MHjwYOzZs8fTXMgEql9//bVl+0knnYQuXbrglltuwbx58/DRRx9hzJgxWLp0aUwdS5YsAWMMU6ZM8dTm6NGjUVRUhOuvvx75+fk488wzHcstX74cJ5xwgmO0OnHgQEI6QRAEQRAEQRAEQRAEQRyknHHGGUZS0iFDhuDOO+/E0KFDcfvttzuWP+GEE/DII4/g4Ycfxvnnn48vv/wSjzzyCG677baYsv/+97/Rvn17XHnllRg9ejTy8vKwePFiHHnkkUaZAQMG4JtvvkHTpk1x00034cwzz8S1116LTz/91CI8d+/eHZ9++imaNGmCUaNG4W9/+xuOO+44TJo0ydM4+/btiw4dOmD+/PmW7aqq4r333kPXrl0xduxYjBw5EsFgEE888URMHZWVlQAQd2HCzJlnnom2bduisLAQo0aNgqLESq2VlZVYtGiRkXCVOHBh3Pz8BUEQBEEQBEEQBEEQBEEQhyQdOnRA9+7d8f7778ctt2TJEgwYMABvvPEGhg4dup96l5iHH34Y9957L7Zu3YrMzMykj7/tttswb948/Prrr2lLXDpr1iz8/e9/x5YtWygi/QCHItIJgiAIgiAIgiAIgiAIgjjgue6665CTk4Mnn3wypeMXL16MSZMmpU1ED4fDuP/++zFx4kQS0Q8CfPXdAYIgCIIgCIIgCIIgCIIgiNqSkZGBF198Ed9++21Kx69YsSKt/dmyZQsuu+wy3HzzzWmtl6gfKCKdIAiCIA5irrjiCiNLvaRDhw644oor9ms/Nm3aBMYYZs+e7an8xo0bcf3116Nz587IzMxEo0aNcMwxx+Cf//wntm7dapS74oorwBgzXo0bN0aHDh2M5EjV1dUxdffv399yjPm1Zs0aT/1bt24dGjVqhBEjRsTs27VrFw4//HCcfPLJiEQiCeuaO3cuWrZsiYqKCk9t72/c5qpFixae65g7dy4YY3jmmWdi9n311VdQVRW33HKLZXsoFELXrl1x3333Jd3nSZMmoUePHtA0zbJ9165daNq0Kd55552k6yTqnwP1ftZQcJo/ezK1tWvXYsqUKdi0aZOn4+uTZBLB7d69G/feey969eqFJk2aIBgMokOHDhg9ejRWrVpllJs9e7blPpeRkYHWrVtjwIABmDZtGkpKSmLqnjJliut90sl7105ZWRny8/Ndf2csXboUqqpakva5cbD/PnnwwQfBGMO7777ruP/ss89Gs2bNsG3bNgDi3Niv2dGjRyMYDOKHH36IOf6+++4DYwzvvfeeZfuGDRsQDAaxbNkyT/00I5MWOv1MJWLRokXIysqy/LtHcvrpp+Omm25Kuk6C8MKmTZsS2roA0cSeDcnWRXLqqafihhtuqO9uAAA6duyIu+66K20R7kQ9wwliP7Bs2TI+dOhQ3rp1a+73+3lubi6/6KKL+FdffZVUPZMnT+apXraLFy/mAPjixYtTOt4r/fr14/369fNUtqamhnfp0oVPmzatTvvkhZ07d/KcnBz+9ttv13dXCIJII6NGjeLt27e3bFu1ahVfv379fu3Hb7/9xgHwF154IWHZ9957jzdu3Ji3b9+eP/jgg/zTTz/lixYt4tOnT+fHHXccP+GEE4yyo0aN4pmZmXzZsmV82bJl/LPPPuNz5szh//d//8dVVeXHHHMM37Jli6X+fv368SOOOMI4xvzau3ev5zE9+uijHAB/8803Ldv/7//+j2dkZPCffvopYR179uzhhx9+OH/wwQc9t7u/AcCHDh0aM1crV65Mqp7zzz+fZ2Vl8d9++83YVllZyTt16sS7du3K9+3bZyk/ffp03qpVK15ZWZl0n8vKynjTpk35888/H7NvypQpvFOnTry6ujrpeon65UC8nzUknOZv2bJllnvkG2+84frv5fXr1/NVq1Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1500x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"x_indices = ((x - range[0,0])/bin_width[0]).astype(int)\n",
|
||||
"y_indices = ((y - range[1,0])/bin_width[1]).astype(int)\n",
|
||||
"\n",
|
||||
"CDF_X, CDF_Y = calculate_rosenblatt_lookup_table(H, xedges[:-1], yedges[:-1],ebins[:-1])\n",
|
||||
"\n",
|
||||
"u = CDF_X[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"v = CDF_Y[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
" \n",
|
||||
"plot(x,y,CDF_X,CDF_Y, u,v)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "fe9e3306",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"worked: \n",
|
||||
"could normalize CDF eta_x\n",
|
||||
"construction worked\n",
|
||||
"calculation of conditional CDF etay worked\n",
|
||||
"assignment m_ietay worked\n",
|
||||
"assignment m_ietax worked\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAABcUAAAGPCAYAAACd27gMAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXl8FFX69U91J+mEJARISAiE1QgomwIOKDKAKJgBRAeRF0YBUUZkxFFEkXHcYBQFBnEDXFhkVFxwwQUXQNEfAgJBQUDQgCydlQTSJJDO0n3fP6pv9a3qqurqTkICeb5+MN1Vt+69dbtTSU6dPo/EGGMgCIIgCIIgCIIgCIIgCIIgiAaAra4nQBAEQRAEQRAEQRAEQRAEQRDnChLFCYIgCIIgCIIgCIIgCIIgiAYDieIEQRAEQRAEQRAEQRAEQRBEg4FEcYIgCIIgCIIgCIIgCIIgCKLBQKI4QRAEQRAEQRAEQRAEQRAE0WAgUZwgCIIgCIIgCIIgCIIgCIJoMJAoThAEQRAEQRAEQRAEQRAEQTQYSBQnCIIgCIIgCIIgCIIgCIIgGgwkihMEQRAEQRAEQRAEQRAEQRANBhLFCYIgCIIgCIIgCIIgCKKes3LlSkiSBEmSsGnTpoD9jDGkp6dDkiQMHDjwnM+vtuDnfeTIEUvtJ02ahOuvv752J2WR3377DVFRUdi1a1ddT4XQQKI4QRAEQRAEQRAEQRAEQZwnxMfHY9myZQHbv/vuOxw6dAjx8fF1MKv6wU8//YQ33ngD//nPf+p6KgCAjh074m9/+xvuv//+up4KoYFEcYIgCIIgCIIgCIIgCII4TxgzZgw++OADnD59WrV92bJluPLKK9GmTZs6mlnd88wzz+BPf/oTevfuXddTUbjnnnvw/fffY8uWLXU9FUKARHGCIAiCIAiCIAiCIAiCOE8YO3YsAGD16tXKNpfLhQ8++ACTJk3SPebJJ59Enz590KxZMzRu3Bg9e/bEsmXLwBhTtWvXrh2GDx+Ojz76CN27d0d0dDQ6dOiAF154QdVu06ZNkCQJb775JqZPn44WLVogJiYGAwYMwE8//RQw/s6dO3HDDTegWbNmiI6OxuWXX4733nsvoN22bdvQr18/REdHo2XLlpg1axYqKystrUt+fj4++ugj3HbbbartRvEr/Bz0omjMqKysRHJycsA4AFBcXIyYmBhMnz5d2darVy9ccsklWLp0aUjjELULieIEQRAEQRAEQRAEQRAEcZ7QuHFj3HzzzVi+fLmybfXq1bDZbBgzZozuMUeOHMFdd92F9957Dx9++CH++te/Ytq0aZgzZ05A259//hn33Xcf7r//fnz00Ue46qqr8M9//hMLFiwIaPuvf/0Lhw8fxuuvv47XX38dOTk5GDhwIA4fPqy0+fbbb9GvXz8UFxdj6dKlWLt2LS677DKMGTMGK1euVNrt378fgwcPRnFxMVauXImlS5fip59+shyF8vXXX6OyshKDBg2y1D5cIiMjceutt+q69VevXg23243bb79dtX3gwIH44osvAm5CEHVHRF1PgCAIgiAIgiAIgiAIgiAI60yaNAmDBg3Cvn370KVLFyxfvhyjR482zBNfsWKF8tjr9WLgwIFgjOH555/Ho48+CkmSlP05OTn46aef0KNHDwBARkYGCgoKMGfOHEydOhWNGjVS2jZv3hwfffSRcvzVV1+Niy++GHPnzsVrr70GAJg6dSq6dOmCb775BhERshQ5dOhQFBYW4l//+hfGjx8Pm82G2bNngzGGb775BikpKQCAYcOGoWvXrpbWZOvWrYiJiUHnzp2tLmPY3H777Xjuuefw7rvvYvLkycr2lStXolevXujWrZuqfc+ePbFkyRIcPHjwnMyPCA45xQmCIAiCIAiCIAiCIAjiPGLAgAG46KKLsHz5cvzyyy/YsWOHYXQKAHzzzTe49tprkZCQALvdjsjISDz22GMoKipCQUGBqm2XLl0UQZwzbtw4nD59Grt27QrYLgrqbdu2xVVXXYVvv/0WAJCVlYUDBw7gb3/7GwCgqqpK+feXv/wFubm5OHjwIADZUT548GBFEAcAu91u6H7XkpOTg+bNm6vmU1t069YNvXr1Ut1s+PXXX7F9+3bd1yE5ORkAkJ2dXetzI6xBojhBEARBEARBEARBEARBnEdIkoTbb78db775JpYuXYqOHTuif//+um23b9+OIUOGAABee+01/PDDD9ixYwceeeQRAEBZWZmqfYsWLQL64NuKioosteXt8vPzAQAzZsxAZGSk6t/UqVMBAIWFhUrfZmMHo6ysDNHR0Zba1gSTJk3C1q1bceDAAQCyG9/hcCiZ7yJ8Xtq1JuoOEsUJgiAIgiAIgiAIgiAI4jxj4sSJKCwsxNKlSwMyrEXeeecdREZG4rPPPsMtt9yCq666Cr179zZsn5eXZ7gtMTHRUlveLikpCQAwa9Ys7NixQ/ffZZddpvRtNnYwkpKScPLkyYDtXJAuLy9XbedifLiMHTsWDocDK1euhMfjwf/+9z/ceOONaNq0aUBbPi++HkTdQ6I4QRAEQRAEQRAEQRAEQZxntGrVCg8++CBGjBiBCRMmGLaTJAkRERGw2+3KtrKyMvzvf//Tbb9v3z7s3r1bte3tt99GfHw8evbsqdq+evVqVfHIo0ePYsuWLRg4cCAAoFOnTrj44ouxe/du9O7dW/cfz0EfNGgQNm7cqLjLAcDj8eDdd9+1tB6dO3dGUVERXC6Xanu7du0AAHv27FFt/+STTyz1a0TTpk1x4403YtWqVfjss8+Ql5dnGGFz+PBh2Gw2dOrUqVpjEjUHieIEQRAEQRAEQRAEQRAEcR7yzDPP4OOPP0Zqaqphm2HDhqG0tBTjxo3D+vXr8c4776B///5wOBy67Vu2bIkbbrgBK1aswJdffolbb70V69evx7///W9VkU0AKCgowE033YTPP/8cb7/9Nq699lpER0dj1qxZSptXXnkFGzduxNChQ7F69Wp8//33+PjjjzF37lyMHj1aaffvf/8bAHDNNdfg3Xffxaeffophw4bhzJkzltaCFw/98ccfVduvuOIKdOrUCTNmzMDq1avx5Zdf4q677sLmzZsD+ti0aRMkScITTzxhacxJkyYhNzcX99xzD9LS0nDttdfqttu2bRsuu+wyXRc5UTeQKE4QBEEQFzATJ05UnBGcdu3aYeLEied0HkeOHIEkSVi5cqWl9ocPH8Y999yDjh07IiYmBo0aNUKXLl3w73//W1WcZuLEiZAkSfkXGxuLdu3aKb/Eaz8iCci/LIvHiP/27t1raX4HDx5Eo0aNMG7cuIB9p06dQqtWrdCnTx94PJ6gfa1atQrNmzdHSUmJpbHPNUZrFcpHP1etWgVJkvDqq68G7NuyZQvsdjtmzJih2l5ZWYnOnTvjmWeeCXnOjz76KHr27Amv16vafurUKTRp0gQff/xxyH0Sdc/5ej2rL+itn/aP/v379+OJJ57AkSNHLB1fl4QiWJw+fRpPPfUUevfujcaNG8PhcKBdu3aYNGmSqmDcypUrVde56OhotGjRAoMGDcLcuXMDCtEBwBNPPGF4nXzppZeCzq24uBhpaWmGPzM2b94Mu92uEpeMuNB/nsyfPx+SJBk6O4cOHYpmzZohJycHgPzaaN+zkyZNgsPhwC+//BJw/DPPPANJkvDpp5+qth86dAgOhwNbt261NE8RLq7pfU8FY+PGjYiLi9MtyvfnP/8Z9913X8h9Eg2Xa665RinIOWLECDzyyCO4+eab8fDDD+u2v+yyy7Bw4UL897//xciRI/HDDz9g4cKFeOihhwLaPv3002jbti1uv/12TJo0Campqfj2229x0UUXKW0GDRqE7du3o0mTJrjvvvtw7bXX4u6778aGDRtUInLXrl2xYcMGNG7cGBMmTMDf//53dO/eHY8++qil8+zXrx/atWuHtWvXqrbb7XZ8+umn6Ny5M6ZMmYLx48fD4XDoXqdLS0sBwPQmg8i1116L1q1bw+l0YsKECbDZAqXW0tJSbNy4USk2StQTGEEQBEEQFywTJkxgbdu2VW3btWsXy8rKOqfz+OOPPxgAtmLFiqBtP/30UxYbG8vatm3L5s+fzzZs2MA2btzIFi1axLp3784uu+wype2ECRNYTEwM27p1K9u6dSv75ptv2BtvvMH+3//7f8xut7MuXbqw48ePq/ofMGAA69Chg3KM+O/s2bOWz+m5555jANiaNWtU2//f//t/LDo6mv36669B+zhz5gxr1aoVmz9/vuVxzzUA2M033xywVjt37gypn5EjR7K4uDj2xx9/KNtKS0tZeno669y5MysrK1O1X7RoEUtOTmalpaUhz7m4uJg1adKELV++PGDfE088wdLT01l5eXnI/RJ1y/l4PatP6K3f1q1bVdfI999/nwFg3377bcDxWVlZbNeuXbU8S+sAYI8//njQdllZWaxDhw4sLi6OzZgxg3322Wds06ZNbOXKlewvf/kLA8CKi4sZY4ytWLFCeW23bt3Kvv/+e7ZmzRp23333sYSEBNasWTO2fv16Vf+PP/44A8C+/PLLgOtkXl6epXP56quvGAD21FNPqbafOXOGpaens27dugW9ZjWEnycej4ddffXVrEWLFqyoqEi175VXXmEA2OrVq5Vtjz/+eMB73uVysTZt2rDLL7+cVVRUKNv37NnDoqKi2MSJEwPGvfHGG9mwYcNCOFM/3377LQOg+tkXCoMGDWLjx48P2L5p0yYWGRnJDhw4EFa/BGFG27ZtLb3n+fv7/fffPwezss6CBQtY06ZNQ/q9XuTBBx9kaWlpAb+bVofXX3+dxcbGspMnT9ZYn0T1IVGcOCds3bqV3XzzzaxFixYsMjKSpaSksFGjRrEtW7aE1A//pTMc+AVb75f8mmTAgAFswIABltpWVFSwTp06sblz59bqnKxw8uRJlpCQwD766KO6ngpBXNLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1500x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## check that cpp rosenblatt transforms to uniform distribution \n",
|
||||
"interpolator = Interpolator(xedges[:-1], yedges[:-1], ebins[:-1]) #dummy energy bins\n",
|
||||
"interpolator.rosenblatttransform(H)\n",
|
||||
"CDF_X = interpolator.get_ietax() #lookup table x \n",
|
||||
"CDF_Y = interpolator.get_ietay() #lookup table y \n",
|
||||
"\n",
|
||||
"x_indices = ((x - range[0,0])/bin_width[0]).astype(int)\n",
|
||||
"y_indices = ((y - range[1,0])/bin_width[1]).astype(int)\n",
|
||||
"\n",
|
||||
"u = CDF_X[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"v = CDF_Y[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"\n",
|
||||
"plot(x,y,CDF_X,CDF_Y, u,v)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6df08eb9",
|
||||
"metadata": {},
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a9d088d2",
|
||||
"metadata": {},
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "40a1ffff",
|
||||
"metadata": {},
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "1e0aa9dc",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "base",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.11.11"
|
||||
}
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAHFCAYAAACJnUzWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAA6OdJREFUeJzs/Xl4HNWdLo6/rd60tVuyZK1ty4BsGSxkBzExi4kdSJwhJCgEE25IyALkhvGQDOEmM2HyZU2exzMJYZw7iUmYhBACNyE2mYhfFkDEmInCkiDjRXZsLBvb6la35G5JpW5Jvap+f7RO6dTpc6qqW63Nqvd5eISrq+qcOnXqsy8WWZZlmDBhwoQJE+c4CuZ6AiZMmDBhwsRswGR4JkyYMGFiUcBkeCZMmDBhYlHAZHgmTJgwYWJRwGR4JkyYMGFiUcBkeCZMmDBhYlHAZHgmTJgwYWJRwGR4JkyYMGFiUcBkeCZMmDBhYlHAZHgmcsYbb7yBm266CbW1tXA4HKipqcHWrVvx+uuvZ3WfBx98EBaLJac57N27FxaLBXv37s3peqPYvHkzNm/ebOg8i8UCi8WCgoICuFwuNDY24qabbsLu3bsxMTGRcc3KlSvxuc99Lqv5vPbaa3jwwQcxPDyc1XXsWGT9du/endV9tDA2NoYHH3yQ+06efPJJWCwWnDp1Km/jZYNEIoE1a9bg3/7t3+ZkfBpDQ0MoKyvDb37zm7meyuKBbMJEDvi///f/ygUFBfJll10mP/XUU/Krr74q//znP5cvu+wyuaCgQP7P//xPw/fq7e2VX3/99ZzmIUmS/Prrr8uSJOV0vVFs2rRJ3rRpk6Hzzj//fPn111+XX3/9dfnll1+W/+u//ku+7rrrZADyVVddJQ8PD6uu2bdvn9zT05PVfL7zne/IAOR33303q+vYsV555RUZgLxr166s7qOFs2fPygDkBx54IOO3gYEB+fXXX5ej0WjexssGO3bskKuqquRIJDIn47N48MEH5cbGRjkWi831VBYFTIZnImt0dnbKBQUF8kc+8hE5kUiofkskEvJHPvIRuaCgQO7s7NS8z+jo6ExOM6/IhuGtXbuW+9sTTzwhA5A/8YlPTHs+2TK8sbEx7vHZZnhziUQiIdfX18tf//rX53oqCgKBgGyz2eRnnnlmrqeyKGAyPBNZ47rrrpOtVqvc29vL/f3MmTOy1WqVP/KRjyjHHnjgARmA3NXVJd94441yWVmZXFNTo/qNRjQale+55x65urpaLioqkq+66ir5rbfekhsaGuTPfvazynmEYL/yyivKsc9+9rNySUmJfPz4cfnaa6+VS0pKZI/HI99zzz0ZmsWDDz4ov/e975XLy8tll8slv+c975F//OMfyxMTE6rz8sHwZFmWP/zhD8sWi0U+deqUcox9plQqJX/zm9+UV69eLRcWFsput1u++OKL5R07dqjWi/2PrEFDQ4N83XXXyc8995y8fv162el0yv/yL//CHYus389//nP5K1/5ilxdXS0XFhbK73vf++R9+/YZWoPPfvazckNDgyzLsvzuu+9y50bG/OlPf8pl1D/5yU/klpYW2el0yuXl5fLHPvYx+ciRIxnjGH2vPDz33HMyAPnw4cPKsfHxcXn9+vXyBRdcoNK8/X6/XF1dLW/atElOJpO696Zx9dVXy01NTRl7aGJiQr7gggvkD3/4w6rj1157rXzVVVdlNYaJ3GD68ExkhVQqhVdeeQWXXnopPB4P95zly5ejtbUVe/bsQSqVUv328Y9/HI2Njdi1axd++MMfCsf5/Oc/jx07duDzn/882tvbceONN+KGG24w7LNKJBK4/vrrcc0116C9vR233XYb/uM//gP//u//rjrv1KlT+OIXv4hf/epX+PWvf42Pf/zj+NKXvoRvfvObhsbJFtdffz1kWcaf/vQn4Tnf/va38eCDD+KTn/wkfve73+HZZ5/F7bffrjz7HXfcgS996UsAgF//+td4/fXX8frrr+OSSy5R7rFv3z587Wtfw5e//GW88MILuPHGGzXn9a//+q84efIkfvzjH+PHP/4x+vr6sHnzZpw8eTKr56utrcULL7wAALj99tuVud13333Ca7Zv347bb78da9euxa9//Wt873vfw8GDB3H55Zfj+PHjqnONvlcefve736GqqgoXXXSRcqywsBC/+tWvMDAwgNtuuw0AMDExgU996lOQZRm/+MUvYLVas1qDf/qnf8KxY8fwxz/+UXX8D3/4A06cOIF//Md/VB3fvHkz/vznP2ftjzWRA+aa45pYWAgEAjIA+X/9r/+led7NN98sA5D7+/tlWZ7SSu6///6Mc1kN7/DhwzIARSsh+MUvfqHSFmRZrOEBkH/1q1+prv/whz8sNzU1CeecSqXkRCIhP/zww3JFRYVKQs+XhveHP/xBBiD/+7//u3KM1bo+8pGPyOvXr9ccR8uk2dDQIFutVvnYsWPc33jrd8kll6ie99SpU7LdbpfvuOMO1bPpaXiyrG3SZDW8oaEhuaioKEPrOXPmjOx0OuVbbrlFNU4u75XgwgsvlP/+7/+e+9uzzz4rA5B37Ngh33///XJBQYH80ksv6d6Th1QqJZ9//vlyW1ub6vi1114rX3DBBRmaX0dHhwxA/sMf/pDTeCaMw9TwTMwI5Mk2i2z0pZ6mAQCvvvoqAOATn/iE6vjWrVths9kMjW+xWPDRj35UdaylpQWnT59WHduzZw8+8IEPwO12w2q1wm634/7770coFMLAwIChsbKBbKD95Hvf+14cOHAA27Ztw4svvoiRkZGsx2lpacHq1asNn3/LLbeo3lVDQwOuuOIKvPLKK1mPnQ1ef/11jI+PZ0SpLl++HFdffXWGlmT0vfLQ19eHqqoq7m+f+MQn8A//8A/42te+hm9961v413/9V3zwgx/M7mEmUVBQgLvuugu//e1vcebMGQDAiRMn8MILL2Dbtm0Z3wSZk8/ny2k8E8ZhMjwTWaGyshLFxcV49913Nc87deoUiouLsXTpUtXx2tpa3TFCoRAAoLq6WnXcZrOhoqLC0DyLi4tRWFioOuZ0OhGNRpV//+Uvf8GWLVsAAP/1X/+FP//5z/jrX/+Kb3zjGwCA8fFxQ2NlA0KY6+rqhOfce++9eOSRR/DGG2/g2muvRUVFBa655hq89dZbhscxss40ampquMfIu5gpkPvz5ltXV5cxvpH3KsL4+HjGtTRuu+02JBIJ2Gw2fPnLXzYyfc17FRUVKWb7H/zgBygqKlLMpjTInGZiv5lQw2R4JrKC1WrF+9//frz11lvwer3cc7xeL7q6unD11Vdn+D+M5NsRptbf3686nkwm80qAf/nLX8Jut+O3v/0tPvGJT+CKK67ApZdemrf78/D888/DYrHgfe97n/Acm82Ge+65B/v27cPg4CB+8YtfoLe3Fx/60IcwNjZmaJxs8xoDgQD3GC1gFBYWIhaLZZwXDAazGosGub/f78/4ra+vD5WVlTnfm0VlZSUGBwe5v42OjuLWW2/F6tWrUVRUhDvuuGNaY7ndbnz2s5/Fj3/8YwwODuKnP/0pbrnlFpSVlWWcS+aUz2c1wYfJ8ExkjXvvvReyLGPbtm0ZQSmpVAr/8A//AFmWce+99+Z0f8IMnn32WdXx3bt3I5lM5jZpDiwWC2w2m4opj4+P4+c//3nexqDx05/+FH/4wx/wyU9+EitWrDB0TVlZGbZu3Yp//Md/xODgoJKw7XQ6lfnmA7/4xS9U5tbTp0/jtddeUyXbr1y5Eu+8846K6YVCIbz22muqe2Uzt8svvxxFRUV4+umnVce9Xi/27NmDa665JpfH4WLNmjU4ceIE97c777wTZ86cwa9//Wv85Cc/wfPPP4//+I//mNZ4X/7ylxEMBrF161YMDw/jrrvu4p5HAoPoYBoTMwNjDhETJihceeWV2LFjB+6++25s3LgRd911F1asWIEzZ87gBz/4Ad58803s2LEDV1xxRU73X7t2LT75yU/iu9/9LqxWK66++mocPnwY3/3ud+F2u1FQkB857brrrsOjjz6KW265Bf/7f/9vhEIhPPLIIwrBzhXj4+N44403lP8/efIkfvOb3+C3v/0tNm3apBmdCgAf/ehH0dzcjEsvvRTLli3D6dOnsWPHDjQ0NGDVqlUAgIsvvhgA8L3vfQ+f/exnYbfb0dTUBJfLldOcBwYGcMMNN+ALX/gCJEnCAw88gMLCQpXQcuutt+JHP/oRPv3pT+MLX/gCQqEQvv3tb2PJkiWqe7lcLjQ0NKC9vR3XXHMNli5disrKSqxcuTJj3LKyMtx3333413/9V3zmM5/BJz/5SYRCITz00EMoLCzEAw88kNPz8LB582Y8/PDDGBsbQ3FxsXL8xz/+MZ5++mn89Kc/xdq1a7F27Vrcdddd+Jd/+RdceeWVeO973wsgbaY/77zz8NnPfhZPPvmk7nirV6/G3//93+MPf/gDNm7ciHXr1nHPe+ONN1BRUaG8UxMziDkNmTGxoPH666/LW7dulaurq2WbzSZXVVXJH//4x+XXXnst41wSiXn27FnhbzRIHl5VVZVcWFgoX3bZZfLrr78uu91u+Stf+YpynlYenpFxnnjiCbmpqUl2Op3y+eefL2/fvl3+yU9+khEBmU2UJqj8s5KSEvn888+Xt27dKu/atUtOpVIZ17CRk9/97nflK664Qq6srJQdDoe8YsUK+fbbb1fl7smyLN97771yXV2dXFBQwM3D40ErD+/LX/6yvGzZMtnpdCp5jyx+9rOfyRdeeKFcWFgoX3TRRfKzzz6bEaUpy7L88ssvy+95z3tkp9NpKA/vxz/+sdzS0iI7HA7Z7XbLbW1tqnw5Wc7uvfLQ09MjWywWVZTnwYMH5aKiItWayHJ6/7W2tsorV66Uh4aGZFmW5UOHDskAskpcf/LJJ2UA8i9/+Uvu7xMTE3JDQ4P8pS99yfA9TeQOiywbCBszYWIe4LXXXsOLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## joint conditional distribution \n",
|
||||
"rng = np.random.default_rng(42)\n",
|
||||
"n = 10_000_000 \n",
|
||||
" \n",
|
||||
"## circular distribution x,y dependant \n",
|
||||
"theta = rng.uniform(0, 4*np.pi, size=n)\n",
|
||||
"r = rng.uniform(0.3, 0.4) \n",
|
||||
"#r = 0.3 + 0.4 * (theta / (4*np.pi)) # 螺旋半径从 0.3 到 0.7\n",
|
||||
"x_raw = r * np.cos(theta) + 0.1 * rng.normal(size=n)\n",
|
||||
"y_raw = r * np.sin(theta) + 0.1 * rng.normal(size=n)\n",
|
||||
" \n",
|
||||
"# scale to [0, 1]\n",
|
||||
"x = (x_raw - x_raw.min()) / (x_raw.max() - x_raw.min())\n",
|
||||
"y = (y_raw - y_raw.min()) / (y_raw.max() - y_raw.min())\n",
|
||||
" \n",
|
||||
"plt.scatter(x[:5000], y[:5000], s=1, alpha=0.5)\n",
|
||||
"plt.title(\"Original Distribution (x, y)\")\n",
|
||||
"plt.gca().set_aspect('equal')\n",
|
||||
"plt.xlabel(\"x\")\n",
|
||||
"plt.ylabel(\"y\")\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "a713be97",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def plot(x, y, CDF_X, CDF_Y, u, v): \n",
|
||||
" #plot\n",
|
||||
" fig, axs = plt.subplots(1, 4, figsize=(15, 4))\n",
|
||||
" \n",
|
||||
" axs[0].scatter(x[:5000], y[:5000], s=1, alpha=0.5)\n",
|
||||
" axs[0].set_title(\"Original (x, y)\")\n",
|
||||
" axs[0].set_aspect('equal')\n",
|
||||
" axs[0].set_xlabel(\"x\")\n",
|
||||
" axs[0].set_ylabel(\"y\")\n",
|
||||
" \n",
|
||||
" axs[1].imshow(CDF_X[:,:,0], origin='lower', extent=[0,1,0,1], cmap='viridis')\n",
|
||||
" axs[1].set_title(\"normalized CDF_X (u = F_X(x))\")\n",
|
||||
" axs[1].set_xlabel(\"y\")\n",
|
||||
" axs[1].set_ylabel(\"P_X(x < x')\")\n",
|
||||
"\n",
|
||||
" axs[2].imshow(CDF_Y[:,:,0], origin='lower', extent=[0,1,0,1], cmap='viridis')\n",
|
||||
" axs[2].set_title(\"normalized conditional CDF_Y (u = F_Y|X(y|x))\")\n",
|
||||
" axs[2].set_xlabel(\"P_Y|X(y < y'|x)\")\n",
|
||||
" axs[2].set_ylabel(\"x\")\n",
|
||||
"\n",
|
||||
" axs[3].scatter(u[:5000], v[:5000], s=1, alpha=0.5, color='blue')\n",
|
||||
" axs[3].set_title(\"Mapped (u, v)\")\n",
|
||||
" axs[3].set_aspect('equal')\n",
|
||||
" axs[3].set_xlim(0,1); axs[2].set_ylim(0,1)\n",
|
||||
" axs[3].set_xlabel(\"u\")\n",
|
||||
" axs[3].set_ylabel(\"v\")\n",
|
||||
" \n",
|
||||
" plt.tight_layout()\n",
|
||||
" plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "0e87346a",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"(300, 300, 2)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"#calculate joint probability density function\n",
|
||||
"nBins = 300\n",
|
||||
"range = np.array([[x.min(),x.max()], [y.min(),y.max()]])\n",
|
||||
"bin_width = (range[:,0] - range[:,1])/nBins\n",
|
||||
"H, xedges, yedges = np.histogram2d(x, y, bins=nBins, range=range)\n",
|
||||
"H = H.astype(int) \n",
|
||||
"H = np.repeat(H[:,:,None], 2, axis=2) #add energy axis\n",
|
||||
"\n",
|
||||
"print(H.shape)\n",
|
||||
"ebins = [0,1,2] #dummy energy bins"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "d63e9295",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABdIAAAGWCAYAAAB8cHHOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXl8FEX6xp/qniMBQlCOEAynCIjiQXBVUAFxdRFZdWXlhxeKuiAKq6yK6CIi66LiwXqAFwqsirpeqHhhBF0UlEMFFkUBOQIJkSskQJKZ6fr9UV093T3dMz2TyQXv188409XV1dXVPT3kqbefl3HOOQiCIAiCIAiCIAiCIAiCIAiCcESp6w4QBEEQBEEQBEEQBEEQBEEQRH2GhHSCIAiCIAiCIAiCIAiCIAiCiAMJ6QRBEARBEARBEARBEARBEAQRBxLSCYIgCIIgCIIgCIIgCIIgCCIOJKQTBEEQBEEQBEEQBEEQBEEQRBxISCcIgiAIgiAIgiAIgiAIgiCIOJCQThAEQRAEQRAEQRAEQRAEQRBxICGdIAiCIAiCIAiCIAiCIAiCIOJAQjpBEARBEARBEARBEARBEARBxIGEdIIgCIIgCIIgCIIgCIKo58yePRuMMTDGsHjx4pj1nHN07twZjDH069ev1vtXU8jj3rx5s6f6I0aMwB/+8Iea7ZRHfv75ZwQCAaxataquu0KkARLSCYIgCIIgCIIgCIIgCKKBkJWVhVmzZsWUf/HFF9i4cSOysrLqoFf1g++++w5z5szBP/7xj7ruCgCgS5cuuPLKK3HbbbfVdVeINEBCOkEQBEEQBEEQBEEQBEE0EIYOHYq33noL+/fvt5TPmjULZ555Jtq1a1dHPat7HnzwQfzud79Dr1696rorBrfccgu+/PJLfP3113XdFaKakJBOEARBEARBEARBEARBEA2EYcOGAQDmzZtnlJWWluKtt97CiBEjHLeZPHkyTj/9dBx99NFo2rQpevbsiVmzZoFzbqnXoUMHXHTRRXjnnXdw0kknISMjA506dcITTzxhqbd48WIwxvDyyy9j3LhxaN26NTIzM9G3b1989913MftfsWIF/vjHP+Loo49GRkYGTj31VLzxxhsx9ZYtW4Y+ffogIyMDbdq0wYQJExAKhTyNy86dO/HOO+/g6quvtpS7WcPIY3CyyYlHKBRCq1atYvYDAPv27UNmZibGjRtnlOXn5+P444/HM888k9R+iPoHCekEQRAEQRAEQRAEQRAE0UBo2rQphgwZghdffNEomzdvHhRFwdChQx232bx5M0aOHIk33ngDb7/9Nv70pz9hzJgxmDJlSkzd77//Hrfeeituu+02vPPOO+jduzf++te/4pFHHompe/fdd2PTpk144YUX8MILL2DHjh3o168fNm3aZNRZtGgR+vTpg3379uGZZ57B/Pnzccopp2Do0KGYPXu2UW/dunUYMGAA9u3bh9mzZ+OZZ57Bd99959mm5dNPP0UoFEL//v091U8Vv9+Pq666yvGpgHnz5qGiogLXXXedpbxfv3746KOPYiYuiIaFr647QBAEQRAEQRAEQRAEQRCEd0aMGIH+/fvjf//7H0444QS8+OKL+POf/+zqj/7SSy8ZnzVNQ79+/cA5x7/+9S9MnDgRjDFj/Y4dO/Ddd9/h5JNPBgAMHDgQJSUlmDJlCkaPHo1GjRoZdVu2bIl33nnH2P6ss87Ccccdh6lTp+L5558HAIwePRonnHACPv/8c/h8Qoq84IILsGvXLtx999245pproCgK7r//fnDO8fnnnyMnJwcAMGjQIJx44omexmTp0qXIzMxEt27dvA5jylx33XV4/PHH8frrr+PGG280ymfPno38/Hz06NHDUr9nz56YOXMm1q9fXyv9I2oGikgnCIIgCIIgCIIgCIIgiAZE3759ceyxx+LFF1/EmjVrsHz5cldbFwD4/PPPcd555yE7OxuqqsLv9+Pee+/F7t27UVJSYql7wgknGCK65IorrsD+/fuxatWqmHKzCN++fXv07t0bixYtAgBs2LABP/30E6688koAQDgcNl4XXnghioqKsH79egAicn3AgAGGiA4Aqqq6Rtnb2bFjB1q2bGnpT03Ro0cP5OfnWyYofvzxR3z77beO56FVq1YAgO3bt9d434iag4R0giAIgiAIgiAIgiAIgmhAMMZw3XXX4eWXX8YzzzyDLl264Oyzz3as++233+L8888HADz//PP46quvsHz5ctxzzz0AgEOHDlnqt27dOqYNWbZ7925PdWW9nTt3AgBuv/12+P1+y2v06NEAgF27dhltx9t3Ig4dOoSMjAxPddPBiBEjsHTpUvz0008ARNR/MBg0POzNyH7Zx5poWJCQThAEQRAEQRAEQRAEQRANjGuvvRa7du3CM888E+PJbea1116D3+/HBx98gMsvvxy9e/dGr169XOsXFxe7ljVv3txTXVmvRYsWAIAJEyZg+fLljq9TTjnFaDvevhPRokUL7NmzJ6ZcitiVlZWWcingp8qwYcMQDAYxe/ZsRCIR/Pvf/8Yll1yCo446Kqau7JccD6JhQkI6QRAEQRAEQRAEQRAEQTQwjjnmGNxxxx0YPHgwhg8f7lqPMQafzwdVVY2yQ4cO4d///rdj/f/973/44YcfLGWvvvoqsrKy0LNnT0v5vHnzLAk0t2zZgq+//hr9+vUDAHTt2hXHHXccfvjhB/Tq1cvxJX3d+/fvj4KCAiOKHQAikQhef/11T+PRrVs37N69G6WlpZbyDh06AABWr15tKX/vvfc8tevGUUcdhUsuuQRz587FBx98gOLiYld7nU2bNkFRFHTt2rVa+yTqFhLSCYIgCIIgCIIgCIIgCKIB8uCDD+Ldd99Fbm6ua51BgwahvLwcV1xxBRYuXIjXXnsNZ599NoLBoGP9Nm3a4I9//CNeeuklfPzxx7jqqquwcOFC/P3vf7ckGgWAkpISXHrppViwYAFeffVVnHfeecjIyMCECROMOs8++ywKCgpwwQUXYN68efjyyy/x7rvvYurUqfjzn/9s1Pv73/8OADj33HPx+uuv4/3338egQYNw4MABT2MhE6h+8803lvLTTjsNXbt2xe2334558+bh448/xsiRI7FkyZKYNhYvXgzGGO677z5P+xwxYgSKiopwyy23IC8vD+edd55jvWXLluGUU05xjFYnGg4kpBMEQRAEQRAEQRAEQRDEYcq5555rJCUdPHgw7rnnHgwZMgR33XWXY/1TTjkFjz32GB599FFcfPHF+Oqrr/DYY4/hzjvvjKn7z3/+E+3bt8d1112HESNGIDc3F4sWLcKxxx5r1Onfvz++/fZbNGvWDLfeeivOO+883HTTTfjss88swvOJJ56Izz77DE2bNsXw4cPxl7/8BSeddBImTpzo6Tj79OmDDh06YP78+ZZyVVXx/vvvo1u3bhg1ahSuueYaBINBPPXUUzFtlJeXA0DciQkz5513Htq2bYvCwkIMHz4cihIrtZaXl6OgoMBIuEo0XBg3P39BEARBEARBEARBEARBEMQRSYcOHXDiiSfigw8+iFtv8eLF6N+/P/7zn/9gyJAhtdS7xDz66KN44IEHsH37dmRmZia9/Z133ol58+bhl19+SVvi0lmzZuGvf/0rtm3bRhHpDRyKSCcIgiAIgiAIgiAIgiAIosFz8803Izs7G08//XRK2y9atAgTJ05Mm4geDofx0EMPYcKECSSiHwb46roDBEEQBEEQBEEQBEEQBEEQ1SUjIwP//ve/8d1336W0/fLly9Pan23btuGqq67C3/72t7S2S9QNFJFOEARBEIcx1157rZGlXtKhQwdce+21tdqPzZs3gzGG2bNne6q/adMm3HLLLejSpQsyMzPRqFEjnHDCCfj73/+O7du3G/WuvfZaMMaMV+PGjdGhQwcjOVJlZWVM2/369bNsY36tXbvWU//Wr1+PRo0a4YorrohZt3fvXhxzzDE4/fTTEYlEErY1d+5ctGzZEmVlZZ72Xdu4jVWLFi08tzF37lwwxvDcc8/FrPv666+hqipuv/12S3koFEK3bt3w4IMPJt3niRMnomfPntA0zVK+d+9eNGvWDO+++27SbRJ1T0O9n9UXnMbPnkxt3bp1uO+++7B582ZP29clySSC279/Px544AH06tULTZs2RTAYRIcOHTBixAisWrXKqDd79mzLfS4jIwOtW7dG//79MXXqVJSUlMS0fd9997neJ528d+3s27cPeXl5rr8ZS5YsgaqqlqR9bhzuvyfTpk0DYwzvvfee4/oLLrgARx99NHbs2AFAnBv7NTtixAgEg0GsWbMmZvsHH3wQjDG8//77lvKNGzciGAxi6dKlnvppRiYtdPpOJaKgoABNmjSx/LtHcs455+DWW29Nuk2C8MLmzZsT2roA0cSe9cnWRXLWWWdhzJgxdd0NAEDHjh1x7733pi3CnahjOEHUAkuXLuVDhgzhrVu35n6/n+fk5PDLLruMf/3110m1M2nSJJ7qZbto0SIOgC9atCil7b3St29f3rdvX091q6qqeNeuXfnUqVNrtE9e2LNnD8/OzubvvPNOXXeFIIg0Mnz4cN6+fXtL2apVq/iGDRtqtR+//vorB8BfeumlhHXff/993rhxY96+fXs+bdo0/tlnn/GCggI+ffp0ftJJJ/FTTjnFqDt8+HCemZnJly5dypcuXco///xzPmfOHP5///d/XFVVfsIJJ/Bt27ZZ2u/bty/v1KmTsY35dfDgQc/H9Pjjj3MA/M0337SU/9///R/PyMjgP/74Y8I2Dhw4wI855hg+bdo0z/utbQDwIUOGxIzVihUrkmrn4osv5k2aNOG//vqrUVZeXs47d+7Mu3Xrxg8dOmSpP336dN6qVSteXl6edJ/37dvHmzVrxl988cWYdffddx/v3Lkzr6ysTLpdom5piPez+oTT+C1dutRyj/zPf/7j+u/lDRs28FWrVtVwL70DgE+Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1500x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## check that interpolator transforms to uniform distribution \n",
|
||||
"interpolator = Interpolator(H, xedges[:-1], yedges[:-1], ebins[:-1]) #dummy energy bins\n",
|
||||
"CDF_X = interpolator.get_ietax() #lookup table x \n",
|
||||
"CDF_Y = interpolator.get_ietay() #lookup table y \n",
|
||||
"\n",
|
||||
"x_indices = ((x - range[0,0])/bin_width[0]).astype(int)\n",
|
||||
"y_indices = ((y - range[1,0])/bin_width[1]).astype(int)\n",
|
||||
"\n",
|
||||
"u = CDF_X[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"v = CDF_Y[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"\n",
|
||||
"plot(x,y,CDF_X,CDF_Y, u,v)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "85823394",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"shape X_tab: (301, 300, 2)\n",
|
||||
"shape Y_tab: (300, 301, 2)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABdIAAAGWCAYAAAB8cHHOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXlg1GT+xp83maOFliJXKZZLERDxAlwRVEBcFZFVV5QfqKCoC+KxLh6ILoKyLt6yHuCFAh54rqLihQi6KCgUFRFFAREKLZVCSwu0nZm8vz/evJkkk8xkplNa4PvZHWcmefNeSVP6vN88X8Y55yAIgiAIgiAIgiAIgiAIgiAIwhGlvjtAEARBEARBEARBEARBEARBEA0ZEtIJgiAIgiAIgiAIgiAIgiAIIg4kpBMEQRAEQRAEQRAEQRAEQRBEHEhIJwiCIAiCIAiCIAiCIAiCIIg4kJBOEARBEARBEARBEARBEARBEHEgIZ0gCIIgCIIgCIIgCIIgCIIg4kBCOkEQBEEQBEEQBEEQBEEQBEHEgYR0giAIgiAIgiAIgiAIgiAIgogDCekEQRAEQRAEQRAEQRAEQRAEEQcS0gmCIAiCIAiCIAiCIAiigTN79mwwxsAYw5IlS2L2c87RqVMnMMbQv3///d6/ukKOe9OmTZ7Kjx49Guecc07ddsojv/zyCwKBAFatWlXfXSHSAAnpBEEQBEEQBEEQBEEQBHGAkJ2djVmzZsVs//zzz7FhwwZkZ2fXQ68aBt9++y3mzJmDf/3rX/XdFQBA586dcemll+If//hHfXeFSAMkpBMEQRAEQRAEQRAEQRDEAcKwYcPw1ltvYffu3Zbts2bNwimnnIJ27drVU8/qn/vuuw9/+tOf0KtXr/ruisH111+PL774Al999VV9d4WoJSSkEwRBEARBEARBEARBEMQBwvDhwwEA8+bNM7aVl5fjrbfewujRox2Pufvuu3HyySejWbNmaNKkCXr06IFZs2aBc24p16FDB5x33nl4++23cdxxxyEjIwNHHHEEHnvsMUu5JUuWgDGGl156CePHj0fr1q2RmZmJfv364dtvv41pf+XKlfjLX/6CZs2aISMjAyeeeCJef/31mHLLly9H3759kZGRgTZt2mDixIkIhUKe5mX79u14++23cfnll1u2u1nDyDE42eTEIxQKoVWrVjHtAEBZWRkyMzMxfvx4Y1vPnj1x9NFH46mnnkqqHaLhQUI6QRAEQRAEQRAEQRAEQRwgNGnSBEOHDsXzzz9vbJs3bx4URcGwYcMcj9m0aRPGjBmD119/Hf/973/x17/+FTfccAOmTp0aU/a7777DTTfdhH/84x94++230adPH/z973/HQw89FFP2jjvuwMaNG/Hcc8/hueeew7Zt29C/f39s3LjRKLN48WL07dsXZWVleOqppzB//nyccMIJGDZsGGbPnm2UW7t2LQYOHIiysjLMnj0bTz31FL799lvPNi2ffPIJQqEQBgwY4Kl8qvj9flx22WWOTwXMmzcPVVVVuPLKKy3b+/fvjw8//DBm4YI4sPDVdwcIgiAIgiAIgiAIgiAIgvDO6NGjMWDAAPz444845phj8Pzzz+Piiy929Ud/4YUXjM+apqF///7gnOM///kPJk2aBMaYsX/btm349ttvcfzxxwMABg0ahJKSEkydOhXjxo1Do0aNjLItW7bE22+/bRx/6qmn4qijjsK0adPw7LPPAgDGjRuHY445Bp999hl8PiFFnn322dixYwfuuOMOjBw5Eoqi4J577gHnHJ999hlyc3MBAIMHD0b37t09zcmyZcuQmZmJrl27ep3GlLnyyivx6KOP4rXXXsM111xjbJ89ezZ69uyJY4891lK+R48emDlzJtatW7df+kfUDRSRThAEQRAEQRAEQRAEQRAHEP369cORRx6J559/Hj/88ANWrFjhausCAJ999hnOPPNM5OTkQFVV+P1+3HXXXSgtLUVJSYml7DHHHGOI6JIRI0Zg9+7dWLVqVcx2swjfvn179OnTB4sXLwYArF+/Hj///DMuvfRSAEA4HDZe5557LoqKirBu3ToAInJ94MCBhogOAKqqukbZ29m2bRtatmxp6U9dceyxx6Jnz56WBYqffvoJ33zzjeN5aNWqFQBg69atdd43ou4gIZ0gCIIgCIIgCIIgCIIgDiAYY7jyyivx0ksv4amnnkLnzp1x2mmnOZb95ptvcNZZZwEAnn32WXz55ZdYsWIF7rzzTgDAvn37LOVbt24dU4fcVlpa6qmsLLd9+3YAwC233AK/3295jRs3DgCwY8cOo+54bSdi3759yMjI8FQ2HYwePRrLli3Dzz//DEBE/QeDQcPD3ozsl32uiQMLEtIJgiAIgiAIgiAIgiAI4gDjiiuuwI4dO/DUU0/FeHKbefXVV+H3+/H+++/jkksuQZ8+fdCrVy/X8sXFxa7bmjdv7qmsLNeiRQsAwMSJE7FixQrH1wknnGDUHa/tRLRo0QI7d+6M2S5F7Orqast2KeCnyvDhwxEMBjF79mxEIhG8+OKLuOCCC3DYYYfFlJX9kvNBHJiQkE4QBEEQBEEQBEEQBEEQBxiHH344br31VgwZMgSjRo1yLccYg8/ng6qqxrZ9+/bhxRdfdCz/448/4vvvv7dse+WVV5CdnY0ePXpYts+bN8+SQPP333/HV199hf79+wMAunTpgqOOOgrff/89evXq5fiSvu4DBgzAokWLjCh2AIhEInjttdc8zUfXrl1RWlqK8vJyy/YOHToAAFavXm3Z/u6773qq143DDjsMF1xwAebOnYv3338fxcXFrvY6GzduhKIo6NKlS63aJOoXEtIJgiAIgiAIgiAIgiAI4gDkvvvuwzvvvIO8vDzXMoMHD0ZlZSVGjBiBhQsX4tVXX8Vpp52GYDDoWL5Nmzb4y1/+ghdeeAEfffQRLrvsMixcuBD//Oc/LYlGAaCkpAQXXnghFixYgFdeeQVnnnkmMjIyMHHiRKPM008/jUWLFuHss8/GvHnz8MUXX+Cdd97BtGnTcPHFFxvl/vnPfwIAzjjjDLz22mt47733MHjwYOzZs8fTXMgEql9//bVl+0knnYQuXbrglltuwbx58/DRRx9hzJgxWLp0aUwdS5YsAWMMU6ZM8dTm6NGjUVRUhOuvvx75+fk488wzHcstX74cJ5xwgmO0OnHgQEI6QRAEQRAEQRAEQRAEQRyknHHGGUZS0iFDhuDOO+/E0KFDcfvttzuWP+GEE/DII4/g4Ycfxvnnn48vv/wSjzzyCG677baYsv/+97/Rvn17XHnllRg9ejTy8vKwePFiHHnkkUaZAQMG4JtvvkHTpk1x00034cwzz8S1116LTz/91CI8d+/eHZ9++imaNGmCUaNG4W9/+xuOO+44TJo0ydM4+/btiw4dOmD+/PmW7aqq4r333kPXrl0xduxYjBw5EsFgEE888URMHZWVlQAQd2HCzJlnnom2bduisLAQo0aNgqLESq2VlZVYtGiRkXCVOHBh3Pz8BUEQBEEQBEEQBEEQBEEQhyQdOnRA9+7d8f7778ctt2TJEgwYMABvvPEGhg4dup96l5iHH34Y9957L7Zu3YrMzMykj7/tttswb948/Prrr2lLXDpr1iz8/e9/x5YtWygi/QCHItIJgiAIgiAIgiAIgiAIgjjgue6665CTk4Mnn3wypeMXL16MSZMmpU1ED4fDuP/++zFx4kQS0Q8CfPXdAYIgCIIgCIIgCIIgCIIgiNqSkZGBF198Ed9++21Kx69YsSKt/dmyZQsuu+wy3HzzzWmtl6gfKCKdIAiCIA5irrjiCiNLvaRDhw644oor9ms/Nm3aBMYYZs+e7an8xo0bcf3116Nz587IzMxEo0aNcMwxx+Cf//wntm7dapS74oorwBgzXo0bN0aHDh2M5EjV1dUxdffv399yjPm1Zs0aT/1bt24dGjVqhBEjRsTs27VrFw4//HCcfPLJiEQiCeuaO3cuWrZsiYqKCk9t72/c5qpFixae65g7dy4YY3jmmWdi9n311VdQVRW33HKLZXsoFELXrl1x3333Jd3nSZMmoUePHtA0zbJ9165daNq0Kd55552k6yTqnwP1ftZQcJo/ezK1tWvXYsqUKdi0aZOn4+uTZBLB7d69G/feey969eqFJk2aIBgMokOHDhg9ejRWrVpllJs9e7blPpeRkYHWrVtjwIABmDZtGkpKSmLqnjJliut90sl7105ZWRny8/Ndf2csXboUqqpakva5cbD/PnnwwQfBGMO7777ruP/ss89Gs2bNsG3bNgDi3Niv2dGjRyMYDOKHH36IOf6+++4DYwzvvfeeZfuGDRsQDAaxbNkyT/00I5MWOv1MJWLRokXIysqy/LtHcvrpp+Omm25Kuk6C8MKmTZsS2roA0cSeDcnWRXLqqafihhtuqO9uAAA6duyIu+66K20R7kQ9wwliP7Bs2TI+dOhQ3rp1a+73+3lubi6/6KKL+FdffZVUPZMnT+apXraLFy/mAPjixYtTOt4r/fr14/369fNUtqamhnfp0oVPmzatTvvkhZ07d/KcnBz+9ttv13dXCIJII6NGjeLt27e3bFu1ahVfv379fu3Hb7/9xgHwF154IWHZ9957jzdu3Ji3b9+eP/jgg/zTTz/lixYt4tOnT+fHHXccP+GEE4yyo0aN4pmZmXzZsmV82bJl/LPPPuNz5szh//d//8dVVeXHHHMM37Jli6X+fv368SOOOMI4xvzau3ev5zE9+uijHAB/8803Ldv/7//+j2dkZPCffvopYR179uzhhx9+OH/wwQc9t7u/AcCHDh0aM1crV65Mqp7zzz+fZ2Vl8d9++83YVllZyTt16sS7du3K9+3bZyk/ffp03qpVK15ZWZl0n8vKynjTpk35888/H7NvypQpvFOnTry6ujrpeon65UC8nzUknOZv2bJllnvkG2+84frv5fXr1/NVq1bVcS+9A4BPnjwLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1500x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"x_indices = ((x - range[0,0])/bin_width[0]).astype(int)\n",
|
||||
"y_indices = ((y - range[1,0])/bin_width[1]).astype(int)\n",
|
||||
"\n",
|
||||
"CDF_X, CDF_Y = calculate_rosenblatt_lookup_table(H, xedges[:-1], yedges[:-1],ebins[:-1])\n",
|
||||
"\n",
|
||||
"u = CDF_X[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"v = CDF_Y[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
" \n",
|
||||
"plot(x,y,CDF_X,CDF_Y, u,v)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "fe9e3306",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"worked: \n",
|
||||
"could normalize CDF eta_x\n",
|
||||
"construction worked\n",
|
||||
"calculation of conditional CDF etay worked\n",
|
||||
"assignment m_ietay worked\n",
|
||||
"assignment m_ietax worked\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABcUAAAGPCAYAAACd27gMAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXl8FFX69U91J+mEJARISAiE1QgomwIOKDKAKJgBRAeRF0YBUUZkxFFEkXHcYBQFBnEDXFhkVFxwwQUXQNEfAgJBQUDQgCydlQTSJJDO0n3fP6pv9a3qqurqTkICeb5+MN1Vt+69dbtTSU6dPo/EGGMgCIIgCIIgCIIgCIIgCIIgiAaAra4nQBAEQRAEQRAEQRAEQRAEQRDnChLFCYIgCIIgCIIgCIIgCIIgiAYDieIEQRAEQRAEQRAEQRAEQRBEg4FEcYIgCIIgCIIgCIIgCIIgCKLBQKI4QRAEQRAEQRAEQRAEQRAE0WAgUZwgCIIgCIIgCIIgCIIgCIJoMJAoThAEQRAEQRAEQRAEQRAEQTQYSBQnCIIgCIIgCIIgCIIgCIIgGgwkihMEQRAEQRAEQRAEQRAEQRANBhLFCYIgCIIgCIIgCIIgCKKes3LlSkiSBEmSsGnTpoD9jDGkp6dDkiQMHDjwnM+vtuDnfeTIEUvtJ02ahOuvv752J2WR3377DVFRUdi1a1ddT4XQQKI4QRAEQRAEQRAEQRAEQZwnxMfHY9myZQHbv/vuOxw6dAjx8fF1MKv6wU8//YQ33ngD//nPf+p6KgCAjh074m9/+xvuv//+up4KoYFEcYIgCIIgCIIgCIIgCII4TxgzZgw++OADnD59WrV92bJluPLKK9GmTZs6mlnd88wzz+BPf/oTevfuXddTUbjnnnvw/fffY8uWLXU9FUKARHGCIAiCIAiCIAiCIAiCOE8YO3YsAGD16tXKNpfLhQ8++ACTJk3SPebJJ59Enz590KxZMzRu3Bg9e/bEsmXLwBhTtWvXrh2GDx+Ojz76CN27d0d0dDQ6dOiAF154QdVu06ZNkCQJb775JqZPn44WLVogJiYGAwYMwE8//RQw/s6dO3HDDTegWbNmiI6OxuWXX4733nsvoN22bdvQr18/REdHo2XLlpg1axYqKystrUt+fj4++ugj3HbbbartRvEr/Bz0omjMqKysRHJycsA4AFBcXIyYmBhMnz5d2darVy9ccsklWLp0aUjjELULieIEQRAEQRAEQRAEQRAEcZ7QuHFj3HzzzVi+fLmybfXq1bDZbBgzZozuMUeOHMFdd92F9957Dx9++CH++te/Ytq0aZgzZ05A259//hn33Xcf7r//fnz00Ue46qqr8M9//hMLFiwIaPuvf/0Lhw8fxuuvv47XX38dOTk5GDhwIA4fPqy0+fbbb9GvXz8UFxdj6dKlWLt2LS677DKMGTMGK1euVNrt378fgwcPRnFxMVauXImlS5fip59+shyF8vXXX6OyshKDBg2y1D5cIiMjceutt+q69VevXg23243bb79dtX3gwIH44osvAm5CEHVHRF1PgCAIgiAIgiAIgiAIgiAI60yaNAmDBg3Cvn370KVLFyxfvhyjR482zBNfsWKF8tjr9WLgwIFgjOH555/Ho48+CkmSlP05OTn46aef0KNHDwBARkYGCgoKMGfOHEydOhWNGjVS2jZv3hwfffSRcvzVV1+Niy++GHPnzsVrr70GAJg6dSq6dOmCb775BhERshQ5dOhQFBYW4l//+hfGjx8Pm82G2bNngzGGb775BikpKQCAYcOGoWvXrpbWZOvWrYiJiUHnzp2tLmPY3H777Xjuuefw7rvvYvLkycr2lStXolevXujWrZuqfc+ePbFkyRIcPHjwnMyPCA45xQmCIAiCIAiCIAiCIAjiPGLAgAG46KKLsHz5cvzyyy/YsWOHYXQKAHzzzTe49tprkZCQALvdjsjISDz22GMoKipCQUGBqm2XLl0UQZwzbtw4nD59Grt27QrYLgrqbdu2xVVXXYVvv/0WAJCVlYUDBw7gb3/7GwCgqqpK+feXv/wFubm5OHjwIADZUT548GBFEAcAu91u6H7XkpOTg+bNm6vmU1t069YNvXr1Ut1s+PXXX7F9+3bd1yE5ORkAkJ2dXetzI6xBojhBEARBEARBEARBEARBnEdIkoTbb78db775JpYuXYqOHTuif//+um23b9+OIUOGAABee+01/PDDD9ixYwceeeQRAEBZWZmqfYsWLQL64NuKioosteXt8vPzAQAzZsxAZGSk6t/UqVMBAIWFhUrfZmMHo6ysDNHR0Zba1gSTJk3C1q1bceDAAQCyG9/hcCiZ7yJ8Xtq1JuoOEsUJgiAIgiAIgiAIgiAI4jxj4sSJKCwsxNKlSwMyrEXeeecdREZG4rPPPsMtt9yCq666Cr179zZsn5eXZ7gtMTHRUlveLikpCQAwa9Ys7NixQ/ffZZddpvRtNnYwkpKScPLkyYDtXJAuLy9XbedifLiMHTsWDocDK1euhMfjwf/+9z/ceOONaNq0aUBbPi++HkTdQ6I4QRAEQRAEQRAEQRAEQZxntGrVCg8++CBGjBiBCRMmGLaTJAkRERGw2+3KtrKyMvzvf//Tbb9v3z7s3r1bte3tt99GfHw8evbsqdq+evVqVfHIo0ePYsuWLRg4cCAAoFOnTrj44ouxe/du9O7dW/cfz0EfNGgQNm7cqLjLAcDj8eDdd9+1tB6dO3dGUVERXC6Xanu7du0AAHv27FFt/+STTyz1a0TTpk1x4403YtWqVfjss8+Ql5dnGGFz+PBh2Gw2dOrUqVpjEjUHieIEQRAEQRAEQRAEQRAEcR7yzDPP4OOPP0Zqaqphm2HDhqG0tBTjxo3D+vXr8c4776B///5wOBy67Vu2bIkbbrgBK1aswJdffolbb70V69evx7///W9VkU0AKCgowE033YTPP/8cb7/9Nq699lpER0dj1qxZSptXXnkFGzduxNChQ7F69Wp8//33+PjjjzF37lyMHj1aaffvf/8bAHDNNdfg3Xffxaeffophw4bhzJkzltaCFw/98ccfVduvuOIKdOrUCTNmzMDq1avx5Zdf4q677sLmzZsD+ti0aRMkScITTzxhacxJkyYhNzcX99xzD9LS0nDttdfqttu2bRsuu+wyXRc5UTeQKE4QBEEQFzATJ05UnBGcdu3aYeLEied0HkeOHIEkSVi5cqWl9ocPH8Y999yDjh07IiYmBo0aNUKXLl3w73//W1WcZuLEiZAkSfkXGxuLdu3aKb/Eaz8iCci/LIvHiP/27t1raX4HDx5Eo0aNMG7cuIB9p06dQqtWrdCnTx94PJ6gfa1atQrNmzdHSUmJpbHPNUZrFcpHP1etWgVJkvDqq68G7NuyZQvsdjtmzJih2l5ZWYnOnTvjmWeeCXnOjz76KHr27Amv16vafurUKTRp0gQff/xxyH0Sdc/5ej2rL+itn/aP/v379+OJJ57AkSNHLB1fl4QiWJw+fRpPPfUUevfujcaNG8PhcKBdu3aYNGmSqmDcypUrVde56OhotGjRAoMGDcLcuXMDCtEBwBNPPGF4nXzppZeCzq24uBhpaWmGPzM2b94Mu92uEpeMuNB/nsyfPx+SJBk6O4cOHYpmzZohJycHgPzaaN+zkyZNgsPhwC+//BJw/DPPPANJkvDpp5+qth86dAgOhwNbt261NE8RLq7pfU8FY+PGjYiLi9MtyvfnP/8Z9913X8h9Eg2Xa665RinIOWLECDzyyCO4+eab8fDDD+u2v+yyy7Bw4UL897//xciRI/HDDz9g4cKFeOihhwLaPv3002jbti1uv/12TJo0Campqfj2229x0UUXKW0GDRqE7du3o0mTJrjvvvtw7bXX4u6778aGDRtUInLXrl2xYcMGNG7cGBMmTMDf//53dO/eHY8++qil8+zXrx/atWuHtWvXqrbb7XZ8+umn6Ny5M6ZMmYLx48fD4XDoXqdLS0sBwPQmg8i1116L1q1bw+l0YsKECbDZAqXW0tJSbNy4USk2StQTGEEQBEEQFywTJkxgbdu2VW3btWsXy8rKOqfz+OOPPxgAtmLFiqBtP/30UxYbG8vatm3L5s+fzzZs2MA2btzIFi1axLp3784uu+wype2ECRNYTEwM27p1K9u6dSv75ptv2BtvvMH+3//7f8xut7MuXbqw48ePq/ofMGAA69Chg3KM+O/s2bOWz+m5555jANiaNWtU2//f//t/LDo6mv36669B+zhz5gxr1aoVmz9/vuVxzzUA2M033xywVjt37gypn5EjR7K4uDj2xx9/KNtKS0tZeno669y5MysrK1O1X7RoEUtOTmalpaUhz7m4uJg1adKELV++PGDfE088wdLT01l5eXnI/RJ1y/l4PatP6K3f1q1bVdfI999/nwFg3377bcDxWVlZbNeuXbU8S+sAYI8//njQdllZWaxDhw4sLi6OzZgxg3322Wds06ZNbOXKlewvf/kLA8CKi4sZY4ytWLFCeW23bt3Kvv/+e7ZmzRp23333sYSEBNasWTO2fv16Vf+PP/44A8C+/PLLgOtkXl6epXP56quvGAD21FNPqbafOXOGpaens27dugW9ZjWEnycej4ddffXVrEWLFqyoqEi175VXXmEA2OrVq5Vtjz/+eMB73uVysTZt2rDLL7+cVVRUKNv37NnDoqKi2MSJEwPGvfHGG9mwYcNCOFM/3377LQOg+tkXCoMGDWLjx48P2L5p0yYWGRnJDhw4EFa/BGFG27ZtLb3n+fv7/fffPwezss6CBQtY06ZNQ/q9XuTBBx9kaWlpAb+bVofXX3+dxcbGspMnT9ZYn0T1IVGcOCds3bqV3XzzzaxFixYsMjKSpaSksFGjRrEtW7aE1A//pTMc+AVb75f8mmTAgAFswIABltpWVFSwTp06sblz59bqnKxw8uRJlpCQwD766KO6ngpBXNCE+8tZuOiJIHWLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1500x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## check that cpp rosenblatt transforms to uniform distribution \n",
|
||||
"interpolator = Interpolator(xedges[:-1], yedges[:-1], ebins[:-1]) #dummy energy bins\n",
|
||||
"interpolator.rosenblatttransform(H)\n",
|
||||
"CDF_X = interpolator.get_ietax() #lookup table x \n",
|
||||
"CDF_Y = interpolator.get_ietay() #lookup table y \n",
|
||||
"\n",
|
||||
"x_indices = ((x - range[0,0])/bin_width[0]).astype(int)\n",
|
||||
"y_indices = ((y - range[1,0])/bin_width[1]).astype(int)\n",
|
||||
"\n",
|
||||
"u = CDF_X[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"v = CDF_Y[x_indices, y_indices, 0] #transformed to uniform coordinates\n",
|
||||
"\n",
|
||||
"plot(x,y,CDF_X,CDF_Y, u,v)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6df08eb9",
|
||||
"metadata": {},
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a9d088d2",
|
||||
"metadata": {},
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "40a1ffff",
|
||||
"metadata": {},
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "1e0aa9dc",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "base",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.11.11"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
File diff suppressed because it is too large.
Load diff
@@ -1,174 +1,175 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "efef8e20-6571-4561-8f0b-6048b57907de",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0, '/home/ferjao_k/sw/aare/build')\n",
|
||||
"from aare import Pol1, Pol2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8449f9e9-59ee-4194-b477-d014b8efc93b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Ground truth\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "35772e2c-37c6-4986-a9c0-7dca4a034827",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"x = np.linspace(0, 10, 50)\n",
|
||||
" \n",
|
||||
"# Pol1: y = 3.0 + 1.5*x\n",
|
||||
"true_pol1 = [3.0, 1.5]\n",
|
||||
"y_pol1 = true_pol1[0] + true_pol1[1] * x\n",
|
||||
"y_pol1_noisy = y_pol1 + np.random.default_rng(42).normal(0, 0.5, x.size)\n",
|
||||
" \n",
|
||||
"# Pol2: y = 2.0 - 0.5*x + 0.3*x^2\n",
|
||||
"true_pol2 = [2.0, -0.5, 0.3]\n",
|
||||
"y_pol2 = true_pol2[0] + true_pol2[1] * x + true_pol2[2] * x**2\n",
|
||||
"y_pol2_noisy = y_pol2 + np.random.default_rng(7).normal(0, 1.0, x.size)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a318ab93-fe18-4798-9ddb-6ffc2c2e3875",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Fit with error estimation"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "06c8fddb-56d9-4f84-8b21-4ffd8b7bd26f",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"=== Pol1 ===\n",
|
||||
" True: p0=3.0000 p1=1.5000\n",
|
||||
" Fitted: p0=2.9325 p1=1.5226\n",
|
||||
" Chi2: [7.00870822]\n",
|
||||
"\n",
|
||||
"=== Pol2 ===\n",
|
||||
" True: p0=2.0000 p1=-0.5000 p2=0.3000\n",
|
||||
" Fitted: p0=2.1391 p1=-0.8442 p2=0.3383\n",
|
||||
" Chi2: [34.11224393]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"m1 = Pol1()\n",
|
||||
"m2 = Pol2()\n",
|
||||
" \n",
|
||||
"res1 = m1.fit(x, y_pol1_noisy)\n",
|
||||
"res2 = m2.fit(x, y_pol2_noisy)\n",
|
||||
"\n",
|
||||
"p1 = res1['par']\n",
|
||||
"print(\"=== Pol1 ===\")\n",
|
||||
"print(f\" True: p0={true_pol1[0]:.4f} p1={true_pol1[1]:.4f}\")\n",
|
||||
"print(f\" Fitted: p0={p1[0]:.4f} p1={p1[1]:.4f}\")\n",
|
||||
"print(f\" Chi2: {res1['chi2']}\")\n",
|
||||
"print()\n",
|
||||
"\n",
|
||||
"p2 = res2['par']\n",
|
||||
"print(\"=== Pol2 ===\")\n",
|
||||
"print(f\" True: p0={true_pol2[0]:.4f} p1={true_pol2[1]:.4f} p2={true_pol2[2]:.4f}\")\n",
|
||||
"print(f\" Fitted: p0={p2[0]:.4f} p1={p2[1]:.4f} p2={p2[2]:.4f}\")\n",
|
||||
"print(f\" Chi2: {res2['chi2']}\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7e185144-3b72-41f4-80d0-3c3504e48bba",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Plot"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "9a0b3292-7c57-4483-aecc-7ef683ff0b62",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAHqCAYAAADVi/1VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADH/0lEQVR4nOzde3zO9f/H8cfn2uHaeYwdMHNIq2zoQAmFvo3kVDmmRDoocqyUjhQp5ZBIB/0ckkqKjpgUOuAbUbHChGnMNFw72Obars/vj4uLfZ3Zdl3bnvfbza193p/39ble117W3l6f9/v9MUzTNBERERERERERESlFFncHICIiIiIiIiIiFY+KUiIiIiIiIiIiUupUlBIRERERERERkVKnopSIiIiIiIiIiJQ6FaVERERERERERKTUqSglIiIiIiIiIiKlTkUpEREREREREREpdSpKiYiIiIiIiIhIqVNRSkRERERERERESp2KUiLi0WbNmoVhGK4/3t7eREdHc++995Kamnpe12rVqhWtWrUq0jZnzhx69uzJZZddhsVioXbt2sUXvIiIiIgblOT4ae/evTzzzDNcf/31VK1alZCQEK655hreeecdCgsLi/mTiEh55+3uAEREzsXMmTO5/PLLyc3NZdWqVYwbN46VK1fyxx9/EBgYeMHXff/990lLS+Paa6/F4XBgt9uLMWoRERER9ymJ8dP69euZM2cO99xzD88++yw+Pj4sXryYhx9+mDVr1vB///d/xfwpRKQ8U1FKRMqE+Ph4GjduDEDr1q0pLCzkxRdfZNGiRdx1110XfN2lS5disTgnjXbo0IFNmzYVS7wiIiIi7lYS46fmzZuzfft2fHx8XG0JCQkcOXKEadOmMXr0aGrWrFks8YtI+afleyJSJjVt2hSAXbt2kZeXx8iRI6lTpw6+vr7UqFGDgQMHcujQobNe51hBSkRERKS8K47xU+XKlYsUpI659tprAfjnn3+KPW4RKb80U0pEyqTk5GQAwsPDue2221i+fDkjR47khhtu4Pfff+f5559n9erVrF69GqvV6uZoRURERNyvJMdP3333Hd7e3sTGxpZE6CJSTqkoJSJlQmFhIQUFBeTl5bFy5UrGjBlDcHAwISEhLF26lPHjx/P4448DzinkNWvWpEePHsyZM4cHHnjAzdGLiIiIlL7SGj8lJiby/vvvM2TIEKpUqVJSH0dEyiGtWxGRMqFp06b4+PgQHBxMhw4diIqKYvHixfz6668A9O3bt0j/bt26ERgYyPLly90QrYiIiIj7lcb46ddff6V79+40bdqUcePGFWf4IlIBaKaUiJQJc+bM4YorrsDb25vIyEiqVasGOJ8q4+3tTXh4eJH+hmEQFRVFRkaGO8IVERERcbuSHj9t2LCBhIQELr30Ur755httmSAi501FKREpE6644grX02NOVKVKFQoKCti/f3+RgZVpmqSlpdGkSZPSDFNERETEY5Tk+GnDhg3cfPPN1KpVi8TEREJDQ4s1dhGpGLR8T0TKtP/85z8AzJ07t0j7p59+Sk5Ojuu8iIiIiDhd7Php48aN3HzzzURHR7Ns2TIqV65cYrGKSPmmmVIiUqYlJCTQtm1bnnjiCTIzM2nevLnr6TFXXXUVvXv3PuPrk5KSSEpKAiAtLY3Dhw+zYMECAOrXr0/9+vVL/DOIiIiIlKaLGT9t2bKFm2++GYCxY8eybds2tm3b5jp/ySWXnLQsUETkdFSUEpEyzTAMFi1axKhRo5g5cyZjx46latWq9O7dm5deeumsexvMnz+f0aNHF2nr1q0bAM8//zyjRo0qqdBFRERE3OJixk+rV6927TnVsWPHk87PnDnzpA3URUROxzBN03R3ECIiIiIiIiIiUrFoTykRERERERERESl1KkqJiIiIiIiIiEipU1FKRERERERERERKnYpSIiIiIiIiIiJS6lSUEhERERERERGRUqeilIiIiIiIiIiIlDpvdwfwvxwOB3v27CE4OBjDMNwdjoiIiFRgpmmSlZVFcHAwISEhHj020RhKREREPMWxMVT16tWxWE4/H8rjilJ79uyhZs2a7g5DREREpAibzUZISIi7wzgtjaFERETE0+zevZvo6OjTnve4olRwcDDgDLykBn4Oh4P9+/cTHh5+xoqdlA7lw7MoH55F+fAsyodnKY18ZGZmUrNmTXbv3u0ao1yI6dOnM336dHbu3AlAXFwczz33HO3atQOcdxNHjx7NO++8w8GDB7nuuuuYNm0acXFx5/weGkNVPMqHZ1E+PIvy4VmUD89SmmOos42fPK4odWy6eUhISIkOqPLy8ggJCdEPhAdQPjyL8uFZlA/Ponx4ltLMx8Uu3YuOjubll1+mXr16AMyePZvOnTuzYcMG4uLiGD9+PBMnTmTWrFnExsYyZswYEhIS2LJlyzkXwzSGqniUD8+ifHgW5cOzKB+epTTzcbbxk/42iIiIiJSwjh07cuuttxIbG0tsbCxjx44lKCiINWvWYJomkydP5umnn+aOO+4gPj6e2bNnc/jwYebNm+fu0EVERERKjIpSIiIiIqWosLCQjz76iJycHK6//np27NhBWloabdq0cfWxWq20bNmSn3/+2Y2RioiIiJQsj1u+JyIiIlIe/fHHH1x//fXk5eURFBTEwoULqV+/vqvwFBkZWaR/ZGQku3btOu318vPzyc/Pdx1nZmYCzin5DoejBD6B89qmaZbY9eX8KB+eRfnwLMqHZ1E+PEtp5ONcr11mi1KFhYXY7fYLeq3D4cBut5OXl1eh17P6+Pjg5eXl7jBEREQqhMsuu4yNGzdy6NAhPv30U/r06cPKlStd5/93zwXTNM+4D8O4ceMYPXr0Se379+8nLy/vtK9zOBwUFhZewCc4/nhnu91+UXtslWVeXl4eM350OBzYbDZM0/SYmCoy5cOzKB+eRfnwLKWRj6ysrHPqV+aKUqZpkpaWxqFDhy7qGg6Hg6ysrAo7oDqmUqVKREVFVfjvg4iISEnz9fV1bXTeuHFjfvnlF15//XWeeOIJANLS0qhWrZqrf3p6+kmzp040cuRIhg8f7jo+9pSb8PDwU250bpom+/btu6gx1LHrZGdnX9Q1yrpKlSoRGRnp9vGTw+HAMAw9zcpDKB+eRfnwLMqHZymNfPj5+Z1TvzJXlDpWkIqIiCAgIOCCBgOmaVJQUIC3t7fbBxPuYpomhw8fJj09HaDIIFhERERKnmma5OfnU6dOHaKioli2bBlXXXUVAEeOHGHlypW88sorp3291WrFarWe1G6xWE45wNy7dy82m43IyEiNoS7QieMnwzA8YvxkGMZpcy6lT/nwLMqHZ1E+PEtJ5+Ncr1umilKFhYWuglSVKlUu+DoVfUB1jL+/P+C8ExsREaGlfCIiIiXkqaeeol27dtSsWZOsrCw++ugjVqxYwZIlSzAMg6FDh/LSSy9x6aWXcumll/LSSy8REBBAr169iuX9NYYqPho/iYiIFJ8yVZQ6todUQECAmyMpP459L+12uwZVIiIiJWTfvn307t2bvXv3EhoaSsOGDVmyZAkJCQkAjBgxgtzcXAYMGMDBgwe57rrrSExMJDg4uFjeX2Oo4qXxk4iISPEoU0WpYyrqnbmSoO+liIhIyXvvvffOeN4wDEaNGsWoUaNKNA793i8e+j6KiIgUDy3mFBERERERERGRUqeilJu1atWKoUOHujsMERERkTJD4ycREZHyQUWpMmTFihUYhnHRj3IWERERqSg0fhIREfFcKkqJiIiIiIiIiEipq7BFqc17Mpm4bCtDPtrAhMQtbEq1lfh75uTkcM899xAUFES1atWYMGFCkfNz586lcePGBAcHExUVRa9evUhPTwdg586dtG7dGoDKlStjGAZ9+/YFYMmSJbRo0YJKlSpRpUoVOnTowPbt20v884iIiEjFsnmPjcnLkxn68UaNn0REROSiVcii1OY9NqZ8l8zq7f+SmVvA6u0ZTFq2tcQHVo8//jjff/89CxcuJDExkRUrVrB+/XrX+SNHjvDiiy/y22+/sWjRInbs2OEaONWsWZNPP/0UgC1btrB3715ef/11wDlYGz58OL/88gvLly/HYrFw++2343A4SvTziIiIuJtx8KC7Q6gwNqXamLRsG2v+PkBmrl3jJxERkTJoU6qNicu28vziHUwshd/jZ+Pt1nd3k6Wb92HLLeDSyCAMw0JkiJXk9GwSk/YRXyO0RN4zOzub9957jzlz5pCQkADA7NmziY6OdvXp16+f6+u6desyZcoUrr32WrKzswkKCiIsLAyAiIgIKlWq5OrbpUuXIu/13nvvERERQVJSEvHx8SXyeURERNzuyy8J79MH3noLevZ0dzTl3tLNadhyj3BJeCBeXhYiQ9D4SUREpAxx3mDaii33CL6GgzV/Z/Dn3iyGJcSW2O/ysznvmVKrVq2iY8eOVK9eHcMwWLRoUZHz2dnZPPLII0RHR+Pv788VV1zB9OnTiyveYpFy4DBBft4YhgGAYRgE+/mwKyOnxN5z+/btHDlyhOuvv97VFhYWxmWXXeY63rBhA507d6ZWrVoEBwfTqlUrZ7wpKWe9dq9evahbty4hISHUqVPnnF4nIiJSJhUUwJNPYrntNiw2G8YDD8CWLe6OqtxLOXCYYD8fjZ9ERETKKOcNJjv1woOICPKlXngQtlw7iUn73BbTeRelcnJyaNSoEVOnTj3l+WHDhrFkyRLmzp3Ln3/+ybBhwxg0aBCff/75RQdbXGLCAsjOK8A0TQBM0yQrz07tKoEl9p7H3ut0cnJyaNOmDUFBQcydO5dffvmFhQsXAs5p6WfSsWNHMjIyePfdd1m7di1r1649p9eJiIh4mk2pNiYkbjn9no9paXDzzfDKK66mbY2uJ8nhX8qRVjwxYQFk5dk1fhIRESmj3HGD6WzOuyjVrl07xowZwx133HHK86tXr6ZPnz60atWK2rVr8+CDD9KoUSPWrVt30cELine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1200x500 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n",
|
||||
" \n",
|
||||
"# Pol1\n",
|
||||
"ax1.plot(x, y_pol1_noisy, 'o', ms=4, alpha=0.6, label='data')\n",
|
||||
"ax1.plot(x, y_pol1, 'k--', label='truth')\n",
|
||||
"ax1.plot(x, m1(x, res1['par']), 'r-', lw=2, label='fit')\n",
|
||||
"ax1.set_title('Pol1')\n",
|
||||
"ax1.legend()\n",
|
||||
"ax1.grid(alpha=0.3)\n",
|
||||
" \n",
|
||||
"# Pol2\n",
|
||||
"ax2.plot(x, y_pol2_noisy, 'o', ms=4, alpha=0.6, label='data')\n",
|
||||
"ax2.plot(x, y_pol2, 'k--', label='truth')\n",
|
||||
"ax2.plot(x, m2(x, res2['par']), 'r-', lw=2, label='fit')\n",
|
||||
"ax2.set_title('Pol2')\n",
|
||||
"ax2.legend()\n",
|
||||
"ax2.grid(alpha=0.3)\n",
|
||||
" \n",
|
||||
"fig.tight_layout()\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e694c18a-a203-4132-a472-c93c4d3ce853",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.12.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "efef8e20-6571-4561-8f0b-6048b57907de",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0, '/home/ferjao_k/sw/aare/build')\n",
|
||||
"from aare import Pol1, Pol2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8449f9e9-59ee-4194-b477-d014b8efc93b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Ground truth\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "35772e2c-37c6-4986-a9c0-7dca4a034827",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"x = np.linspace(0, 10, 50)\n",
|
||||
" \n",
|
||||
"# Pol1: y = 3.0 + 1.5*x\n",
|
||||
"true_pol1 = [3.0, 1.5]\n",
|
||||
"y_pol1 = true_pol1[0] + true_pol1[1] * x\n",
|
||||
"y_pol1_noisy = y_pol1 + np.random.default_rng(42).normal(0, 0.5, x.size)\n",
|
||||
" \n",
|
||||
"# Pol2: y = 2.0 - 0.5*x + 0.3*x^2\n",
|
||||
"true_pol2 = [2.0, -0.5, 0.3]\n",
|
||||
"y_pol2 = true_pol2[0] + true_pol2[1] * x + true_pol2[2] * x**2\n",
|
||||
"y_pol2_noisy = y_pol2 + np.random.default_rng(7).normal(0, 1.0, x.size)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a318ab93-fe18-4798-9ddb-6ffc2c2e3875",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Fit with error estimation"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "06c8fddb-56d9-4f84-8b21-4ffd8b7bd26f",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"=== Pol1 ===\n",
|
||||
" True: p0=3.0000 p1=1.5000\n",
|
||||
" Fitted: p0=2.9325 p1=1.5226\n",
|
||||
" Chi2: [7.00870822]\n",
|
||||
"\n",
|
||||
"=== Pol2 ===\n",
|
||||
" True: p0=2.0000 p1=-0.5000 p2=0.3000\n",
|
||||
" Fitted: p0=2.1391 p1=-0.8442 p2=0.3383\n",
|
||||
" Chi2: [34.11224393]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"m1 = Pol1()\n",
|
||||
"m2 = Pol2()\n",
|
||||
" \n",
|
||||
"res1 = m1.fit(x, y_pol1_noisy)\n",
|
||||
"res2 = m2.fit(x, y_pol2_noisy)\n",
|
||||
"\n",
|
||||
"p1 = res1['par']\n",
|
||||
"print(\"=== Pol1 ===\")\n",
|
||||
"print(f\" True: p0={true_pol1[0]:.4f} p1={true_pol1[1]:.4f}\")\n",
|
||||
"print(f\" Fitted: p0={p1[0]:.4f} p1={p1[1]:.4f}\")\n",
|
||||
"print(f\" Chi2: {res1['chi2']}\")\n",
|
||||
"print()\n",
|
||||
"\n",
|
||||
"p2 = res2['par']\n",
|
||||
"print(\"=== Pol2 ===\")\n",
|
||||
"print(f\" True: p0={true_pol2[0]:.4f} p1={true_pol2[1]:.4f} p2={true_pol2[2]:.4f}\")\n",
|
||||
"print(f\" Fitted: p0={p2[0]:.4f} p1={p2[1]:.4f} p2={p2[2]:.4f}\")\n",
|
||||
"print(f\" Chi2: {res2['chi2']}\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7e185144-3b72-41f4-80d0-3c3504e48bba",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Plot"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "9a0b3292-7c57-4483-aecc-7ef683ff0b62",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAABKUAAAHqCAYAAADVi/1VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADH/0lEQVR4nOzde3zO9f/H8cfn2uHaeYwdMHNIq2zoQAmFvo3kVDmmRDoocqyUjhQp5ZBIB/0ckkqKjpgUOuAbUbHChGnMNFw72Obars/vj4uLfZ3Zdl3bnvfbza193p/39ble117W3l6f9/v9MUzTNBERERERERERESlFFncHICIiIiIiIiIiFY+KUiIiIiIiIiIiUupUlBIRERERERERkVKnopSIiIiIiIiIiJQ6FaVERERERERERKTUqSglIiIiIiIiIiKlTkUpEREREREREREpdSpKiYiIiIiIiIhIqVNRSkRERERERERESp2KUiLi0WbNmoVhGK4/3t7eREdHc++995Kamnpe12rVqhWtWrUq0jZnzhx69uzJZZddhsVioXbt2sUXvIiIiIgblOT4ae/evTzzzDNcf/31VK1alZCQEK655hreeecdCgsLi/mTiEh55+3uAEREzsXMmTO5/PLLyc3NZdWqVYwbN46VK1fyxx9/EBgYeMHXff/990lLS+Paa6/F4XBgt9uLMWoRERER9ymJ8dP69euZM2cO99xzD88++yw+Pj4sXryYhx9+mDVr1vB///d/xfwpRKQ8U1FKRMqE+Ph4GjduDEDr1q0pLCzkxRdfZNGiRdx1110XfN2lS5disTgnjXbo0IFNmzYVS7wiIiIi7lYS46fmzZuzfft2fHx8XG0JCQkcOXKEadOmMXr0aGrWrFks8YtI+afleyJSJjVt2hSAXbt2kZeXx8iRI6lTpw6+vr7UqFGDgQMHcujQobNe51hBSkRERKS8K47xU+XKlYsUpI659tprAfjnn3+KPW4RKb80U0pEyqTk5GQAwsPDue2221i+fDkjR47khhtu4Pfff+f5559n9erVrF69GqvV6uZoRURERNyvJMdP3333Hd7e3sTGxpZE6CJSTqkoJSJlQmFhIQUFBeTl5bFy5UrGjBlDcHAwISEhLF26lPHjx/P4448DzinkNWvWpEePHsyZM4cHHnjAzdGLiIiIlL7SGj8lJiby/vvvM2TIEKpUqVJSH0dEyiGtWxGRMqFp06b4+PgQHBxMhw4diIqKYvHixfz6668A9O3bt0j/bt26ERgYyPLly90QrYiIiIj7lcb46ddff6V79+40bdqUcePGFWf4IlIBaKaUiJQJc+bM4YorrsDb25vIyEiqVasGOJ8q4+3tTXh4eJH+hmEQFRVFRkaGO8IVERERcbuSHj9t2LCBhIQELr30Ur755httmSAi501FKREpE6644grX02NOVKVKFQoKCti/f3+RgZVpmqSlpdGkSZPSDFNERETEY5Tk+GnDhg3cfPPN1KpVi8TEREJDQ4s1dhGpGLR8T0TKtP/85z8AzJ07t0j7p59+Sk5Ojuu8iIiIiDhd7Php48aN3HzzzURHR7Ns2TIqV65cYrGKSPmmmVIiUqYlJCTQtm1bnnjiCTIzM2nevLnr6TFXXXUVvXv3PuPrk5KSSEpKAiAtLY3Dhw+zYMECAOrXr0/9+vVL/DOIiIiIlKaLGT9t2bKFm2++GYCxY8eybds2tm3b5jp/ySWXnLQsUETkdFSUEpEyzTAMFi1axKhRo5g5cyZjx46latWq9O7dm5deeumsexvMnz+f0aNHF2nr1q0bAM8//zyjRo0qqdBFRERE3OJixk+rV6927TnVsWPHk87PnDnzpA3URUROxzBN03R3ECIiIiIiIiIiUrFoTykRERERERERESl1KkqJiIiIiIiIiEipU1FKRERERERERERKnYpSIiIiIiIiIiJS6lSUEhERERERERGRUqeilIiIiIiIiIiIlDpvdwfwvxwOB3v27CE4OBjDMNwdjoiIiFRgpmmSlZVFcHAwISEhHj020RhKREREPMWxMVT16tWxWE4/H8rjilJ79uyhZs2a7g5DREREpAibzUZISIi7wzgtjaFERETE0+zevZvo6OjTnve4olRwcDDgDLykBn4Oh4P9+/cTHh5+xoqdlA7lw7MoH55F+fAsyodnKY18ZGZmUrNmTXbv3u0ao1yI6dOnM336dHbu3AlAXFwczz33HO3atQOcdxNHjx7NO++8w8GDB7nuuuuYNm0acXFx5/weGkNVPMqHZ1E+PIvy4VmUD89SmmOos42fPK4odWy6eUhISIkOqPLy8ggJCdEPhAdQPjyL8uFZlA/Ponx4ltLMx8Uu3YuOjubll1+mXr16AMyePZvOnTuzYcMG4uLiGD9+PBMnTmTWrFnExsYyZswYEhIS2LJlyzkXwzSGqniUD8+ifHgW5cOzKB+epTTzcbbxk/42iIiIiJSwjh07cuuttxIbG0tsbCxjx44lKCiINWvWYJomkydP5umnn+aOO+4gPj6e2bNnc/jwYebNm+fu0EVERERKjIpSIiIiIqWosLCQjz76iJycHK6//np27NhBWloabdq0cfWxWq20bNmSn3/+2Y2RioiIiJQsj1u+JyIiIlIe/fHHH1x//fXk5eURFBTEwoULqV+/vqvwFBkZWaR/ZGQku3btOu318vPzyc/Pdx1nZmYCzin5DoejBD6B89qmaZbY9eX8KB+eRfnwLMqHZ1E+PEtp5ONcr11mi1KFhYXY7fYLeq3D4cBut5OXl1eh17P6+Pjg5eXl7jBEREQqhMsuu4yNGzdy6NAhPv30U/r06cPKlStd5/93zwXTNM+4D8O4ceMYPXr0Se379+8nLy/vtK9zOBwUFhZewCc4/nhnu91+UXtslWVeXl4eM350OBzYbDZM0/SYmCoy5cOzKB+eRfnwLKWRj6ysrHPqV+aKUqZpkpaWxqFDhy7qGg6Hg6ysrAo7oDqmUqVKREVFVfjvg4iISEnz9fV1bXTeuHFjfvnlF15//XWeeOIJANLS0qhWrZqrf3p6+kmzp040cuRIhg8f7jo+9pSb8PDwU250bpom+/btu6gx1LHrZGdnX9Q1yrpKlSoRGRnp9vGTw+HAMAw9zcpDKB+eRfnwLMqHZymNfPj5+Z1TvzJXlDpWkIqIiCAgIOCCBgOmaVJQUIC3t7fbBxPuYpomhw8fJj09HaDIIFhERERKnmma5OfnU6dOHaKioli2bBlXXXUVAEeOHGHlypW88sorp3291WrFarWe1G6xWE45wNy7dy82m43IyEiNoS7QieMnwzA8YvxkGMZpcy6lT/nwLMqHZ1E+PEtJ5+Ncr1umilKFhYWuglSVKlUu+DoVfUB1jL+/P+C8ExsREaGlfCIiIiXkqaeeol27dtSsWZOsrCw++ugjVqxYwZIlSzAMg6FDh/LSSy9x6aWXcumll/LSSy8REBBAr169iuX9NYYqPho/iYiIFJ8yVZQ6todUQECAmyMpP459L+12uwZVIiIiJWTfvn307t2bvXv3EhoaSsOGDVmyZAkJCQkAjBgxgtzcXAYMGMDBgwe57rrrSExMJDg4uFjeX2Oo4qXxk4iISPEoU0WpYyrqnbmSoO+liIhIyXvvvffOeN4wDEaNGsWoUaNKNA793i8e+j6KiIgUDy3mFBERERERERGRUqeilJu1atWKoUOHujsMERERkTJD4ycREZHyQUWpMmTFihUYhnHRj3IWERERqSg0fhIREfFcKkqJiIiIiIiIiEipq7BFqc17Mpm4bCtDPtrAhMQtbEq1lfh75uTkcM899xAUFES1atWYMGFCkfNz586lcePGBAcHExUVRa9evUhPTwdg586dtG7dGoDKlStjGAZ9+/YFYMmSJbRo0YJKlSpRpUoVOnTowPbt20v884iIiEjFsnmPjcnLkxn68UaNn0REROSiVcii1OY9NqZ8l8zq7f+SmVvA6u0ZTFq2tcQHVo8//jjff/89CxcuJDExkRUrVrB+/XrX+SNHjvDiiy/y22+/sWjRInbs2OEaONWsWZNPP/0UgC1btrB3715ef/11wDlYGz58OL/88gvLly/HYrFw++2343A4SvTziIiIuJtx8KC7Q6gwNqXamLRsG2v+PkBmrl3jJxERkTJoU6qNicu28vziHUwshd/jZ+Pt1nd3k6Wb92HLLeDSyCAMw0JkiJXk9GwSk/YRXyO0RN4zOzub9957jzlz5pCQkADA7NmziY6OdvXp16+f6+u6desyZcoUrr32WrKzswkKCiIsLAyAiIgIKlWq5OrbpUuXIu/13nvvERERQVJSEvHx8SXyeURERNzuyy8J79MH3noLevZ0dzTl3tLNadhyj3BJeCBeXhYiQ9D4SUREpAxx3mDaii33CL6GgzV/Z/Dn3iyGJcSW2O/ysznvmVKrVq2iY8eOVK9eHcMwWLRoUZHz2dnZPPLII0RHR+Pv788VV1zB9OnTiyveYpFy4DBBft4YhgGAYRgE+/mwKyOnxN5z+/btHDlyhOuvv97VFhYWxmWXXeY63rBhA507d6ZWrVoEBwfTqlUrZ7wpKWe9dq9evahbty4hISHUqVPnnF4nIiJSJhUUwJNPYrntNiw2G8YDD8CWLe6OqtxLOXCYYD8fjZ9ERETKKOcNJjv1woOICPKlXngQtlw7iUn73BbTeRelcnJyaNSoEVOnTj3l+WHDhrFkyRLmzp3Ln3/+ybBhwxg0aBCff/75RQdbXGLCAsjOK8A0TQBM0yQrz07tKoEl9p7H3ut0cnJyaNOmDUFBQcydO5dffvmFhQsXAs5p6WfSsWNHMjIyePfdd1m7di1r1649p9eJiIh4mk2pNiYkbjn9no9paXDzzfDKK66mbY2uJ8nhX8qRVjwxYQFk5dk1fhIRESmj3HGD6WzOuyjVrl07xowZwx133HHK86tXr6ZPnz60atWK2rVr8+CDD9KoUSLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1200x500 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n",
|
||||
" \n",
|
||||
"# Pol1\n",
|
||||
"ax1.plot(x, y_pol1_noisy, 'o', ms=4, alpha=0.6, label='data')\n",
|
||||
"ax1.plot(x, y_pol1, 'k--', label='truth')\n",
|
||||
"ax1.plot(x, m1(x, res1['par']), 'r-', lw=2, label='fit')\n",
|
||||
"ax1.set_title('Pol1')\n",
|
||||
"ax1.legend()\n",
|
||||
"ax1.grid(alpha=0.3)\n",
|
||||
" \n",
|
||||
"# Pol2\n",
|
||||
"ax2.plot(x, y_pol2_noisy, 'o', ms=4, alpha=0.6, label='data')\n",
|
||||
"ax2.plot(x, y_pol2, 'k--', label='truth')\n",
|
||||
"ax2.plot(x, m2(x, res2['par']), 'r-', lw=2, label='fit')\n",
|
||||
"ax2.set_title('Pol2')\n",
|
||||
"ax2.legend()\n",
|
||||
"ax2.grid(alpha=0.3)\n",
|
||||
" \n",
|
||||
"fig.tight_layout()\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e694c18a-a203-4132-a472-c93c4d3ce853",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.12.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
@@ -1,395 +1,397 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "efef8e20-6571-4561-8f0b-6048b57907de",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import time\n",
|
||||
"import random\n",
|
||||
"import numpy as np\n",
|
||||
"np.set_printoptions(suppress=True, precision=6)\n",
|
||||
"\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from scipy.special import erf\n",
|
||||
"\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0, '/home/ferjao_k/sw/aare/build')\n",
|
||||
"\n",
|
||||
"from aare import fit_scurve, fit_scurve2 # lmfit\n",
|
||||
"from aare import RisingScurve, FallingScurve, fit # minuit2 (object based API)\n",
|
||||
"\n",
|
||||
"from pprint import pprint"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8449f9e9-59ee-4194-b477-d014b8efc93b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Model curves"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "35772e2c-37c6-4986-a9c0-7dca4a034827",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def scurve(x, p): # rising Scurve\n",
|
||||
" p0, p1, p2, p3, p4, p5 = p\n",
|
||||
" return (p0 + p1*x) + 0.5 *(1 + erf((x-p2) / (np.sqrt(2)*p3))) * (p4 + p5*(x-p2))\n",
|
||||
"\n",
|
||||
"def scurve2(x, p): #falling Scurve\n",
|
||||
" p0, p1, p2, p3, p4, p5 = p\n",
|
||||
" return (p0 + p1*x) + 0.5 *(1 - erf((x-p2) / (np.sqrt(2)*p3))) * (p4 + p5*(x-p2))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a318ab93-fe18-4798-9ddb-6ffc2c2e3875",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Generate data (1D)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "06c8fddb-56d9-4f84-8b21-4ffd8b7bd26f",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"x = np.linspace(0,120, 121)\n",
|
||||
"\n",
|
||||
"rng = np.random.default_rng(42)\n",
|
||||
"0\n",
|
||||
"p_true_rising = np.array([100.0, 0.25, 60.0, 6.0, 120.0, 1.0])\n",
|
||||
"p_true_falling = np.array([100.0, 0.25, 60.0, 6.0, 120.0, -1.0])\n",
|
||||
"\n",
|
||||
"y_true_rising = scurve(x, p_true_rising)\n",
|
||||
"y_true_falling = scurve2(x, p_true_falling)\n",
|
||||
"\n",
|
||||
"noise_sigma = 4\n",
|
||||
"\n",
|
||||
"y_rising = y_true_rising + rng.normal(0, noise_sigma, size=x.shape)\n",
|
||||
"# y_err_r = np.full_like(x, noise_sigma)\n",
|
||||
"\n",
|
||||
"y_falling = y_true_falling + rng.normal(0, noise_sigma, size=x.shape)\n",
|
||||
"# y_err_f = np.full_like(x, noise_sigma)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7e185144-3b72-41f4-80d0-3c3504e48bba",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Plot synthetic data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "9a0b3292-7c57-4483-aecc-7ef683ff0b62",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAF0CAYAAADcnkHwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACSjElEQVR4nOzdeXwTZf4H8M8kbdP7phcNpSBHS7lBCogUQS4BURBUUFFg/SmwssDqIiqIIh67HosrqwsCUhBPQAVB7kPkKmdpuXvSC3qlV9I0md8fITFpkp5p0+Pzfr2izWQyeSYkmfnO8zzfryCKoggiIiIiIiIisgmJvRtARERERERE1JIw0CYiIiIiIiKyIQbaRERERERERDbEQJuIiIiIiIjIhhhoExEREREREdkQA20iIiIiIiIiG2KgTURERERERGRDDLSJiIiIiIiIbIiBNhEREREREZENMdAmsmL9+vUQBMFwc3BwQHBwMB5//HFcu3bNbP2YmBjExMTU6jUOHjwIQRBw8OBB2zS6lhITE/HUU0+hQ4cOcHZ2hr+/P/r06YO5c+dCoVDYpU1ERESNofJx3vi2aNGiGm8nOTkZgiBg/fr1ZttOTk42LJsxYwbat29vux2oBVEUsWXLFgwZMgQBAQFwdnZGaGgoRo0ahTVr1tilTUQtnYO9G0DU1K1btw5du3aFUqnE77//jhUrVuDAgQO4fPkyfHx8DOt99tlntd52nz598McffyAyMtKWTa6Rs2fPYvDgwYiIiMAbb7yB9u3b486dOzh//jy2bNmCRYsWwdPTs9HbRURE1Jj0x3ljISEhNn+d119/HS+99JLNt1sTixcvxnvvvYfZs2fj73//Ozw8PJCSkoL9+/dj+/btmDVrll3aRdSSMdAmqkZUVBT69esHQNdrrdFosHTpUmzbtg3PPvusYb26BMuenp6Ijo62WVtr4+OPP4ZEIsHBgwfh4eFhWD558mS89dZbEEXRLu0yJooilEolXFxc7N0UIiJqoYyP8w2pY8eODf4alpSVleHjjz/G008/jS+++MLksRkzZkCr1dqlXZWVlZXB2dkZgiDYuylENsGh40S1pD8YZ2dnmyy3NHR89erV6NmzJ9zd3eHh4YGuXbvi1VdfNTxuaej4jBkz4O7ujuvXr2Ps2LFwd3eHXC7HwoULoVKpTLafnp6OyZMnw8PDA97e3pg2bRpOnTplNoTNktzcXHh6esLd3d3i4zU90J04cQLjx4+Hn58fnJ2d0bFjR8yfP99kfywNlVu2bJnZawiCgLlz5+K///0vIiIiIJPJsGbNGgQEBOCpp54y20ZBQQFcXFywYMECwzKFQoFFixYhPDwcTk5OaNu2LebPn4+SkpIa7Q8REREAXL9+Hc8++yw6deoEV1dXtG3bFuPHj8fFixfrtD1Lx0P9cW/jxo2IiIiAq6srevbsiV9++cXs+du3b0ePHj0gk8nQoUMHfPLJJxaPpZWVlJRApVIhODjY4uMSSc3CAa1Wi1WrVqFXr15wcXGBt7c3oqOj8dNPP5nsz7Jly8ye2759e8yYMcNwXz+0/rfffsNzzz2HNm3awNXVFd988w0EQcC+ffvMtrF69WoIgoALFy4Ylp0+fRoTJkyAr68vnJ2d0bt3b3z77bc12h+ihsYebaJaSkpKAgB07ty5yvW2bNmCF198EfPmzcM///lPSCQSXL9+HQkJCdW+hlqtxoQJEzBz5kwsXLgQhw8fxltvvQUvLy+88cYbAHQHzmHDhiEvLw/vvfce7rnnHuzatQtTp06t0X4MHDgQO3bswLRp0/D888/j3nvvrXXP8e7duzF+/HhERETgww8/RLt27ZCcnIzffvutVtsxtm3bNhw5cgRvvPEGgoKCEBAQgKSkJPz3v//Ff/7zH5Ph7F9//TWUSqVhZEFpaSmGDh2K9PR0vPrqq+jRowcuXbqEN954AxcvXsTevXt5pZyIiExoNBpUVFSYLHNwcEBGRgb8/Pzw7rvvok2bNsjLy8OGDRswYMAAnD17Fl26dLHJ6+/YsQOnTp3C8uXL4e7ujvfffx+PPPIIrly5gg4dOgAAdu3ahUcffRT3338/vvnmG1RUVOCf//yn2UV/S/z9/XHPPffgs88+Q0BAAMaOHYsuXbrU+ng4Y8YMxMbGYubMmVi+fDmcnJxw5swZk3notfXcc8/hoYcewsaNG1FSUoJx48YhICAA69atw/Dhw03WXb9+Pfr06YMePXoAAA4cOIDRo0djwIAB+O9//wsvLy9s2bIFU6dORWlpqUlgT2QXIhFZtG7dOhGAePz4cVGtVotFRUXirl27xKCgIPH+++8X1Wq1yfpDhw4Vhw4darg/d+5c0dvbu8rXOHDggAhAPHDggGHZM888IwIQv/32W5N1x44dK3bp0sVw/z//+Y8IQPz1119N1nv++edFAOK6deuqfG2lUilOnDhRBCACEKVSqdi7d29xyZIlYk5OTpXP1evYsaPYsWNHsayszOo6zzzzjBgWFma2fOnSpWLlnyAAopeXl5iXl2ey/MKFCyIA8YsvvjBZfu+994p9+/Y13F+5cqUokUjEU6dOmaz3/fffiwDEnTt31mi/iIio5dMf5y3dKh/jRVEUKyoqxPLycrFTp07i3/72N8PypKQks+OufttJSUmGZZaOhwDEwMBAUaFQGJZlZWWJEolEXLlypWFZ//79RblcLqpUKsOyoqIi0c/Pz+xYasnJkyfFdu3aGfbPw8NDHDdunPjVV1+JWq222ucfPnxYBCAuWbKkyvUAiEuXLjVbHhYWJj7zzDOG+/r35+mnnzZbd8GCBaKLi4tYUFBgWJaQkCACEFetWmVY1rVrV7F3795m/1bjxo0Tg4ODRY1GU+1+ETUkDh0nqkZ0dDQcHR3h4eGB0aNHw8fHB9u3b4eDQ9UDQu69914UFBTgiSeewPbt23Hnzp0av6YgCBg/frzJsh49eiAlJcVw/9ChQ4Y2GXviiSdq9BoymQxbt25FQkICPvroIzz++OO4ffs2VqxYgYiICFy5cgWAbp50RUWFyQ0Arl69ihs3bmDmzJlwdnau8b5V54EHHjBJMgcA3bt3R9++fbFu3TrDssTERJw8eRLPPfecYdkvv/yCqKgo9OrVy6S9o0aNsmt2dyIiarq++uornDp1yuTm4OCAiooKvPPOO4iMjISTkxMcHBzg5OSEa9euITEx0WavP2zYMJNcKYGBgQgICDAc80tKSnD69GlMnDgRTk5OhvXc3d3NzhWs6d+/P65fv45du3bh1VdfxcCBA7Fv3z48/fTTmDBhgiEvi753X3/Tz9/+9ddfAQBz5syxyT7rTZo0yWzZc889h7KyMnzzzTeGZevWrYNMJsOTTz4JQDes//Lly5g2bRoAmLR57NixyMzMNJzHENkLA22iaugPwPv378fzzz+PxMTEGgWzTz31FL788kukpKRg0qRJCAgIwIABA7Bnz55qn+vq6moWvMpkMiiVSsP93NxcBAYGmj3X0rKqREREYP78+YiNjUVqaio+/PBD5Obm4vXXXwcAbNiwAY6OjiY3ALh9+zYAIDQ0tFavVx1rc8iee+45/PHHH7h8+TKAPw+6xv8W2dnZuHDhgll7PTw8IIpirS52EBFR6xAREYF+/fqZ3ABgwYIFeP311zFx4kT8/PPPOHHiBE6dOoWePXuirKzMZq/v5+dntkwmkxleIz8/H6Io1vuY7+joiFGjRmHFihXYvXs30tLSEBMTg19++cUQSA8fPtzk+Km/mH379m1IpVIEBQXVZRetsnTM79atG/r372+4uK7RaBAbG4uHH34Yvr6+AP7Mk7No0SKzY/6LL74IADzmk91xjjZRNfQHYEB31Vmj0WDNmjX4/vvvMXny5Cqf++yzz+LZZ59FSUkJDh8+jKVLl2LcuHG4evUqwsLC6tUuPz8/nDx50mx5VlZWnbcpCAL+9re/Yfny5YiPjwcAjB8/HqdOnTJbt02bNgB0Cdmq4uzsbJbEDbB+ALQ2Z+yJJ57AggULsH79eqxYsQIbN27ExIkTTXq//f394eLigi+//NLiNvz9/atsKxERkV5sbCyefvppvPPOOybL79y5A29v70Zrh4+PDwRBsDgfuz7HfD8/P8yfPx8HDx5EfHw8xo4di88//xxFRUWGdfTHzTZt2kCj0SArK8vqBXFAd4HA0jE/NzfX4vrWjvnPPvssXnzxRSQmJuLmzZvIzMw0qfSib9fixYvx6KOPWtyGrebQE9UVe7SJaun999+Hj48P3njjjRqXxHBzc8OYMWOwZMkSlJeX49KlS/Vux9ChQ1FUVGS4Cq23ZcuWGj0/MzPT4vKMjAwoFApDDVE/Pz+LV/o7d+6Mjh074ssvv7R4UNVr3749cnJyTE4QysvLsXv37hq1U8/HxwcTJ07EV199hV9++QVZWVkmw8YBYNy4cbhx44bFNvfr189i9nMiIiJLBEGATCYzWbZjxw7cunWrUdvh5uaGfv36Ydu2bSgvLzcsLy4utpidvDK1Wm010NUPgdcf87t06WLxuDlmzBgAuszfVWnfvr1JVnAA2L9/P4qLi6ttp7EnnngCzs7OWL9+PdavX4+2bdti5MiRhse7dOmCTp064fz58xaP9/369TMZjk9kD+zRJqolHx8fLF68GC+//DI2b96M6dOnW1xv9uzZcHFxweDBgxEcHIysrCysXLkSXl5e6N+/f73b8cwzz+Cjjz7C9OnT8fbbb+Oee+7Br7/+aghgqyvX8Ze//AUFBQWYNGkSoqKiIJVKcfnyZXz00UeQSCR45ZVXqm3Df/7zH4wfPx7R0dH429/+hnbt2iE1NRW7d+/Gpk2bAABTp07FG2+8gccffxx///vfoVQq8e9//xsajabW+/zcc8/hm2++wdy5cxEaGooRI0aYPD5//nz88MMPuP/++/G3v/0NPXr0gFarRWpqKn777TcsXLgQAwYMqPXrEhFR6zNu3DisX78Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1200x400 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(1, 2, figsize=(12,4), sharex=True)\n",
|
||||
"\n",
|
||||
"ax[0].plot(x, y_true_rising, label=\"true\")\n",
|
||||
"ax[0].scatter(x, y_rising, s=12, label=\"noisy\")\n",
|
||||
"ax[0].set_title(\"Rising S-curve\")\n",
|
||||
"ax[0].legend()\n",
|
||||
"\n",
|
||||
"ax[1].plot(x, y_true_falling, label=\"true\")\n",
|
||||
"ax[1].scatter(x, y_falling, s=12, label=\"noisy\")\n",
|
||||
"ax[1].set_title(\"Falling S-curve\")\n",
|
||||
"ax[1].legend()\n",
|
||||
"\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "95c826d1-9223-4fb2-9ada-d9fc1acad2fb",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Fit the Rising S-curve"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"id": "37d2fda4-40b8-4ea2-8f2c-6c3f2384d158",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Paramter list of the Scurve:\n",
|
||||
"['p0', 'p1', 'mu', 'sigma', 'A', 'C']\n",
|
||||
"6 \n",
|
||||
"\n",
|
||||
"== Tuned fit settings ==\n",
|
||||
"max_calls : 500\n",
|
||||
"tolerance : 0.01\n",
|
||||
"compute_erros: True\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"== Results ==\n",
|
||||
"True rising params: [100. 0.25 60. 6. 120. 1. ]\n",
|
||||
"lmfit rising result: [ 99.467831 0.287017 60.078405 5.715461 117.257939 0.967234]\n",
|
||||
"minuit_grad (obj api ) rising result:\n",
|
||||
"{'par': array([100. , 0.269137, 60.050193, 5.729932, 117.777635,\n",
|
||||
" 0.984073]),\n",
|
||||
" 'par_err': array([1. , 0.054406, 0.747753, 0.909864, 7.887499, 0.174058]),\n",
|
||||
" 'chi2': array([7.351743])}\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"res_r_lmfit = fit_scurve(x, y_rising)\n",
|
||||
"\n",
|
||||
"model_r = RisingScurve()\n",
|
||||
"print(\"Paramter list of the Scurve:\")\n",
|
||||
"# print(model_r.GetParNames())\n",
|
||||
"print(model_r.par_names)\n",
|
||||
"print(model_r.n_par, '\\n')\n",
|
||||
"\n",
|
||||
"model_r.FixParameter('p0', 100) # Fixed p0 = 100 (optimizer does not touch the value)\n",
|
||||
"model_r.SetParameter('A', 100) # Start value A = 100\n",
|
||||
"\n",
|
||||
"print(\"== Tuned fit settings ==\")\n",
|
||||
"model_r.compute_errors = True\n",
|
||||
"model_r.max_calls = 500\n",
|
||||
"model_r.tolerance = 0.01\n",
|
||||
"print(f\"max_calls : {model_r.max_calls}\")\n",
|
||||
"print(f\"tolerance : {model_r.tolerance}\")\n",
|
||||
"print(f\"compute_erros: {model_r.compute_errors}\")\n",
|
||||
"print(\"\\n\")\n",
|
||||
"print(\"== Results ==\")\n",
|
||||
"res_r_munuit_obj = model_r.fit(x, y_rising, np.sqrt(y_rising))\n",
|
||||
"\n",
|
||||
"print(\"True rising params: \", p_true_rising)\n",
|
||||
"print(\"lmfit rising result: \", res_r_lmfit)\n",
|
||||
"print(\"minuit_grad (obj api ) rising result:\")\n",
|
||||
"pprint(res_r_munuit_obj, sort_dicts=False)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5ebf1328-954a-43ce-953b-c402fd7998f9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Fit the Falling S-curve"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "5ce22d2d-1adb-47d9-bf36-6cb535fa8b77",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"True falling params: [100. 0.25 60. 6. 120. -1. ]\n",
|
||||
"lmfit falling result: [101.73448 0.225786 60.251971 6.092915 118.476735 -0.998677]\n",
|
||||
"minuit_grad (obj api ) falling result\n",
|
||||
"{'par': array([101.741754, 0.225723, 60.251604, 6.092887, 118.473268,\n",
|
||||
" -0.998565]),\n",
|
||||
" 'chi2': array([1865.728797])}\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"res_f_lmfit = fit_scurve2(x, y_falling)\n",
|
||||
"\n",
|
||||
"model_f = FallingScurve()\n",
|
||||
"model_f.SetParameter(2, 20) # set start p2 = 20\n",
|
||||
"model_f.SetParameter(4, 90) # set start p4 = 90\n",
|
||||
"model_f.SetParLimits(4, 80, 130)\n",
|
||||
"model_f.max_calls = 150 # this is limits the number of calls done by the minimizer (increasing the `max_calls` may help if the fit fails) \n",
|
||||
"res_f_munuit_obj = model_f.fit(x, y_falling)\n",
|
||||
"\n",
|
||||
"print(\"True falling params: \", p_true_falling)\n",
|
||||
"print(\"lmfit falling result: \", res_f_lmfit)\n",
|
||||
"print(\"minuit_grad (obj api ) falling result\")\n",
|
||||
"pprint(res_f_munuit_obj, sort_dicts=False)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"id": "fe72503a-fc77-4b17-9b17-b98055e2fc87",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAF0CAYAAADcnkHwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADV1UlEQVR4nOzdeZxN5R/A8c+562x39hVjLMmuiEIKJSJCREkIqSwltIhsJSU/8lO0iCGJUpKSfUk/KaFkJ8Y6YzD7cuduz++Pa64ZZqzDDL7v1+u+uOeee85z7tyZc77neZ7vV1NKKYQQQgghhBBCCFEkdMXdACGEEEIIIYQQ4mYigbYQQgghhBBCCFGEJNAWQgghhBBCCCGKkATaQgghhBBCCCFEEZJAWwghhBBCCCGEKEISaAshhBBCCCGEEEVIAm0hhBBCCCGEEKIISaAthBBCCCGEEEIUIQm0hRBCCCGEEEKIIiSBtrhpxcbGomma52EwGIiKiuKJJ55g3759563fpEkTmjRpcln7WLt2LZqmsXbt2qJp9GXatWsXTz/9NBUqVMDLy4vQ0FDq1KlD//79SUtLK5Y2Fbfhw4dTtmxZDAYDgYGBwPk/26ysLEaNGlVsPzchhBBX79zzfN7HkCFDLnk7cXFxaJpGbGzseduOi4vzLOvRowflypUrugO4DEop5s2bx3333Ud4eDheXl6UKVOGFi1aMH369GJpU3GLi4vjkUceITg4GE3TGDhwYIE/yw0bNjBq1ChSUlKKra3i1mQo7gYIca3NnDmTKlWqYLVa+d///sfYsWNZs2YNu3fvJigoyLPe1KlTL3vbderU4bfffqNatWpF2eRLsnXrVu69916qVq3KiBEjKFeuHKdOneLvv/9m3rx5DBkyBH9//+veruK0aNEixo4dy7Bhw2jZsiVmsxk4/2eblZXF6NGjAS775ooQQoiSJfc8n1epUqWKfD9vvvkmL730UpFv91IMHTqU9957j2effZZXXnkFi8XCoUOHWL16NYsWLaJ3797F0q7i9PLLL/P7778zY8YMIiMjiYqKIjIykt9++42KFSt61tuwYQOjR4+mR48enhvwQlwPEmiLm16NGjWoW7cu4A6qnE4nI0eO5Pvvv+eZZ57xrHclwbK/vz/169cvsrZejg8++ACdTsfatWuxWCye5R07duStt95CKVUs7cpLKYXVasXb2/u67G/79u0AvPjii4SHh3uWF8eNECGEENdH3vP8tZQ3eLuesrOz+eCDD+jWrRuffvppvtd69OiBy+UqlnadKzs7Gy8vLzRNuy772759O3fffTft2rXLt7y4rsuEOJcMHRe3nNyT8YkTJ/ItL2jo+LRp07jjjjvw8/PDYrFQpUoV3njjDc/rBQ0d79GjB35+fuzfv59WrVrh5+dHdHQ0gwcPJicnJ9/2jx49SseOHbFYLAQGBvLUU0+xadOm84Y9FeT06dP4+/vj5+dX4OuXeqL7/fffadOmDSEhIXh5eVGxYkUGDhyY73gKGio3atSo8/ahaRr9+/fn448/pmrVqpjNZqZPn054eDhPP/30edtISUnB29ubQYMGeZalpaUxZMgQypcvj8lkonTp0gwcOJDMzMwLHke5cuUYPnw4ABEREWiaxqhRo4D8P9u4uDjCwsIAGD16tGeYYY8ePQA4efIkffr0ITo6GrPZTFhYGPfeey8rV6684P6FEEKULPv37+eZZ56hUqVK+Pj4ULp0adq0acM///xzRdsr6HyYe9774osvqFq1Kj4+Ptxxxx38+OOP571/0aJF1KpVC7PZTIUKFZg8eXKB59JzZWZmkpOTQ1RUVIGv63SXdjnvcrmYMmUKd955J97e3gQGBlK/fn1++OGHfMeTe+7Mq1y5cp7zJJwdWr98+XJ69uxJWFgYPj4+zJ8/H03TWLVq1XnbmDZtGpqmsW3bNs+yP//8k0cffZTg4GC8vLyoXbs2X3/99QWPI/faa//+/fz888+e83hcXNx5Q8dHjRrFK6+8AkD58uU96+Zet61evZomTZoQEhKCt7c3ZcuWpUOHDmRlZV3SZyrEhUiPtrjlHDx4EIDbb7/9guvNmzePvn37MmDAACZMmIBOp2P//v3s3Lnzovuw2+08+uij9OrVi8GDB/PLL7/w1ltvERAQwIgRIwD3ibNp06YkJSXx3nvvcdttt7F06VI6d+58ScfRoEEDfvrpJ5566imee+457r777svuOV62bBlt2rShatWqTJw4kbJlyxIXF8fy5csvazt5ff/996xfv54RI0YQGRlJeHg4Bw8e5OOPP+ajjz7KN5z9q6++wmq1ekYWZGVl0bhxY44ePcobb7xBrVq12LFjByNGjOCff/5h5cqVhV6QLFy4kI8++ojPP/+cpUuXEhAQQJkyZc5bLyoqiqVLl/Lwww/Tq1cvz3C73OD76aefZsuWLYwdO5bbb7+dlJQUtmzZwunTp6/4MxFCCHHtOJ1OHA5HvmUGg4Hjx48TEhLCu+++S1hYGElJScyaNYt77rmHrVu3Urly5SLZ/08//cSmTZsYM2YMfn5+jB8/nvbt27Nnzx4qVKgAwNKlS3nssce4//77mT9/Pg6HgwkTJpx3078goaGh3HbbbUydOpXw8HBatWpF5cqVL7vnuEePHsyZM4devXoxZswYTCYTW7ZsyTcP/XL17NmTRx55hC+++ILMzExat25NeHg4M2fO5MEHH8y3bmxsLHXq1KFWrVoArFmzhocffph77rmHjz/+mICAAObNm0fnzp3JysrKF9jnlTttr3379lSsWJEJEyYA7vN7fHx8vnV79+5NUlISU6ZM4bvvvvPcrKhWrZpnjvd9993HjBkzCAwM5NixYyxduhSbzYaPj88Vfy5CAKCEuEnNnDlTAWrjxo3Kbrer9PR0tXTpUhUZGanuv/9+Zbfb863fuHFj1bhxY8/z/v37q8DAwAvuY82aNQpQa9as8Szr3r27AtTXX3+db91WrVqpypUre55/9NFHClA///xzvvWee+45BaiZM2decN9Wq1W1a9dOAQpQer1e1a5dWw0bNkwlJiZe8L25KlasqCpWrKiys7MLXad79+4qJibmvOUjR45U5/4JAVRAQIBKSkrKt3zbtm0KUJ9++mm+5Xfffbe66667PM/HjRundDqd2rRpU771FixYoAC1ZMmSCx5PbptOnjyZb/m5P9uTJ08qQI0cOfK8bfj5+amBAwdecD9CCCGKX+55vqDHued4pZRyOBzKZrOpSpUqqZdfftmz/ODBg+edd3O3ffDgQc+ygs6HgIqIiFBpaWmeZQkJCUqn06lx48Z5ltWrV09FR0ernJwcz7L09HQVEhJy3rm0IH/88YcqW7as5/gsFotq3bq1mj17tnK5XBd9/y+//KIANWzYsAuuV9i5MSYmRnXv3t3zPPfz6dat23nrDho0SHl7e6uUlBTPsp07dypATZkyxbOsSpUqqnbt2uf9rFq3bq2ioqKU0+m8YFtjYmLUI488km9ZQT/L999//7yfpVJnry3++uuvC+5HiCslQ8fFTa9+/foYjUYsFgsPP/wwQUFBLFq0CIPhwgM67r77blJSUnjyySdZtGgRp06duuR9appGmzZt8i2rVasWhw4d8jxft26dp015Pfnkk5e0D7PZzMKFC9m5cyeTJk3iiSee4OTJk4wdO5aqVauyZ88ewD1P2uFw5HsA7N27l3///ZdevXrh5eV1ycd2MQ888EC+JHMANWvW5K677mLmzJmeZbt27eKPP/6gZ8+enmU//vgjNWrU4M4778zX3hYtWly37O533303sbGxvP3222zcuBG73X7N9ymEEOLKzZ49m02bNuV7GAwGHA4H77zzDtWqVcNkMmEwGDCZTOzbt49du3YV2f6bNm2aL1dKREQE4eHhnnN+ZmYmf/75J+3atcNkMnnW8/PzO+9aoTD16tVj//79LF26lDfeeIMGDRqwatUqunXrxqOPPurJy5Lbu5/7yJ2//fPPPwPQr1+/IjnmXB06dDhvWc+ePcnOzmb+/PmeZTNnzsRsNtOlSxfAPax/9+7dPPXUUwD52tyqVSvi4+M91zHXyp133onJZKJPnz7MmjWLAwcOXNP9iVuPBNrippd7Al69ejXPPfccu3btuqRg9umnn2bGjBkcOnSIDh06EB4ezj333MOKFSsu+l4fH5/zglez2YzVavU8P336NBEREee9t6BlF1K1alUGDhzInDlzOHz4MBMnTuT06dO8+eabAMyaNQuj0ZjvAe65yECBw6uvRmFzyHr27Mlvv/3G7t27gbMn3bw/ixMnTrBt27bz2muxWFBKXdbNjis1f/58unfvzvTp02nQoAHBwcF069aNhISEa75vIYQQl69q1arUrVs33wNg0KBBvPnmm7Rr147Fixfz+++/s2nTJu644w6ys7OLbP8hISHnLTObzZ59JCcno5S66nO+0WikRYsWjB07lmXLlnHkyBGaNGnCjz/+6AmkH3zwwXznz9yb2SdPnkSv1xMZGXklh1iogs751atXp169ep6b606nkzlz5tC2bVuCg4OBs3lyhgwZct45v2/fvgDX/JxfsWJFVq5cSXh4OP369aNixYpUrFiRyZMnX9P9iluHzNEWN73cEzC47zo7nU6mT5/OggUL6Nix4wXf+8wzz/DMM8+QmZnJL7/8wsiRI2ndujV79+4lJibmqtoVEhLCH3/8cd7yqwnoNE3j5ZdfZsyYMZ4M3G3atGHTpk3nrZs7J/no0aMX3KaXl9d5Sdyg8BNgYXPGnnzySQYNGkRsbCxjx47liy++oF27dvl6v0NDQ/H29mbGjBkFbiM0NPSCbS0KoaGhfPDBB3zwwQccPnyYH374gddff53ExESWLl16zfcvhBCiaMyZM4du3brxzjvv5Ft+6tSp61rmKSgoCE3TCpyPfTXn/JCQEAYOHMjatWvZvn07rVq14pNPPiE9Pd2zTu55MywsDKfLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1200x400 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(1, 2, figsize=(12,4))\n",
|
||||
"\n",
|
||||
"# rising\n",
|
||||
"ax[0].plot(x, y_true_rising, label=\"true\")\n",
|
||||
"ax[0].scatter(x, y_rising, s=12, label=\"data\")\n",
|
||||
"ax[0].scatter(x, scurve(x, res_r_lmfit[:6]), s=20, marker=\"x\", label=\"lmfit\")\n",
|
||||
"ax[0].plot(x, model_r(x, res_r_munuit_obj['par']), linewidth=1, color=\"green\", label=\"minuit\")\n",
|
||||
"ax[0].set_title(\"Rising S-curve fits\")\n",
|
||||
"ax[0].legend()\n",
|
||||
"\n",
|
||||
"# falling\n",
|
||||
"ax[1].plot(x, y_true_falling, label=\"true\")\n",
|
||||
"ax[1].scatter(x, y_falling, s=12, label=\"data\")\n",
|
||||
"ax[1].scatter(x, scurve2(x, res_f_lmfit[:6]), s=20, marker=\"x\", label=\"lmfit\")\n",
|
||||
"ax[1].plot(x, model_f(x, res_f_munuit_obj['par']), linewidth=1, color=\"green\", label=\"minuit\")\n",
|
||||
"ax[1].set_title(\"Falling S-curve fits\")\n",
|
||||
"ax[1].legend()\n",
|
||||
"\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "2c42ad91-fd59-47ba-b90d-8300f97eb27a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Quick error check"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"id": "3ba40e9a-0c64-4ef6-ba90-b1116f99b862",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Rising abs error lmfit : [0.532169 0.037017 0.078405 0.284539 2.742061 0.032766]\n",
|
||||
"Rising abs error minuit_grad: [0. 0.019137 0.050193 0.270068 2.222365 0.015927]\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Falling abs error lmfit : [1.73448 0.024214 0.251971 0.092915 1.523265 0.001323]\n",
|
||||
"Falling abs error minuit_grad: [1.741754 0.024277 0.251604 0.092887 1.526732 0.001435]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(\"Rising abs error lmfit : \", np.abs(res_r_lmfit[:6] - p_true_rising))\n",
|
||||
"print(\"Rising abs error minuit_grad: \", np.abs(res_r_munuit_obj['par'] - p_true_rising))\n",
|
||||
"print(\"\\n\")\n",
|
||||
"print(\"Falling abs error lmfit : \", np.abs(res_f_lmfit[:6] - p_true_falling))\n",
|
||||
"print(\"Falling abs error minuit_grad: \", np.abs(res_f_munuit_obj['par'] - p_true_falling))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c9bbabfd-e24e-43ec-8622-983b3757d7a9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Benchmark"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"id": "dc0b52b6-9fc3-453d-b6eb-e325a2b8342e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Rising lmfit : 0.334 ms\n",
|
||||
"Rising minuit2 : 0.239 ms\n",
|
||||
"\n",
|
||||
"Falling lmfit : 0.294 ms\n",
|
||||
"Falling minuit2 : 0.214 ms\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def bench(fn, n_repeats=200):\n",
|
||||
" # warmup\n",
|
||||
" for _ in range(3):\n",
|
||||
" fn()\n",
|
||||
"\n",
|
||||
" t0 = time.perf_counter()\n",
|
||||
" for _ in range(n_repeats):\n",
|
||||
" res = fn()\n",
|
||||
" t1 = time.perf_counter()\n",
|
||||
" return res, (t1 - t0) / n_repeats\n",
|
||||
"\n",
|
||||
"model_rising = RisingScurve()\n",
|
||||
"model_falling = FallingScurve()\n",
|
||||
"\n",
|
||||
"res_r_lmfit, t_r_lmfit = bench(lambda: fit_scurve(x, y_rising))\n",
|
||||
"res_r_minuit, t_r_minuit_grad = bench(lambda: model_rising.fit(x, y_rising))\n",
|
||||
"\n",
|
||||
"res_f_lmfit, t_f_lmfit = bench(lambda: fit_scurve2(x, y_falling))\n",
|
||||
"res_f_minuit, t_f_minuit_grad = bench(lambda: model_falling.fit(x, y_falling))\n",
|
||||
"\n",
|
||||
"print(f\"Rising lmfit : {1e3*t_r_lmfit:.3f} ms\")\n",
|
||||
"print(f\"Rising minuit2 : {1e3*t_r_minuit_grad:.3f} ms\")\n",
|
||||
"print()\n",
|
||||
"print(f\"Falling lmfit : {1e3*t_f_lmfit:.3f} ms\")\n",
|
||||
"print(f\"Falling minuit2 : {1e3*t_f_minuit_grad:.3f} ms\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e694c18a-a203-4132-a472-c93c4d3ce853",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.12.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "efef8e20-6571-4561-8f0b-6048b57907de",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import time\n",
|
||||
"import random\n",
|
||||
"import numpy as np\n",
|
||||
"np.set_printoptions(suppress=True, precision=6)\n",
|
||||
"\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from scipy.special import erf\n",
|
||||
"\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0, '/home/ferjao_k/sw/aare/build')\n",
|
||||
"\n",
|
||||
"from aare import fit_scurve, fit_scurve2 # lmfit\n",
|
||||
"from aare import RisingScurve, FallingScurve, fit # minuit2 (object based API)\n",
|
||||
"\n",
|
||||
"from pprint import pprint"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8449f9e9-59ee-4194-b477-d014b8efc93b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Model curves"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "35772e2c-37c6-4986-a9c0-7dca4a034827",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def scurve(x, p): # rising Scurve\n",
|
||||
" p0, p1, p2, p3, p4, p5 = p\n",
|
||||
" return (p0 + p1*x) + 0.5 *(1 + erf((x-p2) / (np.sqrt(2)*p3))) * (p4 + p5*(x-p2))\n",
|
||||
"\n",
|
||||
"def scurve2(x, p): #falling Scurve\n",
|
||||
" p0, p1, p2, p3, p4, p5 = p\n",
|
||||
" return (p0 + p1*x) + 0.5 *(1 - erf((x-p2) / (np.sqrt(2)*p3))) * (p4 + p5*(x-p2))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a318ab93-fe18-4798-9ddb-6ffc2c2e3875",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Generate data (1D)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "06c8fddb-56d9-4f84-8b21-4ffd8b7bd26f",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"x = np.linspace(0,120, 121)\n",
|
||||
"\n",
|
||||
"rng = np.random.default_rng(42)\n",
|
||||
"0\n",
|
||||
"p_true_rising = np.array([100.0, 0.25, 60.0, 6.0, 120.0, 1.0])\n",
|
||||
"p_true_falling = np.array([100.0, 0.25, 60.0, 6.0, 120.0, -1.0])\n",
|
||||
"\n",
|
||||
"y_true_rising = scurve(x, p_true_rising)\n",
|
||||
"y_true_falling = scurve2(x, p_true_falling)\n",
|
||||
"\n",
|
||||
"noise_sigma = 4\n",
|
||||
"\n",
|
||||
"y_rising = y_true_rising + rng.normal(0, noise_sigma, size=x.shape)\n",
|
||||
"# y_err_r = np.full_like(x, noise_sigma)\n",
|
||||
"\n",
|
||||
"y_falling = y_true_falling + rng.normal(0, noise_sigma, size=x.shape)\n",
|
||||
"# y_err_f = np.full_like(x, noise_sigma)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7e185144-3b72-41f4-80d0-3c3504e48bba",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Plot synthetic data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "9a0b3292-7c57-4483-aecc-7ef683ff0b62",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAAA9oAAAF0CAYAAADcnkHwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACSjElEQVR4nOzdeXwTZf4H8M8kbdP7phcNpSBHS7lBCogUQS4BURBUUFFg/SmwssDqIiqIIh67HosrqwsCUhBPQAVB7kPkKmdpuXvSC3qlV9I0md8fITFpkp5p0+Pzfr2izWQyeSYkmfnO8zzfryCKoggiIiIiIiIisgmJvRtARERERERE1JIw0CYiIiIiIiKyIQbaRERERERERDbEQJuIiIiIiIjIhhhoExEREREREdkQA20iIiIiIiIiG2KgTURERERERGRDDLSJiIiIiIiIbIiBNhEREREREZENMdAmsmL9+vUQBMFwc3BwQHBwMB5//HFcu3bNbP2YmBjExMTU6jUOHjwIQRBw8OBB2zS6lhITE/HUU0+hQ4cOcHZ2hr+/P/r06YO5c+dCoVDYpU1ERESNofJx3vi2aNGiGm8nOTkZgiBg/fr1ZttOTk42LJsxYwbat29vux2oBVEUsWXLFgwZMgQBAQFwdnZGaGgoRo0ahTVr1tilTUQtnYO9G0DU1K1btw5du3aFUqnE77//jhUrVuDAgQO4fPkyfHx8DOt99tlntd52nz598McffyAyMtKWTa6Rs2fPYvDgwYiIiMAbb7yB9u3b486dOzh//jy2bNmCRYsWwdPTs9HbRURE1Jj0x3ljISEhNn+d119/HS+99JLNt1sTixcvxnvvvYfZs2fj73//Ozw8PJCSkoL9+/dj+/btmDVrll3aRdSSMdAmqkZUVBT69esHQNdrrdFosHTpUmzbtg3PPvusYb26BMuenp6Ijo62WVtr4+OPP4ZEIsHBgwfh4eFhWD558mS89dZbEEXRLu0yJooilEolXFxc7N0UIiJqoYyP8w2pY8eODf4alpSVleHjjz/G008/jS+++MLksRkzZkCr1dqlXZWVlZXB2dkZgiDYuylENsGh40S1pD8YZ2dnmyy3NHR89erV6NmzJ9zd3eHh4YGuXbvi1VdfNTxuaej4jBkz4O7ujuvXr2Ps2LFwd3eHXC7HwoULoVKpTLafnp6OyZMnw8PDA97e3pg2bRpOnTplNoTNktzcXHh6esLd3d3i4zU90J04cQLjx4+Hn58fnJ2d0bFjR8yfP99kfywNlVu2bJnZawiCgLlz5+K///0vIiIiIJPJsGbNGgQEBOCpp54y20ZBQQFcXFywYMECwzKFQoFFixYhPDwcTk5OaNu2LebPn4+SkpIa7Q8REREAXL9+Hc8++yw6deoEV1dXtG3bFuPHj8fFixfrtD1Lx0P9cW/jxo2IiIiAq6srevbsiV9++cXs+du3b0ePHj0gk8nQoUMHfPLJJxaPpZWVlJRApVIhODjY4uMSSc3CAa1Wi1WrVqFXr15wcXGBt7c3oqOj8dNPP5nsz7Jly8ye2759e8yYMcNwXz+0/rfffsNzzz2HNm3awNXVFd988w0EQcC+ffvMtrF69WoIgoALFy4Ylp0+fRoTJkyAr68vnJ2d0bt3b3z77bc12h+ihsYebaJaSkpKAgB07ty5yvW2bNmCF198EfPmzcM///lPSCQSXL9+HQkJCdW+hlqtxoQJEzBz5kwsXLgQhw8fxltvvQUvLy+88cYbAHQHzmHDhiEvLw/vvfce7rnnHuzatQtTp06t0X4MHDgQO3bswLRp0/D888/j3nvvrXXP8e7duzF+/HhERETgww8/RLt27ZCcnIzffvutVtsxtm3bNhw5cgRvvPEGgoKCEBAQgKSkJPz3v//Ff/7zH5Ph7F9//TWUSqVhZEFpaSmGDh2K9PR0vPrqq+jRowcuXbqEN954AxcvXsTevXt5pZyIiExoNBpUVFSYLHNwcEBGRgb8/Pzw7rvvok2bNsjLy8OGDRswYMAAnD17Fl26dLHJ6+/YsQOnTp3C8uXL4e7ujvfffx+PPPIIrly5gg4dOgAAdu3ahUcffRT3338/vvnmG1RUVOCf//yn2UV/S/z9/XHPPffgs88+Q0BAAMaOHYsuXbrU+ng4Y8YMxMbGYubMmVi+fDmcnJxw5swZk3notfXcc8/hoYcewsaNG1FSUoJx48YhICAA69atw/Dhw03WXb9+Pfr06YMePXoAAA4cOIDRo0djwIAB+O9//wsvLy9s2bIFU6dORWlpqUlgT2QXIhFZtG7dOhGAePz4cVGtVotFRUXirl27xKCgIPH+++8X1Wq1yfpDhw4Vhw4darg/d+5c0dvbu8rXOHDggAhAPHDggGHZM888IwIQv/32W5N1x44dK3bp0sVw/z//+Y8IQPz1119N1nv++edFAOK6deuqfG2lUilOnDhRBCACEKVSqdi7d29xyZIlYk5OTpXP1evYsaPYsWNHsayszOo6zzzzjBgWFma2fOnSpWLlnyAAopeXl5iXl2ey/MKFCyIA8YsvvjBZfu+994p9+/Y13F+5cqUokUjEU6dOmaz3/fffiwDEnTt31mi/iIio5dMf5y3dKh/jRVEUKyoqxPLycrFTp07i3/72N8PypKQks+OufttJSUmGZZaOhwDEwMBAUaFQGJZlZWWJEolEXLlypWFZ//79RblcLqpUKsOyoqIi0c/Pz+xYasnJkyfFdu3aGfbPw8NDHDdunPjVV1+JWq222ucfPnxYBCAuWbKkyvUAiEuXLjVbHhYWJj7zzDOG+/r35+mnnzZbd8GCBaKLi4tYUFBgWJaQkCACEFetWmVY1rVrV7F3795m/1bjxo0Tg4ODRY1GU+1+ETUkDh0nqkZ0dDQcHR3h4eGB0aNHw8fHB9u3b4eDQ9UDQu69914UFBTgiSeewPbt23Hnzp0av6YgCBg/frzJsh49eiAlJcVw/9ChQ4Y2GXviiSdq9BoymQxbt25FQkICPvroIzz++OO4ffs2VqxYgYiICFy5cgWAbp50RUWFyQ0Arl69ihs3bmDmzJlwdnau8b5V54EHHjBJMgcA3bt3R9++fbFu3TrDssTERJw8eRLPPfecYdkvv/yCqKgo9OrVy6S9o0aNsmt2dyIiarq++uornDp1yuTm4OCAiooKvPPOO4iMjISTkxMcHBzg5OSEa9euITEx0WavP2zYMJNcKYGBgQgICDAc80tKSnD69GlMnDgRTk5OhvXc3d3NzhWs6d+/P65fv45du3bh1VdfxcCBA7Fv3z48/fTTmDBhgiEvi753X3/Tz9/+9ddfAQBz5syxyT7rTZo0yWzZc889h7KyMnzzzTeGZevWrYNMJsOTTz4JQDes//Lly5g2bRoAmLR57NixyMzMNJzHENkLA22iaugPwPv378fzzz+PxMTEGgWzTz31FL788kukpKRg0qRJCAgIwIABA7Bnz55qn+vq6moWvMpkMiiVSsP93NxcBAYGmj3X0rKqREREYP78+YiNjUVqaio+/PBD5Obm4vXXXwcAbNiwAY6OjiY3ALh9+zYAIDQ0tFavVx1rc8iee+45/PHHH7h8+TKAPw+6xv8W2dnZuHDhgll7PTw8IIpirS52EBFR6xAREYF+/fqZ3ABgwYIFeP311zFx4kT8/PPPOHHiBE6dOoWePXuirKzMZq/v5+dntkwmkxleIz8/H6Io1vuY7+joiFGjRmHFihXYvXs30tLSEBMTg19++cUQSA8fPtzk+Km/mH379m1IpVIEBQXVZRetsnTM79atG/r372+4uK7RaBAbG4uHH34Yvr6+AP7Mk7No0SKzY/6LL74IADzmk91xjjZRNfQHYEB31Vmj0WDNmjX4/vvvMXny5Cqf++yzz+LZZ59FSUkJDh8+jKVLl2LcuHG4evUqwsLC6tUuPz8/nDx50mx5VlZWnbcpCAL+9re/Yfny5YiPjwcAjB8/HqdOnTJbt02bNgB0Cdmq4uzsbJbEDbB+ALQ2Z+yJJ57AggULsH79eqxYsQIbN27ExIkTTXq//f394eLigi+//NLiNvz9/atsKxERkV5sbCyefvppvPPOOybL79y5A29v70Zrh4+PDwRBsDgfuz7HfD8/P8yfPx8HDx5EfHw8xo4di88//xxFRUWGdfTHzTZt2kCj0SArK8vqBXFAd4HA0jE/NzfX4vrWjvnPPvssXnzxRSQmJuLmzZvIzMw0qfSib9fixYvx6KOPWtyGrebQE9UVe7SJaun999+Hj48P3njjjRqXxHBzc8OYMWOwZMkSlJeX49KlS/Vux9ChQ1FUVGS4Cq23ZcuWGj0/MzPT4vKMjAwoFApDDVE/Pz+LV/o7d+6Mjh074ssvv7R4UNVr3749cnJyTE4QysvLsXv37hq1U8/HxwcTJ07EV199hV9++QVZWVkmw8YBYNy4cbhx44bFNvfr189i9nMiIiJLBEGATCYzWbZjxw7cunWrUdvh5uaGfv36Ydu2bSgvLzcsLy4utpidvDK1Wm010NUPgdcf87t06WLxuDlmzBgAuszfVWnfvr1JVnAA2L9/P4qLi6ttp7EnnngCzs7OWL9+PdavX4+2bdti5MiRhse7dOmCTp064fz58xaP9/369TMZjk9kD+zRJqolHx8fLF68GC+//DI2b96M6dOnW1xv9uzZcHFxweDBgxEcHIysrCysXLkSXl5e6N+/f73b8cwzz+Cjjz7C9OnT8fbbb+Oee+7Br7/+aghgqyvX8Ze//AUFBQWYNGkSoqKiIJVKcfnyZXz00UeQSCR45ZVXqm3Df/7zH4wfPx7R0dH429/+hnbt2iE1NRW7d+/Gpk2bAABTp07FG2+8gccffxx///vfoVQq8e9//xsajabW+/zcc8/hm2++wdy5cxEaGooRI0aYPD5//nz88MMPuP/++/G3v/0NPXr0gFarRWpqKn777TcsXLgQAwYMqPXrEhFR6zLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1200x400 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(1, 2, figsize=(12,4), sharex=True)\n",
|
||||
"\n",
|
||||
"ax[0].plot(x, y_true_rising, label=\"true\")\n",
|
||||
"ax[0].scatter(x, y_rising, s=12, label=\"noisy\")\n",
|
||||
"ax[0].set_title(\"Rising S-curve\")\n",
|
||||
"ax[0].legend()\n",
|
||||
"\n",
|
||||
"ax[1].plot(x, y_true_falling, label=\"true\")\n",
|
||||
"ax[1].scatter(x, y_falling, s=12, label=\"noisy\")\n",
|
||||
"ax[1].set_title(\"Falling S-curve\")\n",
|
||||
"ax[1].legend()\n",
|
||||
"\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "95c826d1-9223-4fb2-9ada-d9fc1acad2fb",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Fit the Rising S-curve"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"id": "37d2fda4-40b8-4ea2-8f2c-6c3f2384d158",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Paramter list of the Scurve:\n",
|
||||
"['p0', 'p1', 'mu', 'sigma', 'A', 'C']\n",
|
||||
"6 \n",
|
||||
"\n",
|
||||
"== Tuned fit settings ==\n",
|
||||
"max_calls : 500\n",
|
||||
"tolerance : 0.01\n",
|
||||
"compute_erros: True\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"== Results ==\n",
|
||||
"True rising params: [100. 0.25 60. 6. 120. 1. ]\n",
|
||||
"lmfit rising result: [ 99.467831 0.287017 60.078405 5.715461 117.257939 0.967234]\n",
|
||||
"minuit_grad (obj api ) rising result:\n",
|
||||
"{'par': array([100. , 0.269137, 60.050193, 5.729932, 117.777635,\n",
|
||||
" 0.984073]),\n",
|
||||
" 'par_err': array([1. , 0.054406, 0.747753, 0.909864, 7.887499, 0.174058]),\n",
|
||||
" 'chi2': array([7.351743])}\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"res_r_lmfit = fit_scurve(x, y_rising)\n",
|
||||
"\n",
|
||||
"model_r = RisingScurve()\n",
|
||||
"print(\"Paramter list of the Scurve:\")\n",
|
||||
"# print(model_r.GetParNames())\n",
|
||||
"print(model_r.par_names)\n",
|
||||
"print(model_r.n_par, '\\n')\n",
|
||||
"\n",
|
||||
"model_r.FixParameter('p0', 100) # Fixed p0 = 100 (optimizer does not touch the value)\n",
|
||||
"model_r.SetParameter('A', 100) # Start value A = 100\n",
|
||||
"\n",
|
||||
"print(\"== Tuned fit settings ==\")\n",
|
||||
"model_r.compute_errors = True\n",
|
||||
"model_r.max_calls = 500\n",
|
||||
"model_r.tolerance = 0.01\n",
|
||||
"print(f\"max_calls : {model_r.max_calls}\")\n",
|
||||
"print(f\"tolerance : {model_r.tolerance}\")\n",
|
||||
"print(f\"compute_erros: {model_r.compute_errors}\")\n",
|
||||
"print(\"\\n\")\n",
|
||||
"print(\"== Results ==\")\n",
|
||||
"res_r_munuit_obj = model_r.fit(x, y_rising, np.sqrt(y_rising))\n",
|
||||
"\n",
|
||||
"print(\"True rising params: \", p_true_rising)\n",
|
||||
"print(\"lmfit rising result: \", res_r_lmfit)\n",
|
||||
"print(\"minuit_grad (obj api ) rising result:\")\n",
|
||||
"pprint(res_r_munuit_obj, sort_dicts=False)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5ebf1328-954a-43ce-953b-c402fd7998f9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Fit the Falling S-curve"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "5ce22d2d-1adb-47d9-bf36-6cb535fa8b77",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"True falling params: [100. 0.25 60. 6. 120. -1. ]\n",
|
||||
"lmfit falling result: [101.73448 0.225786 60.251971 6.092915 118.476735 -0.998677]\n",
|
||||
"minuit_grad (obj api ) falling result\n",
|
||||
"{'par': array([101.741754, 0.225723, 60.251604, 6.092887, 118.473268,\n",
|
||||
" -0.998565]),\n",
|
||||
" 'chi2': array([1865.728797])}\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"res_f_lmfit = fit_scurve2(x, y_falling)\n",
|
||||
"\n",
|
||||
"model_f = FallingScurve()\n",
|
||||
"model_f.SetParameter(2, 20) # set start p2 = 20\n",
|
||||
"model_f.SetParameter(4, 90) # set start p4 = 90\n",
|
||||
"model_f.SetParLimits(4, 80, 130)\n",
|
||||
"model_f.max_calls = 150 # this is limits the number of calls done by the minimizer (increasing the `max_calls` may help if the fit fails) \n",
|
||||
"res_f_munuit_obj = model_f.fit(x, y_falling)\n",
|
||||
"\n",
|
||||
"print(\"True falling params: \", p_true_falling)\n",
|
||||
"print(\"lmfit falling result: \", res_f_lmfit)\n",
|
||||
"print(\"minuit_grad (obj api ) falling result\")\n",
|
||||
"pprint(res_f_munuit_obj, sort_dicts=False)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"id": "fe72503a-fc77-4b17-9b17-b98055e2fc87",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAAA9oAAAF0CAYAAADcnkHwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADV1UlEQVR4nOzdeZxN5R/A8c+562x39hVjLMmuiEIKJSJCREkIqSwltIhsJSU/8lO0iCGJUpKSfUk/KaFkJ8Y6YzD7cuduz++Pa64ZZqzDDL7v1+u+uOeee85z7tyZc77neZ7vV1NKKYQQQgghhBBCCFEkdMXdACGEEEIIIYQQ4mYigbYQQgghhBBCCFGEJNAWQgghhBBCCCGKkATaQgghhBBCCCFEEZJAWwghhBBCCCGEKEISaAshhBBCCCGEEEVIAm0hhBBCCCGEEKIISaAthBBCCCGEEEIUIQm0hRBCCCGEEEKIIiSBtrhpxcbGomma52EwGIiKiuKJJ55g3759563fpEkTmjRpcln7WLt2LZqmsXbt2qJp9GXatWsXTz/9NBUqVMDLy4vQ0FDq1KlD//79SUtLK5Y2Fbfhw4dTtmxZDAYDgYGBwPk/26ysLEaNGlVsPzchhBBX79zzfN7HkCFDLnk7cXFxaJpGbGzseduOi4vzLOvRowflypUrugO4DEop5s2bx3333Ud4eDheXl6UKVOGFi1aMH369GJpU3GLi4vjkUceITg4GE3TGDhwYIE/yw0bNjBq1ChSUlKKra3i1mQo7gYIca3NnDmTKlWqYLVa+d///sfYsWNZs2YNu3fvJigoyLPe1KlTL3vbderU4bfffqNatWpF2eRLsnXrVu69916qVq3KiBEjKFeuHKdOneLvv/9m3rx5DBkyBH9//+veruK0aNEixo4dy7Bhw2jZsiVmsxk4/2eblZXF6NGjAS775ooQQoiSJfc8n1epUqWKfD9vvvkmL730UpFv91IMHTqU9957j2effZZXXnkFi8XCoUOHWL16NYsWLaJ3797F0q7i9PLLL/P7778zY8YMIiMjiYqKIjIykt9++42KFSt61tuwYQOjR4+mR48enhvwQlwPEmiLm16NGjWoW7cu4A6qnE4nI0eO5Pvvv+eZZ57xrHclwbK/vz/169cvsrZejg8++ACdTsfatWuxWCye5R07duStt95CKVUs7cpLKYXVasXb2/u67G/79u0AvPjii4SHh3uWF8eNECGEENdH3vP8tZQ3eLuesrOz+eCDD+jWrRuffvppvtd69OiBy+UqlnadKzs7Gy8vLzRNuy772759O3fffTft2rXLt7y4rsuEOJcMHRe3nNyT8YkTJ/ItL2jo+LRp07jjjjvw8/PDYrFQpUoV3njjDc/rBQ0d79GjB35+fuzfv59WrVrh5+dHdHQ0gwcPJicnJ9/2jx49SseOHbFYLAQGBvLUU0+xadOm84Y9FeT06dP4+/vj5+dX4OuXeqL7/fffadOmDSEhIXh5eVGxYkUGDhyY73gKGio3atSo8/ahaRr9+/fn448/pmrVqpjNZqZPn054eDhPP/30edtISUnB29ubQYMGeZalpaUxZMgQypcvj8lkonTp0gwcOJDMzMwLHke5cuUYPnw4ABEREWiaxqhRo4D8P9u4uDjCwsIAGD16tGeYYY8ePQA4efIkffr0ITo6GrPZTFhYGPfeey8rV6684P6FEEKULPv37+eZZ56hUqVK+Pj4ULp0adq0acM///xzRdsr6HyYe9774osvqFq1Kj4+Ptxxxx38+OOP571/0aJF1KpVC7PZTIUKFZg8eXKB59JzZWZmkpOTQ1RUVIGv63SXdjnvcrmYMmUKd955J97e3gQGBlK/fn1++OGHfMeTe+7Mq1y5cp7zJJwdWr98+XJ69uxJWFgYPj4+zJ8/H03TWLVq1XnbmDZtGpqmsW3bNs+yP//8k0cffZTg4GC8vLyoXbs2X3/99QWPI/faa//+/fz888+e83hcXNx5Q8dHjRrFK6+8AkD58uU96+Zet61evZomTZoQEhKCt7c3ZcuWpUOHDmRlZV3SZyrEhUiPtrjlHDx4EIDbb7/9guvNmzePvn37MmDAACZMmIBOp2P//v3s3Lnzovuw2+08+uij9OrVi8GDB/PLL7/w1ltvERAQwIgRIwD3ibNp06YkJSXx3nvvcdttt7F06VI6d+58ScfRoEEDfvrpJ5566imee+457r777svuOV62bBlt2rShatWqTJw4kbJlyxIXF8fy5csvazt5ff/996xfv54RI0YQGRlJeHg4Bw8e5OOPP+ajjz7KN5z9q6++wmq1ekYWZGVl0bhxY44ePcobb7xBrVq12LFjByNGjOCff/5h5cqVhV6QLFy4kI8++ojPP/+cpUuXEhAQQJkyZc5bLyoqiqVLl/Lwww/Tq1cvz3C73OD76aefZsuWLYwdO5bbb7+dlJQUtmzZwunTp6/4MxFCCHHtOJ1OHA5HvmUGg4Hjx48TEhLCu+++S1hYGElJScyaNYt77rmHrVu3Urly5SLZ/08//cSmTZsYM2YMfn5+jB8/nvbt27Nnzx4qVKgAwNKlS3nssce4//77mT9/Pg6HgwkTJpx3078goaGh3HbbbUydOpXw8HBatWpF5cqVL7vnuEePHsyZM4devXoxZswYTCYTW7ZsyTcP/XL17NmTRx55hC+++ILMzExat25NeHg4M2fO5MEHH8y3bmxsLHXq1KFWrVoArFmzhocffph77rmHjz/+mICAAObNm0fnzp3JysrKF9jnlTttr3379lSsWJEJEyYA7vN7fHx8vnV79+5NUlISU6ZM4bvvvvPcrKhWrZpnjvd9993HjBkzCAwM5NixYyxduhSbzYaPj88Vfy5CAKCEuEnNnDlTAWrjxo3Kbrer9PR0tXTpUhUZGanuv/9+Zbfb863fuHFj1bhxY8/z/v37q8DAwAvuY82aNQpQa9as8Szr3r27AtTXX3+db91WrVqpypUre55/9NFHClA///xzvvWee+45BaiZM2decN9Wq1W1a9dOAQpQer1e1a5dWw0bNkwlJiZe8L25KlasqCpWrKiys7MLXad79+4qJibmvOUjR45U5/4JAVRAQIBKSkrKt3zbtm0KUJ9++mm+5Xfffbe66667PM/HjRundDqd2rRpU771FixYoAC1ZMmSCx5PbptOnjyZb/m5P9uTJ08qQI0cOfK8bfj5+amBAwdecD9CCCGKX+55vqDHued4pZRyOBzKZrOpSpUqqZdfftmz/ODBg+edd3O3ffDgQc+ygs6HgIqIiFBpaWmeZQkJCUqn06lx48Z5ltWrV09FR0ernJwcz7L09HQVEhJy3rm0IH/88YcqW7as5/gsFotq3bq1mj17tnK5XBd9/y+//KIANWzYsAuuV9i5MSYmRnXv3t3zPPfz6dat23nrDho0SHl7e6uUlBTPsp07dypATZkyxbOsSpUqqnbt2uf9rFq3bq2ioqKU0+m8YFtjYmLUI488km9ZQT/L999//7yfpVJnry3++uuvC+5HiCslQ8fFTa9+/foYjUYsFgsPP/wwQUFBLFq0CIPhwgM67r77blJSUnjyySdZtGgRp06duuR9appGmzZt8i2rVasWhw4d8jxft26dp015Pfnkk5e0D7PZzMKFC9m5cyeTJk3iiSee4OTJk4wdO5aqVauyZ88ewD1P2uFw5HsA7N27l3///ZdevXrh5eV1ycd2MQ888EC+JHMANWvW5K677mLmzJmeZbt27eKPP/6gZ8+enmU//vgjNWrU4M4778zX3hYtWly37O533303sbGxvP3222zcuBG73X7N9ymEEOLKzZ49m02bNuV7GAwGHA4H77zzDtWqVcNkMmEwGDCZTOzbt49du3YV2f6bNm2aL1dKREQE4eHhnnN+ZmYmf/75J+3atcNkMnnW8/PzO+9aoTD16tVj//79LF26lDfeeIMGDRqwatUqunXrxqOPPurJy5Lbu5/7yJ2//fPPPwPQr1+/IjnmXB06dDhvWc+ePcnOzmb+/PmeZTNnzsRsNtOlSxfAPax/9+7dPPXUUwD52tyqVSvi4+M91zHXyp133onJZKJPnz7MmjWLAwcOXNP9iVuPBNrippd7Al69ejXPPfccu3btuqRg9umnn2bGjBkcOnSIDh06EB4ezj333MOKFSsu+l4fH5/zglez2YzVavU8P336NBEREee9t6BlF1K1alUGDhzInDlzOHz4MBMnTuT06dO8+eabAMyaNQuj0ZjvAe65yECBw6uvRmFzyHr27Mlvv/3G7t27gbMn3bw/ixMnTrBt27bz2muxWFBKXdbNjis1f/58unfvzvTp02nQoAHBwcF069aNhISEa75vIYQQl69q1arUrVs33wNg0KBBvPnmm7Rr147Fixfz+++/s2nTJu644w6ys7OLbP8hISHnLTObzZ59JCcno5S66nO+0WikRYsWjB07lmXLlnHkyBGaNGnCjz/+6AmkH3zwwXznz9yb2SdPnkSv1xMZGXklh1iogs751atXp169ep6b606nkzlz5tC2bVuCg4OBs3lyhgwZct45v2/fvgDX/JxfsWJFVq5cSXh4OP369aNixYpUrFiRyZMnX9P9iluHzNEWN73cEzC47zo7nU6mT5/OggUL6Nix4wXf+8wzz/DMM8+QmZnJL7/8wsiRI2ndujV79+4lJibmqtoVEhLCH3/8cd7yqwnoNE3j5ZdfZsyYMZ4M3G3atGHTpk3nrZs7J/no0aMX3KaXl9d5Sdyg8BNgYXPGnnzySQYNGkRsbCxjx47liy++oF27dvl6v0NDQ/H29mbGjBkFbiM0NPSCbS0KoaGhfPDBB3zwwQccPnyYH374gddff53ExESWLl16zfcvhBCiaMyZM4du3brxzjvv5Ft+6tSp61rmKSgoCE3TCpyPfTXn/JCQEAYOHMjatWvZvn07rVq14pNPPiE9Pd2zTuLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1200x400 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(1, 2, figsize=(12,4))\n",
|
||||
"\n",
|
||||
"# rising\n",
|
||||
"ax[0].plot(x, y_true_rising, label=\"true\")\n",
|
||||
"ax[0].scatter(x, y_rising, s=12, label=\"data\")\n",
|
||||
"ax[0].scatter(x, scurve(x, res_r_lmfit[:6]), s=20, marker=\"x\", label=\"lmfit\")\n",
|
||||
"ax[0].plot(x, model_r(x, res_r_munuit_obj['par']), linewidth=1, color=\"green\", label=\"minuit\")\n",
|
||||
"ax[0].set_title(\"Rising S-curve fits\")\n",
|
||||
"ax[0].legend()\n",
|
||||
"\n",
|
||||
"# falling\n",
|
||||
"ax[1].plot(x, y_true_falling, label=\"true\")\n",
|
||||
"ax[1].scatter(x, y_falling, s=12, label=\"data\")\n",
|
||||
"ax[1].scatter(x, scurve2(x, res_f_lmfit[:6]), s=20, marker=\"x\", label=\"lmfit\")\n",
|
||||
"ax[1].plot(x, model_f(x, res_f_munuit_obj['par']), linewidth=1, color=\"green\", label=\"minuit\")\n",
|
||||
"ax[1].set_title(\"Falling S-curve fits\")\n",
|
||||
"ax[1].legend()\n",
|
||||
"\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "2c42ad91-fd59-47ba-b90d-8300f97eb27a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Quick error check"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"id": "3ba40e9a-0c64-4ef6-ba90-b1116f99b862",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Rising abs error lmfit : [0.532169 0.037017 0.078405 0.284539 2.742061 0.032766]\n",
|
||||
"Rising abs error minuit_grad: [0. 0.019137 0.050193 0.270068 2.222365 0.015927]\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Falling abs error lmfit : [1.73448 0.024214 0.251971 0.092915 1.523265 0.001323]\n",
|
||||
"Falling abs error minuit_grad: [1.741754 0.024277 0.251604 0.092887 1.526732 0.001435]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(\"Rising abs error lmfit : \", np.abs(res_r_lmfit[:6] - p_true_rising))\n",
|
||||
"print(\"Rising abs error minuit_grad: \", np.abs(res_r_munuit_obj['par'] - p_true_rising))\n",
|
||||
"print(\"\\n\")\n",
|
||||
"print(\"Falling abs error lmfit : \", np.abs(res_f_lmfit[:6] - p_true_falling))\n",
|
||||
"print(\"Falling abs error minuit_grad: \", np.abs(res_f_munuit_obj['par'] - p_true_falling))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c9bbabfd-e24e-43ec-8622-983b3757d7a9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Benchmark"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"id": "dc0b52b6-9fc3-453d-b6eb-e325a2b8342e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Rising lmfit : 0.334 ms\n",
|
||||
"Rising minuit2 : 0.239 ms\n",
|
||||
"\n",
|
||||
"Falling lmfit : 0.294 ms\n",
|
||||
"Falling minuit2 : 0.214 ms\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def bench(fn, n_repeats=200):\n",
|
||||
" # warmup\n",
|
||||
" for _ in range(3):\n",
|
||||
" fn()\n",
|
||||
"\n",
|
||||
" t0 = time.perf_counter()\n",
|
||||
" for _ in range(n_repeats):\n",
|
||||
" res = fn()\n",
|
||||
" t1 = time.perf_counter()\n",
|
||||
" return res, (t1 - t0) / n_repeats\n",
|
||||
"\n",
|
||||
"model_rising = RisingScurve()\n",
|
||||
"model_falling = FallingScurve()\n",
|
||||
"\n",
|
||||
"res_r_lmfit, t_r_lmfit = bench(lambda: fit_scurve(x, y_rising))\n",
|
||||
"res_r_minuit, t_r_minuit_grad = bench(lambda: model_rising.fit(x, y_rising))\n",
|
||||
"\n",
|
||||
"res_f_lmfit, t_f_lmfit = bench(lambda: fit_scurve2(x, y_falling))\n",
|
||||
"res_f_minuit, t_f_minuit_grad = bench(lambda: model_falling.fit(x, y_falling))\n",
|
||||
"\n",
|
||||
"print(f\"Rising lmfit : {1e3*t_r_lmfit:.3f} ms\")\n",
|
||||
"print(f\"Rising minuit2 : {1e3*t_r_minuit_grad:.3f} ms\")\n",
|
||||
"print()\n",
|
||||
"print(f\"Falling lmfit : {1e3*t_f_lmfit:.3f} ms\")\n",
|
||||
"print(f\"Falling minuit2 : {1e3*t_f_minuit_grad:.3f} ms\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e694c18a-a203-4132-a472-c93c4d3ce853",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.12.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
File diff suppressed because it is too large.
Load diff
File diff suppressed because it is too large.
Load diff
Reference in new issue
Block a user