diff --git a/.clang-format b/.clang-format index 0a9d5a9..c744508 100644 --- a/.clang-format +++ b/.clang-format @@ -1,6 +1,6 @@ --- Language: Cpp -Standard: c++23 +Standard: Latest BasedOnStyle: LLVM # Indentation @@ -31,8 +31,7 @@ AllowShortIfStatementsOnASingleLine: Never AllowShortLoopsOnASingleLine: false # Pointers and references -PointerAlignment: Right -DereferenceBindsToType: true +PointerAlignment: Left ReferenceAlignment: Pointer # Line breaks @@ -50,7 +49,7 @@ IncludeCategories: Priority: 3 # Constructor initialization -ConstructorInitializerAllOnOneLineIfFitsOnLine: false +PackConstructorInitializers: CurrentLine ConstructorInitializerIndentWidth: 4 # Penalty adjustments diff --git a/.clang-tidy b/.clang-tidy index e3aec7f..cfcf01e 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,5 +1,5 @@ --- -Checks: '-*,readability-*,performance-*,bugprone-*,modernize-*' +Checks: '-*,readability-*,performance-*,bugprone-*,modernize-*,-modernize-use-trailing-return-type' HeaderFilterRegex: 'include/(stochfit|levmar)/.*\.h' WarningsAsErrors: '' FormatStyle: file diff --git a/docs/index.html b/docs/index.html new file mode 100644 index 0000000..b410477 --- /dev/null +++ b/docs/index.html @@ -0,0 +1,416 @@ + + + + + +StochFit — Stochastic Reflectivity Fitting + + + + + +
+
+
+ C++20 + Electron + React + TypeScript + CMake + vcpkg + Windows · macOS · Linux +
+

StochFit

+

Stochastic methods for fitting specular X-ray and neutron reflectivity data — now a modern, cross-platform, open-source rewrite.

+
+ View on GitHub + Download latest release + Report a bug +
+
+
+ +
+ +
+

What is StochFit

+

StochFit fits specular X-ray and neutron reflectivity curves using simulated-annealing-based stochastic optimization, with both model-independent (electron density profile) and model-dependent (box model) fitting. It provides an interactive GUI for building a layer model, running the fit live, and inspecting the resulting reflectivity and density profile.

+
+ +
+

A ground-up modernization

+

The original StochFit (2008–2013) was a Windows-only .NET application. The current codebase is a full rewrite: a native C++20 computation core, a cross-platform Electron/React GUI, and a reproducible CMake + vcpkg build — while preserving the same underlying stochastic fitting algorithms described in the original paper.

+
+
+

Rebuilt computation core

+

Simulated annealing (Greedy / SA / STUN policies), unified Parratt recursion, and electron density profile generation reimplemented in modern, OpenMP-parallel C++20 with optional GPU acceleration hooks.

+
+
+

Cross-platform GUI

+

Electron + React + TypeScript replaces the old WinForms UI, driven by a FlatBuffers-serialized FFI bridge into the native library — runs on Windows, macOS, and Linux.

+
+
+

Standardized build

+

CMake 3.21+ with named presets and vcpkg-managed dependencies (LAPACK, levmar, FlatBuffers) replace hand-maintained Visual Studio project files and vendored DLLs.

+
+
+

Real test coverage

+

GTest-based regression tests for the reflectivity math and fitting harness, plus a standalone console utility (mirefl) for manual spot-checks — the original shipped with none.

+
+
+
+ +
+

Notable fixes since the last public release

+ +

Looking for the legacy .NET releases (1.6.5 / 1.7.0)? They're still archived on SourceForge, but are unmaintained — the current, actively developed releases are on GitHub.

+
+ +
+

Architecture

+
Electron GUI (React / TypeScript) + ↕ IPC (FlatBuffers) +Koffi FFI (Node.js ↔ stochfit.dll / .so / .dylib) + ↕ C-style exports +stochfit_core (static library, C++20) + ├─ StochFitHarness — top-level orchestrator, worker thread + ├─ Anneal<Policy> — Greedy / SimulatedAnnealing / STUN + ├─ ParrattReflectivity — recursive reflectivity, OpenMP + Q-smearing + ├─ CEDP — electron density profile generation + └─ levmardll — Levenberg-Marquardt post-refinement
+
+ +
+

Building from source

+

Builds use CMake 3.21+ with named presets; dependencies are managed by vcpkg and bootstrapped automatically on first configure.

+ 
# Windows
+cmake --preset windows
+cmake --build --preset windows --target stochfit_shared
+
+# macOS / Linux
+cmake --preset default
+cmake --build --preset default --target stochfit_shared
+
+# GUI (Electron)
+cd gui
+npm install
+npm start        # dev mode with hot reload
+npm run make      # package a distributable installer
+ + + + + +
PlatformRequirements
WindowsVisual Studio 2022 (C++ workload), CMake 3.21+, Node.js 24+
macOSXcode 16+ command-line tools, macOS 14+ (required for std::jthread)
LinuxCMake 3.21+, a C++20 toolchain, Node.js 24+
+

Full build instructions, presets, and architecture notes are in the README.

+
+ +
+

Citing StochFit

+
+ S. M. Danauskas, D. Li, M. Meron, B. Lin and K. Y. C. Lee
+ Stochastic fitting of specular X-ray reflectivity data using StochFit.
+ J. Appl. Cryst. (2008). 41, 1187–1193. +
+
+ +
+

Get involved

+

Issues, feature requests, and pull requests are welcome on GitHub — that's now the home for all development, downloads, and bug tracking (this SourceForge page previously pointed to its own forums, which are retired).

+
+ github.com/diehard2/stochfit + Issues & bug reports + Discussions +
+
+ +
+ + + + + diff --git a/include/levmar/BoxLayerBuild.h b/include/levmar/BoxLayerBuild.h index 52e660e..19b7763 100644 --- a/include/levmar/BoxLayerBuild.h +++ b/include/levmar/BoxLayerBuild.h @@ -1,20 +1,22 @@ #pragma once -#include "Settings.h" -#include "stochfit/LayerStack.h" #include #include #include +#include "Settings.h" +#include "stochfit/LayerStack.h" + // Compact representation of a box-model layer stack. // All rho values are pre-multiplied by 2 to match the m_DEDP convention used // by ParrattReflectivity (rho[i] = 2 * real_SLD_i in reduced units). // length_mult[i] = {0, -2*length[i]} — precomputed for the Parratt ak term. // sigma_sq[i] = -2*sigma[i]^2 — precomputed for the Nevot-Croce factor. -struct BoxLayers { - std::vector> rho; // size = boxes + 2 - std::vector> length_mult; // size = boxes + 2 - std::vector sigma_sq; // size = boxes + 1 (per interface) +struct BoxLayers +{ + std::vector> rho; // size = boxes + 2 + std::vector> length_mult; // size = boxes + 2 + std::vector sigma_sq; // size = boxes + 1 (per interface) double normfactor = 1.0; // Returns a LayerStack view over this object. @@ -27,7 +29,4 @@ struct BoxLayers { // one_sigma=true: single global roughness p[0]; stride is 2 per box (length, rho_frac). // one_sigma=false: per-interface roughness; stride is 3 per box (length, rho_frac, sigma). // In both cases the last parameter p.back() is the normalization factor. -void BuildBoxLayers(const BoxReflSettings& rs, - std::span p, - bool one_sigma, - BoxLayers& out); +void BuildBoxLayers(const BoxReflSettings& rs, std::span p, bool one_sigma, BoxLayers& out); diff --git a/include/levmar/LevMardll.h b/include/levmar/LevMardll.h index 44c95f4..1d9ca92 100644 --- a/include/levmar/LevMardll.h +++ b/include/levmar/LevMardll.h @@ -24,11 +24,12 @@ // Output: result table — returns bytes written into outBuf. #pragma once -#include "platform.h" #include -extern "C" EXPORT int FastReflfit (const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); +#include "platform.h" + +extern "C" EXPORT int FastReflfit(const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); extern "C" EXPORT int FastReflGenerate(const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); -extern "C" EXPORT int Rhofit (const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); -extern "C" EXPORT int RhoGenerate (const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); +extern "C" EXPORT int Rhofit(const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); +extern "C" EXPORT int RhoGenerate(const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); extern "C" EXPORT int StochFitBoxModel(const uint8_t* inBuf, int inLen, uint8_t* outBuf, int maxLen); diff --git a/include/levmar/RhoCalc.h b/include/levmar/RhoCalc.h index 93cd61c..06f6fab 100644 --- a/include/levmar/RhoCalc.h +++ b/include/levmar/RhoCalc.h @@ -1,58 +1,58 @@ -/* +/* * Copyright (C) 2008 Stephen Danauskas - * + * * This Program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2, or (at your option) * any later version. - * + * * This Program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. - * + * * You should have received a copy of the GNU General Public License * along with GNU Make; see the file COPYING. If not, write to - * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. * http://www.gnu.org/copyleft/gpl.html * */ #pragma once -#include "Settings.h" #include +#include "Settings.h" + class RhoCalc { -private: - vector distarray; - vector rhoarray; - vector rougharray; - vector MIRho; - vector ZIncrement; - - vector m_LengthArray; - vector m_RhoArray; - vector m_SigmaArray; - - bool onesigma; - double SubSLD; - double m_dSupSLD; - int boxnumber; - int Zlength; - - - void Rhocalculate(double SubRough,double Zoffset); - -public: - //Member functions - ~RhoCalc(); + private: + vector distarray; + vector rhoarray; + vector rougharray; + vector MIRho; + vector ZIncrement; + + vector m_LengthArray; + vector m_RhoArray; + vector m_SigmaArray; + + bool onesigma; + double SubSLD; + double m_dSupSLD; + int boxnumber; + int Zlength; + + void Rhocalculate(double SubRough, double Zoffset); + + public: + // Member functions + ~RhoCalc(); void init(const BoxReflSettings& InitStruct); - void mkdensityboxmodel(std::span p); + void mkdensityboxmodel(std::span p); void mkdensity(std::span p); - static void objective(double *p, double *x, int m, int n, void *data); + static void objective(double* p, double* x, int m, int n, void* data); - vector nk; - vector nkb; + vector nk; + vector nkb; }; \ No newline at end of file diff --git a/include/levmar/Settings.h b/include/levmar/Settings.h index e7efda6..10066cf 100644 --- a/include/levmar/Settings.h +++ b/include/levmar/Settings.h @@ -3,7 +3,7 @@ struct BoxReflSettings { - std::string Directory; + std::string Directory; vector Q; vector Refl; vector ReflError; @@ -11,20 +11,20 @@ struct BoxReflSettings vector UL; vector LL; vector ParamPercs; - int QPoints = 0; - bool OneSigma = false; - double SubSLD = 0; - double SupSLD = 0; - int Boxes = 0; + int QPoints = 0; + bool OneSigma = false; + double SubSLD = 0; + double SupSLD = 0; + int Boxes = 0; double Wavelength = 0; double QSpread = 0; - bool ImpNorm = false; - int FitFunc = 0; + bool ImpNorm = false; + int FitFunc = 0; int LowQOffset = 0; int HighQOffset = 0; int Iterations = 0; // EDP Specific Settings vector MIEDP; vector ZIncrement; - int ZLength = 0; + int ZLength = 0; }; diff --git a/include/platform.h b/include/platform.h index 187e9f2..729285a 100644 --- a/include/platform.h +++ b/include/platform.h @@ -10,29 +10,30 @@ // are either stubbed out or replaced with C++23 equivalents. // ── Standard C/C++ includes ──────────────────────────────────────────────── +#include +#include +#include +#include +#include +#include #include #include #include -#include #include -#include -#include -#include -#include #include -#include #include #include +#include #include -#include -#include -#include -#include #include -#include +#include #include -#include +#include +#include +#include #include +#include + #include // std::jthread is NOT used — replaced with std::thread + std::atomic @@ -43,7 +44,7 @@ // ── OMP thread limit ──────────────────────────────────────────────────────── #ifndef MAX_OMP_THREADS -# define MAX_OMP_THREADS 8 + #define MAX_OMP_THREADS 8 #endif // ── Debug flags ────────────────────────────────────────────────────────────── @@ -51,14 +52,16 @@ inline constexpr bool kSingleProcDebug = false; // ── Export macro ──────────────────────────────────────────────────────────── #if defined(_MSC_VER) -# define EXPORT __declspec(dllexport) + #define EXPORT __declspec(dllexport) #else -# define EXPORT __attribute__((visibility("default"))) + #define EXPORT __attribute__((visibility("default"))) #endif // ── MessageBox replacement ────────────────────────────────────────────────── -inline void platform_error(const char* msg) { std::cerr << msg << std::endl; } +inline void platform_error(const char* msg) +{ + std::cerr << msg << std::endl; +} // ── Convenience namespaces (matching existing codebase convention) ────────── using namespace std; - diff --git a/include/stochfit/Algorithm.h b/include/stochfit/Algorithm.h index eec81fc..74fd841 100644 --- a/include/stochfit/Algorithm.h +++ b/include/stochfit/Algorithm.h @@ -1,20 +1,26 @@ #pragma once enum class SaAlgorithm : int { - Greedy = 0, + Greedy = 0, Simulated = 1, - Stun = 2, + Stun = 2, }; -inline constexpr SaAlgorithm AlgorithmFromInt(int v) { +inline constexpr SaAlgorithm AlgorithmFromInt(int v) +{ switch (v) { - case 0: return SaAlgorithm::Greedy; - case 1: return SaAlgorithm::Simulated; - case 2: return SaAlgorithm::Stun; - default: return SaAlgorithm::Greedy; + case 0: + return SaAlgorithm::Greedy; + case 1: + return SaAlgorithm::Simulated; + case 2: + return SaAlgorithm::Stun; + default: + return SaAlgorithm::Greedy; } } -inline constexpr int AlgorithmToInt(SaAlgorithm a) { +inline constexpr int AlgorithmToInt(SaAlgorithm a) +{ return static_cast(a); } diff --git a/include/stochfit/Anneal.h b/include/stochfit/Anneal.h index 46fe8a4..3a66390 100644 --- a/include/stochfit/Anneal.h +++ b/include/stochfit/Anneal.h @@ -1,37 +1,39 @@ #pragma once +#include +#include +#include +#include + #include "AnnealPolicies.h" #include "CEDP.h" #include "ParameterStepper.h" #include "ReflectivityObjective.h" #include "UnifiedReflectivity.h" -#include -#include -#include -#include - -struct AnnealDeps { +struct AnnealDeps +{ std::span yi, eyi; - std::span reflBuf; + std::span reflBuf; bool impNorm = false; }; -template -class Anneal { -public: +template class Anneal +{ + public: template - Anneal(CEDP& edp, - ParrattReflectivity& parratt, - const ReflectivityObjective& objective, - ParameterStepper& stepper, - const ParamVector& initParams, - AnnealDeps deps, - PolicyArgs&&... policyArgs) - : m_edp(&edp), m_parratt(&parratt), m_objective(&objective), - m_stepper(&stepper), m_tempParams(initParams), m_deps(std::move(deps)), + Anneal(CEDP& edp, ParrattReflectivity& parratt, const ReflectivityObjective& objective, ParameterStepper& stepper, + const ParamVector& initParams, AnnealDeps deps, PolicyArgs&&... policyArgs) + : m_edp(&edp), + m_parratt(&parratt), + m_objective(&objective), + m_stepper(&stepper), + m_tempParams(initParams), + m_deps(std::move(deps)), m_rng(std::random_device{}()), - m_policy(std::forward(policyArgs)...) {} + m_policy(std::forward(policyArgs)...) + { + } void InitEnergy(ParamVector& params); @@ -44,50 +46,75 @@ class Anneal { void ComputeSharedRefl(); bool EvaluateAndAccept(ParamVector& params); - double GetTemperature() const { return m_policy.GetTemperature(); } - void SetTemperature(double t) { m_policy.SetTemperature(t); } - double GetRawTemperature() const { return m_policy.GetRawTemperature(); } - double GetAverageFSTUN() const { return m_policy.GetAverageFSTUN(); } - void SetAverageFSTUN(double f) { m_policy.SetAverageFSTUN(f); } - double GetLowestEnergy() const { return m_bestEnergy; } - double GetCurrentEnergy() const { return m_currentEnergy; } - double GetLastChiSquare() const { return m_lastChiSquare; } - -private: - CEDP* m_edp; - ParrattReflectivity* m_parratt; + double GetTemperature() const + { + return m_policy.GetTemperature(); + } + void SetTemperature(double t) + { + m_policy.SetTemperature(t); + } + double GetRawTemperature() const + { + return m_policy.GetRawTemperature(); + } + double GetAverageFSTUN() const + { + return m_policy.GetAverageFSTUN(); + } + void SetAverageFSTUN(double f) + { + m_policy.SetAverageFSTUN(f); + } + double GetLowestEnergy() const + { + return m_bestEnergy; + } + double GetCurrentEnergy() const + { + return m_currentEnergy; + } + double GetLastChiSquare() const + { + return m_lastChiSquare; + } + + private: + CEDP* m_edp; + ParrattReflectivity* m_parratt; const ReflectivityObjective* m_objective; - ParameterStepper* m_stepper; - ParamVector m_tempParams; - AnnealDeps m_deps; - std::mt19937 m_rng; - double m_bestEnergy = std::numeric_limits::max(); + ParameterStepper* m_stepper; + ParamVector m_tempParams; + AnnealDeps m_deps; + std::mt19937 m_rng; + double m_bestEnergy = std::numeric_limits::max(); double m_currentEnergy = std::numeric_limits::max(); double m_lastChiSquare = 0.0; [[no_unique_address]] Policy m_policy; - void ComputeModel(ParamVector& p) { + void ComputeModel(ParamVector& p) + { auto result = m_parratt->CalculateReflectivity(*m_edp); std::ranges::copy(result, m_deps.reflBuf.begin()); if (m_deps.impNorm) { - for (auto& v : m_deps.reflBuf) v *= p.GetImpNorm(); + for (auto& v : m_deps.reflBuf) + v *= p.GetImpNorm(); } } }; -template -void Anneal::InitEnergy(ParamVector& params) { +template void Anneal::InitEnergy(ParamVector& params) +{ m_edp->GenerateEDP(params); ComputeModel(params); - m_bestEnergy = m_currentEnergy = - m_objective->Evaluate(m_deps.reflBuf, m_deps.yi, m_deps.eyi); + m_bestEnergy = m_currentEnergy = m_objective->Evaluate(m_deps.reflBuf, m_deps.yi, m_deps.eyi); m_lastChiSquare = ComputeChiSquare(m_deps.reflBuf, m_deps.yi, m_deps.eyi); } // ── Cooperative methods for persistent OMP parallel regions ────────────────── -template -void Anneal::PrepareCandidate(ParamVector& params) { +template void Anneal::PrepareCandidate(ParamVector& params) +{ #pragma omp single { m_tempParams = params; @@ -100,20 +127,21 @@ void Anneal::PrepareCandidate(ParamVector& params) { // implicit barrier: all threads see completed EDP } -template -void Anneal::ComputeSharedRefl() { +template void Anneal::ComputeSharedRefl() +{ auto result = m_parratt->CalculateReflectivityCooperative(*m_edp); #pragma omp single { std::ranges::copy(result, m_deps.reflBuf.begin()); if (m_deps.impNorm) { - for (auto& v : m_deps.reflBuf) v *= m_tempParams.GetImpNorm(); + for (auto& v : m_deps.reflBuf) + v *= m_tempParams.GetImpNorm(); } } } -template -bool Anneal::EvaluateAndAccept(ParamVector& params) { +template bool Anneal::EvaluateAndAccept(ParamVector& params) +{ const double candE = m_objective->Evaluate(m_deps.reflBuf, m_deps.yi, m_deps.eyi); if (candE < m_bestEnergy) { @@ -132,4 +160,3 @@ bool Anneal::EvaluateAndAccept(ParamVector& params) { return false; } - diff --git a/include/stochfit/AnnealPolicies.h b/include/stochfit/AnnealPolicies.h index d3a99d3..3e68fdd 100644 --- a/include/stochfit/AnnealPolicies.h +++ b/include/stochfit/AnnealPolicies.h @@ -6,71 +6,122 @@ // ── Greedy ──────────────────────────────────────────────────────────────────── // Accepts any move that lowers the energy. No temperature schedule. -struct GreedyPolicy { - bool Accept(double curE, double candE, double /*bestE*/, std::mt19937& /*rng*/) const { +struct GreedyPolicy +{ + bool Accept(double curE, double candE, double /*bestE*/, std::mt19937& /*rng*/) const + { return candE < curE; } - double GetTemperature() const { return 0.0; } - void SetTemperature(double) {} - double GetRawTemperature() const { return 0.0; } - double GetAverageFSTUN() const { return 0.0; } - void SetAverageFSTUN(double) {} + double GetTemperature() const + { + return 0.0; + } + void SetTemperature(double) + { + } + double GetRawTemperature() const + { + return 0.0; + } + double GetAverageFSTUN() const + { + return 0.0; + } + void SetAverageFSTUN(double) + { + } }; // ── Simulated Annealing ─────────────────────────────────────────────────────── // Standard Metropolis acceptance with geometric cooling. // m_dTemp stores β = 1/T. ProbCalc = exp(-β·ΔE)·100. -struct SimulatedPolicy { +struct SimulatedPolicy +{ SimulatedPolicy(double initTemp, double slope, int platIter); bool Accept(double curE, double candE, double /*bestE*/, std::mt19937& rng); - double GetTemperature() const { return 1.0 / m_dTemp; } - void SetTemperature(double t) { m_dTemp = t; } // t is raw β (from session) - double GetRawTemperature() const { return m_dTemp; } // β, for session save - double GetAverageFSTUN() const { return 0.0; } - void SetAverageFSTUN(double) {} + double GetTemperature() const + { + return 1.0 / m_dTemp; + } + void SetTemperature(double t) + { + m_dTemp = t; + } // t is raw β (from session) + double GetRawTemperature() const + { + return m_dTemp; + } // β, for session save + double GetAverageFSTUN() const + { + return 0.0; + } + void SetAverageFSTUN(double) + { + } -private: + private: double m_dTemp; double m_slope; - int m_platIter; - int m_iter = 0; + int m_platIter; + int m_iter = 0; - void Schedule(); - double ProbCalc(double deltaE) const { return std::exp(-m_dTemp * deltaE) * 100.0; } + void Schedule(); + double ProbCalc(double deltaE) const + { + return std::exp(-m_dTemp * deltaE) * 100.0; + } }; // ── STUN (Stochastic Tunneling) ─────────────────────────────────────────────── // Transforms the energy landscape via fSTUN to tunnel through barriers. // Supports adaptive temperature control via a sliding window average. // m_dTemp stores β = 1/T, same convention as SimulatedPolicy. -struct StunPolicy { - StunPolicy(double initTemp, double slope, int platIter, - double gamma, int stunFunc, int tempIter, bool adaptive); +struct StunPolicy +{ + StunPolicy(double initTemp, double slope, int platIter, double gamma, int stunFunc, int tempIter, bool adaptive); bool Accept(double curE, double candE, double bestE, std::mt19937& rng); - double GetTemperature() const { return 1.0 / m_dTemp; } - void SetTemperature(double t) { m_dTemp = t; } - double GetRawTemperature() const { return m_dTemp; } - double GetAverageFSTUN() const { return m_averageFSTUN; } - void SetAverageFSTUN(double f) { m_averageFSTUN = f; } + double GetTemperature() const + { + return 1.0 / m_dTemp; + } + void SetTemperature(double t) + { + m_dTemp = t; + } + double GetRawTemperature() const + { + return m_dTemp; + } + double GetAverageFSTUN() const + { + return m_averageFSTUN; + } + void SetAverageFSTUN(double f) + { + m_averageFSTUN = f; + } -private: + private: double m_dTemp; double m_slope; - int m_platIter; - int m_iter = 0; + int m_platIter; + int m_iter = 0; double m_gamma; double m_averageFSTUN; - int m_stunFunc; - int m_tempIter; - bool m_adaptive; + int m_stunFunc; + int m_tempIter; + bool m_adaptive; std::deque m_qWindow; double fSTUN(double val, double bestE) const; - void AdjustTemp(double averageSTUN); - void Schedule(); - double ProbCalc(double deltaE) const { return std::exp(-m_dTemp * deltaE) * 100.0; } + void AdjustTemp(double averageSTUN); + void Schedule(); + double ProbCalc(double deltaE) const + { + return std::exp(-m_dTemp * deltaE) * 100.0; + } }; diff --git a/include/stochfit/CEDP.h b/include/stochfit/CEDP.h index a72aefd..34006b1 100644 --- a/include/stochfit/CEDP.h +++ b/include/stochfit/CEDP.h @@ -9,75 +9,79 @@ #include "LayerStack.h" #include "ParamVector.h" -class CEDP { -private: - vector m_fDistArray; - vector m_fRhoArray; - vector m_fImagRhoArray; - vector m_fEDSpacingArray; - - // Precomputed per-layer length multiplier {0, -2*dz} for BuildLayerStack. - vector> m_length_mult; - - double m_dRho; - double m_dLambda; - double m_dDz0; - double m_dBeta; - double m_dBeta_Sup; - double m_dBeta_Sub; - double m_dWaveConstant; - - int m_iLayers; - - bool m_bUseSurfAbs; - - // Cached flat-region offsets; updated at end of GenerateEDP. - mutable int m_supOff = 0; - mutable int m_subOff = 0; - - template - void BuildEDP(ParamVector &g); - template - void FillBoxArraysImpl(ParamVector &g); - -public: - CEDP() = default; - explicit CEDP(const ReflSettings &s) { Init(s); } - - void Init(const ReflSettings &InitStruct); - int GetLayerCount() const { return m_iLayers; } - - // Standalone EDP build — includes serial setup, parallel EDP loop, and offset caching. - void GenerateEDP(ParamVector &g); - - // Cooperative path (SA persistent-team use): - // 1. Call FillBoxArrays(g) inside your omp single block before the parallel section. - // 2. Call GenerateEDPCooperative(g) — runs only the parallel omp for. - // 3. Call BuildLayerStackFull() inside the next omp single — computes offsets fresh. - void FillBoxArrays(ParamVector &g); // serial setup only, no omp pragma - void GenerateEDPCooperative(ParamVector &g); // omp for only, no serial omp singles - - int Get_EDPPointCount() const; - bool Get_UseABS() const; - double Get_FilmAbs() const; - double Get_FilmAbsInput() - const; // returns m_dBeta / m_dWaveConstant — inverse of Set_FilmAbs - double Get_Dz() const; - double Get_LeftOffset() const; // superphase padding = 6 * RoughnessMax - double Get_WaveConstant() const; - void Set_FilmAbs(double absorption); - std::pair GetOffSets() const; - - // Returns a zero-copy view of m_DEDP suitable for ParrattReflectivity. - // Valid only after a GenerateEDP() call. No roughness (σ=0), sup/sub offsets - // from the most recent GenerateEDP, transparent flag from m_bUseSurfAbs. - LayerStack BuildLayerStack() const; - - // Like BuildLayerStack but computes sup/sub offsets fresh from m_DEDP — use - // in the cooperative SA path where GenerateEDPCooperative skipped GetOffSets. - // Const (m_supOff/m_subOff are mutable); call from within an omp single. - LayerStack BuildLayerStackFull() const; - - vector> m_EDP; - vector> m_DEDP; +class CEDP +{ + private: + vector m_fDistArray; + vector m_fRhoArray; + vector m_fImagRhoArray; + vector m_fEDSpacingArray; + + // Precomputed per-layer length multiplier {0, -2*dz} for BuildLayerStack. + vector> m_length_mult; + + double m_dRho; + double m_dLambda; + double m_dDz0; + double m_dBeta; + double m_dBeta_Sup; + double m_dBeta_Sub; + double m_dWaveConstant; + + int m_iLayers; + + bool m_bUseSurfAbs; + + // Cached flat-region offsets; updated at end of GenerateEDP. + mutable int m_supOff = 0; + mutable int m_subOff = 0; + + template void BuildEDP(ParamVector& g); + template void FillBoxArraysImpl(ParamVector& g); + + public: + CEDP() = default; + explicit CEDP(const ReflSettings& s) + { + Init(s); + } + + void Init(const ReflSettings& InitStruct); + int GetLayerCount() const + { + return m_iLayers; + } + + // Standalone EDP build — includes serial setup, parallel EDP loop, and offset caching. + void GenerateEDP(ParamVector& g); + + // Cooperative path (SA persistent-team use): + // 1. Call FillBoxArrays(g) inside your omp single block before the parallel section. + // 2. Call GenerateEDPCooperative(g) — runs only the parallel omp for. + // 3. Call BuildLayerStackFull() inside the next omp single — computes offsets fresh. + void FillBoxArrays(ParamVector& g); // serial setup only, no omp pragma + void GenerateEDPCooperative(ParamVector& g); // omp for only, no serial omp singles + + int Get_EDPPointCount() const; + bool Get_UseABS() const; + double Get_FilmAbs() const; + double Get_FilmAbsInput() const; // returns m_dBeta / m_dWaveConstant — inverse of Set_FilmAbs + double Get_Dz() const; + double Get_LeftOffset() const; // superphase padding = 6 * RoughnessMax + double Get_WaveConstant() const; + void Set_FilmAbs(double absorption); + std::pair GetOffSets() const; + + // Returns a zero-copy view of m_DEDP suitable for ParrattReflectivity. + // Valid only after a GenerateEDP() call. No roughness (σ=0), sup/sub offsets + // from the most recent GenerateEDP, transparent flag from m_bUseSurfAbs. + LayerStack BuildLayerStack() const; + + // Like BuildLayerStack but computes sup/sub offsets fresh from m_DEDP — use + // in the cooperative SA path where GenerateEDPCooperative skipped GetOffSets. + // Const (m_supOff/m_subOff are mutable); call from within an omp single. + LayerStack BuildLayerStackFull() const; + + vector> m_EDP; + vector> m_DEDP; }; diff --git a/include/stochfit/LayerStack.h b/include/stochfit/LayerStack.h index 6d18f3a..3a77f5f 100644 --- a/include/stochfit/LayerStack.h +++ b/include/stochfit/LayerStack.h @@ -13,12 +13,13 @@ // Parratt loop pays only one std::exp() call per interface without extra ops. // sigma_sq[i] = -2·sigma[i]² — precomputed for the same reason; empty span when // has_roughness is false (Nevot-Croce term is skipped at compile time). -struct LayerStack { +struct LayerStack +{ std::span> rho; // 2·EDP density per interface std::span> length_mult; // {0, -2·length[i]} - std::span sigma_sq; // -2·σ[i]²; empty when !has_roughness - int sup_offset = 0; - int sub_offset = 0; - bool transparent = true; // imag(rho) all zero → real-only fast path eligible - bool has_roughness = false; // any σ != 0 → Nevot-Croce multiplied into Fresnel rj + std::span sigma_sq; // -2·σ[i]²; empty when !has_roughness + int sup_offset = 0; + int sub_offset = 0; + bool transparent = true; // imag(rho) all zero → real-only fast path eligible + bool has_roughness = false; // any σ != 0 → Nevot-Croce multiplied into Fresnel rj }; diff --git a/include/stochfit/ParamVector.h b/include/stochfit/ParamVector.h index 7de6b8d..c2835e5 100644 --- a/include/stochfit/ParamVector.h +++ b/include/stochfit/ParamVector.h @@ -27,64 +27,89 @@ // on all mutations. All public API uses double. #include -#include "platform.h" + #include "SettingsStruct.h" +#include "platform.h" class ParamVector { -public: + public: explicit ParamVector(const ReflSettings&); // EDP grid (supphase, box1..boxN, subphase) — used by CEDP. - int RealParamsSize() const { return m_boxes + 2; } - double GetRealParams(int i) const; - std::span RealParams() const; + int RealParamsSize() const + { + return m_boxes + 2; + } + double GetRealParams(int i) const; + std::span RealParams() const; // Mutable SA search space (per-box SLDs + optional roughness/surfabs/impnorm). - int ParamCount() const { return m_paramCount; } - int BoxCount() const { return m_boxes; } + int ParamCount() const + { + return m_paramCount; + } + int BoxCount() const + { + return m_boxes; + } double GetMutatableParameter(int i) const; - void SetMutatableParameter(int i, double val); // clamps to bounds + void SetMutatableParameter(int i, double val); // clamps to bounds // Named parameter accessors. double GetRoughness() const; - void SetRoughness(double rough); // no-op if roughness is fixed - double GetImpNorm() const; - void SetImpNorm(double norm); // no-op if impnorm disabled - double GetSurfAbs() const; - void SetSurfAbs(double surfabs); // no-op if surfabs disabled + void SetRoughness(double rough); // no-op if roughness is fixed + double GetImpNorm() const; + void SetImpNorm(double norm); // no-op if impnorm disabled + double GetSurfAbs() const; + void SetSurfAbs(double surfabs); // no-op if surfabs disabled // Phase-endpoint values (not part of the SA mutable space). - void SetSubphase(double subval) { m_edpValues[m_boxes + 1] = subval; } - void SetSupphase(double supval) { m_edpValues[0] = supval; } + void SetSubphase(double subval) + { + m_edpValues[m_boxes + 1] = subval; + } + void SetSupphase(double supval) + { + m_edpValues[0] = supval; + } // Bounds. double GetUpperBounds(int index) const; double GetLowerBounds(int index) const; - void UpdateBoundaries(); + void UpdateBoundaries(); // Feature flags. - bool IsRoughnessFixed() const { return m_fixRoughness; } - bool IsImpNormFixed() const { return m_fixImpNorm; } - bool UsesSurfAbs() const { return m_useSurfAbs; } + bool IsRoughnessFixed() const + { + return m_fixRoughness; + } + bool IsImpNormFixed() const + { + return m_fixImpNorm; + } + bool UsesSurfAbs() const + { + return m_useSurfAbs; + } -private: - std::vector m_edpValues; // [supphase, box1..boxN, subphase] + private: + std::vector m_edpValues; // [supphase, box1..boxN, subphase] std::vector m_mutableParams; std::vector m_high; std::vector m_low; - int m_boxes = 0; - int m_paramCount = 0; - int m_roughnessIdx = -1; - int m_surfAbsIdx = -1; - int m_impNormIdx = -1; - double m_roughness = 0.0; // fixed roughness value (when m_fixRoughness) + int m_boxes = 0; + int m_paramCount = 0; + int m_roughnessIdx = -1; + int m_surfAbsIdx = -1; + int m_impNormIdx = -1; + double m_roughness = 0.0; // fixed roughness value (when m_fixRoughness) double m_roughnessMax = 8.0; - bool m_fixRoughness = false; - bool m_useSurfAbs = false; - bool m_fixImpNorm = false; - bool m_xrOnly = false; + bool m_fixRoughness = false; + bool m_useSurfAbs = false; + bool m_fixImpNorm = false; + bool m_xrOnly = false; void SetBounds(double lowrough, double highrough, double highimp, double highabs); }; diff --git a/include/stochfit/ParameterStepper.h b/include/stochfit/ParameterStepper.h index 2412615..b3b14e1 100644 --- a/include/stochfit/ParameterStepper.h +++ b/include/stochfit/ParameterStepper.h @@ -1,21 +1,24 @@ #pragma once -#include "ParamVector.h" #include -class ParameterStepper { -public: - struct Config { - int sigmaSearch; - int absSearch; - int normSearch; +#include "ParamVector.h" + +class ParameterStepper +{ + public: + struct Config + { + int sigmaSearch; + int absSearch; + int normSearch; double stepSize; }; explicit ParameterStepper(Config cfg); void Step(ParamVector& params); -private: - Config m_cfg; + private: + Config m_cfg; std::mt19937 m_rng; }; diff --git a/include/stochfit/QSmear.h b/include/stochfit/QSmear.h index fd6e600..97c292f 100644 --- a/include/stochfit/QSmear.h +++ b/include/stochfit/QSmear.h @@ -1,46 +1,43 @@ #pragma once -#include "platform.h" #include #include +#include "platform.h" + namespace QSmear { -inline constexpr int Points = 13; +inline constexpr int Points = 13; inline constexpr double WeightsSum = 6.211; // Fills qspreadsinthetai[Points*i .. Points*i+12] and qspreadsinsquaredthetai // for each data point. sinthetai must already be filled (size: datapoints). // Pass a non-empty qerror span for per-point Q errors; empty span uses a // constant fractional spread (qspread). Center point is at index 6. -inline void BuildArrays( - double lambda, double qspread, - std::span sinthetai, - std::span qerror, - std::span qspreadsinthetai, - std::span qspreadsinsquaredthetai) +inline void BuildArrays(double lambda, double qspread, std::span sinthetai, std::span qerror, + std::span qspreadsinthetai, std::span qspreadsinsquaredthetai) { - for (int i = 0; i < (int)sinthetai.size(); i++) { - const double s = sinthetai[i]; + for (int i = 0; i < (int) sinthetai.size(); i++) { + const double s = sinthetai[i]; const double ds = qerror.empty() ? s * qspread : (lambda / (4.0 * std::numbers::pi)) * qerror[i]; - qspreadsinthetai[13*i] = s - 1.2*ds; - qspreadsinthetai[13*i+1] = s - 1.0*ds; - qspreadsinthetai[13*i+2] = s - 0.8*ds; - qspreadsinthetai[13*i+3] = s - 0.6*ds; - qspreadsinthetai[13*i+4] = s - 0.4*ds; - qspreadsinthetai[13*i+5] = s - 0.2*ds; - qspreadsinthetai[13*i+6] = s; - qspreadsinthetai[13*i+7] = s + 0.2*ds; - qspreadsinthetai[13*i+8] = s + 0.4*ds; - qspreadsinthetai[13*i+9] = s + 0.6*ds; - qspreadsinthetai[13*i+10] = s + 0.8*ds; - qspreadsinthetai[13*i+11] = s + 1.0*ds; - qspreadsinthetai[13*i+12] = s + 1.2*ds; - if (qspreadsinthetai[13*i] < 0.0) + qspreadsinthetai[13 * i] = s - 1.2 * ds; + qspreadsinthetai[13 * i + 1] = s - 1.0 * ds; + qspreadsinthetai[13 * i + 2] = s - 0.8 * ds; + qspreadsinthetai[13 * i + 3] = s - 0.6 * ds; + qspreadsinthetai[13 * i + 4] = s - 0.4 * ds; + qspreadsinthetai[13 * i + 5] = s - 0.2 * ds; + qspreadsinthetai[13 * i + 6] = s; + qspreadsinthetai[13 * i + 7] = s + 0.2 * ds; + qspreadsinthetai[13 * i + 8] = s + 0.4 * ds; + qspreadsinthetai[13 * i + 9] = s + 0.6 * ds; + qspreadsinthetai[13 * i + 10] = s + 0.8 * ds; + qspreadsinthetai[13 * i + 11] = s + 1.0 * ds; + qspreadsinthetai[13 * i + 12] = s + 1.2 * ds; + if (qspreadsinthetai[13 * i] < 0.0) platform_error("Error in QSpread please contact the author - the program will now crash :("); } - for (int l = 0; l < 13 * (int)sinthetai.size(); l++) + for (int l = 0; l < 13 * (int) sinthetai.size(); l++) qspreadsinsquaredthetai[l] = qspreadsinthetai[l] * qspreadsinthetai[l]; } @@ -49,22 +46,22 @@ inline void BuildArrays( // refl: size datapoints (output) inline void Apply(std::span qspreadreflpt, std::span refl) { - for (int i = 0; i < (int)refl.size(); i++) { - double c = 0.056 * qspreadreflpt[13*i]; - c += 0.135 * qspreadreflpt[13*i+1]; - c += 0.278 * qspreadreflpt[13*i+2]; - c += 0.487 * qspreadreflpt[13*i+3]; - c += 0.726 * qspreadreflpt[13*i+4]; - c += 0.923 * qspreadreflpt[13*i+5]; - c += qspreadreflpt[13*i+6]; - c += 0.923 * qspreadreflpt[13*i+7]; - c += 0.726 * qspreadreflpt[13*i+8]; - c += 0.487 * qspreadreflpt[13*i+9]; - c += 0.278 * qspreadreflpt[13*i+10]; - c += 0.135 * qspreadreflpt[13*i+11]; - c += 0.056 * qspreadreflpt[13*i+12]; + for (int i = 0; i < (int) refl.size(); i++) { + double c = 0.056 * qspreadreflpt[13 * i]; + c += 0.135 * qspreadreflpt[13 * i + 1]; + c += 0.278 * qspreadreflpt[13 * i + 2]; + c += 0.487 * qspreadreflpt[13 * i + 3]; + c += 0.726 * qspreadreflpt[13 * i + 4]; + c += 0.923 * qspreadreflpt[13 * i + 5]; + c += qspreadreflpt[13 * i + 6]; + c += 0.923 * qspreadreflpt[13 * i + 7]; + c += 0.726 * qspreadreflpt[13 * i + 8]; + c += 0.487 * qspreadreflpt[13 * i + 9]; + c += 0.278 * qspreadreflpt[13 * i + 10]; + c += 0.135 * qspreadreflpt[13 * i + 11]; + c += 0.056 * qspreadreflpt[13 * i + 12]; refl[i] = c / WeightsSum; } } -} // namespace QSmear +} // namespace QSmear diff --git a/include/stochfit/ReflectivityObjective.h b/include/stochfit/ReflectivityObjective.h index 312c91d..9b3227a 100644 --- a/include/stochfit/ReflectivityObjective.h +++ b/include/stochfit/ReflectivityObjective.h @@ -2,36 +2,31 @@ #include -class ReflectivityObjective { -public: +class ReflectivityObjective +{ + public: enum class Type : int { - LogDiff = 0, - InvRatio = 1, - LogDiffErr = 2, + LogDiff = 0, + InvRatio = 1, + LogDiffErr = 2, InvRatioErr = 3, }; - explicit ReflectivityObjective(Type t) : m_type(t) {} + explicit ReflectivityObjective(Type t) : m_type(t) + { + } // Returns sum(residual[i]²) / (n+1). - double Evaluate(std::span model, - std::span yi, - std::span eyi) const; + double Evaluate(std::span model, std::span yi, std::span eyi) const; // Fills residuals[i] for i in [low_q_offset, n - high_q_offset). // Indices outside that range are zero-filled (already set by caller or here). // Non-finite values are clamped to ±1e6. - void FillResiduals(std::span model, - std::span yi, - std::span eyi, - std::span residuals, - int low_q_offset = 0, - int high_q_offset = 0) const; + void FillResiduals(std::span model, std::span yi, std::span eyi, std::span residuals, + int low_q_offset = 0, int high_q_offset = 0) const; -private: + private: Type m_type; }; -double ComputeChiSquare(std::span model, - std::span yi, - std::span eyi); +double ComputeChiSquare(std::span model, std::span yi, std::span eyi); diff --git a/include/stochfit/SettingsStruct.h b/include/stochfit/SettingsStruct.h index c962f0d..e096bfe 100644 --- a/include/stochfit/SettingsStruct.h +++ b/include/stochfit/SettingsStruct.h @@ -4,61 +4,65 @@ // ReflSettings: passed by const reference into all C++ classes. // StochRunState: optional resume state passed to StochFit constructor (nullptr = fresh start). -#include "platform.h" #include #include +#include "platform.h" + // Forward declarations for FlatBuffer-based constructors (implemented in SettingsStruct.cpp). -namespace StochFitProto { struct ReflSettings; struct StochRunState; } +namespace StochFitProto { +struct ReflSettings; +struct StochRunState; +} // namespace StochFitProto struct ReflSettings { - std::string Directory; - std::vector Q; - std::vector Refl; - std::vector ReflError; - std::vector QError; - // Q.size() replaces the old QPoints field - double SubSLD; - double FilmSLD; - double SupSLD; - int Boxes; - double FilmAbs; - double SubAbs; - double SupAbs; - double Wavelength; - bool UseSurfAbs; - double QErr; - double Forcesig; - double RoughnessMax = 8.0; // upper bound for roughness search; determines EDP padding (6× this value) - bool XRonly; - int Resolution; - double FilmLength; - bool Impnorm; - int Objectivefunction; - double Paramtemp; + std::string Directory; + std::vector Q; + std::vector Refl; + std::vector ReflError; + std::vector QError; + // Q.size() replaces the old QPoints field + double SubSLD; + double FilmSLD; + double SupSLD; + int Boxes; + double FilmAbs; + double SubAbs; + double SupAbs; + double Wavelength; + bool UseSurfAbs; + double QErr; + double Forcesig; + double RoughnessMax = 8.0; // upper bound for roughness search; determines EDP padding (6× this value) + bool XRonly; + int Resolution; + double FilmLength; + bool Impnorm; + int Objectivefunction; + double Paramtemp; - //Annealing parameters - int Sigmasearch; - int NormalizationSearchPerc; - int AbsorptionSearchPerc; - int Algorithm; - double Inittemp; - int Platiter; - double Slope; - double Gamma; - int STUNfunc; - bool Adaptive; - int Tempiter; - int STUNdeciter; - double Gammadec; + // Annealing parameters + int Sigmasearch; + int NormalizationSearchPerc; + int AbsorptionSearchPerc; + int Algorithm; + double Inittemp; + int Platiter; + double Slope; + double Gamma; + int STUNfunc; + bool Adaptive; + int Tempiter; + int STUNdeciter; + double Gammadec; - int CritEdgeOffset; - int HighQOffset; - int Iterations; // SA iteration count; also used by LevMar + int CritEdgeOffset; + int HighQOffset; + int Iterations; // SA iteration count; also used by LevMar - ReflSettings() = default; - explicit ReflSettings(const StochFitProto::ReflSettings*); + ReflSettings() = default; + explicit ReflSettings(const StochFitProto::ReflSettings*); }; // Session state for resuming a previous SA run. @@ -67,19 +71,20 @@ struct ReflSettings // filmAbsInput = m_dBeta / m_dWaveConstant (inverse of Set_FilmAbs multiplication). // temperature = raw m_dTemp (not 1/m_dTemp — passed directly to SetTemp()). // surfAbs saved independently so it is not baked into filmAbsInput across save/load cycles. -struct StochRunState { - double roughness; - double filmAbsInput; // pre-multiplication value: Set_FilmAbs(filmAbsInput) → m_dBeta = filmAbsInput * WC - double surfAbs; // params->getSurfAbs() — saved independently - double temperature; // raw m_dTemp (passed directly to SetTemp()) - double impNorm; - double avgfSTUN; - double bestSolution; - double chiSquare; - double goodnessOfFit; - int iteration; - std::vector edValues; // Boxes+2 doubles: [supphase, box1..boxN, subphase] +struct StochRunState +{ + double roughness; + double filmAbsInput; // pre-multiplication value: Set_FilmAbs(filmAbsInput) → m_dBeta = filmAbsInput * WC + double surfAbs; // params->getSurfAbs() — saved independently + double temperature; // raw m_dTemp (passed directly to SetTemp()) + double impNorm; + double avgfSTUN; + double bestSolution; + double chiSquare; + double goodnessOfFit; + int iteration; + std::vector edValues; // Boxes+2 doubles: [supphase, box1..boxN, subphase] - StochRunState() = default; - explicit StochRunState(const StochFitProto::StochRunState*); + StochRunState() = default; + explicit StochRunState(const StochFitProto::StochRunState*); }; diff --git a/include/stochfit/StochFitDll.h b/include/stochfit/StochFitDll.h index 8ecd303..b45adc4 100644 --- a/include/stochfit/StochFitDll.h +++ b/include/stochfit/StochFitDll.h @@ -29,13 +29,14 @@ // GetRunState() → Destroy() Cancel (no save): Init() → Start() → Cancel() #pragma once -#include "platform.h" #include +#include "platform.h" + // SA lifecycle -extern "C" EXPORT void Init(const uint8_t *buf, - int len); // StochFitProto::InitRequest -extern "C" EXPORT const char *GetInitError(); +extern "C" EXPORT void Init(const uint8_t* buf, + int len); // StochFitProto::InitRequest +extern "C" EXPORT const char* GetInitError(); extern "C" EXPORT void Start(int iterations); extern "C" EXPORT void Stop(); extern "C" EXPORT void Destroy(); @@ -43,8 +44,7 @@ extern "C" EXPORT void Cancel(); // SA polling — returns bytes written into outBuf; outBuf must be pre-allocated // by caller -extern "C" EXPORT int GetData(uint8_t *outBuf, int maxLen); // → GetDataResult -extern "C" EXPORT int -GetRunState(uint8_t *outBuf, - int maxLen); // → GetRunStateResult (call after Stop()) -extern "C" EXPORT int SAParams(uint8_t *outBuf, int maxLen); // → SAParamsResult +extern "C" EXPORT int GetData(uint8_t* outBuf, int maxLen); // → GetDataResult +extern "C" EXPORT int GetRunState(uint8_t* outBuf, + int maxLen); // → GetRunStateResult (call after Stop()) +extern "C" EXPORT int SAParams(uint8_t* outBuf, int maxLen); // → SAParamsResult diff --git a/include/stochfit/StochFitHarness.h b/include/stochfit/StochFitHarness.h index a3e1110..5d2dcfb 100644 --- a/include/stochfit/StochFitHarness.h +++ b/include/stochfit/StochFitHarness.h @@ -27,6 +27,8 @@ // On stop: call Stop() to block until worker exits, then GetRunState() to read // the raw internal values, then Destroy() to clean up the object. +#include + #include #include #include @@ -42,92 +44,103 @@ #include "ParameterStepper.h" #include "ReflectivityObjective.h" #include "UnifiedReflectivity.h" -#include -using AnnealVariant = std::variant, - Anneal, Anneal>; - -struct DataSnapshot { - std::vector Q; - std::vector refl; - std::vector z; - std::vector rho; - double roughness = 0.0; - double chiSquare = 0.0; - double goodnessOfFit = 0.0; - bool isFinished = false; - int iteration = 0; -}; +using AnnealVariant = std::variant, Anneal, Anneal>; -class StochFit { -public: - StochFit(const ReflSettings &InitStruct, - const std::unique_ptr &state = {}); - ~StochFit(); - int Start(int iterations); - int Cancel(); - void Stop(); - DataSnapshot GetData(); - StochRunState GetRunState(); - tl::expected GetInitError() const { return m_initError; } - - const ReflSettings& Settings() const { return m_initStruct; } - int GetDataCount() const { return m_datapoints; } - - // GetTemperature() = 1/β (display value, same as old Get_Temp()) - // GetRawTemperature() = β (for session save, same as old Get_RawTemp()) - // SetTemperature(β) = set β directly (session restore) - double GetTemperature() const; - double GetRawTemperature() const; - void SetTemperature(double t); - double GetLowestEnergy() const; - double GetAverageFSTUN() const; - void SetAverageFSTUN(double f); - -private: - DataSnapshot GetCurrentState(); - int Processing(); - - tl::expected m_initError; - - std::thread m_thread; - std::atomic m_stop_requested; - std::atomic m_icurrentiteration{0}; - - string m_Directory; - int m_itotaliterations = 0; - - // m_initStruct must be declared before m_parratt, which holds a const ref to it. - ReflSettings m_initStruct; - - // Measured data sliced by CritEdgeOffset/HighQOffset, owned by StochFit. - int m_datapoints = 0; - std::vector m_xi; // Q values - std::vector m_yi; // reflectivity - std::vector m_eyi; // reflectivity errors - - CEDP m_cEDP; // SA scoring — worker thread only (via annealer) - CEDP m_displayEDP; // display EDP — main thread only (GetCurrentState) - ParamVector params; // SA state — worker thread only - - // Snapshot of the latest accepted solution, readable from the main thread. - // Guarded by m_displayMutex; copy only the minimum inside the lock. - struct DisplayState { - explicit DisplayState(const ReflSettings &s) : params(s) {} - ParamVector params; - std::vector refl; // SA-computed reflectivity at accepted solution +struct DataSnapshot +{ + std::vector Q; + std::vector refl; + std::vector z; + std::vector rho; + double roughness = 0.0; double chiSquare = 0.0; - double goF = 0.0; - }; - mutable std::mutex m_displayMutex; - DisplayState m_displayState; - - ParrattReflectivity m_parratt; - ReflectivityObjective m_objective; - ParameterStepper m_stepper; - - std::optional m_annealer; + double goodnessOfFit = 0.0; + bool isFinished = false; + int iteration = 0; +}; - // SA scratch buffer (one entry per measured Q point) - std::vector m_saReflBuf; +class StochFit +{ + public: + StochFit(const ReflSettings& InitStruct, const std::unique_ptr& state = {}); + ~StochFit(); + int Start(int iterations); + int Cancel(); + void Stop(); + DataSnapshot GetData(); + StochRunState GetRunState(); + tl::expected GetInitError() const + { + return m_initError; + } + + const ReflSettings& Settings() const + { + return m_initStruct; + } + int GetDataCount() const + { + return m_datapoints; + } + + // GetTemperature() = 1/β (display value, same as old Get_Temp()) + // GetRawTemperature() = β (for session save, same as old Get_RawTemp()) + // SetTemperature(β) = set β directly (session restore) + double GetTemperature() const; + double GetRawTemperature() const; + void SetTemperature(double t); + double GetLowestEnergy() const; + double GetAverageFSTUN() const; + void SetAverageFSTUN(double f); + + private: + DataSnapshot GetCurrentState(); + int Processing(); + + tl::expected m_initError; + + std::thread m_thread; + std::atomic m_stop_requested; + std::atomic m_icurrentiteration{0}; + + string m_Directory; + int m_itotaliterations = 0; + + // m_initStruct must be declared before m_parratt, which holds a const ref to it. + ReflSettings m_initStruct; + + // Measured data sliced by CritEdgeOffset/HighQOffset, owned by StochFit. + int m_datapoints = 0; + std::vector m_xi; // Q values + std::vector m_yi; // reflectivity + std::vector m_eyi; // reflectivity errors + + CEDP m_cEDP; // SA scoring — worker thread only (via annealer) + CEDP m_displayEDP; // display EDP — main thread only (GetCurrentState) + ParamVector params; // SA state — worker thread only + + // Snapshot of the latest accepted solution, readable from the main thread. + // Guarded by m_displayMutex; copy only the minimum inside the lock. + struct DisplayState + { + explicit DisplayState(const ReflSettings& s) : params(s) + { + } + ParamVector params; + std::vector refl; // SA-computed reflectivity at accepted solution + double chiSquare = 0.0; + double goF = 0.0; + }; + mutable std::mutex m_displayMutex; + DisplayState m_displayState; + + ParrattReflectivity m_parratt; + ReflectivityObjective m_objective; + ParameterStepper m_stepper; + + std::optional m_annealer; + + // SA scratch buffer (one entry per measured Q point) + std::vector m_saReflBuf; }; diff --git a/include/stochfit/UnifiedReflectivity.h b/include/stochfit/UnifiedReflectivity.h index e5a1742..fe4ad32 100644 --- a/include/stochfit/UnifiedReflectivity.h +++ b/include/stochfit/UnifiedReflectivity.h @@ -4,69 +4,69 @@ #include "LayerStack.h" #include "QSmear.h" -struct ReflConstants { - double k0; - double sup_sld; - double indexsup; - double indexsupsquared; - bool qsmear_enabled = false; - std::vector sinthetai, sinsquaredthetai; +struct ReflConstants +{ + double k0; + double sup_sld; + double indexsup; + double indexsupsquared; + bool qsmear_enabled = false; + std::vector sinthetai, sinsquaredthetai; - explicit ReflConstants(const ReflSettings &s); + explicit ReflConstants(const ReflSettings& s); }; -template struct WaveScratch { - std::vector kk, ak, rj, Rj; - int edp_points = 0; - int num_threads = 0; +template struct WaveScratch +{ + std::vector kk, ak, rj, Rj; + int edp_points = 0; + int num_threads = 0; - void resize(int new_edp_points, int new_num_threads) { - if (new_edp_points == edp_points && new_num_threads == num_threads) - return; - const int n = new_edp_points * new_num_threads; - kk.resize(n); - ak.resize(n); - rj.resize(n); - Rj.resize(n); - edp_points = new_edp_points; - num_threads = new_num_threads; - } + void resize(int new_edp_points, int new_num_threads) + { + if (new_edp_points == edp_points && new_num_threads == num_threads) + return; + const int n = new_edp_points * new_num_threads; + kk.resize(n); + ak.resize(n); + rj.resize(n); + Rj.resize(n); + edp_points = new_edp_points; + num_threads = new_num_threads; + } }; -class ParrattReflectivity { -public: - explicit ParrattReflectivity(const ReflSettings &settings); - ParrattReflectivity(const ReflSettings &settings, int n_layers); +class ParrattReflectivity +{ + public: + explicit ParrattReflectivity(const ReflSettings& settings); + ParrattReflectivity(const ReflSettings& settings, int n_layers); - // Standalone entry points — creates its own OMP parallel team. - // Safe to call from any context (no enclosing parallel required). - auto CalculateReflectivity(const LayerStack &ls) -> std::span; - auto CalculateReflectivity(const CEDP &EDP) -> std::span; + // Standalone entry points — creates its own OMP parallel team. + // Safe to call from any context (no enclosing parallel required). + auto CalculateReflectivity(const LayerStack& ls) -> std::span; + auto CalculateReflectivity(const CEDP& EDP) -> std::span; - // Cooperative entry point — must be called by ALL threads of an enclosing - // OMP parallel region. Shares the caller's thread team (no fork/join). - // BuildLayerStack runs in omp single internally. - auto CalculateReflectivityCooperative(const CEDP &EDP) -> std::span; + // Cooperative entry point — must be called by ALL threads of an enclosing + // OMP parallel region. Shares the caller's thread team (no fork/join). + // BuildLayerStack runs in omp single internally. + auto CalculateReflectivityCooperative(const CEDP& EDP) -> std::span; -private: - template - void ReflectivityCalcCoreImpl(const LayerStack &ls, - std::span sinthetai, - std::span sin2thetai, - std::span out, - int q_end = std::numeric_limits::max()); + private: + template + void ReflectivityCalcCoreImpl(const LayerStack& ls, std::span sinthetai, std::span sin2thetai, + std::span out, int q_end = std::numeric_limits::max()); - void ReflectivityCalc(const LayerStack &ls, - int q_end = std::numeric_limits::max()); - void TransparentReflectivityCalc(const LayerStack &ls); + void ReflectivityCalc(const LayerStack& ls, int q_end = std::numeric_limits::max()); + void TransparentReflectivityCalc(const LayerStack& ls); - const ReflConstants m_consts; - WaveScratch> m_complex; - WaveScratch m_real; + const ReflConstants m_consts; + WaveScratch> m_complex; + WaveScratch m_real; - bool m_qsmear_enabled = false; - std::vector m_refl_out; // Parratt output (sinthetai.size(): N or 13·N) - std::vector m_refl_smeared; // Smeared result (N); only used when qsmear on + bool m_qsmear_enabled = false; + std::vector m_refl_out; // Parratt output (sinthetai.size(): N or 13·N) + std::vector m_refl_smeared; // Smeared result (N); only used when qsmear on - LayerStack m_cooperative_ls; // written by omp single in CalculateReflectivityCooperative + LayerStack m_cooperative_ls; // written by omp single in CalculateReflectivityCooperative }; diff --git a/src/levmardll/BoxLayerBuild.cpp b/src/levmardll/BoxLayerBuild.cpp index 87ee467..6b1a789 100644 --- a/src/levmardll/BoxLayerBuild.cpp +++ b/src/levmardll/BoxLayerBuild.cpp @@ -1,27 +1,26 @@ #include "BoxLayerBuild.h" -#include + #include +#include -LayerStack BoxLayers::View(int boxes) const noexcept { +LayerStack BoxLayers::View(int boxes) const noexcept +{ LayerStack ls; - ls.rho = rho; - ls.length_mult = length_mult; - ls.sigma_sq = sigma_sq; - ls.sup_offset = 0; - ls.sub_offset = boxes; // last film-layer interface (substrate boundary is at boxes+1) - ls.has_roughness = true; // box model always applies Nevot-Croce + ls.rho = rho; + ls.length_mult = length_mult; + ls.sigma_sq = sigma_sq; + ls.sup_offset = 0; + ls.sub_offset = boxes; // last film-layer interface (substrate boundary is at boxes+1) + ls.has_roughness = true; // box model always applies Nevot-Croce // transparent = true when all imag parts are zero (no absorption) - ls.transparent = std::ranges::all_of(rho, [](const auto& c){ return c.imag() == 0.0; }); + ls.transparent = std::ranges::all_of(rho, [](const auto& c) { return c.imag() == 0.0; }); return ls; } -void BuildBoxLayers(const BoxReflSettings& rs, - std::span p, - bool one_sigma, - BoxLayers& out) { +void BuildBoxLayers(const BoxReflSettings& rs, std::span p, bool one_sigma, BoxLayers& out) +{ const int nl = rs.Boxes + 2; - const double rhofactor = 1e-6 * rs.Wavelength * rs.Wavelength / - (2.0 * std::numbers::pi); + const double rhofactor = 1e-6 * rs.Wavelength * rs.Wavelength / (2.0 * std::numbers::pi); // Factor 2 applied here so rho[i] matches the m_DEDP convention used by // ParrattReflectivity (density_profile[i] = 2 * n_i in reduced units). const double rho2 = 2.0 * rhofactor; @@ -31,28 +30,29 @@ void BuildBoxLayers(const BoxReflSettings& rs, out.sigma_sq.resize(nl - 1); // one per interface: 0..boxes // Superstrate - out.rho[0] = rs.SupSLD * rho2; + out.rho[0] = rs.SupSLD * rho2; out.length_mult[0] = {0.0, 0.0}; // Film boxes for (int i = 1; i <= rs.Boxes; ++i) { double sigma; if (one_sigma) { - out.length_mult[i] = {0.0, -2.0 * p[2*(i-1)+1]}; - out.rho[i] = p[2*(i-1)+2] * rs.SubSLD * rho2; - sigma = p[0]; + out.length_mult[i] = {0.0, -2.0 * p[2 * (i - 1) + 1]}; + out.rho[i] = p[2 * (i - 1) + 2] * rs.SubSLD * rho2; + sigma = p[0]; } else { - out.length_mult[i] = {0.0, -2.0 * p[3*(i-1)+1]}; - out.rho[i] = p[3*(i-1)+2] * rs.SubSLD * rho2; - sigma = std::fabs(p[3*(i-1)+3]); - if (sigma < 1e-8) sigma = 1e-8; + out.length_mult[i] = {0.0, -2.0 * p[3 * (i - 1) + 1]}; + out.rho[i] = p[3 * (i - 1) + 2] * rs.SubSLD * rho2; + sigma = std::fabs(p[3 * (i - 1) + 3]); + if (sigma < 1e-8) + sigma = 1e-8; } - out.sigma_sq[i-1] = -2.0 * sigma * sigma; + out.sigma_sq[i - 1] = -2.0 * sigma * sigma; } // Substrate - out.rho[nl-1] = rs.SubSLD * rho2; - out.length_mult[nl-1] = {0.0, 0.0}; + out.rho[nl - 1] = rs.SubSLD * rho2; + out.length_mult[nl - 1] = {0.0, 0.0}; // Interface between last film box and substrate uses global roughness p[0] out.sigma_sq[rs.Boxes] = -2.0 * p[0] * p[0]; diff --git a/src/levmardll/LevMardll.cpp b/src/levmardll/LevMardll.cpp index 47ed37f..47c637b 100644 --- a/src/levmardll/LevMardll.cpp +++ b/src/levmardll/LevMardll.cpp @@ -18,477 +18,445 @@ */ #include "LevMardll.h" -#include "BoxLayerBuild.h" -#include "RhoCalc.h" -#include "Settings.h" -#include "SettingsBridge.h" -#include "generated/stochfit_generated.h" -#include "platform.h" -#include "stochfit/ReflectivityObjective.h" -#include "stochfit/UnifiedReflectivity.h" + #include +#include #include #include #include -#include #include #include #include +#include "BoxLayerBuild.h" +#include "RhoCalc.h" +#include "Settings.h" +#include "SettingsBridge.h" +#include "generated/stochfit_generated.h" +#include "platform.h" +#include "stochfit/ReflectivityObjective.h" +#include "stochfit/UnifiedReflectivity.h" + // ── Helpers // ─────────────────────────────────────────────────────────────────── -static std::vector copy_fbs_vec(const flatbuffers::Vector *v) { - if (!v || v->size() == 0) - return {}; - return std::vector(v->data(), v->data() + v->size()); +static std::vector copy_fbs_vec(const flatbuffers::Vector* v) +{ + if (!v || v->size() == 0) + return {}; + return std::vector(v->data(), v->data() + v->size()); } -static void fill_box_settings(const StochFitProto::BoxReflSettings *s, - BoxReflSettings &rs) { - rs.Directory = s->directory() ? s->directory()->str() : ""; - rs.Q = copy_fbs_vec(s->q()); - rs.Refl = copy_fbs_vec(s->refl()); - rs.ReflError = copy_fbs_vec(s->refl_error()); - rs.QError = copy_fbs_vec(s->q_error()); - rs.UL = copy_fbs_vec(s->ul()); - rs.LL = copy_fbs_vec(s->ll()); - rs.ParamPercs = copy_fbs_vec(s->param_percs()); - rs.MIEDP = copy_fbs_vec(s->miedp()); - rs.ZIncrement = copy_fbs_vec(s->z_increment()); - rs.QPoints = s->q_points(); - rs.OneSigma = s->one_sigma(); - rs.SubSLD = s->sub_sld(); - rs.SupSLD = s->sup_sld(); - rs.Boxes = s->boxes(); - rs.Wavelength = s->wavelength(); - rs.QSpread = s->q_spread(); - rs.ImpNorm = s->imp_norm(); - rs.FitFunc = s->fit_func(); - rs.LowQOffset = s->low_q_offset(); - rs.HighQOffset = s->high_q_offset(); - rs.Iterations = s->iterations(); - rs.ZLength = s->z_length(); +static void fill_box_settings(const StochFitProto::BoxReflSettings* s, BoxReflSettings& rs) +{ + rs.Directory = s->directory() ? s->directory()->str() : ""; + rs.Q = copy_fbs_vec(s->q()); + rs.Refl = copy_fbs_vec(s->refl()); + rs.ReflError = copy_fbs_vec(s->refl_error()); + rs.QError = copy_fbs_vec(s->q_error()); + rs.UL = copy_fbs_vec(s->ul()); + rs.LL = copy_fbs_vec(s->ll()); + rs.ParamPercs = copy_fbs_vec(s->param_percs()); + rs.MIEDP = copy_fbs_vec(s->miedp()); + rs.ZIncrement = copy_fbs_vec(s->z_increment()); + rs.QPoints = s->q_points(); + rs.OneSigma = s->one_sigma(); + rs.SubSLD = s->sub_sld(); + rs.SupSLD = s->sup_sld(); + rs.Boxes = s->boxes(); + rs.Wavelength = s->wavelength(); + rs.QSpread = s->q_spread(); + rs.ImpNorm = s->imp_norm(); + rs.FitFunc = s->fit_func(); + rs.LowQOffset = s->low_q_offset(); + rs.HighQOffset = s->high_q_offset(); + rs.Iterations = s->iterations(); + rs.ZLength = s->z_length(); } -template -static int finish_into(flatbuffers::FlatBufferBuilder &fbb, - flatbuffers::Offset root, uint8_t *outBuf, - int maxLen) { - fbb.Finish(root); - int written = static_cast(fbb.GetSize()); - if (written > maxLen) - return -1; - std::memcpy(outBuf, fbb.GetBufferPointer(), written); - return written; +template static int finish_into(flatbuffers::FlatBufferBuilder& fbb, flatbuffers::Offset root, uint8_t* outBuf, int maxLen) +{ + fbb.Finish(root); + int written = static_cast(fbb.GetSize()); + if (written > maxLen) + return -1; + std::memcpy(outBuf, fbb.GetBufferPointer(), written); + return written; } // ── LevmarReflTask // ──────────────────────────────────────────────────────────── All state needed // by the levmar residual callback. Stack-allocated per FFI call. -struct LevmarReflTask { - ParrattReflectivity parratt; - ReflectivityObjective objective; - BoxLayers layers; - std::span Realrefl; - std::span Realreflerrors; - const BoxReflSettings *rs = nullptr; +struct LevmarReflTask +{ + ParrattReflectivity parratt; + ReflectivityObjective objective; + BoxLayers layers; + std::span Realrefl; + std::span Realreflerrors; + const BoxReflSettings* rs = nullptr; }; // Static C-style callback required by levmar's void* interface. -static void LevmarReflResidual(double *p, double *x, int m, int n, void *data) { - auto *self = static_cast(data); - BuildBoxLayers(*self->rs, std::span(p, m), self->rs->OneSigma, - self->layers); - auto refl = - self->parratt.CalculateReflectivity(self->layers.View(self->rs->Boxes)); - if (self->rs->ImpNorm) { - const double nf = self->layers.normfactor; - for (auto &r : refl) - r *= nf; - } - self->objective.FillResiduals(refl, self->Realrefl, self->Realreflerrors, - std::span(x, n), self->rs->LowQOffset, - self->rs->HighQOffset); +static void LevmarReflResidual(double* p, double* x, int m, int n, void* data) +{ + auto* self = static_cast(data); + BuildBoxLayers(*self->rs, std::span(p, m), self->rs->OneSigma, self->layers); + auto refl = self->parratt.CalculateReflectivity(self->layers.View(self->rs->Boxes)); + if (self->rs->ImpNorm) { + const double nf = self->layers.normfactor; + for (auto& r : refl) + r *= nf; + } + self->objective.FillResiduals(refl, self->Realrefl, self->Realreflerrors, std::span(x, n), self->rs->LowQOffset, + self->rs->HighQOffset); } // ── FastReflfit // ─────────────────────────────────────────────────────────────── -extern "C" EXPORT int FastReflfit(const uint8_t *inBuf, int /*inLen*/, - uint8_t *outBuf, int maxLen) { - auto *req = flatbuffers::GetRoot(inBuf); - - BoxReflSettings rs{}; - fill_box_settings(req->settings(), rs); - - std::vector params = copy_fbs_vec(req->parameters()); - int paramsize = static_cast(params.size()); - - double opts[LM_OPTS_SZ] = {LM_INIT_MU, 1e-15, 1e-15, 1e-20, -LM_DIFF_DELTA}; - - LevmarReflTask task{ - ParrattReflectivity(ToReflSettings(rs)), - ReflectivityObjective(ReflectivityObjective::Type{rs.FitFunc}), - {}, - rs.Refl, - rs.ReflError, - &rs}; - - std::vector xvec(rs.QPoints, 0.0); - std::vector work(LM_DIF_WORKSZ(paramsize, rs.QPoints) + - paramsize * rs.QPoints); - auto covar = std::span(work).subspan(LM_DIF_WORKSZ(paramsize, rs.QPoints)); - std::vector info(LM_INFO_SZ, 0.0); - - if (rs.UL.empty()) - dlevmar_dif(LevmarReflResidual, params.data(), xvec.data(), paramsize, - rs.QPoints, 1000, opts, info.data(), work.data(), covar.data(), - (void *)(&task)); - else - dlevmar_bc_dif(LevmarReflResidual, params.data(), xvec.data(), paramsize, - rs.QPoints, rs.LL.data(), rs.UL.data(), nullptr, 1000, opts, - info.data(), work.data(), covar.data(), (void *)(&task)); - - std::vector covarOut(paramsize); - for (int i = 0; i < paramsize; i++) - covarOut[i] = std::sqrt(covar[i * (paramsize + 1)]); - - flatbuffers::FlatBufferBuilder fbb(4096); - auto result = StochFitProto::CreateReflFitResult( - fbb, fbb.CreateVector(params), fbb.CreateVector(covarOut), - fbb.CreateVector(info)); - return finish_into(fbb, result, outBuf, maxLen); +extern "C" EXPORT int FastReflfit(const uint8_t* inBuf, int /*inLen*/, uint8_t* outBuf, int maxLen) +{ + auto* req = flatbuffers::GetRoot(inBuf); + + BoxReflSettings rs{}; + fill_box_settings(req->settings(), rs); + + std::vector params = copy_fbs_vec(req->parameters()); + int paramsize = static_cast(params.size()); + + double opts[LM_OPTS_SZ] = {LM_INIT_MU, 1e-15, 1e-15, 1e-20, -LM_DIFF_DELTA}; + + LevmarReflTask task{ParrattReflectivity(ToReflSettings(rs)), + ReflectivityObjective(ReflectivityObjective::Type{rs.FitFunc}), + {}, + rs.Refl, + rs.ReflError, + &rs}; + + std::vector xvec(rs.QPoints, 0.0); + std::vector work(LM_DIF_WORKSZ(paramsize, rs.QPoints) + paramsize * rs.QPoints); + auto covar = std::span(work).subspan(LM_DIF_WORKSZ(paramsize, rs.QPoints)); + std::vector info(LM_INFO_SZ, 0.0); + + if (rs.UL.empty()) + dlevmar_dif(LevmarReflResidual, params.data(), xvec.data(), paramsize, rs.QPoints, 1000, opts, info.data(), work.data(), + covar.data(), (void*) (&task)); + else + dlevmar_bc_dif(LevmarReflResidual, params.data(), xvec.data(), paramsize, rs.QPoints, rs.LL.data(), rs.UL.data(), nullptr, 1000, + opts, info.data(), work.data(), covar.data(), (void*) (&task)); + + std::vector covarOut(paramsize); + for (int i = 0; i < paramsize; i++) + covarOut[i] = std::sqrt(covar[i * (paramsize + 1)]); + + flatbuffers::FlatBufferBuilder fbb(4096); + auto result = StochFitProto::CreateReflFitResult(fbb, fbb.CreateVector(params), fbb.CreateVector(covarOut), fbb.CreateVector(info)); + return finish_into(fbb, result, outBuf, maxLen); } // ── FastReflGenerate // ────────────────────────────────────────────────────────── -extern "C" EXPORT int FastReflGenerate(const uint8_t *inBuf, int /*inLen*/, - uint8_t *outBuf, int maxLen) { - auto *req = flatbuffers::GetRoot(inBuf); +extern "C" EXPORT int FastReflGenerate(const uint8_t* inBuf, int /*inLen*/, uint8_t* outBuf, int maxLen) +{ + auto* req = flatbuffers::GetRoot(inBuf); - BoxReflSettings rs{}; - fill_box_settings(req->settings(), rs); + BoxReflSettings rs{}; + fill_box_settings(req->settings(), rs); - std::vector params = copy_fbs_vec(req->parameters()); + std::vector params = copy_fbs_vec(req->parameters()); - ParrattReflectivity parratt(ToReflSettings(rs)); - BoxLayers layers; - BuildBoxLayers(rs, params, rs.OneSigma, layers); - auto refl = parratt.CalculateReflectivity(layers.View(rs.Boxes)); - if (rs.ImpNorm) { - const double nf = layers.normfactor; - for (auto &r : refl) - r *= nf; - } + ParrattReflectivity parratt(ToReflSettings(rs)); + BoxLayers layers; + BuildBoxLayers(rs, params, rs.OneSigma, layers); + auto refl = parratt.CalculateReflectivity(layers.View(rs.Boxes)); + if (rs.ImpNorm) { + const double nf = layers.normfactor; + for (auto& r : refl) + r *= nf; + } - std::vector reflOut(refl.begin(), refl.begin() + rs.QPoints); + std::vector reflOut(refl.begin(), refl.begin() + rs.QPoints); - flatbuffers::FlatBufferBuilder fbb(4096); - auto result = - StochFitProto::CreateReflGenerateResult(fbb, fbb.CreateVector(reflOut)); - return finish_into(fbb, result, outBuf, maxLen); + flatbuffers::FlatBufferBuilder fbb(4096); + auto result = StochFitProto::CreateReflGenerateResult(fbb, fbb.CreateVector(reflOut)); + return finish_into(fbb, result, outBuf, maxLen); } // ── Rhofit ─────────────────────────────────────────────────────────────────── -extern "C" EXPORT int Rhofit(const uint8_t *inBuf, int /*inLen*/, - uint8_t *outBuf, int maxLen) { - auto *req = flatbuffers::GetRoot(inBuf); - - BoxReflSettings rs{}; - fill_box_settings(req->settings(), rs); - - std::vector params = copy_fbs_vec(req->parameters()); - int paramsize = static_cast(params.size()); - - double opts[LM_OPTS_SZ] = {LM_INIT_MU, 1e-15, 1e-15, 1e-20, -LM_DIFF_DELTA}; - - RhoCalc Rho; - Rho.init(rs); - - std::vector xvec(rs.ZLength, 0.0); - std::vector work(LM_DIF_WORKSZ(paramsize, rs.ZLength) + - paramsize * rs.ZLength); - auto covar = std::span(work).subspan(LM_DIF_WORKSZ(paramsize, rs.ZLength)); - std::vector info(LM_INFO_SZ, 0.0); - - if (rs.UL.empty()) - dlevmar_dif(RhoCalc::objective, params.data(), xvec.data(), paramsize, - rs.ZLength, 1000, opts, info.data(), work.data(), covar.data(), - (void *)(&Rho)); - else - dlevmar_bc_dif(RhoCalc::objective, params.data(), xvec.data(), paramsize, - rs.ZLength, rs.LL.data(), rs.UL.data(), nullptr, 1000, opts, - info.data(), work.data(), covar.data(), (void *)(&Rho)); - - std::vector covarOut(paramsize); - for (int i = 0; i < paramsize; i++) - covarOut[i] = std::sqrt(covar[i * (paramsize + 1)]); - - flatbuffers::FlatBufferBuilder fbb(4096); - auto result = StochFitProto::CreateRhoFitResult(fbb, fbb.CreateVector(params), - fbb.CreateVector(covarOut), - fbb.CreateVector(info)); - return finish_into(fbb, result, outBuf, maxLen); +extern "C" EXPORT int Rhofit(const uint8_t* inBuf, int /*inLen*/, uint8_t* outBuf, int maxLen) +{ + auto* req = flatbuffers::GetRoot(inBuf); + + BoxReflSettings rs{}; + fill_box_settings(req->settings(), rs); + + std::vector params = copy_fbs_vec(req->parameters()); + int paramsize = static_cast(params.size()); + + double opts[LM_OPTS_SZ] = {LM_INIT_MU, 1e-15, 1e-15, 1e-20, -LM_DIFF_DELTA}; + + RhoCalc Rho; + Rho.init(rs); + + std::vector xvec(rs.ZLength, 0.0); + std::vector work(LM_DIF_WORKSZ(paramsize, rs.ZLength) + paramsize * rs.ZLength); + auto covar = std::span(work).subspan(LM_DIF_WORKSZ(paramsize, rs.ZLength)); + std::vector info(LM_INFO_SZ, 0.0); + + if (rs.UL.empty()) + dlevmar_dif(RhoCalc::objective, params.data(), xvec.data(), paramsize, rs.ZLength, 1000, opts, info.data(), work.data(), + covar.data(), (void*) (&Rho)); + else + dlevmar_bc_dif(RhoCalc::objective, params.data(), xvec.data(), paramsize, rs.ZLength, rs.LL.data(), rs.UL.data(), nullptr, 1000, + opts, info.data(), work.data(), covar.data(), (void*) (&Rho)); + + std::vector covarOut(paramsize); + for (int i = 0; i < paramsize; i++) + covarOut[i] = std::sqrt(covar[i * (paramsize + 1)]); + + flatbuffers::FlatBufferBuilder fbb(4096); + auto result = StochFitProto::CreateRhoFitResult(fbb, fbb.CreateVector(params), fbb.CreateVector(covarOut), fbb.CreateVector(info)); + return finish_into(fbb, result, outBuf, maxLen); } // ── RhoGenerate // ─────────────────────────────────────────────────────────────── -extern "C" EXPORT int RhoGenerate(const uint8_t *inBuf, int /*inLen*/, - uint8_t *outBuf, int maxLen) { - auto *req = flatbuffers::GetRoot(inBuf); +extern "C" EXPORT int RhoGenerate(const uint8_t* inBuf, int /*inLen*/, uint8_t* outBuf, int maxLen) +{ + auto* req = flatbuffers::GetRoot(inBuf); - BoxReflSettings rs{}; - fill_box_settings(req->settings(), rs); + BoxReflSettings rs{}; + fill_box_settings(req->settings(), rs); - std::vector params = copy_fbs_vec(req->parameters()); + std::vector params = copy_fbs_vec(req->parameters()); - RhoCalc Rho; - Rho.init(rs); - Rho.mkdensity(params); - Rho.mkdensityboxmodel(params); + RhoCalc Rho; + Rho.init(rs); + Rho.mkdensity(params); + Rho.mkdensityboxmodel(params); - std::vector ed(Rho.nk.begin(), Rho.nk.begin() + rs.ZLength); - std::vector boxED(Rho.nkb.begin(), Rho.nkb.begin() + rs.ZLength); + std::vector ed(Rho.nk.begin(), Rho.nk.begin() + rs.ZLength); + std::vector boxED(Rho.nkb.begin(), Rho.nkb.begin() + rs.ZLength); - flatbuffers::FlatBufferBuilder fbb(8192); - auto result = StochFitProto::CreateRhoGenerateResult( - fbb, fbb.CreateVector(ed), fbb.CreateVector(boxED)); - return finish_into(fbb, result, outBuf, maxLen); + flatbuffers::FlatBufferBuilder fbb(8192); + auto result = StochFitProto::CreateRhoGenerateResult(fbb, fbb.CreateVector(ed), fbb.CreateVector(boxED)); + return finish_into(fbb, result, outBuf, maxLen); } // ── StochFitBoxModel helpers // ────────────────────────────────────────────────── -struct BoxSolution { - std::vector params; - std::vector covar; // per-parameter sigma: sqrt(|diag(covariance)|) - std::array info{}; - double score = 1e300; +struct BoxSolution +{ + std::vector params; + std::vector covar; // per-parameter sigma: sqrt(|diag(covariance)|) + std::array info{}; + double score = 1e300; - bool operator<(const BoxSolution &o) const { return score < o.score; } + bool operator<(const BoxSolution& o) const + { + return score < o.score; + } }; -static bool IsReasonable(const BoxSolution &sol, const BoxReflSettings &rs, - double cutoff) { - BoxLayers layers; - BuildBoxLayers(rs, sol.params, rs.OneSigma, layers); - for (int i = 1; i <= rs.Boxes; ++i) { - if (layers.length_mult[i].imag() > 0.0) { - return false; - } - if (layers.rho[i].real() < 0.0) { - return false; +static bool IsReasonable(const BoxSolution& sol, const BoxReflSettings& rs, double cutoff) +{ + BoxLayers layers; + BuildBoxLayers(rs, sol.params, rs.OneSigma, layers); + for (int i = 1; i <= rs.Boxes; ++i) { + if (layers.length_mult[i].imag() > 0.0) { + return false; + } + if (layers.rho[i].real() < 0.0) { + return false; + } } - } - if (rs.OneSigma && cutoff > 0.0) { - for (int i = 0; i < (int)sol.params.size(); ++i) { - if (sol.covar[i] > cutoff * std::fabs(sol.params[i])) - return false; + if (rs.OneSigma && cutoff > 0.0) { + for (int i = 0; i < (int) sol.params.size(); ++i) { + if (sol.covar[i] > cutoff * std::fabs(sol.params[i])) + return false; + } } - } - return true; + return true; } -static bool ApproxEqual(const BoxSolution &a, const BoxSolution &b, - double tol = 0.005) { - if (std::fabs(a.score / b.score - 1.0) > tol) - return false; - for (size_t i = 0; i < a.params.size(); ++i) { - if (b.params[i] != 0.0 && std::fabs(a.params[i] / b.params[i] - 1.0) > tol) - return false; - } - return true; +static bool ApproxEqual(const BoxSolution& a, const BoxSolution& b, double tol = 0.005) +{ + if (std::fabs(a.score / b.score - 1.0) > tol) + return false; + for (size_t i = 0; i < a.params.size(); ++i) { + if (b.params[i] != 0.0 && std::fabs(a.params[i] / b.params[i] - 1.0) > tol) + return false; + } + return true; } // ── StochFitBoxModel // ────────────────────────────────────────────────────────── -extern "C" EXPORT int StochFitBoxModel(const uint8_t *inBuf, int /*inLen*/, - uint8_t *outBuf, int maxLen) { - auto *req = flatbuffers::GetRoot(inBuf); - - BoxReflSettings rs{}; - fill_box_settings(req->settings(), rs); - - std::vector params = copy_fbs_vec(req->parameters()); - int paramsize = static_cast(params.size()); - int QSize = rs.QPoints; - - double opts[LM_OPTS_SZ] = {LM_INIT_MU, 1e-15, 1e-15, 1e-20, -LM_DIFF_DELTA}; - - // Compute initial chi-square using the selected FitFunc (same metric as the - // fits). - double bestchisquare = 0; - { - LevmarReflTask tmp{.parratt = ParrattReflectivity(ToReflSettings(rs)), - .objective = ReflectivityObjective( - ReflectivityObjective::Type{rs.FitFunc}), - .layers = {}, - .Realrefl = rs.Refl, - .Realreflerrors = rs.ReflError, - .rs = &rs}; - BuildBoxLayers(rs, params, rs.OneSigma, tmp.layers); - auto r = tmp.parratt.CalculateReflectivity(tmp.layers.View(rs.Boxes)); - if (rs.ImpNorm) - for (auto &rv : r) - rv *= tmp.layers.normfactor; - std::vector residuals(QSize, 0.0); - tmp.objective.FillResiduals(r, rs.Refl, rs.ReflError, residuals, - rs.LowQOffset, rs.HighQOffset); - for (double res : residuals) - bestchisquare += res * res; - } - - BoxSolution original; - original.params = params; - original.covar.assign(paramsize, 0.0); - original.info[1] = bestchisquare; - original.score = bestchisquare; - - std::vector temp; - temp.reserve(6000); - - omp_set_num_threads(omp_get_num_procs()); +extern "C" EXPORT int StochFitBoxModel(const uint8_t* inBuf, int /*inLen*/, uint8_t* outBuf, int maxLen) +{ + auto* req = flatbuffers::GetRoot(inBuf); + + BoxReflSettings rs{}; + fill_box_settings(req->settings(), rs); + + std::vector params = copy_fbs_vec(req->parameters()); + int paramsize = static_cast(params.size()); + int QSize = rs.QPoints; + + double opts[LM_OPTS_SZ] = {LM_INIT_MU, 1e-15, 1e-15, 1e-20, -LM_DIFF_DELTA}; + + // Compute initial chi-square using the selected FitFunc (same metric as the + // fits). + double bestchisquare = 0; + { + LevmarReflTask tmp{.parratt = ParrattReflectivity(ToReflSettings(rs)), + .objective = ReflectivityObjective(ReflectivityObjective::Type{rs.FitFunc}), + .layers = {}, + .Realrefl = rs.Refl, + .Realreflerrors = rs.ReflError, + .rs = &rs}; + BuildBoxLayers(rs, params, rs.OneSigma, tmp.layers); + auto r = tmp.parratt.CalculateReflectivity(tmp.layers.View(rs.Boxes)); + if (rs.ImpNorm) + for (auto& rv : r) + rv *= tmp.layers.normfactor; + std::vector residuals(QSize, 0.0); + tmp.objective.FillResiduals(r, rs.Refl, rs.ReflError, residuals, rs.LowQOffset, rs.HighQOffset); + for (double res : residuals) + bestchisquare += res * res; + } + + BoxSolution original; + original.params = params; + original.covar.assign(paramsize, 0.0); + original.info[1] = bestchisquare; + original.score = bestchisquare; + + std::vector temp; + temp.reserve(6000); + + omp_set_num_threads(omp_get_num_procs()); #pragma omp parallel - { - LevmarReflTask task{.parratt = ParrattReflectivity(ToReflSettings(rs)), - .objective = ReflectivityObjective( - ReflectivityObjective::Type{rs.FitFunc}), - .layers = {}, - .Realrefl = rs.Refl, - .Realreflerrors = rs.ReflError, - .rs = &rs}; - - std::mt19937 randgen(std::random_device{}() + omp_get_thread_num()); - auto IRandom = [&](double max, double min) { - return std::uniform_real_distribution(min, max)(randgen); - }; - - BoxSolution localanswer; - localanswer.params.resize(paramsize); - localanswer.covar.resize(paramsize); - std::vector locparameters(paramsize); - std::vector localvec; - localvec.reserve(1000); - - std::vector locinfo(LM_INFO_SZ); - std::vector work(LM_DIF_WORKSZ(paramsize, QSize) + - paramsize * QSize); - auto covar = std::span(work).subspan(LM_DIF_WORKSZ(paramsize, QSize)); - std::vector xvec(QSize, 0.0); + { + LevmarReflTask task{.parratt = ParrattReflectivity(ToReflSettings(rs)), + .objective = ReflectivityObjective(ReflectivityObjective::Type{rs.FitFunc}), + .layers = {}, + .Realrefl = rs.Refl, + .Realreflerrors = rs.ReflError, + .rs = &rs}; + + std::mt19937 randgen(std::random_device{}() + omp_get_thread_num()); + auto IRandom = [&](double max, double min) { return std::uniform_real_distribution(min, max)(randgen); }; + + BoxSolution localanswer; + localanswer.params.resize(paramsize); + localanswer.covar.resize(paramsize); + std::vector locparameters(paramsize); + std::vector localvec; + localvec.reserve(1000); + + std::vector locinfo(LM_INFO_SZ); + std::vector work(LM_DIF_WORKSZ(paramsize, QSize) + paramsize * QSize); + auto covar = std::span(work).subspan(LM_DIF_WORKSZ(paramsize, QSize)); + std::vector xvec(QSize, 0.0); #pragma omp for schedule(runtime) - for (int i = 0; i < rs.Iterations; i++) { - locparameters[0] = - IRandom(params[0] * rs.ParamPercs[4], params[0] * rs.ParamPercs[5]); - for (int k = 0; k < rs.Boxes; k++) { - if (rs.OneSigma) { - locparameters[2 * k + 1] = - IRandom(params[2 * k + 1] * rs.ParamPercs[0], - params[2 * k + 1] * rs.ParamPercs[1]); - locparameters[2 * k + 2] = - IRandom(params[2 * k + 2] * rs.ParamPercs[2], - params[2 * k + 2] * rs.ParamPercs[3]); - } else { - locparameters[3 * k + 1] = - IRandom(params[3 * k + 1] * rs.ParamPercs[0], - params[3 * k + 1] * rs.ParamPercs[1]); - locparameters[3 * k + 2] = - IRandom(params[3 * k + 2] * rs.ParamPercs[2], - params[3 * k + 2] * rs.ParamPercs[3]); - locparameters[3 * k + 3] = - IRandom(params[3 * k + 3] * rs.ParamPercs[4], - params[3 * k + 3] * rs.ParamPercs[5]); - } - } - locparameters[paramsize - 1] = params[paramsize - 1]; - - if (rs.UL.empty()) { - dlevmar_dif(LevmarReflResidual, locparameters.data(), xvec.data(), - paramsize, QSize, 500, opts, locinfo.data(), work.data(), - covar.data(), (void *)(&task)); - } else { - dlevmar_bc_dif(LevmarReflResidual, locparameters.data(), xvec.data(), - paramsize, QSize, rs.LL.data(), rs.UL.data(), nullptr, - 500, opts, locinfo.data(), work.data(), covar.data(), - (void *)(&task)); - } - - localanswer.params = locparameters; - for (int j = 0; j < paramsize; ++j) { - localanswer.covar[j] = std::sqrt(std::fabs(covar[j * (paramsize + 1)])); - } - std::ranges::copy(locinfo, localanswer.info.begin()); - localanswer.score = locinfo[1]; - - if (locinfo[1] < bestchisquare && - IsReasonable(localanswer, rs, rs.ParamPercs[6])) { - bool unique = true; - for (const auto &v : localvec) { - if (ApproxEqual(localanswer, v)) { - unique = false; - break; - } + for (int i = 0; i < rs.Iterations; i++) { + locparameters[0] = IRandom(params[0] * rs.ParamPercs[4], params[0] * rs.ParamPercs[5]); + for (int k = 0; k < rs.Boxes; k++) { + if (rs.OneSigma) { + locparameters[2 * k + 1] = IRandom(params[2 * k + 1] * rs.ParamPercs[0], params[2 * k + 1] * rs.ParamPercs[1]); + locparameters[2 * k + 2] = IRandom(params[2 * k + 2] * rs.ParamPercs[2], params[2 * k + 2] * rs.ParamPercs[3]); + } else { + locparameters[3 * k + 1] = IRandom(params[3 * k + 1] * rs.ParamPercs[0], params[3 * k + 1] * rs.ParamPercs[1]); + locparameters[3 * k + 2] = IRandom(params[3 * k + 2] * rs.ParamPercs[2], params[3 * k + 2] * rs.ParamPercs[3]); + locparameters[3 * k + 3] = IRandom(params[3 * k + 3] * rs.ParamPercs[4], params[3 * k + 3] * rs.ParamPercs[5]); + } + } + locparameters[paramsize - 1] = params[paramsize - 1]; + + if (rs.UL.empty()) { + dlevmar_dif(LevmarReflResidual, locparameters.data(), xvec.data(), paramsize, QSize, 500, opts, locinfo.data(), work.data(), + covar.data(), (void*) (&task)); + } else { + dlevmar_bc_dif(LevmarReflResidual, locparameters.data(), xvec.data(), paramsize, QSize, rs.LL.data(), rs.UL.data(), nullptr, + 500, opts, locinfo.data(), work.data(), covar.data(), (void*) (&task)); + } + + localanswer.params = locparameters; + for (int j = 0; j < paramsize; ++j) { + localanswer.covar[j] = std::sqrt(std::fabs(covar[j * (paramsize + 1)])); + } + std::ranges::copy(locinfo, localanswer.info.begin()); + localanswer.score = locinfo[1]; + + if (locinfo[1] < bestchisquare && IsReasonable(localanswer, rs, rs.ParamPercs[6])) { + bool unique = true; + for (const auto& v : localvec) { + if (ApproxEqual(localanswer, v)) { + unique = false; + break; + } + } + if (unique) + localvec.push_back(localanswer); + } } - if (unique) - localvec.push_back(localanswer); - } - } #pragma omp critical(AddVecs) - { - for (const auto &v : localvec) - temp.push_back(v); + { + for (const auto& v : localvec) + temp.push_back(v); + } } - } - - temp.push_back(original); - - std::vector allsolutions; - allsolutions.reserve(6000); - int tempsize = static_cast(temp.size()); - allsolutions.push_back(temp[0]); - for (int i = 1; i < tempsize; i++) { - int sz = static_cast(allsolutions.size()); - for (int j = 0; j < sz; j++) { - if (ApproxEqual(temp[i], allsolutions[j])) - break; - if (j == sz - 1) - allsolutions.push_back(temp[i]); + + temp.push_back(original); + + std::vector allsolutions; + allsolutions.reserve(6000); + int tempsize = static_cast(temp.size()); + allsolutions.push_back(temp[0]); + for (int i = 1; i < tempsize; i++) { + int sz = static_cast(allsolutions.size()); + for (int j = 0; j < sz; j++) { + if (ApproxEqual(temp[i], allsolutions[j])) + break; + if (j == sz - 1) + allsolutions.push_back(temp[i]); + } } - } - - if (!allsolutions.empty()) - std::sort(allsolutions.begin(), allsolutions.end()); - - int n = static_cast(std::min(allsolutions.size(), 999)); - std::vector outParams(params); - std::vector covarArray(static_cast(n * paramsize)); - std::vector infoOut(static_cast(n * LM_INFO_SZ)); - std::vector paramArray(static_cast(n * paramsize)); - std::vector chiSquareArray(n); - - for (int i = 0; i < n; i++) { - for (int j = 0; j < paramsize; j++) { - paramArray[i * paramsize + j] = allsolutions[i].params[j]; - covarArray[i * paramsize + j] = allsolutions[i].covar[j]; + + if (!allsolutions.empty()) + std::sort(allsolutions.begin(), allsolutions.end()); + + int n = static_cast(std::min(allsolutions.size(), 999)); + std::vector outParams(params); + std::vector covarArray(static_cast(n * paramsize)); + std::vector infoOut(static_cast(n * LM_INFO_SZ)); + std::vector paramArray(static_cast(n * paramsize)); + std::vector chiSquareArray(n); + + for (int i = 0; i < n; i++) { + for (int j = 0; j < paramsize; j++) { + paramArray[i * paramsize + j] = allsolutions[i].params[j]; + covarArray[i * paramsize + j] = allsolutions[i].covar[j]; + } + std::copy(allsolutions[i].info.begin(), allsolutions[i].info.end(), infoOut.begin() + i * LM_INFO_SZ); + chiSquareArray[i] = allsolutions[i].score; } - std::copy(allsolutions[i].info.begin(), allsolutions[i].info.end(), - infoOut.begin() + i * LM_INFO_SZ); - chiSquareArray[i] = allsolutions[i].score; - } - - flatbuffers::FlatBufferBuilder fbb(static_cast(maxLen)); - auto result = StochFitProto::CreateBoxStochFitResult( - fbb, fbb.CreateVector(outParams), fbb.CreateVector(covarArray), - fbb.CreateVector(infoOut), fbb.CreateVector(paramArray), - fbb.CreateVector(chiSquareArray), n); - return finish_into(fbb, result, outBuf, maxLen); + + flatbuffers::FlatBufferBuilder fbb(static_cast(maxLen)); + auto result = + StochFitProto::CreateBoxStochFitResult(fbb, fbb.CreateVector(outParams), fbb.CreateVector(covarArray), fbb.CreateVector(infoOut), + fbb.CreateVector(paramArray), fbb.CreateVector(chiSquareArray), n); + return finish_into(fbb, result, outBuf, maxLen); } diff --git a/src/levmardll/RhoCalc.cpp b/src/levmardll/RhoCalc.cpp index b293c20..9ff76e5 100644 --- a/src/levmardll/RhoCalc.cpp +++ b/src/levmardll/RhoCalc.cpp @@ -18,37 +18,38 @@ * */ -#include "platform.h" #include "RhoCalc.h" + #include "Settings.h" +#include "platform.h" void RhoCalc::init(const BoxReflSettings& InitStruct) { - onesigma = InitStruct.OneSigma; - boxnumber = InitStruct.Boxes; - ZIncrement= InitStruct.ZIncrement; - Zlength = InitStruct.ZLength; - MIRho = InitStruct.MIEDP; - SubSLD = InitStruct.SubSLD; - m_dSupSLD = InitStruct.SupSLD; - - nk.resize(Zlength); - nkb.resize(Zlength); - - distarray.resize(boxnumber+1); - rhoarray.resize(boxnumber+1); - rougharray.resize(boxnumber+1); - - m_LengthArray.resize(boxnumber); - m_RhoArray.resize(boxnumber); - m_SigmaArray.resize(boxnumber); + onesigma = InitStruct.OneSigma; + boxnumber = InitStruct.Boxes; + ZIncrement = InitStruct.ZIncrement; + Zlength = InitStruct.ZLength; + MIRho = InitStruct.MIEDP; + SubSLD = InitStruct.SubSLD; + m_dSupSLD = InitStruct.SupSLD; + + nk.resize(Zlength); + nkb.resize(Zlength); + + distarray.resize(boxnumber + 1); + rhoarray.resize(boxnumber + 1); + rougharray.resize(boxnumber + 1); + + m_LengthArray.resize(boxnumber); + m_RhoArray.resize(boxnumber); + m_SigmaArray.resize(boxnumber); } RhoCalc::~RhoCalc() = default; void RhoCalc::objective(double* par, double* x, int m, int n, void* data) { - RhoCalc* rhoinst = (RhoCalc*)data; + RhoCalc* rhoinst = (RhoCalc*) data; rhoinst->mkdensity(std::span(par, m)); for (int i = 0; i < rhoinst->Zlength; ++i) @@ -57,100 +58,83 @@ void RhoCalc::objective(double* par, double* x, int m, int n, void* data) void RhoCalc::Rhocalculate(double SubRough, double Zoffset) { - const double SuperphaseSLD = m_dSupSLD; - const double sqrt2 = sqrt(2.0); - double dist = 0; + const double SuperphaseSLD = m_dSupSLD; + const double sqrt2 = sqrt(2.0); + double dist = 0; - for (int i = 0; i <= boxnumber; i++) - { + for (int i = 0; i <= boxnumber; i++) { double deltarho, thick, roughness; - if (i == 0) - { - deltarho = m_RhoArray[0] * SubSLD - SuperphaseSLD; - thick = 0; + if (i == 0) { + deltarho = m_RhoArray[0] * SubSLD - SuperphaseSLD; + thick = 0; roughness = m_SigmaArray[0]; - } - else if (i == boxnumber) - { - deltarho = SubSLD - m_RhoArray[i-1] * SubSLD; + } else if (i == boxnumber) { + deltarho = SubSLD - m_RhoArray[i - 1] * SubSLD; roughness = SubRough; - thick = m_LengthArray[i-1]; - } - else - { - deltarho = (m_RhoArray[i] - m_RhoArray[i-1]) * SubSLD; - thick = m_LengthArray[i-1]; + thick = m_LengthArray[i - 1]; + } else { + deltarho = (m_RhoArray[i] - m_RhoArray[i - 1]) * SubSLD; + thick = m_LengthArray[i - 1]; roughness = m_SigmaArray[i]; } - dist += thick; - distarray[i] = dist; - rhoarray[i] = deltarho / 2.0; - rougharray[i] = roughness * sqrt2; - } - - #pragma omp parallel for - for (int j = 0; j < Zlength; j++) - { - double summ = SuperphaseSLD; - for (int i = 0; i <= boxnumber; i++) - summ += rhoarray[i] * (1.0 + erf((ZIncrement[j] - distarray[i] - Zoffset) / rougharray[i])); - - if (SubRough != 1e-16) - nk[j] = summ / SubSLD; - else - nkb[j] = summ / SubSLD; - } + dist += thick; + distarray[i] = dist; + rhoarray[i] = deltarho / 2.0; + rougharray[i] = roughness * sqrt2; + } + +#pragma omp parallel for + for (int j = 0; j < Zlength; j++) { + double summ = SuperphaseSLD; + for (int i = 0; i <= boxnumber; i++) + summ += rhoarray[i] * (1.0 + erf((ZIncrement[j] - distarray[i] - Zoffset) / rougharray[i])); + + if (SubRough != 1e-16) + nk[j] = summ / SubSLD; + else + nkb[j] = summ / SubSLD; + } } void RhoCalc::mkdensityboxmodel(std::span p) { - constexpr double SubRough = 1e-16; - const double ZOffset = p[1]; - - if (onesigma) - { - for (int i = 0; i < boxnumber; i++) - { - m_LengthArray[i] = p[2*i+2]; - m_RhoArray[i] = p[2*i+3]; - m_SigmaArray[i] = 1e-16; - } - } - else - { - for (int i = 0; i < boxnumber; i++) - { - m_LengthArray[i] = p[3*i+2]; - m_RhoArray[i] = p[3*i+3]; - m_SigmaArray[i] = 1e-16; - } - } - Rhocalculate(SubRough, ZOffset); + constexpr double SubRough = 1e-16; + const double ZOffset = p[1]; + + if (onesigma) { + for (int i = 0; i < boxnumber; i++) { + m_LengthArray[i] = p[2 * i + 2]; + m_RhoArray[i] = p[2 * i + 3]; + m_SigmaArray[i] = 1e-16; + } + } else { + for (int i = 0; i < boxnumber; i++) { + m_LengthArray[i] = p[3 * i + 2]; + m_RhoArray[i] = p[3 * i + 3]; + m_SigmaArray[i] = 1e-16; + } + } + Rhocalculate(SubRough, ZOffset); } void RhoCalc::mkdensity(std::span p) { - const double SubRough = p[0]; - const double ZOffset = p[1]; - - if (onesigma) - { - for (int i = 0; i < boxnumber; i++) - { - m_LengthArray[i] = p[2*i+2]; - m_RhoArray[i] = p[2*i+3]; - m_SigmaArray[i] = p[0]; - } - } - else - { - for (int i = 0; i < boxnumber; i++) - { - m_LengthArray[i] = p[3*i+2]; - m_RhoArray[i] = p[3*i+3]; - m_SigmaArray[i] = p[3*i+4]; - } - } - Rhocalculate(SubRough, ZOffset); + const double SubRough = p[0]; + const double ZOffset = p[1]; + + if (onesigma) { + for (int i = 0; i < boxnumber; i++) { + m_LengthArray[i] = p[2 * i + 2]; + m_RhoArray[i] = p[2 * i + 3]; + m_SigmaArray[i] = p[0]; + } + } else { + for (int i = 0; i < boxnumber; i++) { + m_LengthArray[i] = p[3 * i + 2]; + m_RhoArray[i] = p[3 * i + 3]; + m_SigmaArray[i] = p[3 * i + 4]; + } + } + Rhocalculate(SubRough, ZOffset); } diff --git a/src/levmardll/SettingsBridge.cpp b/src/levmardll/SettingsBridge.cpp index 6d432dc..0ff5d25 100644 --- a/src/levmardll/SettingsBridge.cpp +++ b/src/levmardll/SettingsBridge.cpp @@ -1,18 +1,19 @@ #include "SettingsBridge.h" + #include -ReflSettings ToReflSettings(const BoxReflSettings& rs) { +ReflSettings ToReflSettings(const BoxReflSettings& rs) +{ ReflSettings out; - out.Q = rs.Q; - out.QError = rs.QError; + out.Q = rs.Q; + out.QError = rs.QError; out.Wavelength = rs.Wavelength; // ReflConstants uses: sup_sld = SupSLD; indexsup = 1 - SupSLD/2. // The Parratt formula requires SupSLD = 2 * raw_SLD * rhofactor where // rhofactor = 1e-6 * lambda^2 / (2*pi). BoxReflSettings stores the raw // SLD in 10^-6 Å^-2, so convert here. - const double rhofactor = 1e-6 * rs.Wavelength * rs.Wavelength / - (2.0 * std::numbers::pi); + const double rhofactor = 1e-6 * rs.Wavelength * rs.Wavelength / (2.0 * std::numbers::pi); out.SupSLD = 2.0 * rs.SupSLD * rhofactor; - out.QErr = rs.QSpread; // both are in percentage units + out.QErr = rs.QSpread; // both are in percentage units return out; } diff --git a/src/stochfitdll/AnnealPolicies.cpp b/src/stochfitdll/AnnealPolicies.cpp index 3814657..178fd87 100644 --- a/src/stochfitdll/AnnealPolicies.cpp +++ b/src/stochfitdll/AnnealPolicies.cpp @@ -1,12 +1,18 @@ #include "AnnealPolicies.h" + #include // ── SimulatedPolicy ─────────────────────────────────────────────────────────── SimulatedPolicy::SimulatedPolicy(double initTemp, double slope, int platIter) - : m_dTemp(1.0 / initTemp), m_slope(slope), m_platIter(platIter) {} + : m_dTemp(1.0 / initTemp), + m_slope(slope), + m_platIter(platIter) +{ +} -bool SimulatedPolicy::Accept(double curE, double candE, double, std::mt19937& rng) { +bool SimulatedPolicy::Accept(double curE, double candE, double, std::mt19937& rng) +{ const double deltaE = candE - curE; ++m_iter; if (m_iter % m_platIter == 0) @@ -14,21 +20,28 @@ bool SimulatedPolicy::Accept(double curE, double candE, double, std::mt19937& rn return std::uniform_real_distribution(0.0, 100.0)(rng) < ProbCalc(deltaE); } -void SimulatedPolicy::Schedule() { +void SimulatedPolicy::Schedule() +{ if (m_dTemp > 1e-30) m_dTemp /= m_slope; } // ── StunPolicy ──────────────────────────────────────────────────────────────── -StunPolicy::StunPolicy(double initTemp, double slope, int platIter, - double gamma, int stunFunc, int tempIter, bool adaptive) +StunPolicy::StunPolicy(double initTemp, double slope, int platIter, double gamma, int stunFunc, int tempIter, bool adaptive) : m_dTemp(adaptive ? 10.0 : 1.0 / initTemp), - m_slope(slope), m_platIter(platIter), - m_gamma(gamma), m_averageFSTUN(initTemp), - m_stunFunc(stunFunc), m_tempIter(tempIter), m_adaptive(adaptive) {} + m_slope(slope), + m_platIter(platIter), + m_gamma(gamma), + m_averageFSTUN(initTemp), + m_stunFunc(stunFunc), + m_tempIter(tempIter), + m_adaptive(adaptive) +{ +} -bool StunPolicy::Accept(double curE, double candE, double bestE, std::mt19937& rng) { +bool StunPolicy::Accept(double curE, double candE, double bestE, std::mt19937& rng) +{ const double deltaE = fSTUN(candE, bestE) - fSTUN(curE, bestE); if (m_adaptive) { @@ -40,7 +53,8 @@ bool StunPolicy::Accept(double curE, double candE, double bestE, std::mt19937& r // Recompute rolling average once the window is full. double windowAvg = 0.0; if (static_cast(m_qWindow.size()) >= m_platIter) { - for (double v : m_qWindow) windowAvg += v; + for (double v : m_qWindow) + windowAvg += v; windowAvg /= static_cast(m_qWindow.size()); } @@ -57,28 +71,34 @@ bool StunPolicy::Accept(double curE, double candE, double bestE, std::mt19937& r return std::uniform_real_distribution(0.0, 100.0)(rng) < ProbCalc(deltaE); } -double StunPolicy::fSTUN(double val, double bestE) const { +double StunPolicy::fSTUN(double val, double bestE) const +{ const double x = m_gamma * (val - bestE); switch (m_stunFunc) { - case 0: return -std::exp(-x); - case 1: return std::sinh(x) - 1.0; - default: return std::log(x + std::sqrt(x * x + 1.0)) - 1.0; // asinh - 1 + case 0: + return -std::exp(-x); + case 1: + return std::sinh(x) - 1.0; + default: + return std::log(x + std::sqrt(x * x + 1.0)) - 1.0; // asinh - 1 } } -void StunPolicy::AdjustTemp(double averageSTUN) { +void StunPolicy::AdjustTemp(double averageSTUN) +{ // Wenzel & Hamacher (PRL 82:3003, 1999): reduce β when avg fSTUN exceeds // threshold (tunneling phase), increase β otherwise (local-search phase). if (averageSTUN > m_averageFSTUN) { if (m_dTemp > 1e-200) - m_dTemp *= m_slope; // β↓, T↑ — keep tunneling + m_dTemp *= m_slope; // β↓, T↑ — keep tunneling } else { if (m_dTemp < 1e200) - m_dTemp /= m_slope; // β↑, T↓ — settle into basin + m_dTemp /= m_slope; // β↑, T↓ — settle into basin } } -void StunPolicy::Schedule() { +void StunPolicy::Schedule() +{ if (m_dTemp > 1e-30) m_dTemp /= m_slope; } diff --git a/src/stochfitdll/CEDP.cpp b/src/stochfitdll/CEDP.cpp index e784b0f..6733dc8 100644 --- a/src/stochfitdll/CEDP.cpp +++ b/src/stochfitdll/CEDP.cpp @@ -1,94 +1,100 @@ #include "CEDP.h" + #include "platform.h" -std::pair CEDP::GetOffSets() const { - const double front = m_DEDP.front().real(); - const double back = m_DEDP.back().real(); - - const int sup_offset = static_cast( - std::ranges::find_if(m_DEDP, [&](const auto &v) { - return v.real() != front; - }) - m_DEDP.begin()); - - const int sub_offset = static_cast(m_DEDP.size()) - 1 - - static_cast( - std::ranges::find_if(m_DEDP | std::views::reverse, [&](const auto &v) { - return v.real() != back; - }) - m_DEDP.rbegin()); - - return {sup_offset, sub_offset}; -} - -void CEDP::Init(const ReflSettings &InitStruct) { - int resolution = InitStruct.Resolution > 0 ? InitStruct.Resolution : 3; - m_dDz0 = 1.0 / resolution; - m_dLambda = InitStruct.Wavelength; - m_bUseSurfAbs = InitStruct.UseSurfAbs; - m_dWaveConstant = m_dLambda * m_dLambda / (2.0 * std::numbers::pi); - m_dRho = InitStruct.FilmSLD * 1e-6 * m_dWaveConstant; - // Padding sized to 6× the roughness upper bound: the erf tail falls below - // 1e-9 at 6σ. Scales automatically with RoughnessMax so high-roughness - // films (e.g. gold at 15 Å roughness → 90 Å padding) work without changes. - // FilmSlack adds 7 Å past the last box so the substrate erf tail converges. - const double leftOffset = 6.0 * InitStruct.RoughnessMax; - const double substrateOffset = 6.0 * InitStruct.RoughnessMax; - const int FilmSlack = 7; - - m_iLayers = static_cast(leftOffset + InitStruct.FilmLength + FilmSlack + - substrateOffset); - m_iLayers *= resolution; - - if (InitStruct.UseSurfAbs != 0) { - m_dBeta = InitStruct.FilmAbs * m_dWaveConstant; - m_dBeta_Sub = InitStruct.SubAbs * m_dWaveConstant; - m_dBeta_Sup = InitStruct.SupAbs * m_dWaveConstant; - } else { - m_dBeta = m_dBeta_Sub = m_dBeta_Sup = 0; - } - - m_EDP.resize(m_iLayers); - m_DEDP.resize(m_iLayers); - m_fEDSpacingArray.resize(m_iLayers); - - m_fDistArray.resize(InitStruct.Boxes + 2); - m_fRhoArray.resize(InitStruct.Boxes + 2); - m_fImagRhoArray.resize(InitStruct.Boxes + 2); - - // Precompute per-layer length multiplier for BuildLayerStack. - m_length_mult.assign(m_iLayers, std::complex{0.0, -2.0 * m_dDz0}); - - for (int i = 0; i < m_iLayers; i++) { - m_fEDSpacingArray[i] = i * m_dDz0 - leftOffset; - } - - for (int k = 0; k < InitStruct.Boxes + 2; k++) { - m_fDistArray[k] = - k * (InitStruct.FilmLength + (double)FilmSlack) / InitStruct.Boxes; - } -} - -void CEDP::GenerateEDP(ParamVector &g) { - // Standalone path: wrap serial setup in its own single so standalone callers - // (display, InitEnergy, LevMar) work correctly outside a persistent OMP team. +std::pair CEDP::GetOffSets() const +{ + const double front = m_DEDP.front().real(); + const double back = m_DEDP.back().real(); + + const int sup_offset = + static_cast(std::ranges::find_if(m_DEDP, [&](const auto& v) { return v.real() != front; }) - m_DEDP.begin()); + + const int sub_offset = + static_cast(m_DEDP.size()) - 1 - + static_cast(std::ranges::find_if(m_DEDP | std::views::reverse, [&](const auto& v) { return v.real() != back; }) - + m_DEDP.rbegin()); + + return {sup_offset, sub_offset}; +} + +void CEDP::Init(const ReflSettings& InitStruct) +{ + int resolution = InitStruct.Resolution > 0 ? InitStruct.Resolution : 3; + m_dDz0 = 1.0 / resolution; + m_dLambda = InitStruct.Wavelength; + m_bUseSurfAbs = InitStruct.UseSurfAbs; + m_dWaveConstant = m_dLambda * m_dLambda / (2.0 * std::numbers::pi); + m_dRho = InitStruct.FilmSLD * 1e-6 * m_dWaveConstant; + // Padding sized to 6× the roughness upper bound: the erf tail falls below + // 1e-9 at 6σ. Scales automatically with RoughnessMax so high-roughness + // films (e.g. gold at 15 Å roughness → 90 Å padding) work without changes. + // FilmSlack adds 7 Å past the last box so the substrate erf tail converges. + const double leftOffset = 6.0 * InitStruct.RoughnessMax; + const double substrateOffset = 6.0 * InitStruct.RoughnessMax; + const int FilmSlack = 7; + + m_iLayers = static_cast(leftOffset + InitStruct.FilmLength + FilmSlack + substrateOffset); + m_iLayers *= resolution; + + if (InitStruct.UseSurfAbs != 0) { + m_dBeta = InitStruct.FilmAbs * m_dWaveConstant; + m_dBeta_Sub = InitStruct.SubAbs * m_dWaveConstant; + m_dBeta_Sup = InitStruct.SupAbs * m_dWaveConstant; + } else { + m_dBeta = m_dBeta_Sub = m_dBeta_Sup = 0; + } + + m_EDP.resize(m_iLayers); + m_DEDP.resize(m_iLayers); + m_fEDSpacingArray.resize(m_iLayers); + + m_fDistArray.resize(InitStruct.Boxes + 2); + m_fRhoArray.resize(InitStruct.Boxes + 2); + m_fImagRhoArray.resize(InitStruct.Boxes + 2); + + // Precompute per-layer length multiplier for BuildLayerStack. + m_length_mult.assign(m_iLayers, std::complex{0.0, -2.0 * m_dDz0}); + + for (int i = 0; i < m_iLayers; i++) { + m_fEDSpacingArray[i] = i * m_dDz0 - leftOffset; + } + + for (int k = 0; k < InitStruct.Boxes + 2; k++) { + m_fDistArray[k] = k * (InitStruct.FilmLength + (double) FilmSlack) / InitStruct.Boxes; + } +} + +void CEDP::GenerateEDP(ParamVector& g) +{ + // Standalone path: wrap serial setup in its own single so standalone callers + // (display, InitEnergy, LevMar) work correctly outside a persistent OMP team. #pragma omp single - { FillBoxArrays(g); } + { + FillBoxArrays(g); + } - if (!m_bUseSurfAbs) BuildEDP(g); - else BuildEDP(g); + if (!m_bUseSurfAbs) + BuildEDP(g); + else + BuildEDP(g); #pragma omp single - { - auto [sup, sub] = GetOffSets(); - m_supOff = sup; - m_subOff = sub; - } + { + auto [sup, sub] = GetOffSets(); + m_supOff = sup; + m_subOff = sub; + } } // Cooperative path: FillBoxArrays was already called in the caller's omp single // (merged with Step). GetOffSets is called later in BuildLayerStackFull. -void CEDP::GenerateEDPCooperative(ParamVector &g) { - if (!m_bUseSurfAbs) BuildEDP(g); - else BuildEDP(g); +void CEDP::GenerateEDPCooperative(ParamVector& g) +{ + if (!m_bUseSurfAbs) + BuildEDP(g); + else + BuildEDP(g); } // The code for the ED calculation section is loosely based on the electron @@ -99,136 +105,148 @@ void CEDP::GenerateEDPCooperative(ParamVector &g) { // films, the absorbance is negligible. // Fills m_fRhoArray (and m_fImagRhoArray for absorbing films). No omp pragma — // designed to be called from within the caller's omp single block. -void CEDP::FillBoxArrays(ParamVector &g) { - if (!m_bUseSurfAbs) FillBoxArraysImpl(g); - else FillBoxArraysImpl(g); -} - -template -void CEDP::FillBoxArraysImpl(ParamVector &g) { - const int refllayers = g.RealParamsSize() - 1; - - if constexpr (!Absorbing) - m_EDP[0].imag(0.0); - - for (int k = 0; k < refllayers; k++) { - m_fRhoArray[k] = - m_dRho * (g.GetRealParams(k + 1) - g.GetRealParams(k)) * 0.5; - - if constexpr (Absorbing) { - if (k == 0) { - m_fImagRhoArray[k] = - (m_dBeta * g.GetSurfAbs() * g.GetRealParams(k + 1) / - g.GetRealParams(refllayers) - - m_dBeta_Sup) / - 2.0; - } else if (k == refllayers - 1) { - m_fImagRhoArray[k] = - (m_dBeta_Sub - m_dBeta * g.GetSurfAbs() * g.GetRealParams(k) / - g.GetRealParams(refllayers)) / - 2.0; - } else { - m_fImagRhoArray[k] = - (m_dBeta * g.GetSurfAbs() * g.GetRealParams(k + 1) / - g.GetRealParams(refllayers) - - (m_dBeta * g.GetSurfAbs() * g.GetRealParams(k) / - g.GetRealParams(refllayers)) / - 2.0); - } - } - } +void CEDP::FillBoxArrays(ParamVector& g) +{ + if (!m_bUseSurfAbs) + FillBoxArraysImpl(g); + else + FillBoxArraysImpl(g); } -template -void CEDP::BuildEDP(ParamVector &g) { - const int reflpoints = m_iLayers; - const int refllayers = g.RealParamsSize() - 1; - const double supersld = g.GetRealParams(0) * m_dRho; +template void CEDP::FillBoxArraysImpl(ParamVector& g) +{ + const int refllayers = g.RealParamsSize() - 1; - double roughness = g.GetRoughness(); - if (roughness < 1e-6) roughness = 1e-6; - roughness = 1.0 / (roughness * std::sqrt(2.0)); + if constexpr (!Absorbing) + m_EDP[0].imag(0.0); - // m_fRhoArray was filled by FillBoxArrays (in an omp single with implicit barrier) - // before this function is called. No serial setup here. + for (int k = 0; k < refllayers; k++) { + m_fRhoArray[k] = m_dRho * (g.GetRealParams(k + 1) - g.GetRealParams(k)) * 0.5; + + if constexpr (Absorbing) { + if (k == 0) { + m_fImagRhoArray[k] = (m_dBeta * g.GetSurfAbs() * g.GetRealParams(k + 1) / g.GetRealParams(refllayers) - m_dBeta_Sup) / 2.0; + } else if (k == refllayers - 1) { + m_fImagRhoArray[k] = (m_dBeta_Sub - m_dBeta * g.GetSurfAbs() * g.GetRealParams(k) / g.GetRealParams(refllayers)) / 2.0; + } else { + m_fImagRhoArray[k] = (m_dBeta * g.GetSurfAbs() * g.GetRealParams(k + 1) / g.GetRealParams(refllayers) - + (m_dBeta * g.GetSurfAbs() * g.GetRealParams(k) / g.GetRealParams(refllayers)) / 2.0); + } + } + } +} - double dist; -#pragma omp for private(dist) - for (int i = 0; i < reflpoints; i++) { - if constexpr (Absorbing) - m_EDP[i] = std::complex(supersld, m_dBeta); - else - m_EDP[i].real(supersld); +template void CEDP::BuildEDP(ParamVector& g) +{ + const int reflpoints = m_iLayers; + const int refllayers = g.RealParamsSize() - 1; + const double supersld = g.GetRealParams(0) * m_dRho; - for (int k = 0; k < refllayers; k++) { - dist = (m_fEDSpacingArray[i] - m_fDistArray[k]) * roughness; + double roughness = g.GetRoughness(); + if (roughness < 1e-6) + roughness = 1e-6; + roughness = 1.0 / (roughness * std::sqrt(2.0)); - if (dist > 6.0) { - if constexpr (Absorbing) - m_EDP[i] += std::complex(m_fRhoArray[k] * 2.0, - m_fImagRhoArray[k] * 2.0); - else - m_EDP[i] += m_fRhoArray[k] * 2.0; - } else if (dist > -6.0) { - const double erf_val = 1.0 + std::erf(dist); + // m_fRhoArray was filled by FillBoxArrays (in an omp single with implicit barrier) + // before this function is called. No serial setup here. + + double dist; +#pragma omp for private(dist) + for (int i = 0; i < reflpoints; i++) { if constexpr (Absorbing) - m_EDP[i] += std::complex(m_fRhoArray[k] * erf_val, - m_fImagRhoArray[k] * erf_val); + m_EDP[i] = std::complex(supersld, m_dBeta); else - m_EDP[i] += m_fRhoArray[k] * erf_val; - } + m_EDP[i].real(supersld); + + for (int k = 0; k < refllayers; k++) { + dist = (m_fEDSpacingArray[i] - m_fDistArray[k]) * roughness; + + if (dist > 6.0) { + if constexpr (Absorbing) + m_EDP[i] += std::complex(m_fRhoArray[k] * 2.0, m_fImagRhoArray[k] * 2.0); + else + m_EDP[i] += m_fRhoArray[k] * 2.0; + } else if (dist > -6.0) { + const double erf_val = 1.0 + std::erf(dist); + if constexpr (Absorbing) + m_EDP[i] += std::complex(m_fRhoArray[k] * erf_val, m_fImagRhoArray[k] * erf_val); + else + m_EDP[i] += m_fRhoArray[k] * erf_val; + } + } + + m_DEDP[i] = 2.0 * m_EDP[i]; } - - m_DEDP[i] = 2.0 * m_EDP[i]; - } } -LayerStack CEDP::BuildLayerStack() const { - LayerStack ls; - ls.rho = m_DEDP; - ls.length_mult = m_length_mult; - ls.sup_offset = m_supOff; - ls.sub_offset = m_subOff; - ls.transparent = !m_bUseSurfAbs; - ls.has_roughness = false; - return ls; +LayerStack CEDP::BuildLayerStack() const +{ + LayerStack ls; + ls.rho = m_DEDP; + ls.length_mult = m_length_mult; + ls.sup_offset = m_supOff; + ls.sub_offset = m_subOff; + ls.transparent = !m_bUseSurfAbs; + ls.has_roughness = false; + return ls; } // Cooperative variant: computes sup/sub offsets fresh (skips the cached values // that GenerateEDPCooperative doesn't update). m_supOff/m_subOff are mutable so // this can be called on a const CEDP& and still keeps the cache in sync. -LayerStack CEDP::BuildLayerStackFull() const { - auto [sup, sub] = GetOffSets(); - m_supOff = sup; - m_subOff = sub; - LayerStack ls; - ls.rho = m_DEDP; - ls.length_mult = m_length_mult; - ls.sup_offset = sup; - ls.sub_offset = sub; - ls.transparent = !m_bUseSurfAbs; - ls.has_roughness = false; - return ls; +LayerStack CEDP::BuildLayerStackFull() const +{ + auto [sup, sub] = GetOffSets(); + m_supOff = sup; + m_subOff = sub; + LayerStack ls; + ls.rho = m_DEDP; + ls.length_mult = m_length_mult; + ls.sup_offset = sup; + ls.sub_offset = sub; + ls.transparent = !m_bUseSurfAbs; + ls.has_roughness = false; + return ls; } -int CEDP::Get_EDPPointCount() const { return m_iLayers; } +int CEDP::Get_EDPPointCount() const +{ + return m_iLayers; +} -bool CEDP::Get_UseABS() const { return m_bUseSurfAbs; } +bool CEDP::Get_UseABS() const +{ + return m_bUseSurfAbs; +} -double CEDP::Get_Dz() const { return m_dDz0; } +double CEDP::Get_Dz() const +{ + return m_dDz0; +} -double CEDP::Get_LeftOffset() const { - return m_fEDSpacingArray.empty() ? 0.0 : -m_fEDSpacingArray[0]; +double CEDP::Get_LeftOffset() const +{ + return m_fEDSpacingArray.empty() ? 0.0 : -m_fEDSpacingArray[0]; } -double CEDP::Get_FilmAbs() const { return m_dBeta; } +double CEDP::Get_FilmAbs() const +{ + return m_dBeta; +} -double CEDP::Get_FilmAbsInput() const { - // Returns the value that, when passed to Set_FilmAbs(), reproduces m_dBeta. - // Set_FilmAbs(x) stores x * m_dWaveConstant, so x = m_dBeta / m_dWaveConstant. - return (m_dWaveConstant > 0.0) ? m_dBeta / m_dWaveConstant : 0.0; +double CEDP::Get_FilmAbsInput() const +{ + // Returns the value that, when passed to Set_FilmAbs(), reproduces m_dBeta. + // Set_FilmAbs(x) stores x * m_dWaveConstant, so x = m_dBeta / m_dWaveConstant. + return (m_dWaveConstant > 0.0) ? m_dBeta / m_dWaveConstant : 0.0; } -double CEDP::Get_WaveConstant() const { return m_dWaveConstant; } +double CEDP::Get_WaveConstant() const +{ + return m_dWaveConstant; +} -void CEDP::Set_FilmAbs(double abs) { m_dBeta = abs * m_dWaveConstant; } +void CEDP::Set_FilmAbs(double abs) +{ + m_dBeta = abs * m_dWaveConstant; +} diff --git a/src/stochfitdll/ParamVector.cpp b/src/stochfitdll/ParamVector.cpp index 1389895..bcbc4a9 100644 --- a/src/stochfitdll/ParamVector.cpp +++ b/src/stochfitdll/ParamVector.cpp @@ -18,9 +18,10 @@ * */ -#include "platform.h" #include "ParamVector.h" +#include "platform.h" + ParamVector::ParamVector(const ReflSettings& s) : m_boxes(s.Boxes), m_paramCount(s.Boxes), @@ -37,7 +38,7 @@ ParamVector::ParamVector(const ReflSettings& s) } if (s.Forcesig > 0.0) { m_fixRoughness = true; - m_roughness = s.Forcesig; + m_roughness = s.Forcesig; } else { m_roughnessIdx = m_paramCount++; } @@ -65,15 +66,15 @@ void ParamVector::SetBounds(double lowrough, double highrough, double highimp, d if (!m_fixRoughness) { m_high[m_roughnessIdx] = highrough; - m_low[m_roughnessIdx] = lowrough; + m_low[m_roughnessIdx] = lowrough; } if (m_fixImpNorm) { m_high[m_impNormIdx] = highimp; - m_low[m_impNormIdx] = 0.0; + m_low[m_impNormIdx] = 0.0; } if (m_useSurfAbs) { m_high[m_surfAbsIdx] = highabs; - m_low[m_surfAbsIdx] = 0.0; + m_low[m_surfAbsIdx] = 0.0; } } diff --git a/src/stochfitdll/ParameterStepper.cpp b/src/stochfitdll/ParameterStepper.cpp index f8e1169..79141ee 100644 --- a/src/stochfitdll/ParameterStepper.cpp +++ b/src/stochfitdll/ParameterStepper.cpp @@ -1,37 +1,31 @@ #include "ParameterStepper.h" -ParameterStepper::ParameterStepper(Config cfg) - : m_cfg(cfg), m_rng(std::random_device{}()) {} +ParameterStepper::ParameterStepper(Config cfg) : m_cfg(cfg), m_rng(std::random_device{}()) +{ +} -void ParameterStepper::Step(ParamVector& params) { +void ParameterStepper::Step(ParamVector& params) +{ constexpr double kRoughMult = 5.0 / 3.0; - const int ii = std::uniform_int_distribution(0, params.BoxCount() - 1)(m_rng); + const int ii = std::uniform_int_distribution(0, params.BoxCount() - 1)(m_rng); const int perc = std::uniform_int_distribution(1, 100)(m_rng); const int sigmaTop = m_cfg.sigmaSearch; - const int absTop = sigmaTop + m_cfg.absSearch; - const int normTop = absTop + m_cfg.normSearch; + const int absTop = sigmaTop + m_cfg.absSearch; + const int normTop = absTop + m_cfg.normSearch; if (perc > normTop) { - params.SetMutatableParameter(ii, - std::uniform_real_distribution( - params.GetMutatableParameter(ii) - m_cfg.stepSize, - params.GetMutatableParameter(ii) + m_cfg.stepSize)(m_rng)); + params.SetMutatableParameter(ii, std::uniform_real_distribution(params.GetMutatableParameter(ii) - m_cfg.stepSize, + params.GetMutatableParameter(ii) + m_cfg.stepSize)(m_rng)); } else if (perc <= sigmaTop) { - params.SetRoughness( - std::uniform_real_distribution( - params.GetRoughness() * (1.0 - kRoughMult * m_cfg.stepSize), - params.GetRoughness() * (1.0 + kRoughMult * m_cfg.stepSize))(m_rng)); + params.SetRoughness(std::uniform_real_distribution(params.GetRoughness() * (1.0 - kRoughMult * m_cfg.stepSize), + params.GetRoughness() * (1.0 + kRoughMult * m_cfg.stepSize))(m_rng)); } else if (perc <= absTop) { - params.SetSurfAbs( - std::uniform_real_distribution( - params.GetSurfAbs() * (1.0 - m_cfg.stepSize), - params.GetSurfAbs() * (1.0 + m_cfg.stepSize))(m_rng)); + params.SetSurfAbs(std::uniform_real_distribution(params.GetSurfAbs() * (1.0 - m_cfg.stepSize), + params.GetSurfAbs() * (1.0 + m_cfg.stepSize))(m_rng)); } else { - params.SetImpNorm( - std::uniform_real_distribution( - params.GetImpNorm() * (1.0 - m_cfg.stepSize), - params.GetImpNorm() * (1.0 + m_cfg.stepSize))(m_rng)); + params.SetImpNorm(std::uniform_real_distribution(params.GetImpNorm() * (1.0 - m_cfg.stepSize), + params.GetImpNorm() * (1.0 + m_cfg.stepSize))(m_rng)); } } diff --git a/src/stochfitdll/ReflectivityObjective.cpp b/src/stochfitdll/ReflectivityObjective.cpp index 44b6853..52363e3 100644 --- a/src/stochfitdll/ReflectivityObjective.cpp +++ b/src/stochfitdll/ReflectivityObjective.cpp @@ -1,92 +1,92 @@ #include "ReflectivityObjective.h" + #include -void ReflectivityObjective::FillResiduals(std::span model, - std::span yi, - std::span eyi, - std::span residuals, - int low_q_offset, - int high_q_offset) const { +void ReflectivityObjective::FillResiduals(std::span model, std::span yi, std::span eyi, + std::span residuals, int low_q_offset, int high_q_offset) const +{ const int n = static_cast(model.size()); std::fill(residuals.begin(), residuals.end(), 0.0); const int hi = n - high_q_offset; switch (m_type) { - case Type::LogDiff: - for (int i = low_q_offset; i < hi; ++i) { - double r = std::log(yi[i]) - std::log(model[i]); - residuals[i] = std::isfinite(r) ? r : (r > 0 ? 1e6 : -1e6); - } - break; - case Type::InvRatio: - for (int i = low_q_offset; i < hi; ++i) { - double r = yi[i] / model[i]; - if (r < 1.0) r = 1.0 / r; - double res = 1.0 - r; - residuals[i] = std::isfinite(res) ? res : -1e6; - } - break; - case Type::LogDiffErr: - for (int i = low_q_offset; i < hi; ++i) { - double d = std::log(yi[i]) - std::log(model[i]); - double res = d / std::sqrt(std::fabs(std::log(eyi[i]))); - residuals[i] = std::isfinite(res) ? res : (res > 0 ? 1e6 : -1e6); - } - break; - case Type::InvRatioErr: - for (int i = low_q_offset; i < hi; ++i) { - double r = yi[i] / model[i]; - if (r < 1.0) r = 1.0 / r; - double emap = yi[i] / eyi[i]; - double res = (1.0 - r) * emap; - residuals[i] = std::isfinite(res) ? res : -1e6; - } - break; + case Type::LogDiff: + for (int i = low_q_offset; i < hi; ++i) { + double r = std::log(yi[i]) - std::log(model[i]); + residuals[i] = std::isfinite(r) ? r : (r > 0 ? 1e6 : -1e6); + } + break; + case Type::InvRatio: + for (int i = low_q_offset; i < hi; ++i) { + double r = yi[i] / model[i]; + if (r < 1.0) + r = 1.0 / r; + double res = 1.0 - r; + residuals[i] = std::isfinite(res) ? res : -1e6; + } + break; + case Type::LogDiffErr: + for (int i = low_q_offset; i < hi; ++i) { + double d = std::log(yi[i]) - std::log(model[i]); + double res = d / std::sqrt(std::fabs(std::log(eyi[i]))); + residuals[i] = std::isfinite(res) ? res : (res > 0 ? 1e6 : -1e6); + } + break; + case Type::InvRatioErr: + for (int i = low_q_offset; i < hi; ++i) { + double r = yi[i] / model[i]; + if (r < 1.0) + r = 1.0 / r; + double emap = yi[i] / eyi[i]; + double res = (1.0 - r) * emap; + residuals[i] = std::isfinite(res) ? res : -1e6; + } + break; } } -double ReflectivityObjective::Evaluate(std::span model, - std::span yi, - std::span eyi) const { +double ReflectivityObjective::Evaluate(std::span model, std::span yi, std::span eyi) const +{ const int n = static_cast(model.size()); double score = 0.0; switch (m_type) { - case Type::LogDiff: - for (int i = 0; i < n; ++i) { - double d = std::log(yi[i]) - std::log(model[i]); - score += d * d; - } - break; - case Type::InvRatio: - for (int i = 0; i < n; ++i) { - double r = yi[i] / model[i]; - if (r < 1.0) r = 1.0 / r; - score += (1.0 - r) * (1.0 - r); - } - break; - case Type::LogDiffErr: - for (int i = 0; i < n; ++i) { - double d = std::log(yi[i]) - std::log(model[i]); - score += d * d / std::fabs(std::log(eyi[i])); - } - break; - case Type::InvRatioErr: - for (int i = 0; i < n; ++i) { - double r = yi[i] / model[i]; - if (r < 1.0) r = 1.0 / r; - double emap = (yi[i] / eyi[i]) * (yi[i] / eyi[i]); - score += (1.0 - r) * (1.0 - r) * emap; - } - break; + case Type::LogDiff: + for (int i = 0; i < n; ++i) { + double d = std::log(yi[i]) - std::log(model[i]); + score += d * d; + } + break; + case Type::InvRatio: + for (int i = 0; i < n; ++i) { + double r = yi[i] / model[i]; + if (r < 1.0) + r = 1.0 / r; + score += (1.0 - r) * (1.0 - r); + } + break; + case Type::LogDiffErr: + for (int i = 0; i < n; ++i) { + double d = std::log(yi[i]) - std::log(model[i]); + score += d * d / std::fabs(std::log(eyi[i])); + } + break; + case Type::InvRatioErr: + for (int i = 0; i < n; ++i) { + double r = yi[i] / model[i]; + if (r < 1.0) + r = 1.0 / r; + double emap = (yi[i] / eyi[i]) * (yi[i] / eyi[i]); + score += (1.0 - r) * (1.0 - r) * emap; + } + break; } return score / (n + 1); } -double ComputeChiSquare(std::span model, - std::span yi, - std::span eyi) { +double ComputeChiSquare(std::span model, std::span yi, std::span eyi) +{ const int n = static_cast(model.size()); double chi = 0.0; for (int i = 0; i < n; ++i) { diff --git a/src/stochfitdll/SettingsStruct.cpp b/src/stochfitdll/SettingsStruct.cpp index f3aabaf..41976fd 100644 --- a/src/stochfitdll/SettingsStruct.cpp +++ b/src/stochfitdll/SettingsStruct.cpp @@ -1,66 +1,68 @@ -#include "platform.h" #include "SettingsStruct.h" + #include "generated/stochfit_generated.h" +#include "platform.h" static std::vector copy_fbs_vec(const flatbuffers::Vector* v) { - if (!v || v->size() == 0) return {}; + if (!v || v->size() == 0) + return {}; return std::vector(v->data(), v->data() + v->size()); } ReflSettings::ReflSettings(const StochFitProto::ReflSettings* s) { - Directory = s->directory() ? s->directory()->str() : ""; - Q = copy_fbs_vec(s->q()); - Refl = copy_fbs_vec(s->refl()); - ReflError = copy_fbs_vec(s->refl_error()); - QError = copy_fbs_vec(s->q_error()); - SubSLD = s->sub_sld(); - FilmSLD = s->film_sld(); - SupSLD = s->sup_sld(); - Boxes = s->boxes(); - FilmAbs = s->film_abs(); - SubAbs = s->sub_abs(); - SupAbs = s->sup_abs(); - Wavelength = s->wavelength(); - UseSurfAbs = s->use_surf_abs(); - QErr = s->q_err(); - Forcesig = s->forcesig(); - XRonly = s->xr_only(); - Resolution = s->resolution(); - FilmLength = s->film_length(); - Impnorm = s->impnorm(); - Objectivefunction = s->objectivefunction(); - Paramtemp = s->paramtemp(); - Sigmasearch = s->sigmasearch(); + Directory = s->directory() ? s->directory()->str() : ""; + Q = copy_fbs_vec(s->q()); + Refl = copy_fbs_vec(s->refl()); + ReflError = copy_fbs_vec(s->refl_error()); + QError = copy_fbs_vec(s->q_error()); + SubSLD = s->sub_sld(); + FilmSLD = s->film_sld(); + SupSLD = s->sup_sld(); + Boxes = s->boxes(); + FilmAbs = s->film_abs(); + SubAbs = s->sub_abs(); + SupAbs = s->sup_abs(); + Wavelength = s->wavelength(); + UseSurfAbs = s->use_surf_abs(); + QErr = s->q_err(); + Forcesig = s->forcesig(); + XRonly = s->xr_only(); + Resolution = s->resolution(); + FilmLength = s->film_length(); + Impnorm = s->impnorm(); + Objectivefunction = s->objectivefunction(); + Paramtemp = s->paramtemp(); + Sigmasearch = s->sigmasearch(); NormalizationSearchPerc = s->normalization_search_perc(); AbsorptionSearchPerc = s->absorption_search_perc(); - Algorithm = s->algorithm(); - Inittemp = s->inittemp(); - Platiter = s->platiter(); - Slope = s->slope(); - Gamma = s->gamma(); - STUNfunc = s->stun_func(); - Adaptive = s->adaptive(); - Tempiter = s->tempiter(); - STUNdeciter = s->stun_dec_iter(); - Gammadec = s->gammadec(); - CritEdgeOffset = s->crit_edge_offset(); - HighQOffset = s->high_q_offset(); - Iterations = s->iterations(); + Algorithm = s->algorithm(); + Inittemp = s->inittemp(); + Platiter = s->platiter(); + Slope = s->slope(); + Gamma = s->gamma(); + STUNfunc = s->stun_func(); + Adaptive = s->adaptive(); + Tempiter = s->tempiter(); + STUNdeciter = s->stun_dec_iter(); + Gammadec = s->gammadec(); + CritEdgeOffset = s->crit_edge_offset(); + HighQOffset = s->high_q_offset(); + Iterations = s->iterations(); } StochRunState::StochRunState(const StochFitProto::StochRunState* st) { - roughness = st->roughness(); - filmAbsInput = st->film_abs_input(); - surfAbs = st->surf_abs(); - temperature = st->temperature(); - impNorm = st->imp_norm(); - avgfSTUN = st->avg_f_stun(); - bestSolution = st->best_solution(); - chiSquare = st->chi_square(); + roughness = st->roughness(); + filmAbsInput = st->film_abs_input(); + surfAbs = st->surf_abs(); + temperature = st->temperature(); + impNorm = st->imp_norm(); + avgfSTUN = st->avg_f_stun(); + bestSolution = st->best_solution(); + chiSquare = st->chi_square(); goodnessOfFit = st->goodness_of_fit(); - iteration = st->iteration(); - edValues = copy_fbs_vec(st->ed_values()); + iteration = st->iteration(); + edValues = copy_fbs_vec(st->ed_values()); } diff --git a/src/stochfitdll/StochFitDll.cpp b/src/stochfitdll/StochFitDll.cpp index e32b4b8..4b79255 100644 --- a/src/stochfitdll/StochFitDll.cpp +++ b/src/stochfitdll/StochFitDll.cpp @@ -19,13 +19,15 @@ */ #include "StochFitDll.h" + +#include + +#include + #include "SettingsStruct.h" #include "StochFitHarness.h" #include "generated/stochfit_generated.h" #include "platform.h" -#include - -#include // ── Global state // ────────────────────────────────────────────────────────────── @@ -35,116 +37,120 @@ static std::string g_last_init_error; // ── Helper: finish a FlatBuffer and copy into caller-supplied buffer // ────────── Returns bytes written, or -1 if the buffer is too small. -template -static int finish_into(flatbuffers::FlatBufferBuilder &fbb, - flatbuffers::Offset root, uint8_t *outBuf, - int maxLen) { - fbb.Finish(root); - int written = static_cast(fbb.GetSize()); - if (written > maxLen) - return -1; - std::memcpy(outBuf, fbb.GetBufferPointer(), written); - return written; +template static int finish_into(flatbuffers::FlatBufferBuilder& fbb, flatbuffers::Offset root, uint8_t* outBuf, int maxLen) +{ + fbb.Finish(root); + int written = static_cast(fbb.GetSize()); + if (written > maxLen) + return -1; + std::memcpy(outBuf, fbb.GetBufferPointer(), written); + return written; } // ── Lifecycle // ───────────────────────────────────────────────────────────────── -extern "C" EXPORT void Init(const uint8_t *buf, int len) { - g_last_init_error.clear(); +extern "C" EXPORT void Init(const uint8_t* buf, int len) +{ + g_last_init_error.clear(); - const auto *req = StochFitProto::GetInitRequest(buf); - ReflSettings rs(req->settings()); + const auto* req = StochFitProto::GetInitRequest(buf); + ReflSettings rs(req->settings()); - std::unique_ptr state; - if (req->state() != nullptr) { - state = std::make_unique(req->state()); - } + std::unique_ptr state; + if (req->state() != nullptr) { + state = std::make_unique(req->state()); + } - stochfit = std::make_unique(rs, state); + stochfit = std::make_unique(rs, state); - auto result = stochfit->GetInitError(); - if (!result) { - g_last_init_error = result.error(); - stochfit.reset(); - } + auto result = stochfit->GetInitError(); + if (!result) { + g_last_init_error = result.error(); + stochfit.reset(); + } } -extern "C" EXPORT const char *GetInitError() { - return g_last_init_error.c_str(); +extern "C" EXPORT const char* GetInitError() +{ + return g_last_init_error.c_str(); } -extern "C" EXPORT void Start(int iterations) { - if (stochfit) { - stochfit->Start(iterations); - } +extern "C" EXPORT void Start(int iterations) +{ + if (stochfit) { + stochfit->Start(iterations); + } } -extern "C" EXPORT void Stop() { - if (stochfit) { - stochfit->Stop(); - } +extern "C" EXPORT void Stop() +{ + if (stochfit) { + stochfit->Stop(); + } } -extern "C" EXPORT void Destroy() { stochfit.reset(); } +extern "C" EXPORT void Destroy() +{ + stochfit.reset(); +} -extern "C" EXPORT void Cancel() { - if (stochfit) { - stochfit->Stop(); - } - stochfit.reset(); +extern "C" EXPORT void Cancel() +{ + if (stochfit) { + stochfit->Stop(); + } + stochfit.reset(); } // ── SA polling // ──────────────────────────────────────────────────────────────── -extern "C" EXPORT int GetData(uint8_t *outBuf, int maxLen) { - if (!stochfit) - return -1; +extern "C" EXPORT int GetData(uint8_t* outBuf, int maxLen) +{ + if (!stochfit) + return -1; - const DataSnapshot snap = stochfit->GetData(); + const DataSnapshot snap = stochfit->GetData(); - flatbuffers::FlatBufferBuilder fbb(static_cast(maxLen)); - auto zVec = fbb.CreateVector(snap.z); - auto rhoVec = fbb.CreateVector(snap.rho); - auto qVec = fbb.CreateVector(snap.Q); - auto rfVec = fbb.CreateVector(snap.refl); + flatbuffers::FlatBufferBuilder fbb(static_cast(maxLen)); + auto zVec = fbb.CreateVector(snap.z); + auto rhoVec = fbb.CreateVector(snap.rho); + auto qVec = fbb.CreateVector(snap.Q); + auto rfVec = fbb.CreateVector(snap.refl); - auto result = StochFitProto::CreateGetDataResult( - fbb, zVec, rhoVec, qVec, rfVec, - snap.roughness, snap.chiSquare, snap.goodnessOfFit, - snap.isFinished, snap.iteration); + auto result = StochFitProto::CreateGetDataResult(fbb, zVec, rhoVec, qVec, rfVec, snap.roughness, snap.chiSquare, snap.goodnessOfFit, + snap.isFinished, snap.iteration); - return finish_into(fbb, result, outBuf, maxLen); + return finish_into(fbb, result, outBuf, maxLen); } -extern "C" EXPORT int GetRunState(uint8_t *outBuf, int maxLen) { - if (!stochfit) - return -1; +extern "C" EXPORT int GetRunState(uint8_t* outBuf, int maxLen) +{ + if (!stochfit) + return -1; - const StochRunState s = stochfit->GetRunState(); - flatbuffers::FlatBufferBuilder fbb(4096); - auto edVec = fbb.CreateVector(s.edValues); - auto result = StochFitProto::CreateGetRunStateResult( - fbb, s.roughness, s.filmAbsInput, s.surfAbs, s.temperature, s.impNorm, - s.avgfSTUN, s.bestSolution, s.chiSquare, s.goodnessOfFit, edVec, - 0 /*iteration filled by caller*/); + const StochRunState s = stochfit->GetRunState(); + flatbuffers::FlatBufferBuilder fbb(4096); + auto edVec = fbb.CreateVector(s.edValues); + auto result = + StochFitProto::CreateGetRunStateResult(fbb, s.roughness, s.filmAbsInput, s.surfAbs, s.temperature, s.impNorm, s.avgfSTUN, + s.bestSolution, s.chiSquare, s.goodnessOfFit, edVec, 0 /*iteration filled by caller*/); - return finish_into(fbb, result, outBuf, maxLen); + return finish_into(fbb, result, outBuf, maxLen); } -extern "C" EXPORT int SAParams(uint8_t *outBuf, int maxLen) { - if (!stochfit) - return -1; +extern "C" EXPORT int SAParams(uint8_t* outBuf, int maxLen) +{ + if (!stochfit) + return -1; - int mode = 0; - double temp = stochfit->GetTemperature(); - double lowestEnergy = stochfit->GetLowestEnergy(); - mode = (temp < 1e-20) ? -1 : 1; + int mode = 0; + double temp = stochfit->GetTemperature(); + double lowestEnergy = stochfit->GetLowestEnergy(); + mode = (temp < 1e-20) ? -1 : 1; - flatbuffers::FlatBufferBuilder fbb(64); - auto result = - StochFitProto::CreateSaParamsResult(fbb, lowestEnergy, temp, mode); - return finish_into(fbb, result, outBuf, maxLen); + flatbuffers::FlatBufferBuilder fbb(64); + auto result = StochFitProto::CreateSaParamsResult(fbb, lowestEnergy, temp, mode); + return finish_into(fbb, result, outBuf, maxLen); } - diff --git a/src/stochfitdll/StochFitHarness.cpp b/src/stochfitdll/StochFitHarness.cpp index 128c999..5be4789 100644 --- a/src/stochfitdll/StochFitHarness.cpp +++ b/src/stochfitdll/StochFitHarness.cpp @@ -19,49 +19,56 @@ */ #include "StochFitHarness.h" + +#include + #include "ParamVector.h" #include "platform.h" -#include #ifdef _WIN32 -# define NOMINMAX -# define WIN32_LEAN_AND_MEAN -# include -# include -# include + #define NOMINMAX + #define WIN32_LEAN_AND_MEAN + #include + + #include + #include -static DWORD_PTR GetPCoreMask() { +static DWORD_PTR GetPCoreMask() +{ DWORD bufLen = 0; GetLogicalProcessorInformationEx(RelationProcessorCore, nullptr, &bufLen); std::vector buf(bufLen); auto* info = reinterpret_cast(buf.data()); if (!GetLogicalProcessorInformationEx(RelationProcessorCore, info, &bufLen)) - return ~(DWORD_PTR)0; + return ~(DWORD_PTR) 0; uint8_t maxClass = 0; for (auto* p = info; reinterpret_cast(p) < buf.data() + bufLen; - p = reinterpret_cast( - reinterpret_cast(p) + p->Size)) + p = reinterpret_cast(reinterpret_cast(p) + p->Size)) maxClass = std::max(maxClass, p->Processor.EfficiencyClass); DWORD_PTR mask = 0; bool hybrid = false; for (auto* p = info; reinterpret_cast(p) < buf.data() + bufLen; - p = reinterpret_cast( - reinterpret_cast(p) + p->Size)) { - if (p->Processor.EfficiencyClass < maxClass) hybrid = true; - else mask |= p->Processor.GroupMask[0].Mask; + p = reinterpret_cast(reinterpret_cast(p) + p->Size)) { + if (p->Processor.EfficiencyClass < maxClass) + hybrid = true; + else + mask |= p->Processor.GroupMask[0].Mask; } - return hybrid ? mask : ~(DWORD_PTR)0; + return hybrid ? mask : ~(DWORD_PTR) 0; } -static void PinOMPThreadsToPCores() { +static void PinOMPThreadsToPCores() +{ static std::once_flag pinned; std::call_once(pinned, [] { const DWORD_PTR mask = GetPCoreMask(); const int pCoreCount = static_cast(__popcnt64(mask)); - #pragma omp parallel - { SetThreadAffinityMask(GetCurrentThread(), mask); } + #pragma omp parallel + { + SetThreadAffinityMask(GetCurrentThread(), mask); + } omp_set_num_threads(pCoreCount); }); } @@ -75,321 +82,318 @@ static void PinOMPThreadsToPCores() { // built by ReflConstants — must see the same window as the sliced data, or the // model buffer and reflBuf sizes disagree. Invalid inputs are returned // unchanged; the constructor body reports them via m_initError. -static ReflSettings SliceToFitWindow(ReflSettings s) { - const int qsize = static_cast(s.Q.size()); - const int n = qsize - s.CritEdgeOffset - s.HighQOffset; - if (s.CritEdgeOffset < 0 || s.HighQOffset < 0 || n <= 0 || - static_cast(s.Refl.size()) < qsize || - static_cast(s.ReflError.size()) < qsize) +static ReflSettings SliceToFitWindow(ReflSettings s) +{ + const int qsize = static_cast(s.Q.size()); + const int n = qsize - s.CritEdgeOffset - s.HighQOffset; + if (s.CritEdgeOffset < 0 || s.HighQOffset < 0 || n <= 0 || static_cast(s.Refl.size()) < qsize || + static_cast(s.ReflError.size()) < qsize) + return s; + + const int off = s.CritEdgeOffset; + auto slice = [&](std::vector& v) { + if (static_cast(v.size()) >= off + n) + v.assign(v.begin() + off, v.begin() + off + n); + }; + slice(s.Q); + slice(s.Refl); + slice(s.ReflError); + slice(s.QError); // may be empty (no smearing); left alone in that case + s.CritEdgeOffset = 0; + s.HighQOffset = 0; return s; - - const int off = s.CritEdgeOffset; - auto slice = [&](std::vector &v) { - if (static_cast(v.size()) >= off + n) - v.assign(v.begin() + off, v.begin() + off + n); - }; - slice(s.Q); - slice(s.Refl); - slice(s.ReflError); - slice(s.QError); // may be empty (no smearing); left alone in that case - s.CritEdgeOffset = 0; - s.HighQOffset = 0; - return s; } -StochFit::StochFit(const ReflSettings &InitStruct, - const std::unique_ptr &state) +StochFit::StochFit(const ReflSettings& InitStruct, const std::unique_ptr& state) // m_initStruct first; m_cEDP/m_displayEDP before m_parratt (declaration order). : m_initStruct(SliceToFitWindow(InitStruct)), - m_cEDP(m_initStruct), m_displayEDP(m_initStruct), - params(InitStruct), m_displayState(InitStruct), + m_cEDP(m_initStruct), + m_displayEDP(m_initStruct), + params(InitStruct), + m_displayState(InitStruct), m_parratt(m_initStruct, m_cEDP.GetLayerCount()), m_objective(ReflectivityObjective::Type{InitStruct.Objectivefunction}), m_stepper({.sigmaSearch = InitStruct.Sigmasearch, .absSearch = InitStruct.AbsorptionSearchPerc, .normSearch = InitStruct.NormalizationSearchPerc, - .stepSize = InitStruct.Paramtemp}) { - m_stop_requested = false; - - m_Directory = InitStruct.Directory; - - // m_initStruct was sliced to the fit window by SliceToFitWindow; validate - // against the caller's original settings so invalid inputs are rejected. - const int qsize = static_cast(InitStruct.Q.size()); - m_datapoints = qsize - InitStruct.HighQOffset - InitStruct.CritEdgeOffset; - if (InitStruct.CritEdgeOffset < 0 || InitStruct.HighQOffset < 0 || - m_datapoints <= 0 || - static_cast(InitStruct.Refl.size()) < qsize || - static_cast(InitStruct.ReflError.size()) < qsize) { - m_initError = tl::unexpected(std::string("Invalid Q/Refl/ReflError sizes")); - return; - } - m_xi = m_initStruct.Q; - m_yi = m_initStruct.Refl; - m_eyi = m_initStruct.ReflError; - - - // Apply run state if provided. - // filmAbsInput is the pre-multiplication value: Set_FilmAbs(x) stores x*WC. - // temperature is raw m_dTemp (β) stored directly via SetTemperature. - // surfAbs is saved independently so it is never baked into filmAbsInput. - if (state && - static_cast(state->edValues.size()) == params.RealParamsSize()) { - params.SetRoughness(state->roughness); - params.SetSupphase(state->edValues[0]); - for (int i = 1; i < static_cast(state->edValues.size()) - 1; i++) - params.SetMutatableParameter(i - 1, state->edValues[i]); - params.SetSubphase(state->edValues[state->edValues.size() - 1]); - m_cEDP.Set_FilmAbs(state->filmAbsInput); - params.SetSurfAbs(state->surfAbs); - params.SetImpNorm(state->impNorm); - } - params.UpdateBoundaries(); - m_displayState.params = params; - - // Build SA scratch buffer and deps. - const int nd = m_datapoints; - m_saReflBuf.resize(nd); - AnnealDeps deps; - deps.yi = m_yi; - deps.eyi = m_eyi; - deps.reflBuf = m_saReflBuf; - deps.impNorm = InitStruct.Impnorm; - - // Construct the algorithm-specific annealer variant. - const auto algo = AlgorithmFromInt(InitStruct.Algorithm); - switch (algo) { - case SaAlgorithm::Greedy: - m_annealer.emplace(std::in_place_type>, m_cEDP, - m_parratt, m_objective, m_stepper, params, deps); - break; - case SaAlgorithm::Simulated: - m_annealer.emplace(std::in_place_type>, m_cEDP, - m_parratt, m_objective, m_stepper, params, deps, - InitStruct.Inittemp, InitStruct.Slope, - InitStruct.Platiter); - break; - case SaAlgorithm::Stun: - m_annealer.emplace( - std::in_place_type>, m_cEDP, m_parratt, m_objective, - m_stepper, params, deps, InitStruct.Inittemp, InitStruct.Slope, - InitStruct.Platiter, InitStruct.Gamma, - InitStruct.STUNfunc, InitStruct.Tempiter, InitStruct.Adaptive); - break; - } - - // Apply session temperature/avgfSTUN after the annealer is constructed. - if (state != nullptr && - static_cast(state->edValues.size()) == params.RealParamsSize()) { - std::visit( - [&](auto &a) { - a.SetTemperature(state->temperature); - a.SetAverageFSTUN(state->avgfSTUN); - }, - *m_annealer); - } - - // Compute the initial energy to seed best/current state. - std::visit([&](auto &a) { a.InitEnergy(params); }, *m_annealer); - - m_initError = {}; + .stepSize = InitStruct.Paramtemp}) +{ + m_stop_requested = false; + + m_Directory = InitStruct.Directory; + + // m_initStruct was sliced to the fit window by SliceToFitWindow; validate + // against the caller's original settings so invalid inputs are rejected. + const int qsize = static_cast(InitStruct.Q.size()); + m_datapoints = qsize - InitStruct.HighQOffset - InitStruct.CritEdgeOffset; + if (InitStruct.CritEdgeOffset < 0 || InitStruct.HighQOffset < 0 || m_datapoints <= 0 || + static_cast(InitStruct.Refl.size()) < qsize || static_cast(InitStruct.ReflError.size()) < qsize) { + m_initError = tl::unexpected(std::string("Invalid Q/Refl/ReflError sizes")); + return; + } + m_xi = m_initStruct.Q; + m_yi = m_initStruct.Refl; + m_eyi = m_initStruct.ReflError; + + // Apply run state if provided. + // filmAbsInput is the pre-multiplication value: Set_FilmAbs(x) stores x*WC. + // temperature is raw m_dTemp (β) stored directly via SetTemperature. + // surfAbs is saved independently so it is never baked into filmAbsInput. + if (state && static_cast(state->edValues.size()) == params.RealParamsSize()) { + params.SetRoughness(state->roughness); + params.SetSupphase(state->edValues[0]); + for (int i = 1; i < static_cast(state->edValues.size()) - 1; i++) + params.SetMutatableParameter(i - 1, state->edValues[i]); + params.SetSubphase(state->edValues[state->edValues.size() - 1]); + m_cEDP.Set_FilmAbs(state->filmAbsInput); + params.SetSurfAbs(state->surfAbs); + params.SetImpNorm(state->impNorm); + } + params.UpdateBoundaries(); + m_displayState.params = params; + + // Build SA scratch buffer and deps. + const int nd = m_datapoints; + m_saReflBuf.resize(nd); + AnnealDeps deps; + deps.yi = m_yi; + deps.eyi = m_eyi; + deps.reflBuf = m_saReflBuf; + deps.impNorm = InitStruct.Impnorm; + + // Construct the algorithm-specific annealer variant. + const auto algo = AlgorithmFromInt(InitStruct.Algorithm); + switch (algo) { + case SaAlgorithm::Greedy: + m_annealer.emplace(std::in_place_type>, m_cEDP, m_parratt, m_objective, m_stepper, params, deps); + break; + case SaAlgorithm::Simulated: + m_annealer.emplace(std::in_place_type>, m_cEDP, m_parratt, m_objective, m_stepper, params, deps, + InitStruct.Inittemp, InitStruct.Slope, InitStruct.Platiter); + break; + case SaAlgorithm::Stun: + m_annealer.emplace(std::in_place_type>, m_cEDP, m_parratt, m_objective, m_stepper, params, deps, + InitStruct.Inittemp, InitStruct.Slope, InitStruct.Platiter, InitStruct.Gamma, InitStruct.STUNfunc, + InitStruct.Tempiter, InitStruct.Adaptive); + break; + } + + // Apply session temperature/avgfSTUN after the annealer is constructed. + if (state != nullptr && static_cast(state->edValues.size()) == params.RealParamsSize()) { + std::visit( + [&](auto& a) { + a.SetTemperature(state->temperature); + a.SetAverageFSTUN(state->avgfSTUN); + }, + *m_annealer); + } + + // Compute the initial energy to seed best/current state. + std::visit([&](auto& a) { a.InitEnergy(params); }, *m_annealer); + + m_initError = {}; } -StochFit::~StochFit() { - if (m_thread.joinable()) { - m_stop_requested = true; - m_thread.join(); - } +StochFit::~StochFit() +{ + if (m_thread.joinable()) { + m_stop_requested = true; + m_thread.join(); + } } // ── Harness accessors // ───────────────────────────────────────────────────────── -double StochFit::GetTemperature() const { - if (!m_annealer) - return 0.0; - return std::visit([](const auto &a) { return a.GetTemperature(); }, - *m_annealer); +double StochFit::GetTemperature() const +{ + if (!m_annealer) + return 0.0; + return std::visit([](const auto& a) { return a.GetTemperature(); }, *m_annealer); } -double StochFit::GetRawTemperature() const { - if (!m_annealer) - return 0.0; - return std::visit([](const auto &a) { return a.GetRawTemperature(); }, - *m_annealer); +double StochFit::GetRawTemperature() const +{ + if (!m_annealer) + return 0.0; + return std::visit([](const auto& a) { return a.GetRawTemperature(); }, *m_annealer); } -void StochFit::SetTemperature(double t) { - if (m_annealer) - std::visit([t](auto &a) { a.SetTemperature(t); }, *m_annealer); +void StochFit::SetTemperature(double t) +{ + if (m_annealer) + std::visit([t](auto& a) { a.SetTemperature(t); }, *m_annealer); } -double StochFit::GetLowestEnergy() const { - if (!m_annealer) - return 0.0; - return std::visit([](const auto &a) { return a.GetLowestEnergy(); }, - *m_annealer); +double StochFit::GetLowestEnergy() const +{ + if (!m_annealer) + return 0.0; + return std::visit([](const auto& a) { return a.GetLowestEnergy(); }, *m_annealer); } -double StochFit::GetAverageFSTUN() const { - if (!m_annealer) - return 0.0; - return std::visit([](const auto &a) { return a.GetAverageFSTUN(); }, - *m_annealer); +double StochFit::GetAverageFSTUN() const +{ + if (!m_annealer) + return 0.0; + return std::visit([](const auto& a) { return a.GetAverageFSTUN(); }, *m_annealer); } -void StochFit::SetAverageFSTUN(double f) { - if (m_annealer) - std::visit([f](auto &a) { a.SetAverageFSTUN(f); }, *m_annealer); +void StochFit::SetAverageFSTUN(double f) +{ + if (m_annealer) + std::visit([f](auto& a) { a.SetAverageFSTUN(f); }, *m_annealer); } // ── Main CPU loop // ───────────────────────────────────────────────────────────── -int StochFit::Processing() { - try { - if (!m_annealer) - return -1; +int StochFit::Processing() +{ + try { + if (!m_annealer) + return -1; #ifdef _WIN32 - PinOMPThreadsToPCores(); + PinOMPThreadsToPCores(); #endif - const int nThreads = omp_get_max_threads(); + const int nThreads = omp_get_max_threads(); #pragma omp parallel num_threads(nThreads) - { - for (int isteps = 0; - isteps < m_itotaliterations && !m_stop_requested.load(); ++isteps) { + { + for (int isteps = 0; isteps < m_itotaliterations && !m_stop_requested.load(); ++isteps) { - // All threads: mutate candidate (omp single inside) + build EDP (omp for inside). - std::visit([&](auto &a) { a.PrepareCandidate(params); }, *m_annealer); + // All threads: mutate candidate (omp single inside) + build EDP (omp for inside). + std::visit([&](auto& a) { a.PrepareCandidate(params); }, *m_annealer); - // All threads: cooperative Parratt Q-point distribution (omp for inside). - std::visit([&](auto &a) { a.ComputeSharedRefl(); }, *m_annealer); + // All threads: cooperative Parratt Q-point distribution (omp for inside). + std::visit([&](auto& a) { a.ComputeSharedRefl(); }, *m_annealer); - // Serial: accept/reject + display snapshot. + // Serial: accept/reject + display snapshot. #pragma omp single - { - bool accepted = false; - double gof = 0.0, chi = 0.0; - - accepted = std::visit( - [&](auto &a) { return a.EvaluateAndAccept(params); }, *m_annealer); - if (accepted || isteps == 0) { - gof = std::visit([](const auto &a) { return a.GetCurrentEnergy(); }, - *m_annealer); - chi = std::visit([](const auto &a) { return a.GetLastChiSquare(); }, - *m_annealer); - } - - if (accepted || isteps == 0) { - std::lock_guard lock(m_displayMutex); - m_displayState.params = params; - m_displayState.refl.assign(m_saReflBuf.begin(), m_saReflBuf.end()); - m_displayState.chiSquare = chi; - m_displayState.goF = gof; - } - m_icurrentiteration.store(isteps + 1, std::memory_order_relaxed); + { + bool accepted = false; + double gof = 0.0, chi = 0.0; + + accepted = std::visit([&](auto& a) { return a.EvaluateAndAccept(params); }, *m_annealer); + if (accepted || isteps == 0) { + gof = std::visit([](const auto& a) { return a.GetCurrentEnergy(); }, *m_annealer); + chi = std::visit([](const auto& a) { return a.GetLastChiSquare(); }, *m_annealer); + } + + if (accepted || isteps == 0) { + std::lock_guard lock(m_displayMutex); + m_displayState.params = params; + m_displayState.refl.assign(m_saReflBuf.begin(), m_saReflBuf.end()); + m_displayState.chiSquare = chi; + m_displayState.goF = gof; + } + m_icurrentiteration.store(isteps + 1, std::memory_order_relaxed); + } + // implicit barrier: all threads sync before the next iteration + } } - // implicit barrier: all threads sync before the next iteration - } - } - m_icurrentiteration.store(m_itotaliterations, std::memory_order_relaxed); - return 0; - } catch (const std::exception &ex) { - std::cerr << "[StochFit] Processing() caught exception: " << ex.what() - << std::endl; - return -1; - } catch (...) { - std::cerr << "[StochFit] Processing() caught unknown exception" << std::endl; - return -1; - } + m_icurrentiteration.store(m_itotaliterations, std::memory_order_relaxed); + return 0; + } catch (const std::exception& ex) { + std::cerr << "[StochFit] Processing() caught exception: " << ex.what() << std::endl; + return -1; + } catch (...) { + std::cerr << "[StochFit] Processing() caught unknown exception" << std::endl; + return -1; + } } // ── Display snapshot (main thread only) // ────────────────────────────────────── -DataSnapshot StochFit::GetCurrentState() { - // Copy the minimum inside the lock, then compute outside. - DisplayState snap(m_initStruct); - { - std::lock_guard lock(m_displayMutex); - snap = m_displayState; - } - - m_displayEDP.GenerateEDP(snap.params); - - DataSnapshot out; - out.roughness = snap.params.GetRoughness(); - out.chiSquare = snap.chiSquare; - out.goodnessOfFit = snap.goF; - - // snap.refl is the SA scoring buffer: one value per measured Q point. - const int nd = m_datapoints; - out.Q.assign(m_xi.begin(), m_xi.begin() + nd); - out.refl = std::move(snap.refl); - out.refl.resize(nd); - - const int nl = m_displayEDP.Get_EDPPointCount(); - out.z.resize(nl); - out.rho.resize(nl); - const double dz = m_displayEDP.Get_Dz(); - const double leftOffset = m_displayEDP.Get_LeftOffset(); - const double rhoSub = m_displayEDP.m_EDP[nl - 1].real(); - for (int i = 0; i < nl; i++) { - out.z[i] = i * dz - leftOffset; - out.rho[i] = m_displayEDP.m_EDP[i].real() / rhoSub; - } - - return out; +DataSnapshot StochFit::GetCurrentState() +{ + // Copy the minimum inside the lock, then compute outside. + DisplayState snap(m_initStruct); + { + std::lock_guard lock(m_displayMutex); + snap = m_displayState; + } + + m_displayEDP.GenerateEDP(snap.params); + + DataSnapshot out; + out.roughness = snap.params.GetRoughness(); + out.chiSquare = snap.chiSquare; + out.goodnessOfFit = snap.goF; + + // snap.refl is the SA scoring buffer: one value per measured Q point. + const int nd = m_datapoints; + out.Q.assign(m_xi.begin(), m_xi.begin() + nd); + out.refl = std::move(snap.refl); + out.refl.resize(nd); + + const int nl = m_displayEDP.Get_EDPPointCount(); + out.z.resize(nl); + out.rho.resize(nl); + const double dz = m_displayEDP.Get_Dz(); + const double leftOffset = m_displayEDP.Get_LeftOffset(); + const double rhoSub = m_displayEDP.m_EDP[nl - 1].real(); + for (int i = 0; i < nl; i++) { + out.z[i] = i * dz - leftOffset; + out.rho[i] = m_displayEDP.m_EDP[i].real() / rhoSub; + } + + return out; } // ── Public API // ──────────────────────────────────────────────────────────────── -int StochFit::Start(int iterations) { - m_itotaliterations = iterations; - m_stop_requested = false; - m_thread = std::thread([this] { Processing(); }); - return 0; +int StochFit::Start(int iterations) +{ + m_itotaliterations = iterations; + m_stop_requested = false; + m_thread = std::thread([this] { Processing(); }); + return 0; } -int StochFit::Cancel() { - if (m_thread.joinable()) - m_stop_requested = true; - return 0; +int StochFit::Cancel() +{ + if (m_thread.joinable()) + m_stop_requested = true; + return 0; } -void StochFit::Stop() { - if (m_thread.joinable()) { - m_stop_requested = true; - m_thread.join(); - } +void StochFit::Stop() +{ + if (m_thread.joinable()) { + m_stop_requested = true; + m_thread.join(); + } } -StochRunState StochFit::GetRunState() { - // Called after Stop() — thread is joined, no lock needed. - StochRunState s; - s.roughness = params.GetRoughness(); - s.filmAbsInput = m_cEDP.Get_FilmAbsInput(); - s.surfAbs = params.GetSurfAbs(); - s.temperature = GetRawTemperature(); - s.impNorm = params.GetImpNorm(); - s.avgfSTUN = GetAverageFSTUN(); - s.bestSolution = GetLowestEnergy(); - s.chiSquare = m_displayState.chiSquare; - s.goodnessOfFit = m_displayState.goF; - auto span = params.RealParams(); - s.edValues.assign(span.begin(), span.end()); - return s; +StochRunState StochFit::GetRunState() +{ + // Called after Stop() — thread is joined, no lock needed. + StochRunState s; + s.roughness = params.GetRoughness(); + s.filmAbsInput = m_cEDP.Get_FilmAbsInput(); + s.surfAbs = params.GetSurfAbs(); + s.temperature = GetRawTemperature(); + s.impNorm = params.GetImpNorm(); + s.avgfSTUN = GetAverageFSTUN(); + s.bestSolution = GetLowestEnergy(); + s.chiSquare = m_displayState.chiSquare; + s.goodnessOfFit = m_displayState.goF; + auto span = params.RealParams(); + s.edValues.assign(span.begin(), span.end()); + return s; } -DataSnapshot StochFit::GetData() { - const int iter = m_icurrentiteration.load(std::memory_order_relaxed); - DataSnapshot snap = GetCurrentState(); - snap.iteration = iter; - snap.isFinished = (iter >= m_itotaliterations); - return snap; +DataSnapshot StochFit::GetData() +{ + const int iter = m_icurrentiteration.load(std::memory_order_relaxed); + DataSnapshot snap = GetCurrentState(); + snap.iteration = iter; + snap.isFinished = (iter >= m_itotaliterations); + return snap; } diff --git a/src/stochfitdll/UnifiedReflectivity.cpp b/src/stochfitdll/UnifiedReflectivity.cpp index 194d223..cff6e9c 100644 --- a/src/stochfitdll/UnifiedReflectivity.cpp +++ b/src/stochfitdll/UnifiedReflectivity.cpp @@ -1,12 +1,12 @@ #include "UnifiedReflectivity.h" + #include namespace { -void InitScratchArrays(WaveScratch> &complex_scratch, - WaveScratch &real_scratch, int num_threads, - int n_layers) { - complex_scratch.resize(n_layers, num_threads); - real_scratch.resize(n_layers, num_threads); +void InitScratchArrays(WaveScratch>& complex_scratch, WaveScratch& real_scratch, int num_threads, int n_layers) +{ + complex_scratch.resize(n_layers, num_threads); + real_scratch.resize(n_layers, num_threads); } // Returns the first Q index where no EDP layer is evanescent. Exploits the fact @@ -14,129 +14,120 @@ void InitScratchArrays(WaveScratch> &complex_scratch, // single max-density precompute + binary search replaces the original O(Q·EDP) // scan. Only valid for the measurement Q grid (monotone); never call for // qspread arrays. -auto FindComplexToRealOffset( - std::span sinsquaredthetai, - std::span> density_profile, - double indexsupsquared_real, double sup_sld) -> int { - const double max_density = std::ranges::max( - density_profile | - std::views::transform([](const auto &dp) { return dp.real(); })); - - const auto it = std::ranges::lower_bound( - sinsquaredthetai, max_density, [&](double sin2, double target) { - return indexsupsquared_real * sin2 + sup_sld < target; - }); - return static_cast(it - sinsquaredthetai.begin()); +auto FindComplexToRealOffset(std::span sinsquaredthetai, std::span> density_profile, + double indexsupsquared_real, double sup_sld) -> int +{ + const double max_density = std::ranges::max(density_profile | std::views::transform([](const auto& dp) { return dp.real(); })); + + const auto it = std::ranges::lower_bound(sinsquaredthetai, max_density, + [&](double sin2, double target) { return indexsupsquared_real * sin2 + sup_sld < target; }); + return static_cast(it - sinsquaredthetai.begin()); } -} // namespace - -ReflConstants::ReflConstants(const ReflSettings &s) { - k0 = 2.0 * std::numbers::pi / s.Wavelength; - sup_sld = s.SupSLD; - indexsup = 1.0 - s.SupSLD / 2.0; - indexsupsquared = indexsup * indexsup; - - const double qspread_frac = s.QErr / 100.0; - qsmear_enabled = qspread_frac >= 0.005 && !s.QError.empty(); - - const int meas_n = static_cast(s.Q.size()); - const int grid_n = qsmear_enabled ? meas_n * QSmear::Points : meas_n; - sinthetai.resize(grid_n); - sinsquaredthetai.resize(grid_n); - - if (qsmear_enabled) { - std::vector meas_sin(meas_n); - for (int i = 0; i < meas_n; ++i) - meas_sin[i] = s.Q[i] * s.Wavelength / (4.0 * std::numbers::pi); - QSmear::BuildArrays(s.Wavelength, qspread_frac, meas_sin, s.QError, - sinthetai, sinsquaredthetai); - } else { - for (int i = 0; i < meas_n; ++i) { - sinthetai[i] = s.Q[i] * s.Wavelength / (4.0 * std::numbers::pi); - sinsquaredthetai[i] = sinthetai[i] * sinthetai[i]; +} // namespace + +ReflConstants::ReflConstants(const ReflSettings& s) +{ + k0 = 2.0 * std::numbers::pi / s.Wavelength; + sup_sld = s.SupSLD; + indexsup = 1.0 - s.SupSLD / 2.0; + indexsupsquared = indexsup * indexsup; + + const double qspread_frac = s.QErr / 100.0; + qsmear_enabled = qspread_frac >= 0.005 && !s.QError.empty(); + + const int meas_n = static_cast(s.Q.size()); + const int grid_n = qsmear_enabled ? meas_n * QSmear::Points : meas_n; + sinthetai.resize(grid_n); + sinsquaredthetai.resize(grid_n); + + if (qsmear_enabled) { + std::vector meas_sin(meas_n); + for (int i = 0; i < meas_n; ++i) + meas_sin[i] = s.Q[i] * s.Wavelength / (4.0 * std::numbers::pi); + QSmear::BuildArrays(s.Wavelength, qspread_frac, meas_sin, s.QError, sinthetai, sinsquaredthetai); + } else { + for (int i = 0; i < meas_n; ++i) { + sinthetai[i] = s.Q[i] * s.Wavelength / (4.0 * std::numbers::pi); + sinsquaredthetai[i] = sinthetai[i] * sinthetai[i]; + } } - } } -ParrattReflectivity::ParrattReflectivity(const ReflSettings &settings) - : m_consts(settings), - m_qsmear_enabled(m_consts.qsmear_enabled) { +ParrattReflectivity::ParrattReflectivity(const ReflSettings& settings) : m_consts(settings), m_qsmear_enabled(m_consts.qsmear_enabled) +{ - const int meas_n = static_cast(settings.Q.size()); - m_refl_out.resize(m_consts.sinthetai.size()); - if (m_qsmear_enabled) - m_refl_smeared.resize(meas_n); + const int meas_n = static_cast(settings.Q.size()); + m_refl_out.resize(m_consts.sinthetai.size()); + if (m_qsmear_enabled) + m_refl_smeared.resize(meas_n); } -ParrattReflectivity::ParrattReflectivity(const ReflSettings &settings, - int n_layers) - : ParrattReflectivity(settings) { - InitScratchArrays(m_complex, m_real, omp_get_max_threads(), n_layers); +ParrattReflectivity::ParrattReflectivity(const ReflSettings& settings, int n_layers) : ParrattReflectivity(settings) +{ + InitScratchArrays(m_complex, m_real, omp_get_max_threads(), n_layers); } // Standalone: creates its own OMP team. Safe for levmar/InitEnergy/any serial context. -auto ParrattReflectivity::CalculateReflectivity(const LayerStack &ls) - -> std::span { - const int n_layers = static_cast(ls.rho.size()); - InitScratchArrays(m_complex, m_real, omp_get_max_threads(), n_layers); +auto ParrattReflectivity::CalculateReflectivity(const LayerStack& ls) -> std::span +{ + const int n_layers = static_cast(ls.rho.size()); + InitScratchArrays(m_complex, m_real, omp_get_max_threads(), n_layers); #pragma omp parallel - { - if (m_qsmear_enabled) { - if (ls.has_roughness) - ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, - m_consts.sinsquaredthetai, m_refl_out); - else - ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, - m_consts.sinsquaredthetai, m_refl_out); - } else if (ls.transparent && !ls.has_roughness) { - TransparentReflectivityCalc(ls); - } else { - ReflectivityCalc(ls); + { + if (m_qsmear_enabled) { + if (ls.has_roughness) + ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, m_consts.sinsquaredthetai, m_refl_out); + else + ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, m_consts.sinsquaredthetai, m_refl_out); + } else if (ls.transparent && !ls.has_roughness) { + TransparentReflectivityCalc(ls); + } else { + ReflectivityCalc(ls); + } } - } - if (m_qsmear_enabled) { - QSmear::Apply(m_refl_out, m_refl_smeared); - return m_refl_smeared; - } - return m_refl_out; + if (m_qsmear_enabled) { + QSmear::Apply(m_refl_out, m_refl_smeared); + return m_refl_smeared; + } + return m_refl_out; } // Cooperative: must be called by ALL threads of an enclosing OMP parallel region. // BuildLayerStack and scratch sizing happen in omp single; Q-point work is shared. -auto ParrattReflectivity::CalculateReflectivityCooperative(const CEDP &EDP) - -> std::span { +auto ParrattReflectivity::CalculateReflectivityCooperative(const CEDP& EDP) -> std::span +{ #pragma omp single - { - // BuildLayerStackFull computes GetOffSets inline — merges 2 barriers into 1. - m_cooperative_ls = EDP.BuildLayerStackFull(); - } - // implicit barrier: all threads see m_cooperative_ls - - if (m_qsmear_enabled) { - if (m_cooperative_ls.has_roughness) - ReflectivityCalcCoreImpl(m_cooperative_ls, m_consts.sinthetai, - m_consts.sinsquaredthetai, m_refl_out); - else - ReflectivityCalcCoreImpl(m_cooperative_ls, m_consts.sinthetai, - m_consts.sinsquaredthetai, m_refl_out); + { + // BuildLayerStackFull computes GetOffSets inline — merges 2 barriers into 1. + m_cooperative_ls = EDP.BuildLayerStackFull(); + } + // implicit barrier: all threads see m_cooperative_ls + + if (m_qsmear_enabled) { + if (m_cooperative_ls.has_roughness) + ReflectivityCalcCoreImpl(m_cooperative_ls, m_consts.sinthetai, m_consts.sinsquaredthetai, m_refl_out); + else + ReflectivityCalcCoreImpl(m_cooperative_ls, m_consts.sinthetai, m_consts.sinsquaredthetai, m_refl_out); #pragma omp single - { QSmear::Apply(m_refl_out, m_refl_smeared); } - return m_refl_smeared; - } + { + QSmear::Apply(m_refl_out, m_refl_smeared); + } + return m_refl_smeared; + } - if (m_cooperative_ls.transparent && !m_cooperative_ls.has_roughness) - TransparentReflectivityCalc(m_cooperative_ls); - else - ReflectivityCalc(m_cooperative_ls); + if (m_cooperative_ls.transparent && !m_cooperative_ls.has_roughness) + TransparentReflectivityCalc(m_cooperative_ls); + else + ReflectivityCalc(m_cooperative_ls); - return m_refl_out; + return m_refl_out; } -auto ParrattReflectivity::CalculateReflectivity(const CEDP &EDP) - -> std::span { - return CalculateReflectivity(EDP.BuildLayerStack()); +auto ParrattReflectivity::CalculateReflectivity(const CEDP& EDP) -> std::span +{ + return CalculateReflectivity(EDP.BuildLayerStack()); } // Complex Parratt loop over Q indices [0, q_end). @@ -154,167 +145,152 @@ auto ParrattReflectivity::CalculateReflectivity(const CEDP &EDP) // - HasRoughness=true: Nevot-Croce factor exp(sigma_sq[i]*kk[i]*kk[i+1]) // is multiplied into each Fresnel rj[i]. When false, the compiler omits it. template -void ParrattReflectivity::ReflectivityCalcCoreImpl( - const LayerStack &ls, std::span sinthetai, - std::span sin2thetai, std::span out, int q_end) { - const auto &density_profile = ls.rho; - const int n_layers = static_cast(density_profile.size()); - const int n_q = std::min(q_end, static_cast(sinthetai.size())); +void ParrattReflectivity::ReflectivityCalcCoreImpl(const LayerStack& ls, std::span sinthetai, + std::span sin2thetai, std::span out, int q_end) +{ + const auto& density_profile = ls.rho; + const int n_layers = static_cast(density_profile.size()); + const int n_q = std::min(q_end, static_cast(sinthetai.size())); - const int sup_offset = ls.sup_offset; - const int sub_offset = ls.sub_offset; + const int sup_offset = ls.sup_offset; + const int sub_offset = ls.sub_offset; - // Per-thread scratch setup — each thread in the enclosing parallel team computes - // its own slice. No inner omp parallel needed; caller provides the team. - const int threadnum = omp_get_thread_num(); - const int arrayoffset = threadnum * n_layers; + // Per-thread scratch setup — each thread in the enclosing parallel team computes + // its own slice. No inner omp parallel needed; caller provides the team. + const int threadnum = omp_get_thread_num(); + const int arrayoffset = threadnum * n_layers; - auto *kk_slice = m_complex.kk.data() + arrayoffset; - auto *ak_slice = m_complex.ak.data() + arrayoffset; - auto *rj_slice = m_complex.rj.data() + arrayoffset; - auto *Rj_slice = m_complex.Rj.data() + arrayoffset; + auto* kk_slice = m_complex.kk.data() + arrayoffset; + auto* ak_slice = m_complex.ak.data() + arrayoffset; + auto* rj_slice = m_complex.rj.data() + arrayoffset; + auto* Rj_slice = m_complex.Rj.data() + arrayoffset; - ak_slice[0] = 1.0; - Rj_slice[sub_offset + 1] = 0.0; + ak_slice[0] = 1.0; + Rj_slice[sub_offset + 1] = 0.0; #pragma omp for schedule(runtime) - for (int l = 0; l < n_q; l++) { - kk_slice[0] = m_consts.k0 * m_consts.indexsup * sinthetai[l]; - - // Flat-front wavevector — density uniform at [1, sup_offset-1] - const auto kk1 = - m_consts.k0 * std::sqrt(m_consts.indexsupsquared * sin2thetai[l] - - density_profile[1] + m_consts.sup_sld); - const auto ak1 = std::exp(ls.length_mult[1] * kk1); - - // Flat-back wavevector - const auto kk_sub = - m_consts.k0 * - std::sqrt(m_consts.indexsupsquared * sin2thetai[l] - - density_profile[n_layers - 1] + m_consts.sup_sld); - - // Varying region: kk and ak for [sup_offset, sub_offset] - for (int i = sup_offset; i <= sub_offset; i++) { - kk_slice[i] = - m_consts.k0 * std::sqrt(m_consts.indexsupsquared * sin2thetai[l] - - density_profile[i] + m_consts.sup_sld); - ak_slice[i] = std::exp(ls.length_mult[i] * kk_slice[i]); - } - - // Boundary entries for rj at the flat/varying interfaces. - const int rj_start = sup_offset >= 2 ? sup_offset - 1 : 0; - if (sup_offset >= 2) { - kk_slice[sup_offset - 1] = kk1; - ak_slice[sup_offset - 1] = ak1; - } - kk_slice[sub_offset + 1] = kk_sub; - - for (int i = rj_start; i <= sub_offset; i++) { - rj_slice[i] = - (kk_slice[i] - kk_slice[i + 1]) / (kk_slice[i] + kk_slice[i + 1]); - if constexpr (HasRoughness) - rj_slice[i] *= std::exp(ls.sigma_sq[i] * kk_slice[i] * kk_slice[i + 1]); - } - - for (int i = sub_offset; i >= rj_start; i--) - Rj_slice[i] = ak_slice[i] * (Rj_slice[i + 1] + rj_slice[i]) / - (Rj_slice[i + 1] * rj_slice[i] + 1.0); - - if (sup_offset >= 2) { - // Replace std::pow(ak1, n) with an explicit multiply loop — pow calls - // exp(n*log(z)) which is ~30x slower than n multiplications for small n. - // sup_offset is typically 1-3 so this loop runs 0-1 times in practice. - auto ak_pow = std::complex{1.0, 0.0}; - for (int p = 0; p < sup_offset - 2; ++p) ak_pow *= ak1; - out[l] = std::norm(ak_pow * Rj_slice[sup_offset - 1]); - } else { - out[l] = std::norm(Rj_slice[0]); - } - } + for (int l = 0; l < n_q; l++) { + kk_slice[0] = m_consts.k0 * m_consts.indexsup * sinthetai[l]; + + // Flat-front wavevector — density uniform at [1, sup_offset-1] + const auto kk1 = m_consts.k0 * std::sqrt(m_consts.indexsupsquared * sin2thetai[l] - density_profile[1] + m_consts.sup_sld); + const auto ak1 = std::exp(ls.length_mult[1] * kk1); + + // Flat-back wavevector + const auto kk_sub = + m_consts.k0 * std::sqrt(m_consts.indexsupsquared * sin2thetai[l] - density_profile[n_layers - 1] + m_consts.sup_sld); + + // Varying region: kk and ak for [sup_offset, sub_offset] + for (int i = sup_offset; i <= sub_offset; i++) { + kk_slice[i] = m_consts.k0 * std::sqrt(m_consts.indexsupsquared * sin2thetai[l] - density_profile[i] + m_consts.sup_sld); + ak_slice[i] = std::exp(ls.length_mult[i] * kk_slice[i]); + } + + // Boundary entries for rj at the flat/varying interfaces. + const int rj_start = sup_offset >= 2 ? sup_offset - 1 : 0; + if (sup_offset >= 2) { + kk_slice[sup_offset - 1] = kk1; + ak_slice[sup_offset - 1] = ak1; + } + kk_slice[sub_offset + 1] = kk_sub; + + for (int i = rj_start; i <= sub_offset; i++) { + rj_slice[i] = (kk_slice[i] - kk_slice[i + 1]) / (kk_slice[i] + kk_slice[i + 1]); + if constexpr (HasRoughness) + rj_slice[i] *= std::exp(ls.sigma_sq[i] * kk_slice[i] * kk_slice[i + 1]); + } + + for (int i = sub_offset; i >= rj_start; i--) + Rj_slice[i] = ak_slice[i] * (Rj_slice[i + 1] + rj_slice[i]) / (Rj_slice[i + 1] * rj_slice[i] + 1.0); + + if (sup_offset >= 2) { + // Replace std::pow(ak1, n) with an explicit multiply loop — pow calls + // exp(n*log(z)) which is ~30x slower than n multiplications for small n. + // sup_offset is typically 1-3 so this loop runs 0-1 times in practice. + auto ak_pow = std::complex{1.0, 0.0}; + for (int p = 0; p < sup_offset - 2; ++p) + ak_pow *= ak1; + out[l] = std::norm(ak_pow * Rj_slice[sup_offset - 1]); + } else { + out[l] = std::norm(Rj_slice[0]); + } + } } -void ParrattReflectivity::ReflectivityCalc(const LayerStack &ls, int q_end) { - assert(!m_qsmear_enabled); - if (ls.has_roughness) - ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, - m_consts.sinsquaredthetai, m_refl_out, - q_end); - else - ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, - m_consts.sinsquaredthetai, m_refl_out, - q_end); +void ParrattReflectivity::ReflectivityCalc(const LayerStack& ls, int q_end) +{ + assert(!m_qsmear_enabled); + if (ls.has_roughness) + ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, m_consts.sinsquaredthetai, m_refl_out, q_end); + else + ReflectivityCalcCoreImpl(ls, m_consts.sinthetai, m_consts.sinsquaredthetai, m_refl_out, q_end); } -void ParrattReflectivity::TransparentReflectivityCalc(const LayerStack &ls) { - assert(!m_qsmear_enabled); - const auto &density_profile = ls.rho; - const int n_layers = static_cast(density_profile.size()); - const int n_q = static_cast(m_refl_out.size()); +void ParrattReflectivity::TransparentReflectivityCalc(const LayerStack& ls) +{ + assert(!m_qsmear_enabled); + const auto& density_profile = ls.rho; + const int n_layers = static_cast(density_profile.size()); + const int n_q = static_cast(m_refl_out.size()); - const int sup_offset = ls.sup_offset; - const int sub_offset = ls.sub_offset; - const int rj_start = sup_offset >= 2 ? sup_offset - 1 : 0; + const int sup_offset = ls.sup_offset; + const int sub_offset = ls.sub_offset; + const int rj_start = sup_offset >= 2 ? sup_offset - 1 : 0; - const int complex_to_real_offset = - FindComplexToRealOffset(m_consts.sinsquaredthetai, density_profile, - m_consts.indexsupsquared, m_consts.sup_sld); + const int complex_to_real_offset = + FindComplexToRealOffset(m_consts.sinsquaredthetai, density_profile, m_consts.indexsupsquared, m_consts.sup_sld); - // ls.length_mult[i].imag() = -2*length[i]; uniform for CEDP path. - const double angle_mult_flat = ls.length_mult[1].imag(); + // ls.length_mult[i].imag() = -2*length[i]; uniform for CEDP path. + const double angle_mult_flat = ls.length_mult[1].imag(); - // Both loops share the caller's parallel team so there is only one fork+join - // per SA run instead of two per SA iteration. nowait on the complex loop lets - // threads proceed to the real loop without an extra barrier between them. - const int threadnum = omp_get_thread_num(); - const int arrayoffset = threadnum * n_layers; + // Both loops share the caller's parallel team so there is only one fork+join + // per SA run instead of two per SA iteration. nowait on the complex loop lets + // threads proceed to the real loop without an extra barrier between them. + const int threadnum = omp_get_thread_num(); + const int arrayoffset = threadnum * n_layers; - auto *kk_slice = m_complex.kk.data() + arrayoffset; - auto *ak_slice = m_complex.ak.data() + arrayoffset; - auto *rj_slice = m_complex.rj.data() + arrayoffset; - auto *Rj_slice = m_complex.Rj.data() + arrayoffset; - auto *dkk_slice = m_real.kk.data() + arrayoffset; - auto *drj_slice = m_real.rj.data() + arrayoffset; + auto* kk_slice = m_complex.kk.data() + arrayoffset; + auto* ak_slice = m_complex.ak.data() + arrayoffset; + auto* rj_slice = m_complex.rj.data() + arrayoffset; + auto* Rj_slice = m_complex.Rj.data() + arrayoffset; + auto* dkk_slice = m_real.kk.data() + arrayoffset; + auto* drj_slice = m_real.rj.data() + arrayoffset; - ak_slice[0] = 1.0; - Rj_slice[sub_offset + 1] = 0.0; + ak_slice[0] = 1.0; + Rj_slice[sub_offset + 1] = 0.0; // Complex Q-points: evanescent layers present; use full complex arithmetic. #pragma omp for nowait schedule(runtime) for (int l = 0; l < complex_to_real_offset; l++) { - kk_slice[0] = m_consts.k0 * m_consts.indexsup * m_consts.sinthetai[l]; - - const auto kk1 = - m_consts.k0 * std::sqrt(m_consts.indexsupsquared * m_consts.sinsquaredthetai[l] - - density_profile[1] + m_consts.sup_sld); - const auto ak1 = std::exp(ls.length_mult[1] * kk1); - const auto kk_sub = - m_consts.k0 * std::sqrt(m_consts.indexsupsquared * m_consts.sinsquaredthetai[l] - - density_profile[n_layers - 1] + m_consts.sup_sld); - - for (int i = sup_offset; i <= sub_offset; i++) { - kk_slice[i] = - m_consts.k0 * std::sqrt(m_consts.indexsupsquared * m_consts.sinsquaredthetai[l] - - density_profile[i] + m_consts.sup_sld); - ak_slice[i] = std::exp(ls.length_mult[i] * kk_slice[i]); - } - - if (sup_offset >= 2) { - kk_slice[sup_offset - 1] = kk1; - ak_slice[sup_offset - 1] = ak1; - } - kk_slice[sub_offset + 1] = kk_sub; - - for (int i = rj_start; i <= sub_offset; i++) - rj_slice[i] = (kk_slice[i] - kk_slice[i + 1]) / (kk_slice[i] + kk_slice[i + 1]); - - for (int i = sub_offset; i >= rj_start; i--) - Rj_slice[i] = ak_slice[i] * (Rj_slice[i + 1] + rj_slice[i]) / - (Rj_slice[i + 1] * rj_slice[i] + 1.0); - - // kk1 is always real for the flat superstrate region, so |ak1| == 1 and - // the pow(ak1, n) phase factor drops out of the norm. - m_refl_out[l] = (sup_offset >= 2) ? std::norm(Rj_slice[sup_offset - 1]) - : std::norm(Rj_slice[0]); + kk_slice[0] = m_consts.k0 * m_consts.indexsup * m_consts.sinthetai[l]; + + const auto kk1 = + m_consts.k0 * std::sqrt(m_consts.indexsupsquared * m_consts.sinsquaredthetai[l] - density_profile[1] + m_consts.sup_sld); + const auto ak1 = std::exp(ls.length_mult[1] * kk1); + const auto kk_sub = m_consts.k0 * std::sqrt(m_consts.indexsupsquared * m_consts.sinsquaredthetai[l] - + density_profile[n_layers - 1] + m_consts.sup_sld); + + for (int i = sup_offset; i <= sub_offset; i++) { + kk_slice[i] = + m_consts.k0 * std::sqrt(m_consts.indexsupsquared * m_consts.sinsquaredthetai[l] - density_profile[i] + m_consts.sup_sld); + ak_slice[i] = std::exp(ls.length_mult[i] * kk_slice[i]); + } + + if (sup_offset >= 2) { + kk_slice[sup_offset - 1] = kk1; + ak_slice[sup_offset - 1] = ak1; + } + kk_slice[sub_offset + 1] = kk_sub; + + for (int i = rj_start; i <= sub_offset; i++) + rj_slice[i] = (kk_slice[i] - kk_slice[i + 1]) / (kk_slice[i] + kk_slice[i + 1]); + + for (int i = sub_offset; i >= rj_start; i--) + Rj_slice[i] = ak_slice[i] * (Rj_slice[i + 1] + rj_slice[i]) / (Rj_slice[i + 1] * rj_slice[i] + 1.0); + + // kk1 is always real for the flat superstrate region, so |ak1| == 1 and + // the pow(ak1, n) phase factor drops out of the norm. + m_refl_out[l] = (sup_offset >= 2) ? std::norm(Rj_slice[sup_offset - 1]) : std::norm(Rj_slice[0]); } // Real Q-points: all wavevectors real; exp({0,-2dz}*real_kk) is a pure @@ -322,46 +298,37 @@ void ParrattReflectivity::TransparentReflectivityCalc(const LayerStack &ls) { // call that the general complex overload would perform. #pragma omp for schedule(runtime) for (int l = complex_to_real_offset; l < n_q; l++) { - dkk_slice[0] = m_consts.k0 * m_consts.indexsup * m_consts.sinthetai[l]; - - const auto dkk1 = - m_consts.k0 * - sqrt((m_consts.indexsupsquared * m_consts.sinsquaredthetai[l]) - - density_profile[1].real() + m_consts.sup_sld); - const double angle0 = angle_mult_flat * dkk1; - const auto dak1 = std::complex(std::cos(angle0), std::sin(angle0)); - - const auto dkk_sub = - m_consts.k0 * - sqrt((m_consts.indexsupsquared * m_consts.sinsquaredthetai[l]) - - density_profile.back().real() + m_consts.sup_sld); - - for (int i = sup_offset; i <= sub_offset; i++) { - dkk_slice[i] = - m_consts.k0 * - sqrt((m_consts.indexsupsquared * m_consts.sinsquaredthetai[l]) - - density_profile[i].real() + m_consts.sup_sld); - const double angle_i = ls.length_mult[i].imag() * dkk_slice[i]; - ak_slice[i] = {std::cos(angle_i), std::sin(angle_i)}; - } - - if (sup_offset >= 2) { - dkk_slice[sup_offset - 1] = dkk1; - ak_slice[sup_offset - 1] = dak1; - } - dkk_slice[sub_offset + 1] = dkk_sub; - - for (int i = rj_start; i <= sub_offset; i++) - drj_slice[i] = (dkk_slice[i] - dkk_slice[i + 1]) / - (dkk_slice[i] + dkk_slice[i + 1]); - - for (int i = sub_offset; i >= rj_start; i--) - Rj_slice[i] = ak_slice[i] * (Rj_slice[i + 1] + drj_slice[i]) / - (Rj_slice[i + 1] * drj_slice[i] + 1.0); - - // |dak1| == 1 (pure phase), so the flat-front propagation factor drops - // out of the norm. - m_refl_out[l] = (sup_offset >= 2) ? std::norm(Rj_slice[sup_offset - 1]) - : std::norm(Rj_slice[0]); - } + dkk_slice[0] = m_consts.k0 * m_consts.indexsup * m_consts.sinthetai[l]; + + const auto dkk1 = + m_consts.k0 * sqrt((m_consts.indexsupsquared * m_consts.sinsquaredthetai[l]) - density_profile[1].real() + m_consts.sup_sld); + const double angle0 = angle_mult_flat * dkk1; + const auto dak1 = std::complex(std::cos(angle0), std::sin(angle0)); + + const auto dkk_sub = m_consts.k0 * sqrt((m_consts.indexsupsquared * m_consts.sinsquaredthetai[l]) - density_profile.back().real() + + m_consts.sup_sld); + + for (int i = sup_offset; i <= sub_offset; i++) { + dkk_slice[i] = m_consts.k0 * + sqrt((m_consts.indexsupsquared * m_consts.sinsquaredthetai[l]) - density_profile[i].real() + m_consts.sup_sld); + const double angle_i = ls.length_mult[i].imag() * dkk_slice[i]; + ak_slice[i] = {std::cos(angle_i), std::sin(angle_i)}; + } + + if (sup_offset >= 2) { + dkk_slice[sup_offset - 1] = dkk1; + ak_slice[sup_offset - 1] = dak1; + } + dkk_slice[sub_offset + 1] = dkk_sub; + + for (int i = rj_start; i <= sub_offset; i++) + drj_slice[i] = (dkk_slice[i] - dkk_slice[i + 1]) / (dkk_slice[i] + dkk_slice[i + 1]); + + for (int i = sub_offset; i >= rj_start; i--) + Rj_slice[i] = ak_slice[i] * (Rj_slice[i + 1] + drj_slice[i]) / (Rj_slice[i + 1] * drj_slice[i] + 1.0); + + // |dak1| == 1 (pure phase), so the flat-front propagation factor drops + // out of the norm. + m_refl_out[l] = (sup_offset >= 2) ? std::norm(Rj_slice[sup_offset - 1]) : std::norm(Rj_slice[0]); + } } diff --git a/tests/MIRefl.cpp b/tests/MIRefl.cpp index 61f33cd..4590d7e 100644 --- a/tests/MIRefl.cpp +++ b/tests/MIRefl.cpp @@ -1,87 +1,82 @@ // MIRefl.cpp : Defines the entry point for the console application. -#include "platform.h" #include "CEDP.h" -#include "UnifiedReflectivity.h" #include "ParamVector.h" +#include "UnifiedReflectivity.h" +#include "platform.h" static void FillInitStruct(ReflSettings& Initstruct); // Test Q-range (Angstrom^-1), 100 points from 0.02 to 0.60 -static const double qrange[] = { - 0.020, 0.026, 0.032, 0.038, 0.044, 0.050, 0.056, 0.062, 0.068, 0.074, - 0.080, 0.086, 0.092, 0.098, 0.104, 0.110, 0.116, 0.122, 0.128, 0.134, - 0.140, 0.146, 0.152, 0.158, 0.164, 0.170, 0.176, 0.182, 0.188, 0.194, - 0.200, 0.206, 0.212, 0.218, 0.224, 0.230, 0.236, 0.242, 0.248, 0.254, - 0.260, 0.266, 0.272, 0.278, 0.284, 0.290, 0.296, 0.302, 0.308, 0.314, - 0.320, 0.326, 0.332, 0.338, 0.344, 0.350, 0.356, 0.362, 0.368, 0.374, - 0.380, 0.386, 0.392, 0.398, 0.404, 0.410, 0.416, 0.422, 0.428, 0.434, - 0.440, 0.446, 0.452, 0.458, 0.464, 0.470, 0.476, 0.482, 0.488, 0.494, - 0.500, 0.506, 0.512, 0.518, 0.524, 0.530, 0.536, 0.542, 0.548, 0.554, - 0.560, 0.566, 0.572, 0.578, 0.584, 0.590, 0.594, 0.596, 0.598, 0.600 -}; +static const double qrange[] = {0.020, 0.026, 0.032, 0.038, 0.044, 0.050, 0.056, 0.062, 0.068, 0.074, 0.080, 0.086, 0.092, 0.098, 0.104, + 0.110, 0.116, 0.122, 0.128, 0.134, 0.140, 0.146, 0.152, 0.158, 0.164, 0.170, 0.176, 0.182, 0.188, 0.194, + 0.200, 0.206, 0.212, 0.218, 0.224, 0.230, 0.236, 0.242, 0.248, 0.254, 0.260, 0.266, 0.272, 0.278, 0.284, + 0.290, 0.296, 0.302, 0.308, 0.314, 0.320, 0.326, 0.332, 0.338, 0.344, 0.350, 0.356, 0.362, 0.368, 0.374, + 0.380, 0.386, 0.392, 0.398, 0.404, 0.410, 0.416, 0.422, 0.428, 0.434, 0.440, 0.446, 0.452, 0.458, 0.464, + 0.470, 0.476, 0.482, 0.488, 0.494, 0.500, 0.506, 0.512, 0.518, 0.524, 0.530, 0.536, 0.542, 0.548, 0.554, + 0.560, 0.566, 0.572, 0.578, 0.584, 0.590, 0.594, 0.596, 0.598, 0.600}; int main(int argc, char* argv[]) { - CEDP EDPGen; - ReflSettings InitStruct = {}; - int calculations = 10; + CEDP EDPGen; + ReflSettings InitStruct = {}; + int calculations = 10; - // Normalized SLD values (divided by FilmSLD = 9.38) - // Superphase (air): 0.0, Box1: 8.911, Box2: 13.5072, Subphase: 9.38 - double FilmSLD = 9.38; - double SLD[] = {0.0, 8.911, 13.5072, 9.38}; + // Normalized SLD values (divided by FilmSLD = 9.38) + // Superphase (air): 0.0, Box1: 8.911, Box2: 13.5072, Subphase: 9.38 + double FilmSLD = 9.38; + double SLD[] = {0.0, 8.911, 13.5072, 9.38}; - FillInitStruct(InitStruct); + FillInitStruct(InitStruct); - // Build a ParamVector and populate it with normalized SLD values - ParamVector params(InitStruct); - for(int i = 0; i < InitStruct.Boxes; i++) - params.SetMutatableParameter(i, SLD[i+1] / FilmSLD); - params.SetRoughness(3.15); + // Build a ParamVector and populate it with normalized SLD values + ParamVector params(InitStruct); + for (int i = 0; i < InitStruct.Boxes; i++) + params.SetMutatableParameter(i, SLD[i + 1] / FilmSLD); + params.SetRoughness(3.15); - EDPGen.Init(InitStruct); - ParrattReflectivity Refl(InitStruct); + EDPGen.Init(InitStruct); + ParrattReflectivity Refl(InitStruct); - int t_on = clock(); + int t_on = clock(); - for(int i = 0; i < calculations; i++) - { - EDPGen.GenerateEDP(params); - Refl.CalculateReflectivity(EDPGen); - } + for (int i = 0; i < calculations; i++) { + EDPGen.GenerateEDP(params); + Refl.CalculateReflectivity(EDPGen); + } - int t_off = clock(); + int t_off = clock(); - cout << calculations << " calculations in: " << (((static_cast(t_off - t_on))/(CLOCKS_PER_SEC)))*(1000000) << " microseconds\n\n"; - cout << calculations/(((static_cast(t_off - t_on))/(CLOCKS_PER_SEC))) << " calcuations per second\n\n"; - cout << "1 calculation in: " << (((static_cast(t_off - t_on))/(CLOCKS_PER_SEC)))*(1E6)/calculations << " microseconds\n\n"; + cout << calculations << " calculations in: " << (((static_cast(t_off - t_on)) / (CLOCKS_PER_SEC))) * (1000000) + << " microseconds\n\n"; + cout << calculations / (((static_cast(t_off - t_on)) / (CLOCKS_PER_SEC))) << " calcuations per second\n\n"; + cout << "1 calculation in: " << (((static_cast(t_off - t_on)) / (CLOCKS_PER_SEC))) * (1E6) / calculations << " microseconds\n\n"; - return 0; + return 0; } static void FillInitStruct(ReflSettings& InitStruct) { - InitStruct.Wavelength = 1.24; - InitStruct.QErr = 0; - InitStruct.XRonly = 0; - InitStruct.Impnorm = 0; - InitStruct.Q = std::vector(std::begin(qrange), std::end(qrange)); - InitStruct.CritEdgeOffset = 0; - InitStruct.HighQOffset = 0; - - InitStruct.UseSurfAbs = 0; - InitStruct.SupAbs = 0; - InitStruct.SubAbs = 0; - InitStruct.FilmAbs = 0; - InitStruct.Boxes = 2; - InitStruct.FilmLength = 26; - InitStruct.Resolution = 10; - InitStruct.FilmSLD = 9.38; - InitStruct.SupSLD = 0.0; - InitStruct.SubSLD = 9.38; - InitStruct.Forcesig = 0; - InitStruct.Objectivefunction = 0; - - // Refl, ReflError, QError left empty — only CalculateReflectivity is called + InitStruct.Wavelength = 1.24; + InitStruct.QErr = 0; + InitStruct.XRonly = 0; + InitStruct.Impnorm = 0; + InitStruct.Q = std::vector(std::begin(qrange), std::end(qrange)); + InitStruct.CritEdgeOffset = 0; + InitStruct.HighQOffset = 0; + + InitStruct.UseSurfAbs = 0; + InitStruct.SupAbs = 0; + InitStruct.SubAbs = 0; + InitStruct.FilmAbs = 0; + InitStruct.Boxes = 2; + InitStruct.FilmLength = 26; + InitStruct.Resolution = 10; + InitStruct.FilmSLD = 9.38; + InitStruct.SupSLD = 0.0; + InitStruct.SubSLD = 9.38; + InitStruct.Forcesig = 0; + InitStruct.Objectivefunction = 0; + + // Refl, ReflError, QError left empty — only CalculateReflectivity is called } diff --git a/tests/stochfit_profile.cpp b/tests/stochfit_profile.cpp index ef57400..0c7a6c5 100644 --- a/tests/stochfit_profile.cpp +++ b/tests/stochfit_profile.cpp @@ -12,20 +12,17 @@ #include #ifdef _WIN32 -# define NOMINMAX -# include + #define NOMINMAX + #include #endif #include "StochFitHarness.h" #ifndef RESOURCES_DIR -#define RESOURCES_DIR "../../resources" + #define RESOURCES_DIR "../../resources" #endif -static bool LoadReflData(const std::string& path, - std::vector& q, - std::vector& refl, - std::vector& reflErr) +static bool LoadReflData(const std::string& path, std::vector& q, std::vector& refl, std::vector& reflErr) { std::ifstream f(path); if (!f) { @@ -34,10 +31,12 @@ static bool LoadReflData(const std::string& path, } std::string line; while (std::getline(f, line)) { - if (line.empty() || line[0] == '#') continue; + if (line.empty() || line[0] == '#') + continue; std::istringstream ss(line); double qv, rv, ev; - if (!(ss >> qv >> rv >> ev)) continue; + if (!(ss >> qv >> rv >> ev)) + continue; q.push_back(qv); refl.push_back(rv); reflErr.push_back(ev); @@ -51,10 +50,10 @@ int main(int argc, char* argv[]) SetConsoleOutputCP(CP_UTF8); #endif std::string dataPath = RESOURCES_DIR "/test1refl.txt"; - int iterations = 5'000'000; - int resolution = 3; - bool xrOnly = true; - double targetChi = 0.0; // 0 = disabled + int iterations = 5'000'000; + int resolution = 3; + bool xrOnly = true; + double targetChi = 0.0; // 0 = disabled for (int i = 1; i < argc; ++i) { std::string_view arg = argv[i]; @@ -75,54 +74,55 @@ int main(int argc, char* argv[]) } std::vector q, refl, reflErr; - if (!LoadReflData(dataPath, q, refl, reflErr)) return 1; + if (!LoadReflData(dataPath, q, refl, reflErr)) + return 1; std::cout << "Loaded " << q.size() << " data points from " << dataPath << "\n"; ReflSettings settings{}; - settings.Q = q; - settings.Refl = refl; + settings.Q = q; + settings.Refl = refl; settings.ReflError = reflErr; settings.Directory = "."; // SLD parameters (silicon substrate, X-ray at 1.24 Å) - settings.SubSLD = 9.38; - settings.FilmSLD = 9.38; - settings.SupSLD = 0.0; - settings.SubAbs = 2e-8; - settings.FilmAbs = 1e-14; - settings.SupAbs = 0.0; + settings.SubSLD = 9.38; + settings.FilmSLD = 9.38; + settings.SupSLD = 0.0; + settings.SubAbs = 2e-8; + settings.FilmAbs = 1e-14; + settings.SupAbs = 0.0; settings.UseSurfAbs = false; settings.Wavelength = 1.24; // Box model - settings.Boxes = 40; + settings.Boxes = 40; settings.FilmLength = 25.0; settings.Resolution = resolution; - settings.Forcesig = 0.0; + settings.Forcesig = 0.0; // Q trimming settings.CritEdgeOffset = 0; - settings.HighQOffset = 0; + settings.HighQOffset = 0; // Algorithm — Simulated Annealing with Log(R) objective (matches tuned defaults) - settings.Objectivefunction = 0; - settings.Algorithm = 0; - settings.Paramtemp = 0.03; - settings.Inittemp = 10.0; - settings.Platiter = 4000; - settings.Slope = 0.95; - settings.Gamma = 0.05; - settings.Tempiter = 100; - settings.STUNdeciter = 200000; - settings.Gammadec = 0.85; - settings.STUNfunc = 0; - settings.Adaptive = false; - settings.Sigmasearch = 10; + settings.Objectivefunction = 0; + settings.Algorithm = 0; + settings.Paramtemp = 0.03; + settings.Inittemp = 10.0; + settings.Platiter = 4000; + settings.Slope = 0.95; + settings.Gamma = 0.05; + settings.Tempiter = 100; + settings.STUNdeciter = 200000; + settings.Gammadec = 0.85; + settings.STUNfunc = 0; + settings.Adaptive = false; + settings.Sigmasearch = 10; settings.NormalizationSearchPerc = 0; - settings.AbsorptionSearchPerc = 0; - settings.Impnorm = false; - settings.XRonly = xrOnly; - settings.Iterations = iterations; + settings.AbsorptionSearchPerc = 0; + settings.Impnorm = false; + settings.XRonly = xrOnly; + settings.Iterations = iterations; StochFit harness(settings); if (auto err = harness.GetInitError(); !err) { @@ -132,21 +132,20 @@ int main(int argc, char* argv[]) auto wallStart = std::chrono::steady_clock::now(); harness.Start(iterations); - std::printf("SA started: %d iterations, %d boxes, %zu data points, resolution=%d, xrOnly=%s\n\n", - iterations, settings.Boxes, q.size(), resolution, xrOnly ? "true" : "false"); + std::printf("SA started: %d iterations, %d boxes, %zu data points, resolution=%d, xrOnly=%s\n\n", iterations, settings.Boxes, q.size(), + resolution, xrOnly ? "true" : "false"); while (true) { std::this_thread::sleep_for(std::chrono::seconds(2)); DataSnapshot snap = harness.GetData(); - double elapsed = std::chrono::duration( - std::chrono::steady_clock::now() - wallStart).count(); + double elapsed = std::chrono::duration(std::chrono::steady_clock::now() - wallStart).count(); double pct = 100.0 * snap.iteration / iterations; - double ips = snap.iteration / elapsed; - std::printf("[%5.1f%%] iter=%-9d χ²=%.4e GoF=%.4f rough=%.3f %.0f iter/s\n", - pct, snap.iteration, snap.chiSquare, snap.goodnessOfFit, - snap.roughness, ips); + double ips = snap.iteration / elapsed; + std::printf("[%5.1f%%] iter=%-9d χ²=%.4e GoF=%.4f rough=%.3f %.0f iter/s\n", pct, snap.iteration, snap.chiSquare, + snap.goodnessOfFit, snap.roughness, ips); std::fflush(stdout); - if (snap.isFinished) break; + if (snap.isFinished) + break; if (targetChi > 0.0 && snap.chiSquare > 0.0 && snap.chiSquare <= targetChi) { std::printf(" --> χ² target %.6g reached\n", targetChi); break; @@ -154,8 +153,7 @@ int main(int argc, char* argv[]) } harness.Stop(); - double totalSec = std::chrono::duration( - std::chrono::steady_clock::now() - wallStart).count(); + double totalSec = std::chrono::duration(std::chrono::steady_clock::now() - wallStart).count(); DataSnapshot final = harness.GetData(); std::printf("\n=== Profile Complete ===\n"); diff --git a/tests/test_reflectivity.cpp b/tests/test_reflectivity.cpp index 419fff4..c33a587 100644 --- a/tests/test_reflectivity.cpp +++ b/tests/test_reflectivity.cpp @@ -4,51 +4,47 @@ // Any change that corrupts the Parratt recursion or EDP generation will fail. #include -#include #include #include #include #include +#include + #include #include -#include #include +#include #include // Identical to the qrange in src/mirefl/MIRefl.cpp -static const double qrange[] = { - 0.020, 0.026, 0.032, 0.038, 0.044, 0.050, 0.056, 0.062, 0.068, 0.074, - 0.080, 0.086, 0.092, 0.098, 0.104, 0.110, 0.116, 0.122, 0.128, 0.134, - 0.140, 0.146, 0.152, 0.158, 0.164, 0.170, 0.176, 0.182, 0.188, 0.194, - 0.200, 0.206, 0.212, 0.218, 0.224, 0.230, 0.236, 0.242, 0.248, 0.254, - 0.260, 0.266, 0.272, 0.278, 0.284, 0.290, 0.296, 0.302, 0.308, 0.314, - 0.320, 0.326, 0.332, 0.338, 0.344, 0.350, 0.356, 0.362, 0.368, 0.374, - 0.380, 0.386, 0.392, 0.398, 0.404, 0.410, 0.416, 0.422, 0.428, 0.434, - 0.440, 0.446, 0.452, 0.458, 0.464, 0.470, 0.476, 0.482, 0.488, 0.494, - 0.500, 0.506, 0.512, 0.518, 0.524, 0.530, 0.536, 0.542, 0.548, 0.554, - 0.560, 0.566, 0.572, 0.578, 0.584, 0.590, 0.594, 0.596, 0.598, 0.600 -}; +static const double qrange[] = {0.020, 0.026, 0.032, 0.038, 0.044, 0.050, 0.056, 0.062, 0.068, 0.074, 0.080, 0.086, 0.092, 0.098, 0.104, + 0.110, 0.116, 0.122, 0.128, 0.134, 0.140, 0.146, 0.152, 0.158, 0.164, 0.170, 0.176, 0.182, 0.188, 0.194, + 0.200, 0.206, 0.212, 0.218, 0.224, 0.230, 0.236, 0.242, 0.248, 0.254, 0.260, 0.266, 0.272, 0.278, 0.284, + 0.290, 0.296, 0.302, 0.308, 0.314, 0.320, 0.326, 0.332, 0.338, 0.344, 0.350, 0.356, 0.362, 0.368, 0.374, + 0.380, 0.386, 0.392, 0.398, 0.404, 0.410, 0.416, 0.422, 0.428, 0.434, 0.440, 0.446, 0.452, 0.458, 0.464, + 0.470, 0.476, 0.482, 0.488, 0.494, 0.500, 0.506, 0.512, 0.518, 0.524, 0.530, 0.536, 0.542, 0.548, 0.554, + 0.560, 0.566, 0.572, 0.578, 0.584, 0.590, 0.594, 0.596, 0.598, 0.600}; static void FillInitStruct(ReflSettings& s) { - s.Wavelength = 1.24; - s.QErr = 0; - s.XRonly = false; - s.Impnorm = false; - s.Q = std::vector(std::begin(qrange), std::end(qrange)); - s.CritEdgeOffset = 0; - s.HighQOffset = 0; - s.UseSurfAbs = false; - s.SupAbs = 0; - s.SubAbs = 0; - s.FilmAbs = 0; - s.Boxes = 2; - s.FilmLength = 26; - s.Resolution = 10; - s.FilmSLD = 9.38; - s.SupSLD = 0.0; - s.SubSLD = 9.38; - s.Forcesig = 0; + s.Wavelength = 1.24; + s.QErr = 0; + s.XRonly = false; + s.Impnorm = false; + s.Q = std::vector(std::begin(qrange), std::end(qrange)); + s.CritEdgeOffset = 0; + s.HighQOffset = 0; + s.UseSurfAbs = false; + s.SupAbs = 0; + s.SubAbs = 0; + s.FilmAbs = 0; + s.Boxes = 2; + s.FilmLength = 26; + s.Resolution = 10; + s.FilmSLD = 9.38; + s.SupSLD = 0.0; + s.SubSLD = 9.38; + s.Forcesig = 0; s.Objectivefunction = 0; // Refl, ReflError, QError left empty — only CalculateReflectivity is called } @@ -56,9 +52,7 @@ static void FillInitStruct(ReflSettings& s) // Relative-tolerance check. At a thin-film interference minimum R can be // many orders of magnitude smaller than neighbouring points, so pure // absolute tolerance is unsuitable here. -static void ExpectRelative(double actual, double expected, - double reltol = 1e-3, - const char* label = "") +static void ExpectRelative(double actual, double expected, double reltol = 1e-3, const char* label = "") { double tol = std::fabs(expected) * reltol + 1e-15; EXPECT_NEAR(actual, expected, tol) << label; @@ -68,8 +62,8 @@ static void ExpectRelative(double actual, double expected, // Q < Qc uses the complex path; Q > Qc uses the real-kk fast path. TEST(Reflectivity, TwoBoxLipidFilmParratt) { - const double FilmSLD = 9.38; - const double SLD[] = {0.0, 8.911, 13.5072, 9.38}; // sup, box1, box2, sub + const double FilmSLD = 9.38; + const double SLD[] = {0.0, 8.911, 13.5072, 9.38}; // sup, box1, box2, sub ReflSettings init = {}; FillInitStruct(init); @@ -91,35 +85,35 @@ TEST(Reflectivity, TwoBoxLipidFilmParratt) // ── Q < Qc region (complex Parratt path) ────────────────────────────── // Q=0.020 total external reflection (below critical edge) - ExpectRelative(reflOut[0], 1.0, 1e-6, "Q=0.020"); + ExpectRelative(reflOut[0], 1.0, 1e-6, "Q=0.020"); // Q=0.026 just above critical Q, sharp drop - ExpectRelative(reflOut[1], 9.57934e-02, 1e-3, "Q=0.026"); + ExpectRelative(reflOut[1], 9.57934e-02, 1e-3, "Q=0.026"); // ── Q > Qc region (real Parratt path) ───────────────────────────────── // Q=0.032 continuing fall - ExpectRelative(reflOut[2], 2.96789e-02, 1e-3, "Q=0.032"); + ExpectRelative(reflOut[2], 2.96789e-02, 1e-3, "Q=0.032"); // Q=0.038 - ExpectRelative(reflOut[3], 1.36002e-02, 1e-3, "Q=0.038"); + ExpectRelative(reflOut[3], 1.36002e-02, 1e-3, "Q=0.038"); // Q=0.050 Kiessig fringe region - ExpectRelative(reflOut[5], 4.56059e-03, 1e-3, "Q=0.050"); + ExpectRelative(reflOut[5], 4.56059e-03, 1e-3, "Q=0.050"); // Q=0.080 - ExpectRelative(reflOut[10], 7.70483e-04, 1e-3, "Q=0.080"); + ExpectRelative(reflOut[10], 7.70483e-04, 1e-3, "Q=0.080"); // Q=0.140 rapid fall-off - ExpectRelative(reflOut[20], 3.53317e-05, 1e-3, "Q=0.140"); + ExpectRelative(reflOut[20], 3.53317e-05, 1e-3, "Q=0.140"); // Q=0.188 near thin-film interference minimum — use looser tolerance - ExpectRelative(reflOut[28], 6.14088e-08, 2e-2, "Q=0.188"); + ExpectRelative(reflOut[28], 6.14088e-08, 2e-2, "Q=0.188"); // Q=0.200 recovering after minimum - ExpectRelative(reflOut[30], 2.13742e-07, 1e-3, "Q=0.200"); + ExpectRelative(reflOut[30], 2.13742e-07, 1e-3, "Q=0.200"); // Q=0.260 second fringe maximum - ExpectRelative(reflOut[40], 2.61797e-06, 1e-3, "Q=0.260"); + ExpectRelative(reflOut[40], 2.61797e-06, 1e-3, "Q=0.260"); // Q=0.320 - ExpectRelative(reflOut[50], 9.57848e-07, 1e-3, "Q=0.320"); + ExpectRelative(reflOut[50], 9.57848e-07, 1e-3, "Q=0.320"); // Q=0.440 - ExpectRelative(reflOut[70], 1.04542e-07, 1e-3, "Q=0.440"); + ExpectRelative(reflOut[70], 1.04542e-07, 1e-3, "Q=0.440"); // Q=0.560 high-Q tail - ExpectRelative(reflOut[90], 4.49987e-10, 2e-2, "Q=0.560"); + ExpectRelative(reflOut[90], 4.49987e-10, 2e-2, "Q=0.560"); // Q=0.600 - ExpectRelative(reflOut[99], 1.19181e-09, 2e-2, "Q=0.600"); + ExpectRelative(reflOut[99], 1.19181e-09, 2e-2, "Q=0.600"); } // ── Box model Parratt tests ─────────────────────────────────────────────────── @@ -134,14 +128,13 @@ TEST(Reflectivity, TwoBoxLipidFilmParratt) // // Reference values are derived analytically from the Fresnel/Parratt formulae. -static void MakeBoxReflSettings(ReflSettings &s, double wavelength, - const std::vector &Q) +static void MakeBoxReflSettings(ReflSettings& s, double wavelength, const std::vector& Q) { s = {}; s.Wavelength = wavelength; - s.SupSLD = 0.0; // air (Parratt units: 2*raw*rhofactor = 0 for air) - s.Q = Q; - s.QErr = 0; + s.SupSLD = 0.0; // air (Parratt units: 2*raw*rhofactor = 0 for air) + s.Q = Q; + s.QErr = 0; } // Single-box film where SLD_film = SLD_sub (degenerate: film has same density as @@ -153,32 +146,32 @@ static void MakeBoxReflSettings(ReflSettings &s, double wavelength, // r = 0.000847/0.139153 = 6.087e-3 -> R = 3.703e-5 TEST(Reflectivity, BoxModelNoFilm) { - const double SubSLD = 9.38; - const double lambda = 1.24; + const double SubSLD = 9.38; + const double lambda = 1.24; const double rhofactor = 1e-6 * lambda * lambda / (2.0 * std::numbers::pi); - const double rho2 = 2.0 * rhofactor; + const double rho2 = 2.0 * rhofactor; ReflSettings s; MakeBoxReflSettings(s, lambda, {0.14}); std::vector> rho = { - {0.0, 0.0}, // superstrate (air) - {SubSLD * rho2, 0.0}, // film = Si (no contrast → r12 = 0) - {SubSLD * rho2, 0.0}, // substrate (Si) + {0.0, 0.0}, // superstrate (air) + {SubSLD * rho2, 0.0}, // film = Si (no contrast → r12 = 0) + {SubSLD * rho2, 0.0}, // substrate (Si) }; std::vector> lm = { - {0.0, 0.0}, // superstrate (zero length) - {0.0, -2.0 * 200.0}, // film 200 Å (irrelevant since r12=0) - {0.0, 0.0}, // substrate + {0.0, 0.0}, // superstrate (zero length) + {0.0, -2.0 * 200.0}, // film 200 Å (irrelevant since r12=0) + {0.0, 0.0}, // substrate }; LayerStack ls{}; - ls.rho = rho; - ls.length_mult = lm; - ls.sup_offset = 0; - ls.sub_offset = 1; + ls.rho = rho; + ls.length_mult = lm; + ls.sup_offset = 0; + ls.sub_offset = 1; ls.has_roughness = false; - ls.transparent = true; + ls.transparent = true; ParrattReflectivity parratt(s); auto refl = parratt.CalculateReflectivity(ls); @@ -192,28 +185,27 @@ TEST(Reflectivity, BoxModelNoFilm) // Also checks Q=0.016 (below Qc_Si ≈ 0.0217) where R must be ≈ 1. TEST(Reflectivity, BoxModelSingleFilm) { - const double SubSLD = 9.38; - const double filmSLD = SubSLD / 2.0; // 4.69 - const double d = 200.0; - const double lambda = 1.24; + const double SubSLD = 9.38; + const double filmSLD = SubSLD / 2.0; // 4.69 + const double d = 200.0; + const double lambda = 1.24; const double rhofactor = 1e-6 * lambda * lambda / (2.0 * std::numbers::pi); - const double rho2 = 2.0 * rhofactor; + const double rho2 = 2.0 * rhofactor; // -- Q = 0.016 (below Qc_Si ≈ 0.0217): total reflection ---------------------- { ReflSettings s; MakeBoxReflSettings(s, lambda, {0.016}); - std::vector> rho = { - {0.0, 0.0}, {filmSLD*rho2, 0.0}, {SubSLD*rho2, 0.0} - }; - std::vector> lm = { - {0.0, 0.0}, {0.0, -2.0*d}, {0.0, 0.0} - }; + std::vector> rho = {{0.0, 0.0}, {filmSLD * rho2, 0.0}, {SubSLD * rho2, 0.0}}; + std::vector> lm = {{0.0, 0.0}, {0.0, -2.0 * d}, {0.0, 0.0}}; LayerStack ls{}; - ls.rho = rho; ls.length_mult = lm; - ls.sup_offset = 0; ls.sub_offset = 1; - ls.has_roughness = false; ls.transparent = false; + ls.rho = rho; + ls.length_mult = lm; + ls.sup_offset = 0; + ls.sub_offset = 1; + ls.has_roughness = false; + ls.transparent = false; ParrattReflectivity parratt(s); auto refl = parratt.CalculateReflectivity(ls); @@ -227,16 +219,15 @@ TEST(Reflectivity, BoxModelSingleFilm) ReflSettings s; MakeBoxReflSettings(s, lambda, {0.14}); - std::vector> rho = { - {0.0, 0.0}, {filmSLD*rho2, 0.0}, {SubSLD*rho2, 0.0} - }; - std::vector> lm = { - {0.0, 0.0}, {0.0, -2.0*d}, {0.0, 0.0} - }; + std::vector> rho = {{0.0, 0.0}, {filmSLD * rho2, 0.0}, {SubSLD * rho2, 0.0}}; + std::vector> lm = {{0.0, 0.0}, {0.0, -2.0 * d}, {0.0, 0.0}}; LayerStack ls{}; - ls.rho = rho; ls.length_mult = lm; - ls.sup_offset = 0; ls.sub_offset = 1; - ls.has_roughness = false; ls.transparent = true; + ls.rho = rho; + ls.length_mult = lm; + ls.sup_offset = 0; + ls.sub_offset = 1; + ls.has_roughness = false; + ls.transparent = true; ParrattReflectivity parratt(s); auto refl = parratt.CalculateReflectivity(ls); @@ -257,10 +248,10 @@ TEST(Reflectivity, BoxModelSingleFilm) // misaligned with the data by CritEdgeOffset points. TEST(Harness, QOffsetsSliceFitWindow) { - constexpr int kLowOff = 5; + constexpr int kLowOff = 5; constexpr int kHighOff = 5; - constexpr int kFull = 100; - constexpr int kWindow = kFull - kLowOff - kHighOff; + constexpr int kFull = 100; + constexpr int kWindow = kFull - kLowOff - kHighOff; ReflSettings init = {}; FillInitStruct(init); @@ -286,7 +277,7 @@ TEST(Harness, QOffsetsSliceFitWindow) init.ReflError.assign(kFull, 1.0); std::copy(model.begin(), model.end(), init.Refl.begin() + kLowOff); init.CritEdgeOffset = kLowOff; - init.HighQOffset = kHighOff; + init.HighQOffset = kHighOff; StochFit fit(init); ASSERT_TRUE(fit.GetInitError().has_value()) << fit.GetInitError().error(); @@ -299,6 +290,6 @@ TEST(Harness, QOffsetsSliceFitWindow) const DataSnapshot snap = fit.GetData(); ASSERT_EQ(snap.Q.size(), static_cast(kWindow)); EXPECT_DOUBLE_EQ(snap.Q.front(), qrange[kLowOff]); - EXPECT_DOUBLE_EQ(snap.Q.back(), qrange[kFull - kHighOff - 1]); + EXPECT_DOUBLE_EQ(snap.Q.back(), qrange[kFull - kHighOff - 1]); EXPECT_EQ(snap.refl.size(), static_cast(kWindow)); }