Skip to content

Commit a206d17

Browse files
mgoversmgovers
committed
add scaling functionality to statistics.hpp
Signed-off-by: Martijn Govers <Martijn.Govers@Alliander.com>
1 parent 4ab0e9c commit a206d17

File tree

2 files changed

+289
-0
lines changed

2 files changed

+289
-0
lines changed

‎power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -208,6 +208,52 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
208208
};
209209

210210
namespace statistics {
211+
// Var(s x) ≈ Var(x) * ||s||²
212+
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
213+
requires is_in_list_c<RandVarType<sym>, UniformRealRandVar<sym>, IndependentRealRandVar<sym>,
214+
UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>
215+
inline auto scale(RandVarType<sym> const& var, double scale_factor) {
216+
return RandVarType<sym>{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
217+
}
218+
219+
template <typename RandVarType>
220+
requires is_in_list_c<RandVarType, IndependentRealRandVar<asymmetric_t>, IndependentComplexRandVar<asymmetric_t>>
221+
inline auto scale(RandVarType const& var, RealValue<asymmetric_t> const& scale_factor) {
222+
return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
223+
}
224+
225+
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
226+
requires is_in_list_c<RandVarType<sym>, UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>
227+
inline auto scale(RandVarType<sym> const& var, DoubleComplex const& scale_factor) {
228+
return RandVarType<sym>{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
229+
}
230+
231+
inline auto scale(IndependentComplexRandVar<asymmetric_t> const& var, ComplexValue<asymmetric_t> const& scale_factor) {
232+
return IndependentComplexRandVar<asymmetric_t>{.value = var.value * scale_factor,
233+
.variance = var.variance * abs2(scale_factor)};
234+
}
235+
236+
template <symmetry_tag sym, typename ScaleType>
237+
requires is_in_list_c<ScaleType, RealValue<symmetric_t>, RealValue<asymmetric_t>>
238+
inline auto scale(DecomposedComplexRandVar<sym> const& var, ScaleType const& scale_factor) {
239+
return DecomposedComplexRandVar<sym>{.real_component = scale(var.real_component, scale_factor),
240+
.imag_component = scale(var.imag_component, scale_factor)};
241+
}
242+
243+
template <symmetry_tag sym, typename ScaleType>
244+
requires(std::same_as<ScaleType, ComplexValue<symmetric_t>> ||
245+
(is_asymmetric_v<sym> && std::same_as<ScaleType, ComplexValue<asymmetric_t>>))
246+
inline auto scale(DecomposedComplexRandVar<sym> const& var, ScaleType const& scale_factor) {
247+
ComplexValue<sym> const scaled_value = var.value() * scale_factor;
248+
return DecomposedComplexRandVar<sym>{
249+
.real_component = {.value = real(scaled_value),
250+
.variance = var.real_component.variance * abs2(real(scale_factor)) +
251+
var.imag_component.variance * abs2(imag(scale_factor))},
252+
.imag_component = {.value = imag(scaled_value),
253+
.variance = var.real_component.variance * abs2(imag(scale_factor)) +
254+
var.imag_component.variance * abs2(real(scale_factor))}};
255+
}
256+
211257
// combine multiple random variables of one quantity using Kalman filter
212258
template <std::ranges::view RandVarsView>
213259
requires std::same_as<std::ranges::range_value_t<RandVarsView>,

‎tests/cpp_unit_tests/test_statistics.cpp

Lines changed: 243 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -815,6 +815,249 @@ TEST_CASE("Test statistics") {
815815
}
816816
}
817817

818+
TEST_CASE("Test statistics - scale") {
819+
using statistics::scale;
820+
821+
SUBCASE("UniformRealRandVar<symmetric_t> | IndependentRealRandVar<symmetric_t>") {
822+
auto const check = []<typename RandVar> {
823+
RandVar const var{.value = 1.0, .variance = 2.0};
824+
auto const scaled = scale(var, 3.0);
825+
CHECK(scaled.value == doctest::Approx(3.0));
826+
CHECK(scaled.variance == doctest::Approx(18.0));
827+
};
828+
829+
SUBCASE("UniformRealRandVar<symmetric_t>") { check.template operator()<UniformRealRandVar<symmetric_t>>(); }
830+
SUBCASE("IndependentRealRandVar<symmetric_t>") {
831+
check.template operator()<IndependentRealRandVar<symmetric_t>>();
832+
}
833+
}
834+
SUBCASE("UniformRealRandVar<asymmetric_t>") {
835+
UniformRealRandVar<asymmetric_t> const var{.value = {1.0, 2.0, 3.0}, .variance = 2.0};
836+
SUBCASE("Scalar scale") {
837+
auto const scaled = scale(var, 3.0);
838+
CHECK(scaled.value(0) == doctest::Approx(3.0));
839+
CHECK(scaled.value(1) == doctest::Approx(6.0));
840+
CHECK(scaled.value(2) == doctest::Approx(9.0));
841+
CHECK(scaled.variance == doctest::Approx(18.0));
842+
}
843+
// scaling asymmetrically would promote the UniformRealRandVar to an IndependentRealRandVar, because the
844+
// individual phases scale differently
845+
}
846+
SUBCASE("IndependentRealRandVar<asymmetric_t>") {
847+
auto const var = IndependentRealRandVar<asymmetric_t>{.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}};
848+
SUBCASE("Scalar scale") {
849+
auto const scaled = scale(var, 3.0);
850+
CHECK(scaled.value(0) == doctest::Approx(3.0));
851+
CHECK(scaled.value(1) == doctest::Approx(6.0));
852+
CHECK(scaled.value(2) == doctest::Approx(9.0));
853+
CHECK(scaled.variance(0) == doctest::Approx(18.0));
854+
CHECK(scaled.variance(1) == doctest::Approx(27.0));
855+
CHECK(scaled.variance(2) == doctest::Approx(36.0));
856+
}
857+
SUBCASE("Asymmetric scale") {
858+
auto const scaled = scale(var, RealValue<asymmetric_t>{1.0, 2.0, 3.0});
859+
CHECK(scaled.value(0) == doctest::Approx(1.0));
860+
CHECK(scaled.value(1) == doctest::Approx(4.0));
861+
CHECK(scaled.value(2) == doctest::Approx(9.0));
862+
CHECK(scaled.variance(0) == doctest::Approx(2.0));
863+
CHECK(scaled.variance(1) == doctest::Approx(12.0));
864+
CHECK(scaled.variance(2) == doctest::Approx(36.0));
865+
}
866+
}
867+
SUBCASE("UniformComplexRandVar<symmetric_t> | IndependentComplexRandVar<symmetric_t>") {
868+
auto const check = []<typename RandVar> {
869+
RandVar const var{.value = {1.0, 2.0}, .variance = 2.0};
870+
871+
SUBCASE("Real scalar scale") {
872+
auto const scaled = scale(var, 3.0);
873+
CHECK(real(scaled.value) == doctest::Approx(3.0));
874+
CHECK(imag(scaled.value) == doctest::Approx(6.0));
875+
CHECK(scaled.variance == doctest::Approx(18.0));
876+
}
877+
SUBCASE("Complex scalar scale") {
878+
auto const scaled = scale(var, DoubleComplex{3.0, 4.0});
879+
CHECK(real(scaled.value) == doctest::Approx(3.0 * 1.0 - 4.0 * 2.0));
880+
CHECK(imag(scaled.value) == doctest::Approx(3.0 * 2.0 + 4.0 * 1.0));
881+
CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0)));
882+
}
883+
};
884+
SUBCASE("UniformComplexRandVar<symmetric_t>") {
885+
check.template operator()<UniformComplexRandVar<symmetric_t>>();
886+
}
887+
SUBCASE("IndependentComplexRandVar<symmetric_t>") {
888+
check.template operator()<IndependentComplexRandVar<symmetric_t>>();
889+
}
890+
}
891+
SUBCASE("UniformComplexRandVar<asymmetric_t>") {
892+
UniformComplexRandVar<asymmetric_t> const var{
893+
.value = {RealValue<asymmetric_t>{1.0, 2.0, 3.0}, RealValue<asymmetric_t>{4.0, 5.0, 6.0}}, .variance = 2.0};
894+
SUBCASE("Real scalar scale") {
895+
auto const scaled = scale(var, 3.0);
896+
CHECK(real(scaled.value(0)) == doctest::Approx(3.0));
897+
CHECK(real(scaled.value(1)) == doctest::Approx(6.0));
898+
CHECK(real(scaled.value(2)) == doctest::Approx(9.0));
899+
CHECK(imag(scaled.value(0)) == doctest::Approx(12.0));
900+
CHECK(imag(scaled.value(1)) == doctest::Approx(15.0));
901+
CHECK(imag(scaled.value(2)) == doctest::Approx(18.0));
902+
CHECK(scaled.variance == doctest::Approx(18.0));
903+
}
904+
SUBCASE("Complex scalar scale") {
905+
DoubleComplex const scale_factor{3.0, 4.0};
906+
auto const scaled = scale(var, scale_factor);
907+
for (auto i = 0; i < 3; ++i) {
908+
CHECK(real(scaled.value(i)) ==
909+
real(scale(UniformComplexRandVar<symmetric_t>{.value = var.value(i), .variance = var.variance},
910+
scale_factor)
911+
.value));
912+
CHECK(real(scaled.value(i)) ==
913+
real(scale(UniformComplexRandVar<symmetric_t>{.value = var.value(i), .variance = var.variance},
914+
scale_factor)
915+
.value));
916+
}
917+
CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0)));
918+
}
919+
// scaling asymmetrically would promote the UniformComplexRandVar to an IndependentComplexRandVar, because the
920+
// individual phases scale differently
921+
}
922+
SUBCASE("IndependentComplexRandVar<asymmetric_t>") {
923+
auto const var = IndependentComplexRandVar<asymmetric_t>{
924+
.value = {RealValue<asymmetric_t>{1.0, 2.0, 3.0}, RealValue<asymmetric_t>{4.0, 5.0, 6.0}},
925+
.variance = {2.0, 3.0, 4.0}};
926+
auto const individual_phases = [&var] {
927+
std::array<UniformComplexRandVar<symmetric_t>, 3> result;
928+
for (auto i = 0; i < 3; ++i) {
929+
result[i] = {.value = var.value(i), .variance = var.variance(i)};
930+
}
931+
return result;
932+
}();
933+
SUBCASE("Real scalar scale") {
934+
auto const scale_factor = 3.0;
935+
auto const scaled = scale(var, scale_factor);
936+
for (auto i = 0; i < 3; ++i) {
937+
auto const expected = scale(individual_phases[i], scale_factor);
938+
CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
939+
CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
940+
CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
941+
}
942+
}
943+
SUBCASE("Complex scalar scale") {
944+
DoubleComplex const scale_factor{3.0, 4.0};
945+
auto const scaled = scale(var, scale_factor);
946+
for (auto i = 0; i < 3; ++i) {
947+
auto const expected = scale(individual_phases[i], scale_factor);
948+
CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
949+
CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
950+
CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
951+
}
952+
}
953+
SUBCASE("Real asymmetric scale") {
954+
RealValue<asymmetric_t> const scale_factor{3.0, 4.0, 5.0};
955+
auto const scaled = scale(var, scale_factor);
956+
for (auto i = 0; i < 3; ++i) {
957+
auto const expected = scale(individual_phases[i], scale_factor[i]);
958+
CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
959+
CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
960+
CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
961+
}
962+
}
963+
SUBCASE("Complex asymmetric scale") {
964+
ComplexValue<asymmetric_t> const scale_factor{RealValue<asymmetric_t>{3.0, 4.0, 5.0},
965+
RealValue<asymmetric_t>{6.0, 7.0, 8.0}};
966+
auto const scaled = scale(var, scale_factor);
967+
for (auto i = 0; i < 3; ++i) {
968+
auto const expected = scale(individual_phases[i], scale_factor(i));
969+
CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
970+
CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
971+
CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
972+
}
973+
}
974+
}
975+
SUBCASE("DecomposedComplexRandVar<symmetric_t>") {
976+
DecomposedComplexRandVar<symmetric_t> const var{.real_component = {.value = 1.0, .variance = 2.0},
977+
.imag_component = {.value = 4.0, .variance = 5.0}};
978+
979+
SUBCASE("Real scalar scale") {
980+
auto const scale_factor = 3.0;
981+
auto const scaled = scale(var, scale_factor);
982+
CHECK(scaled.real_component.value == scale(var.real_component, scale_factor).value);
983+
CHECK(scaled.imag_component.value == scale(var.imag_component, scale_factor).value);
984+
CHECK(scaled.real_component.variance == scale(var.real_component, scale_factor).variance);
985+
CHECK(scaled.imag_component.variance == scale(var.imag_component, scale_factor).variance);
986+
}
987+
SUBCASE("Complex scalar scale") {
988+
auto const scale_factor = DoubleComplex{3.0, 4.0};
989+
auto const scaled = scale(var, scale_factor);
990+
CHECK(scaled.real_component.value == real(var.value() * scale_factor));
991+
CHECK(scaled.imag_component.value == imag(var.value() * scale_factor));
992+
CHECK(scaled.real_component.variance ==
993+
real(scale_factor) * real(scale_factor) * var.real_component.variance +
994+
imag(scale_factor) * imag(scale_factor) * var.imag_component.variance);
995+
CHECK(scaled.imag_component.variance ==
996+
real(scale_factor) * real(scale_factor) * var.imag_component.variance +
997+
imag(scale_factor) * imag(scale_factor) * var.real_component.variance);
998+
}
999+
}
1000+
SUBCASE("DecomposedComplexRandVar<asymmetric_t>") {
1001+
DecomposedComplexRandVar<asymmetric_t> const var{
1002+
.real_component = {.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}},
1003+
.imag_component = {.value = {4.0, 5.0, 6.0}, .variance = {5.0, 6.0, 7.0}}};
1004+
1005+
SUBCASE("Real scalar scale") {
1006+
auto const scale_factor = 3.0;
1007+
auto const scaled = scale(var, scale_factor);
1008+
1009+
for (int i = 0; i < 3; ++i) {
1010+
CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i));
1011+
CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i));
1012+
CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i));
1013+
CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i));
1014+
}
1015+
}
1016+
SUBCASE("Real asymmetric scale") {
1017+
auto const scale_factor = RealValue<asymmetric_t>{3.0, 4.0, 5.0};
1018+
auto const scaled = scale(var, scale_factor);
1019+
1020+
for (int i = 0; i < 3; ++i) {
1021+
CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i));
1022+
CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i));
1023+
CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i));
1024+
CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i));
1025+
}
1026+
}
1027+
SUBCASE("Complex scalar scale") {
1028+
auto const scale_factor = DoubleComplex{3.0, 4.0};
1029+
auto const scaled = scale(var, scale_factor);
1030+
1031+
for (int i = 0; i < 3; ++i) {
1032+
CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor));
1033+
CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor));
1034+
CHECK(scaled.real_component.variance(i) ==
1035+
real(scale_factor) * real(scale_factor) * var.real_component.variance(i) +
1036+
imag(scale_factor) * imag(scale_factor) * var.imag_component.variance(i));
1037+
CHECK(scaled.imag_component.variance(i) ==
1038+
real(scale_factor) * real(scale_factor) * var.imag_component.variance(i) +
1039+
imag(scale_factor) * imag(scale_factor) * var.real_component.variance(i));
1040+
}
1041+
}
1042+
SUBCASE("Complex asymmetric scale") {
1043+
auto const scale_factor = ComplexValue<asymmetric_t>{RealValue<asymmetric_t>{3.0, 4.0, 5.0},
1044+
RealValue<asymmetric_t>{6.0, 7.0, 8.0}};
1045+
auto const scaled = scale(var, scale_factor);
1046+
1047+
for (int i = 0; i < 3; ++i) {
1048+
CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor(i)));
1049+
CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor(i)));
1050+
CHECK(scaled.real_component.variance(i) ==
1051+
real(scale_factor(i)) * real(scale_factor(i)) * var.real_component.variance(i) +
1052+
imag(scale_factor(i)) * imag(scale_factor(i)) * var.imag_component.variance(i));
1053+
CHECK(scaled.imag_component.variance(i) ==
1054+
real(scale_factor(i)) * real(scale_factor(i)) * var.imag_component.variance(i) +
1055+
imag(scale_factor(i)) * imag(scale_factor(i)) * var.real_component.variance(i));
1056+
}
1057+
}
1058+
}
1059+
}
1060+
8181061
TEST_CASE("Test statistics - combine") {
8191062
using statistics::combine;
8201063
using std::views::take;

0 commit comments

Comments
 (0)