From c153797c36a96d18dcbb2da1ef8bd192600d6e5f Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 9 Apr 2025 14:57:17 +0200 Subject: [PATCH 1/7] first implementation of ILSE with rescaled variances Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 46 ++++ .../iterative_linear_se_solver.hpp | 78 +++--- tests/cpp_unit_tests/test_statistics.cpp | 243 ++++++++++++++++++ 3 files changed, 338 insertions(+), 29 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index a4fdab969..d581dafdd 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -204,4 +204,50 @@ template struct PolarComplexRandVar { 9.0}}; } }; + +namespace statistics { +template typename RandVarType> + requires is_in_list_c, UniformRealRandVar, IndependentRealRandVar, + UniformComplexRandVar, IndependentComplexRandVar> +inline auto scale(RandVarType const& var, double scale_factor) { + return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)}; +} + +template + requires is_in_list_c, IndependentComplexRandVar> +inline auto scale(RandVarType const& var, RealValue const& scale_factor) { + return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)}; +} + +template typename RandVarType> + requires is_in_list_c, UniformComplexRandVar, IndependentComplexRandVar> +inline auto scale(RandVarType const& var, DoubleComplex const& scale_factor) { + return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)}; +} + +inline auto scale(IndependentComplexRandVar const& var, ComplexValue const& scale_factor) { + return IndependentComplexRandVar{.value = var.value * scale_factor, + .variance = var.variance * abs2(scale_factor)}; +} + +template + requires is_in_list_c, RealValue> +inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { + return DecomposedComplexRandVar{.real_component = scale(var.real_component, scale_factor), + .imag_component = scale(var.imag_component, scale_factor)}; +} + +template + requires is_in_list_c, ComplexValue> +inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { + auto const scaled_value = var.value() * scale_factor; + return DecomposedComplexRandVar{ + .real_component = {.value = real(scaled_value), + .variance = var.real_component.variance * abs2(real(scale_factor)) + + var.imag_component.variance * abs2(imag(scale_factor))}, + .imag_component = {.value = imag(scaled_value), + .variance = var.real_component.variance * abs2(imag(scale_factor)) + + var.imag_component.variance * abs2(real(scale_factor))}}; +} +} // namespace statistics } // namespace power_grid_model diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index f99003a5d..3f783b77a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -98,12 +98,6 @@ template class IterativeLinearSESolver { auto const observability_result = necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()); - // prepare matrix - sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization"); - prepare_matrix(y_bus, measured_values); - // prefactorize - sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation()); - // initialize voltage with initial angle sub_timer = Timer(calculation_info, 2223, "Initialize voltages"); RealValue const mean_angle_shift = measured_values.mean_angle_shift(); @@ -117,8 +111,21 @@ template class IterativeLinearSESolver { if (num_iter++ == max_iter) { throw IterationDiverge{max_iter, max_dev, err_tol}; } + + // prepare matrix + sub_timer = Timer(calculation_info, 2222, "Linearize measurements"); + // get generated (measured/estimated) voltage phasor + // with current result voltage angle + ComplexValueVector linearized_u = linearize_measurements(output.u, measured_values); + + // prepare matrix + sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization"); + prepare_matrix(y_bus, measured_values, linearized_u); + // prefactorize + sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation()); + sub_timer = Timer(calculation_info, 2224, "Calculate rhs"); - prepare_rhs(y_bus, measured_values, output.u); + prepare_rhs(y_bus, measured_values, linearized_u); // solve with prefactorization sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation (pre-factorized)"); sparse_solver_.solve_with_prefactorized_matrix(data_gain_, perm_, x_rhs_, x_rhs_); @@ -162,10 +169,12 @@ template class IterativeLinearSESolver { return ComplexDiagonalTensor{static_cast>(RealValue{1.0} / value)}; } - void prepare_matrix(YBus const& y_bus, MeasuredValues const& measured_value) { + void prepare_matrix(YBus const& y_bus, MeasuredValues const& measured_value, + ComplexValueVector const& linearized_u) { MathModelParam const& param = y_bus.math_model_param(); IdxVector const& row_indptr = y_bus.row_indptr_lu(); IdxVector const& col_indices = y_bus.col_indices_lu(); + std::vector const& branch_bus_idx = y_bus.math_topology().branch_bus_idx; // loop data index, all rows and columns for (Idx row = 0; row != n_bus_; ++row) { @@ -196,10 +205,12 @@ template class IterativeLinearSESolver { if (measured_value.has_shunt(obj)) { // G += (-Ys)^H * (variance^-1) * (-Ys) auto const& shunt_power = measured_value.shunt_power(obj); - block.g() += - dot(hermitian_transpose(param.shunt_param[obj]), - diagonal_inverse(static_cast>(shunt_power).variance), - param.shunt_param[obj]); + auto const shunt_current = + statistics::scale(shunt_power, ComplexValue{1.0 / linearized_u[row]}); + block.g() += dot( + hermitian_transpose(param.shunt_param[obj]), + diagonal_inverse(static_cast>(shunt_current).variance), + param.shunt_param[obj]); } } // branch @@ -213,9 +224,14 @@ template class IterativeLinearSESolver { if (std::invoke(has_branch_[measured_side], measured_value, obj)) { // G += Y{side, b0}^H * (variance^-1) * Y{side, b1} auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj); + // the current needs to be calculated with the voltage of the measured bus side + // NOTE: not the current bus! + Idx const measured_bus = branch_bus_idx[obj][measured_side]; + auto const current = + statistics::scale(power, ComplexValue{1.0 / linearized_u[measured_bus]}); block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]), - diagonal_inverse(static_cast>(power).variance), + diagonal_inverse(static_cast>(current).variance), param.branch_param[obj].value[measured_side * 2 + b1]); } } @@ -231,7 +247,9 @@ template class IterativeLinearSESolver { // assign variance to diagonal of 3x3 tensor, for asym auto const& injection = measured_value.bus_injection(row); block.r() = ComplexTensor{static_cast>( - -(static_cast>(injection).variance))}; + -(static_cast>( + statistics::scale(injection, ComplexValue{1.0 / linearized_u[row]})) + .variance))}; } } // injection measurement not exist @@ -259,12 +277,9 @@ template class IterativeLinearSESolver { } void prepare_rhs(YBus const& y_bus, MeasuredValues const& measured_value, - ComplexValueVector const& current_u) { + ComplexValueVector const& linearized_u) { MathModelParam const& param = y_bus.math_model_param(); - std::vector const branch_bus_idx = y_bus.math_topology().branch_bus_idx; - // get generated (measured/estimated) voltage phasor - // with current result voltage angle - ComplexValueVector u = linearize_measurements(current_u, measured_value); + std::vector const& branch_bus_idx = y_bus.math_topology().branch_bus_idx; // loop all bus to fill rhs for (Idx bus = 0; bus != n_bus_; ++bus) { @@ -275,7 +290,7 @@ template class IterativeLinearSESolver { // fill block with voltage measurement if (measured_value.has_voltage(bus)) { // eta += u / variance - rhs_block.eta() += u[bus] / measured_value.voltage_var(bus); + rhs_block.eta() += linearized_u[bus] / measured_value.voltage_var(bus); } // fill block with branch, shunt measurement, need to convert to current for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx]; @@ -285,12 +300,13 @@ template class IterativeLinearSESolver { // shunt if (type == YBusElementType::shunt) { if (measured_value.has_shunt(obj)) { - PowerSensorCalcParam const& m = measured_value.shunt_power(obj); + PowerSensorCalcParam const& power = measured_value.shunt_power(obj); + auto const current = statistics::scale(power, ComplexValue{1.0 / linearized_u[bus]}); + auto const rot_invariant_current = static_cast>(current); // eta += (-Ys)^H * (variance^-1) * i_shunt rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]), - diagonal_inverse(static_cast>(m).variance), - conj(m.value() / u[bus])); + diagonal_inverse(rot_invariant_current.variance), conj(rot_invariant_current.value)); } } // branch @@ -302,23 +318,27 @@ template class IterativeLinearSESolver { for (IntS const measured_side : std::array{0, 1}) { // has measurement if (std::invoke(has_branch_[measured_side], measured_value, obj)) { - PowerSensorCalcParam const& m = + PowerSensorCalcParam const& power = std::invoke(branch_power_[measured_side], measured_value, obj); // the current needs to be calculated with the voltage of the measured bus side // NOTE: not the current bus! Idx const measured_bus = branch_bus_idx[obj][measured_side]; + auto const current = + statistics::scale(power, ComplexValue{1.0 / linearized_u[measured_bus]}); + auto const rot_invariant_current = static_cast>(current); // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t} - rhs_block.eta() += - dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]), - diagonal_inverse(static_cast>(m).variance), - conj(m.value() / u[measured_bus])); + rhs_block.eta() += dot( + hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]), + diagonal_inverse(rot_invariant_current.variance), conj(rot_invariant_current.value)); } } } } // fill block with injection measurement, need to convert to current if (measured_value.has_bus_injection(bus)) { - rhs_block.tau() = conj(measured_value.bus_injection(bus).value() / u[bus]); + rhs_block.tau() = conj( + statistics::scale(measured_value.bus_injection(bus), ComplexValue{1.0 / linearized_u[bus]}) + .value()); } } } diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index 6a8462f73..d68040482 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -815,6 +815,249 @@ TEST_CASE("Test statistics") { } } +TEST_CASE("Test statistics - scale") { + using statistics::scale; + + SUBCASE("UniformRealRandVar | IndependentRealRandVar") { + auto const check = [] { + RandVar const var{.value = 1.0, .variance = 2.0}; + auto const scaled = scale(var, 3.0); + CHECK(scaled.value == doctest::Approx(3.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + }; + + SUBCASE("UniformRealRandVar") { check.template operator()>(); } + SUBCASE("IndependentRealRandVar") { + check.template operator()>(); + } + } + SUBCASE("UniformRealRandVar") { + UniformRealRandVar const var{.value = {1.0, 2.0, 3.0}, .variance = 2.0}; + SUBCASE("Scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(scaled.value(0) == doctest::Approx(3.0)); + CHECK(scaled.value(1) == doctest::Approx(6.0)); + CHECK(scaled.value(2) == doctest::Approx(9.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + } + // scaling asymmetrically would promote the UniformRealRandVar to an IndependentRealRandVar, because the + // individual phases scale differently + } + SUBCASE("IndependentRealRandVar") { + auto const var = IndependentRealRandVar{.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}}; + SUBCASE("Scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(scaled.value(0) == doctest::Approx(3.0)); + CHECK(scaled.value(1) == doctest::Approx(6.0)); + CHECK(scaled.value(2) == doctest::Approx(9.0)); + CHECK(scaled.variance(0) == doctest::Approx(18.0)); + CHECK(scaled.variance(1) == doctest::Approx(27.0)); + CHECK(scaled.variance(2) == doctest::Approx(36.0)); + } + SUBCASE("Asymmetric scale") { + auto const scaled = scale(var, RealValue{1.0, 2.0, 3.0}); + CHECK(scaled.value(0) == doctest::Approx(1.0)); + CHECK(scaled.value(1) == doctest::Approx(4.0)); + CHECK(scaled.value(2) == doctest::Approx(9.0)); + CHECK(scaled.variance(0) == doctest::Approx(2.0)); + CHECK(scaled.variance(1) == doctest::Approx(12.0)); + CHECK(scaled.variance(2) == doctest::Approx(36.0)); + } + } + SUBCASE("UniformComplexRandVar | IndependentComplexRandVar") { + auto const check = [] { + RandVar const var{.value = {1.0, 2.0}, .variance = 2.0}; + + SUBCASE("Real scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(real(scaled.value) == doctest::Approx(3.0)); + CHECK(imag(scaled.value) == doctest::Approx(6.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + } + SUBCASE("Complex scalar scale") { + auto const scaled = scale(var, DoubleComplex{3.0, 4.0}); + CHECK(real(scaled.value) == doctest::Approx(3.0 * 1.0 - 4.0 * 2.0)); + CHECK(imag(scaled.value) == doctest::Approx(3.0 * 2.0 + 4.0 * 1.0)); + CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0))); + } + }; + SUBCASE("UniformComplexRandVar") { + check.template operator()>(); + } + SUBCASE("IndependentComplexRandVar") { + check.template operator()>(); + } + } + SUBCASE("UniformComplexRandVar") { + UniformComplexRandVar const var{ + .value = {RealValue{1.0, 2.0, 3.0}, RealValue{4.0, 5.0, 6.0}}, .variance = 2.0}; + SUBCASE("Real scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(real(scaled.value(0)) == doctest::Approx(3.0)); + CHECK(real(scaled.value(1)) == doctest::Approx(6.0)); + CHECK(real(scaled.value(2)) == doctest::Approx(9.0)); + CHECK(imag(scaled.value(0)) == doctest::Approx(12.0)); + CHECK(imag(scaled.value(1)) == doctest::Approx(15.0)); + CHECK(imag(scaled.value(2)) == doctest::Approx(18.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + } + SUBCASE("Complex scalar scale") { + DoubleComplex const scale_factor{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + CHECK(real(scaled.value(i)) == + real(scale(UniformComplexRandVar{.value = var.value(i), .variance = var.variance}, + scale_factor) + .value)); + CHECK(real(scaled.value(i)) == + real(scale(UniformComplexRandVar{.value = var.value(i), .variance = var.variance}, + scale_factor) + .value)); + } + CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0))); + } + // scaling asymmetrically would promote the UniformComplexRandVar to an IndependentComplexRandVar, because the + // individual phases scale differently + } + SUBCASE("IndependentComplexRandVar") { + auto const var = IndependentComplexRandVar{ + .value = {RealValue{1.0, 2.0, 3.0}, RealValue{4.0, 5.0, 6.0}}, + .variance = {2.0, 3.0, 4.0}}; + auto const individual_phases = [&var] { + std::array, 3> result; + for (auto i = 0; i < 3; ++i) { + result[i] = {.value = var.value(i), .variance = var.variance(i)}; + } + return result; + }(); + SUBCASE("Real scalar scale") { + auto const scale_factor = 3.0; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + SUBCASE("Complex scalar scale") { + DoubleComplex const scale_factor{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + SUBCASE("Real asymmetric scale") { + RealValue const scale_factor{3.0, 4.0, 5.0}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor[i]); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + SUBCASE("Complex asymmetric scale") { + ComplexValue const scale_factor{RealValue{3.0, 4.0, 5.0}, + RealValue{6.0, 7.0, 8.0}}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor(i)); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + } + SUBCASE("DecomposedComplexRandVar") { + DecomposedComplexRandVar const var{.real_component = {.value = 1.0, .variance = 2.0}, + .imag_component = {.value = 4.0, .variance = 5.0}}; + + SUBCASE("Real scalar scale") { + auto const scale_factor = 3.0; + auto const scaled = scale(var, scale_factor); + CHECK(scaled.real_component.value == scale(var.real_component, scale_factor).value); + CHECK(scaled.imag_component.value == scale(var.imag_component, scale_factor).value); + CHECK(scaled.real_component.variance == scale(var.real_component, scale_factor).variance); + CHECK(scaled.imag_component.variance == scale(var.imag_component, scale_factor).variance); + } + SUBCASE("Complex scalar scale") { + auto const scale_factor = DoubleComplex{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + CHECK(scaled.real_component.value == real(var.value() * scale_factor)); + CHECK(scaled.imag_component.value == imag(var.value() * scale_factor)); + CHECK(scaled.real_component.variance == + real(scale_factor) * real(scale_factor) * var.real_component.variance + + imag(scale_factor) * imag(scale_factor) * var.imag_component.variance); + CHECK(scaled.imag_component.variance == + real(scale_factor) * real(scale_factor) * var.imag_component.variance + + imag(scale_factor) * imag(scale_factor) * var.real_component.variance); + } + } + SUBCASE("DecomposedComplexRandVar") { + DecomposedComplexRandVar const var{ + .real_component = {.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}}, + .imag_component = {.value = {4.0, 5.0, 6.0}, .variance = {5.0, 6.0, 7.0}}}; + + SUBCASE("Real scalar scale") { + auto const scale_factor = 3.0; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i)); + CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i)); + CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i)); + CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i)); + } + } + SUBCASE("Real asymmetric scale") { + auto const scale_factor = RealValue{3.0, 4.0, 5.0}; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i)); + CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i)); + CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i)); + CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i)); + } + } + SUBCASE("Complex scalar scale") { + auto const scale_factor = DoubleComplex{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor)); + CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor)); + CHECK(scaled.real_component.variance(i) == + real(scale_factor) * real(scale_factor) * var.real_component.variance(i) + + imag(scale_factor) * imag(scale_factor) * var.imag_component.variance(i)); + CHECK(scaled.imag_component.variance(i) == + real(scale_factor) * real(scale_factor) * var.imag_component.variance(i) + + imag(scale_factor) * imag(scale_factor) * var.real_component.variance(i)); + } + } + SUBCASE("Complex asymmetric scale") { + auto const scale_factor = ComplexValue{RealValue{3.0, 4.0, 5.0}, + RealValue{6.0, 7.0, 8.0}}; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor(i))); + CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor(i))); + CHECK(scaled.real_component.variance(i) == + real(scale_factor(i)) * real(scale_factor(i)) * var.real_component.variance(i) + + imag(scale_factor(i)) * imag(scale_factor(i)) * var.imag_component.variance(i)); + CHECK(scaled.imag_component.variance(i) == + real(scale_factor(i)) * real(scale_factor(i)) * var.imag_component.variance(i) + + imag(scale_factor(i)) * imag(scale_factor(i)) * var.real_component.variance(i)); + } + } + } +} + TEST_SUITE_END(); } // namespace power_grid_model From f57ccb3d850281572b674cdf5f2fe0b06c47816f Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 15:57:12 +0200 Subject: [PATCH 2/7] fix eigen operations Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 2 +- .../iterative_linear_se_solver.hpp | 39 +++++++++++-------- 2 files changed, 23 insertions(+), 18 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index d581dafdd..73f9ac80a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -240,7 +240,7 @@ inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& sca template requires is_in_list_c, ComplexValue> inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { - auto const scaled_value = var.value() * scale_factor; + ComplexValue const scaled_value = var.value() * scale_factor; return DecomposedComplexRandVar{ .real_component = {.value = real(scaled_value), .variance = var.real_component.variance * abs2(real(scale_factor)) + diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index 3f783b77a..0efc952e0 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -205,12 +205,13 @@ template class IterativeLinearSESolver { if (measured_value.has_shunt(obj)) { // G += (-Ys)^H * (variance^-1) * (-Ys) auto const& shunt_power = measured_value.shunt_power(obj); - auto const shunt_current = + auto const shunt_current_conj = statistics::scale(shunt_power, ComplexValue{1.0 / linearized_u[row]}); - block.g() += dot( - hermitian_transpose(param.shunt_param[obj]), - diagonal_inverse(static_cast>(shunt_current).variance), - param.shunt_param[obj]); + block.g() += + dot(hermitian_transpose(param.shunt_param[obj]), + diagonal_inverse( + static_cast>(shunt_current_conj).variance), + param.shunt_param[obj]); } } // branch @@ -227,11 +228,12 @@ template class IterativeLinearSESolver { // the current needs to be calculated with the voltage of the measured bus side // NOTE: not the current bus! Idx const measured_bus = branch_bus_idx[obj][measured_side]; - auto const current = + auto const current_conj = statistics::scale(power, ComplexValue{1.0 / linearized_u[measured_bus]}); block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]), - diagonal_inverse(static_cast>(current).variance), + diagonal_inverse( + static_cast>(current_conj).variance), param.branch_param[obj].value[measured_side * 2 + b1]); } } @@ -301,12 +303,13 @@ template class IterativeLinearSESolver { if (type == YBusElementType::shunt) { if (measured_value.has_shunt(obj)) { PowerSensorCalcParam const& power = measured_value.shunt_power(obj); - auto const current = statistics::scale(power, ComplexValue{1.0 / linearized_u[bus]}); - auto const rot_invariant_current = static_cast>(current); + auto const current_conj = statistics::scale(power, ComplexValue{1.0 / linearized_u[bus]}); + auto const rot_invariant_current_conj = + static_cast>(current_conj); // eta += (-Ys)^H * (variance^-1) * i_shunt - rhs_block.eta() -= - dot(hermitian_transpose(param.shunt_param[obj]), - diagonal_inverse(rot_invariant_current.variance), conj(rot_invariant_current.value)); + rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]), + diagonal_inverse(rot_invariant_current_conj.variance), + conj(rot_invariant_current_conj.value)); } } // branch @@ -323,13 +326,15 @@ template class IterativeLinearSESolver { // the current needs to be calculated with the voltage of the measured bus side // NOTE: not the current bus! Idx const measured_bus = branch_bus_idx[obj][measured_side]; - auto const current = + auto const current_conj = statistics::scale(power, ComplexValue{1.0 / linearized_u[measured_bus]}); - auto const rot_invariant_current = static_cast>(current); + auto const rot_invariant_current_conj = + static_cast>(current_conj); // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t} - rhs_block.eta() += dot( - hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]), - diagonal_inverse(rot_invariant_current.variance), conj(rot_invariant_current.value)); + rhs_block.eta() += + dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]), + diagonal_inverse(rot_invariant_current_conj.variance), + conj(rot_invariant_current_conj.value)); } } } From 076c6f8125a6f6ef217e8a1b57223e5579b26629 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 16:02:26 +0200 Subject: [PATCH 3/7] update tests to reflect migration to scaled power sensor variance Signed-off-by: Martijn Govers --- tests/data/state_estimation/line-power-sensor/params.json | 2 +- .../state_estimation/single-line-load-il/sym_output.json | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/data/state_estimation/line-power-sensor/params.json b/tests/data/state_estimation/line-power-sensor/params.json index 47b5ca68b..fb888eb3a 100644 --- a/tests/data/state_estimation/line-power-sensor/params.json +++ b/tests/data/state_estimation/line-power-sensor/params.json @@ -3,6 +3,6 @@ "rtol": 1e-03, "atol": { "default": 1e-05, - ".+_residual": 0.1 + ".+_residual": 0.2 } } \ No newline at end of file diff --git a/tests/data/state_estimation/single-line-load-il/sym_output.json b/tests/data/state_estimation/single-line-load-il/sym_output.json index 7e34500c8..c0f22327a 100644 --- a/tests/data/state_estimation/single-line-load-il/sym_output.json +++ b/tests/data/state_estimation/single-line-load-il/sym_output.json @@ -8,8 +8,8 @@ { "id": 21, "energized": 1, - "p_from": 1501363.636462, - "p_to": -1499319.107870 + "p_from": 1500680.11502645, + "p_to": -1498637.44762342 } ], "node": [ @@ -28,14 +28,14 @@ { "id": 41, "energized": 1, - "p": 1501363.636462 + "p": 1500680.11502645 } ], "sym_load": [ { "id": 31, "energized": 1, - "p": 1499319.107870, + "p": 1498637.44762342, "q": 0.0 } ] From e9766949ddd9971f6c6c799e4e35a0ee9f4be4c9 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 16:23:56 +0200 Subject: [PATCH 4/7] minor additional requirement Signed-off-by: Martijn Govers --- .../include/power_grid_model/common/statistics.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index 73f9ac80a..8c5349a4c 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -238,7 +238,8 @@ inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& sca } template - requires is_in_list_c, ComplexValue> + requires std::is_same> || + (is_asymmetric_v && std::same_as>) inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { ComplexValue const scaled_value = var.value() * scale_factor; return DecomposedComplexRandVar{ From f621ab210d3d23e1800c76dc6b23875bec3d4fc4 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 17:16:26 +0200 Subject: [PATCH 5/7] make clang happy Signed-off-by: Martijn Govers --- .../include/power_grid_model/common/statistics.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index 8c5349a4c..643a03c5a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -238,8 +238,8 @@ inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& sca } template - requires std::is_same> || - (is_asymmetric_v && std::same_as>) + requires(std::same_as> || + (is_asymmetric_v && std::same_as>)) inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { ComplexValue const scaled_value = var.value() * scale_factor; return DecomposedComplexRandVar{ From 83b4c48ed1a8d954d938278cb965a6da59fe6b3c Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Fri, 11 Apr 2025 09:11:26 +0200 Subject: [PATCH 6/7] improve naming Signed-off-by: Martijn Govers --- .../iterative_linear_se_solver.hpp | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index 0efc952e0..c9acf74d4 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -116,16 +116,16 @@ template class IterativeLinearSESolver { sub_timer = Timer(calculation_info, 2222, "Linearize measurements"); // get generated (measured/estimated) voltage phasor // with current result voltage angle - ComplexValueVector linearized_u = linearize_measurements(output.u, measured_values); + ComplexValueVector u = linearize_measurements(output.u, measured_values); // prepare matrix sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization"); - prepare_matrix(y_bus, measured_values, linearized_u); + prepare_matrix(y_bus, measured_values, u); // prefactorize sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation()); sub_timer = Timer(calculation_info, 2224, "Calculate rhs"); - prepare_rhs(y_bus, measured_values, linearized_u); + prepare_rhs(y_bus, measured_values, u); // solve with prefactorization sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation (pre-factorized)"); sparse_solver_.solve_with_prefactorized_matrix(data_gain_, perm_, x_rhs_, x_rhs_); @@ -170,7 +170,7 @@ template class IterativeLinearSESolver { } void prepare_matrix(YBus const& y_bus, MeasuredValues const& measured_value, - ComplexValueVector const& linearized_u) { + ComplexValueVector const& u) { MathModelParam const& param = y_bus.math_model_param(); IdxVector const& row_indptr = y_bus.row_indptr_lu(); IdxVector const& col_indices = y_bus.col_indices_lu(); @@ -206,7 +206,7 @@ template class IterativeLinearSESolver { // G += (-Ys)^H * (variance^-1) * (-Ys) auto const& shunt_power = measured_value.shunt_power(obj); auto const shunt_current_conj = - statistics::scale(shunt_power, ComplexValue{1.0 / linearized_u[row]}); + statistics::scale(shunt_power, ComplexValue{1.0 / u[row]}); block.g() += dot(hermitian_transpose(param.shunt_param[obj]), diagonal_inverse( @@ -229,7 +229,7 @@ template class IterativeLinearSESolver { // NOTE: not the current bus! Idx const measured_bus = branch_bus_idx[obj][measured_side]; auto const current_conj = - statistics::scale(power, ComplexValue{1.0 / linearized_u[measured_bus]}); + statistics::scale(power, ComplexValue{1.0 / u[measured_bus]}); block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]), diagonal_inverse( @@ -250,7 +250,7 @@ template class IterativeLinearSESolver { auto const& injection = measured_value.bus_injection(row); block.r() = ComplexTensor{static_cast>( -(static_cast>( - statistics::scale(injection, ComplexValue{1.0 / linearized_u[row]})) + statistics::scale(injection, ComplexValue{1.0 / u[row]})) .variance))}; } } @@ -279,7 +279,7 @@ template class IterativeLinearSESolver { } void prepare_rhs(YBus const& y_bus, MeasuredValues const& measured_value, - ComplexValueVector const& linearized_u) { + ComplexValueVector const& u) { MathModelParam const& param = y_bus.math_model_param(); std::vector const& branch_bus_idx = y_bus.math_topology().branch_bus_idx; @@ -292,7 +292,7 @@ template class IterativeLinearSESolver { // fill block with voltage measurement if (measured_value.has_voltage(bus)) { // eta += u / variance - rhs_block.eta() += linearized_u[bus] / measured_value.voltage_var(bus); + rhs_block.eta() += u[bus] / measured_value.voltage_var(bus); } // fill block with branch, shunt measurement, need to convert to current for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx]; @@ -303,7 +303,7 @@ template class IterativeLinearSESolver { if (type == YBusElementType::shunt) { if (measured_value.has_shunt(obj)) { PowerSensorCalcParam const& power = measured_value.shunt_power(obj); - auto const current_conj = statistics::scale(power, ComplexValue{1.0 / linearized_u[bus]}); + auto const current_conj = statistics::scale(power, ComplexValue{1.0 / u[bus]}); auto const rot_invariant_current_conj = static_cast>(current_conj); // eta += (-Ys)^H * (variance^-1) * i_shunt @@ -327,7 +327,7 @@ template class IterativeLinearSESolver { // NOTE: not the current bus! Idx const measured_bus = branch_bus_idx[obj][measured_side]; auto const current_conj = - statistics::scale(power, ComplexValue{1.0 / linearized_u[measured_bus]}); + statistics::scale(power, ComplexValue{1.0 / u[measured_bus]}); auto const rot_invariant_current_conj = static_cast>(current_conj); // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t} @@ -341,9 +341,8 @@ template class IterativeLinearSESolver { } // fill block with injection measurement, need to convert to current if (measured_value.has_bus_injection(bus)) { - rhs_block.tau() = conj( - statistics::scale(measured_value.bus_injection(bus), ComplexValue{1.0 / linearized_u[bus]}) - .value()); + rhs_block.tau() = + conj(statistics::scale(measured_value.bus_injection(bus), ComplexValue{1.0 / u[bus]}).value()); } } } From a71a4eeb05c8d87614f084b756f22a0f1e012f14 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 14 Apr 2025 15:54:17 +0200 Subject: [PATCH 7/7] add scaling functionality to statistics.hpp Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 45 ++++ tests/cpp_unit_tests/test_statistics.cpp | 243 ++++++++++++++++++ 2 files changed, 288 insertions(+) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index 247a9b676..55fcd1ff3 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -208,6 +208,51 @@ template struct PolarComplexRandVar { }; namespace statistics { +template typename RandVarType> + requires is_in_list_c, UniformRealRandVar, IndependentRealRandVar, + UniformComplexRandVar, IndependentComplexRandVar> +inline auto scale(RandVarType const& var, double scale_factor) { + return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)}; +} + +template + requires is_in_list_c, IndependentComplexRandVar> +inline auto scale(RandVarType const& var, RealValue const& scale_factor) { + return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)}; +} + +template typename RandVarType> + requires is_in_list_c, UniformComplexRandVar, IndependentComplexRandVar> +inline auto scale(RandVarType const& var, DoubleComplex const& scale_factor) { + return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)}; +} + +inline auto scale(IndependentComplexRandVar const& var, ComplexValue const& scale_factor) { + return IndependentComplexRandVar{.value = var.value * scale_factor, + .variance = var.variance * abs2(scale_factor)}; +} + +template + requires is_in_list_c, RealValue> +inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { + return DecomposedComplexRandVar{.real_component = scale(var.real_component, scale_factor), + .imag_component = scale(var.imag_component, scale_factor)}; +} + +template + requires(std::same_as> || + (is_asymmetric_v && std::same_as>)) +inline auto scale(DecomposedComplexRandVar const& var, ScaleType const& scale_factor) { + ComplexValue const scaled_value = var.value() * scale_factor; + return DecomposedComplexRandVar{ + .real_component = {.value = real(scaled_value), + .variance = var.real_component.variance * abs2(real(scale_factor)) + + var.imag_component.variance * abs2(imag(scale_factor))}, + .imag_component = {.value = imag(scaled_value), + .variance = var.real_component.variance * abs2(imag(scale_factor)) + + var.imag_component.variance * abs2(real(scale_factor))}}; +} + // combine multiple random variables of one quantity using Kalman filter template requires std::same_as, diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index 4df368497..504433959 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -815,6 +815,249 @@ TEST_CASE("Test statistics") { } } +TEST_CASE("Test statistics - scale") { + using statistics::scale; + + SUBCASE("UniformRealRandVar | IndependentRealRandVar") { + auto const check = [] { + RandVar const var{.value = 1.0, .variance = 2.0}; + auto const scaled = scale(var, 3.0); + CHECK(scaled.value == doctest::Approx(3.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + }; + + SUBCASE("UniformRealRandVar") { check.template operator()>(); } + SUBCASE("IndependentRealRandVar") { + check.template operator()>(); + } + } + SUBCASE("UniformRealRandVar") { + UniformRealRandVar const var{.value = {1.0, 2.0, 3.0}, .variance = 2.0}; + SUBCASE("Scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(scaled.value(0) == doctest::Approx(3.0)); + CHECK(scaled.value(1) == doctest::Approx(6.0)); + CHECK(scaled.value(2) == doctest::Approx(9.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + } + // scaling asymmetrically would promote the UniformRealRandVar to an IndependentRealRandVar, because the + // individual phases scale differently + } + SUBCASE("IndependentRealRandVar") { + auto const var = IndependentRealRandVar{.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}}; + SUBCASE("Scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(scaled.value(0) == doctest::Approx(3.0)); + CHECK(scaled.value(1) == doctest::Approx(6.0)); + CHECK(scaled.value(2) == doctest::Approx(9.0)); + CHECK(scaled.variance(0) == doctest::Approx(18.0)); + CHECK(scaled.variance(1) == doctest::Approx(27.0)); + CHECK(scaled.variance(2) == doctest::Approx(36.0)); + } + SUBCASE("Asymmetric scale") { + auto const scaled = scale(var, RealValue{1.0, 2.0, 3.0}); + CHECK(scaled.value(0) == doctest::Approx(1.0)); + CHECK(scaled.value(1) == doctest::Approx(4.0)); + CHECK(scaled.value(2) == doctest::Approx(9.0)); + CHECK(scaled.variance(0) == doctest::Approx(2.0)); + CHECK(scaled.variance(1) == doctest::Approx(12.0)); + CHECK(scaled.variance(2) == doctest::Approx(36.0)); + } + } + SUBCASE("UniformComplexRandVar | IndependentComplexRandVar") { + auto const check = [] { + RandVar const var{.value = {1.0, 2.0}, .variance = 2.0}; + + SUBCASE("Real scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(real(scaled.value) == doctest::Approx(3.0)); + CHECK(imag(scaled.value) == doctest::Approx(6.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + } + SUBCASE("Complex scalar scale") { + auto const scaled = scale(var, DoubleComplex{3.0, 4.0}); + CHECK(real(scaled.value) == doctest::Approx(3.0 * 1.0 - 4.0 * 2.0)); + CHECK(imag(scaled.value) == doctest::Approx(3.0 * 2.0 + 4.0 * 1.0)); + CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0))); + } + }; + SUBCASE("UniformComplexRandVar") { + check.template operator()>(); + } + SUBCASE("IndependentComplexRandVar") { + check.template operator()>(); + } + } + SUBCASE("UniformComplexRandVar") { + UniformComplexRandVar const var{ + .value = {RealValue{1.0, 2.0, 3.0}, RealValue{4.0, 5.0, 6.0}}, .variance = 2.0}; + SUBCASE("Real scalar scale") { + auto const scaled = scale(var, 3.0); + CHECK(real(scaled.value(0)) == doctest::Approx(3.0)); + CHECK(real(scaled.value(1)) == doctest::Approx(6.0)); + CHECK(real(scaled.value(2)) == doctest::Approx(9.0)); + CHECK(imag(scaled.value(0)) == doctest::Approx(12.0)); + CHECK(imag(scaled.value(1)) == doctest::Approx(15.0)); + CHECK(imag(scaled.value(2)) == doctest::Approx(18.0)); + CHECK(scaled.variance == doctest::Approx(18.0)); + } + SUBCASE("Complex scalar scale") { + DoubleComplex const scale_factor{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + CHECK(real(scaled.value(i)) == + real(scale(UniformComplexRandVar{.value = var.value(i), .variance = var.variance}, + scale_factor) + .value)); + CHECK(real(scaled.value(i)) == + real(scale(UniformComplexRandVar{.value = var.value(i), .variance = var.variance}, + scale_factor) + .value)); + } + CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0))); + } + // scaling asymmetrically would promote the UniformComplexRandVar to an IndependentComplexRandVar, because the + // individual phases scale differently + } + SUBCASE("IndependentComplexRandVar") { + auto const var = IndependentComplexRandVar{ + .value = {RealValue{1.0, 2.0, 3.0}, RealValue{4.0, 5.0, 6.0}}, + .variance = {2.0, 3.0, 4.0}}; + auto const individual_phases = [&var] { + std::array, 3> result; + for (auto i = 0; i < 3; ++i) { + result[i] = {.value = var.value(i), .variance = var.variance(i)}; + } + return result; + }(); + SUBCASE("Real scalar scale") { + auto const scale_factor = 3.0; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + SUBCASE("Complex scalar scale") { + DoubleComplex const scale_factor{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + SUBCASE("Real asymmetric scale") { + RealValue const scale_factor{3.0, 4.0, 5.0}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor[i]); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + SUBCASE("Complex asymmetric scale") { + ComplexValue const scale_factor{RealValue{3.0, 4.0, 5.0}, + RealValue{6.0, 7.0, 8.0}}; + auto const scaled = scale(var, scale_factor); + for (auto i = 0; i < 3; ++i) { + auto const expected = scale(individual_phases[i], scale_factor(i)); + CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value))); + CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value))); + CHECK(scaled.variance(i) == doctest::Approx(expected.variance)); + } + } + } + SUBCASE("DecomposedComplexRandVar") { + DecomposedComplexRandVar const var{.real_component = {.value = 1.0, .variance = 2.0}, + .imag_component = {.value = 4.0, .variance = 5.0}}; + + SUBCASE("Real scalar scale") { + auto const scale_factor = 3.0; + auto const scaled = scale(var, scale_factor); + CHECK(scaled.real_component.value == scale(var.real_component, scale_factor).value); + CHECK(scaled.imag_component.value == scale(var.imag_component, scale_factor).value); + CHECK(scaled.real_component.variance == scale(var.real_component, scale_factor).variance); + CHECK(scaled.imag_component.variance == scale(var.imag_component, scale_factor).variance); + } + SUBCASE("Complex scalar scale") { + auto const scale_factor = DoubleComplex{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + CHECK(scaled.real_component.value == real(var.value() * scale_factor)); + CHECK(scaled.imag_component.value == imag(var.value() * scale_factor)); + CHECK(scaled.real_component.variance == + real(scale_factor) * real(scale_factor) * var.real_component.variance + + imag(scale_factor) * imag(scale_factor) * var.imag_component.variance); + CHECK(scaled.imag_component.variance == + real(scale_factor) * real(scale_factor) * var.imag_component.variance + + imag(scale_factor) * imag(scale_factor) * var.real_component.variance); + } + } + SUBCASE("DecomposedComplexRandVar") { + DecomposedComplexRandVar const var{ + .real_component = {.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}}, + .imag_component = {.value = {4.0, 5.0, 6.0}, .variance = {5.0, 6.0, 7.0}}}; + + SUBCASE("Real scalar scale") { + auto const scale_factor = 3.0; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i)); + CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i)); + CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i)); + CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i)); + } + } + SUBCASE("Real asymmetric scale") { + auto const scale_factor = RealValue{3.0, 4.0, 5.0}; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i)); + CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i)); + CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i)); + CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i)); + } + } + SUBCASE("Complex scalar scale") { + auto const scale_factor = DoubleComplex{3.0, 4.0}; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor)); + CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor)); + CHECK(scaled.real_component.variance(i) == + real(scale_factor) * real(scale_factor) * var.real_component.variance(i) + + imag(scale_factor) * imag(scale_factor) * var.imag_component.variance(i)); + CHECK(scaled.imag_component.variance(i) == + real(scale_factor) * real(scale_factor) * var.imag_component.variance(i) + + imag(scale_factor) * imag(scale_factor) * var.real_component.variance(i)); + } + } + SUBCASE("Complex asymmetric scale") { + auto const scale_factor = ComplexValue{RealValue{3.0, 4.0, 5.0}, + RealValue{6.0, 7.0, 8.0}}; + auto const scaled = scale(var, scale_factor); + + for (int i = 0; i < 3; ++i) { + CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor(i))); + CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor(i))); + CHECK(scaled.real_component.variance(i) == + real(scale_factor(i)) * real(scale_factor(i)) * var.real_component.variance(i) + + imag(scale_factor(i)) * imag(scale_factor(i)) * var.imag_component.variance(i)); + CHECK(scaled.imag_component.variance(i) == + real(scale_factor(i)) * real(scale_factor(i)) * var.imag_component.variance(i) + + imag(scale_factor(i)) * imag(scale_factor(i)) * var.real_component.variance(i)); + } + } + } +} + TEST_CASE("Test statistics - combine") { using statistics::combine; using std::views::take;