Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
main_optimizer_deeper.cpp
Go to the documentation of this file.
1#include <iostream>
2#include <chrono>
3#include <memory>
4#include <unordered_map>
5#include <vector>
6#include <cassert>
7#include <cmath>
8
9#include "Block.h"
10#include "BlockAccessor.h"
11#include "DependentParameter.h"
12#include "ParamOptimizer.h"
13#include "SourcesView.h"
14
15static std::shared_ptr<Parameter> make_param(const std::string& block, const LhaID& id, double v) {
16 return std::make_shared<Parameter>(ParamId(block, id), v, 0., 0.);
17}
18
19int main() {
20 using clock = std::chrono::steady_clock;
21
22 const int N = 600;
23 const int OPS = 300;
24 const int i1 = 5, i2 = 77, j1 = 123, j2 = 451;
25
26 auto SRC_A = std::make_shared<Block>(); SRC_A->blockname = "SRC_A"; SRC_A->set_scale(1.0);
27 auto SRC_B = std::make_shared<Block>(); SRC_B->blockname = "SRC_B"; SRC_B->set_scale(1.0);
28 auto SRC_C = std::make_shared<Block>(); SRC_C->blockname = "SRC_C"; SRC_C->set_scale(1.0);
29
30 auto SCALE = std::make_shared<Block>(); SCALE->blockname = "SCALE"; SCALE->set_scale(1.0);
31 SCALE->store(LhaID(1), make_param("SCALE", LhaID(1), 1.25)); // scale factor
32
33 for (int i = 0; i < N; ++i) {
34 LhaID id(i);
35 SRC_A->store(id, make_param("SRC_A", id, 1.0*i));
36 SRC_B->store(id, make_param("SRC_B", id, 2.0*i));
37 SRC_C->store(id, make_param("SRC_C", id, 0.5*i));
38 }
39
40 int upd_FUSED_AB = 0, upd_FUSED_BC = 0, upd_SCALED_A = 0, upd_MIXED = 0, upd_REDUCED = 0;
41 int upd_DP1 = 0, upd_DP2 = 0, upd_DP3 = 0, upd_DP4 = 0;
42
43
44 auto FUSED_AB = std::make_shared<DependentBlock>(
45 std::unordered_map<std::string, std::shared_ptr<Block>>{
46 {"SRC_A", SRC_A}, {"SRC_B", SRC_B}
47 },
48 [&upd_FUSED_AB, N](const BlockSrc& srcs, std::shared_ptr<DependentBlock> self) {
49 ++upd_FUSED_AB;
50 for (int i = 0; i < N; ++i) {
51 LhaID id(i);
52 double a = srcs.get_val("SRC_A", id);
53 double b = srcs.get_val("SRC_B", id);
54 double z = a + 2.0*b;
55 if (self->contains(id)) self->assign(id, z);
56 else self->store(id, std::make_shared<Parameter>(ParamId(self->get_name(), id), z, 0., 0.));
57 }
58 }
59 );
60 FUSED_AB->blockname = "FUSED_AB"; FUSED_AB->set_scale(1.0); FUSED_AB->init();
61
62 auto FUSED_BC = std::make_shared<DependentBlock>(
63 std::unordered_map<std::string, std::shared_ptr<Block>>{
64 {"SRC_B", SRC_B}, {"SRC_C", SRC_C}
65 },
66 [&upd_FUSED_BC, N](const BlockSrc& srcs, std::shared_ptr<DependentBlock> self) {
67 ++upd_FUSED_BC;
68 for (int i = 0; i < N; ++i) {
69 LhaID id(i);
70 double b = srcs.get_val("SRC_B", id);
71 double c = srcs.get_val("SRC_C", id);
72 double y = b - c;
73 if (self->contains(id)) self->assign(id, y);
74 else self->store(id, std::make_shared<Parameter>(ParamId(self->get_name(), id), y, 0., 0.));
75 }
76 }
77 );
78 FUSED_BC->blockname = "FUSED_BC"; FUSED_BC->set_scale(1.0); FUSED_BC->init();
79
80 auto SCALED_A = std::make_shared<DependentBlock>(
81 std::unordered_map<std::string, std::shared_ptr<Block>>{
82 {"SRC_A", SRC_A}, {"SCALE", SCALE}
83 },
84 [&upd_SCALED_A, N](const BlockSrc& srcs, std::shared_ptr<DependentBlock> self) {
85 ++upd_SCALED_A;
86 double sf = srcs.get_val("SCALE", 1);
87 for (int i = 0; i < N; ++i) {
88 LhaID id(i);
89 double a = srcs.get_val("SRC_A", id);
90 double s = sf * a;
91 if (self->contains(id)) self->assign(id, s);
92 else self->store(id, std::make_shared<Parameter>(ParamId(self->get_name(), id), s, 0., 0.));
93 }
94 }
95 );
96 SCALED_A->blockname = "SCALED_A"; SCALED_A->set_scale(1.0); SCALED_A->init();
97
98 auto MIXED = std::make_shared<DependentBlock>(
99 std::unordered_map<std::string, std::shared_ptr<Block>>{
100 {"SRC_A", SRC_A}, {"FUSED_AB", FUSED_AB}, {"FUSED_BC", FUSED_BC}
101 },
102 [&upd_MIXED, N](const BlockSrc& srcs, std::shared_ptr<DependentBlock> self) {
103 ++upd_MIXED;
104 for (int i = 0; i < N; ++i) {
105 LhaID id(i);
106 double a = srcs.get_val("SRC_A", id);
107 double z = srcs.get_val("FUSED_AB", id);
108 double y = srcs.get_val("FUSED_BC", id);
109 double m = 0.5*a + z + 0.1*y;
110 if (self->contains(id)) self->assign(id, m);
111 else self->store(id, std::make_shared<Parameter>(ParamId(self->get_name(), id), m, 0., 0.));
112 }
113 }
114 );
115 MIXED->blockname = "MIXED"; MIXED->set_scale(1.0); MIXED->init();
116
117 auto REDUCED = std::make_shared<DependentBlock>(
118 std::unordered_map<std::string, std::shared_ptr<Block>>{
119 {"MIXED", MIXED}
120 },
121 [&upd_REDUCED, N](const BlockSrc& srcs, std::shared_ptr<DependentBlock> self) {
122 ++upd_REDUCED;
123 double sum = 0.0;
124 for (int i = 0; i < N; ++i) sum += srcs.get_val("MIXED", LhaID(i));
125 if (self->contains(LhaID(0))) self->assign(LhaID(0), sum);
126 else self->store(LhaID(0), std::make_shared<Parameter>(ParamId(self->get_name(), LhaID(0)), sum, 0., 0.));
127 }
128 );
129 REDUCED->blockname = "REDUCED"; REDUCED->set_scale(1.0); REDUCED->init();
130
131
132 MIXED->freeze();
133 REDUCED->freeze();
134
135 FUSED_AB->update();
136 FUSED_BC->update();
137 SCALED_A->update();
138
139 MIXED->unfreeze();
140 MIXED->update();
141 REDUCED->unfreeze();
142 REDUCED->update();
143
144 upd_FUSED_AB = upd_FUSED_BC = upd_SCALED_A = upd_MIXED = upd_REDUCED = 0;
145
146 auto DERIVED1 = std::make_shared<Block>(); DERIVED1->blockname = "DERIVED1"; DERIVED1->set_scale(1.0);
147 auto DERIVED2 = std::make_shared<Block>(); DERIVED2->blockname = "DERIVED2"; DERIVED2->set_scale(1.0);
148
149 auto A_i1 = SRC_A->retrieve(LhaID(i1))->get_id();
150 auto A_i2 = SRC_A->retrieve(LhaID(i2))->get_id();
151 auto Z_i1 = FUSED_AB->retrieve(LhaID(i1))->get_id();
152 auto Z_i2 = FUSED_AB->retrieve(LhaID(i2))->get_id();
153 auto M_j1 = MIXED->retrieve(LhaID(j1))->get_id();
154 auto M_j2 = MIXED->retrieve(LhaID(j2))->get_id();
155 auto S_j1 = SCALED_A->retrieve(LhaID(j1))->get_id();
156 auto S_j2 = SCALED_A->retrieve(LhaID(j2))->get_id();
157
158 auto DP1 = std::make_shared<DependentParameter>(
159 ParamId{ParameterType::SM, "DERIVED1", LhaID(1001)},
160 std::unordered_map<ParamId, std::shared_ptr<Parameter>>{
161 {A_i1, SRC_A->retrieve(LhaID(i1))}, {Z_i1, FUSED_AB->retrieve(LhaID(i1))}
162 },
163 [A_i1, Z_i1, &upd_DP1](const ParamSrc& src, std::shared_ptr<DependentParameter> self){
164 ++upd_DP1; self->set_expected(src.get_val(A_i1) + src.get_val(Z_i1));
165 }
166 ); DP1->init(); DERIVED1->store(LhaID(1001), DP1);
167
168 auto DP2 = std::make_shared<DependentParameter>(
169 ParamId{ParameterType::SM, "DERIVED1", LhaID(1002)},
170 std::unordered_map<ParamId, std::shared_ptr<Parameter>>{
171 {A_i2, SRC_A->retrieve(LhaID(i2))}, {Z_i2, FUSED_AB->retrieve(LhaID(i2))}
172 },
173 [A_i2, Z_i2, &upd_DP2](const ParamSrc& src, std::shared_ptr<DependentParameter> self){
174 ++upd_DP2; self->set_expected(src.get_val(A_i2) - src.get_val(Z_i2));
175 }
176 ); DP2->init(); DERIVED1->store(LhaID(1002), DP2);
177
178 auto DP3 = std::make_shared<DependentParameter>(
179 ParamId{ParameterType::SM, "DERIVED2", LhaID(2001)},
180 std::unordered_map<ParamId, std::shared_ptr<Parameter>>{
181 {M_j1, MIXED->retrieve(LhaID(j1))}, {S_j1, SCALED_A->retrieve(LhaID(j1))}
182 },
183 [M_j1, S_j1, &upd_DP3](const ParamSrc& src, std::shared_ptr<DependentParameter> self){
184 ++upd_DP3; self->set_expected(src.get_val(M_j1) + src.get_val(S_j1));
185 }
186 ); DP3->init(); DERIVED2->store(LhaID(2001), DP3);
187
188 auto DP4 = std::make_shared<DependentParameter>(
189 ParamId{ParameterType::SM, "DERIVED2", LhaID(2002)},
190 std::unordered_map<ParamId, std::shared_ptr<Parameter>>{
191 {M_j2, MIXED->retrieve(LhaID(j2))}, {S_j2, SCALED_A->retrieve(LhaID(j2))}
192 },
193 [M_j2, S_j2, &upd_DP4](const ParamSrc& src, std::shared_ptr<DependentParameter> self){
194 ++upd_DP4; self->set_expected(src.get_val(M_j2) - 0.3*src.get_val(S_j2));
195 }
196 ); DP4->init(); DERIVED2->store(LhaID(2002), DP4);
197
198 auto BA1 = std::make_shared<BlockAccessor>();
199 BA1->emplace("SRC_A", SRC_A);
200 BA1->emplace("SRC_B", SRC_B);
201 BA1->emplace("FUSED_AB", FUSED_AB);
202 BA1->emplace("SCALE", SCALE);
203 BA1->emplace("SCALED_A", SCALED_A);
204
205 auto BA2 = std::make_shared<BlockAccessor>();
206 BA2->emplace("SRC_C", SRC_C);
207 BA2->emplace("FUSED_BC", FUSED_BC);
208 BA2->emplace("MIXED", MIXED);
209 BA2->emplace("REDUCED", REDUCED);
210 BA2->emplace("DERIVED1", DERIVED1);
211 BA2->emplace("DERIVED2", DERIVED2);
212
213 ParamOptimizer opt({BA1, BA2});
214
215 upd_FUSED_AB = upd_FUSED_BC = upd_SCALED_A = upd_MIXED = upd_REDUCED = 0;
216 upd_DP1 = upd_DP2 = upd_DP3 = upd_DP4 = 0;
217
218 auto t0 = clock::now();
219 for (int k = 0; k < OPS; ++k) {
220 int i = k % N;
221 SRC_A->assign(LhaID(i), 0.1*k + std::sin(0.001*k));
222 SRC_B->assign(LhaID(i), 0.2*k + std::cos(0.001*k));
223 SRC_C->assign(LhaID(i), 0.05*k);
224
225 if ((k % 1000) == 0) SCALE->assign(LhaID(1), 1.0 + 0.001*k);
226
227 if ((k % 5000) == 0) {
228 LhaID nid(N + (k/5000));
229 SRC_A->store(nid, make_param("SRC_A", nid, 42.0 + k*1e-3));
230 SRC_A->notifyObservers(); // propager l’ajout d’ID
231 }
232 }
233 auto t1 = clock::now();
234 auto naive_ms = std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count();
235
236 std::cout << "[NAIF] "
237 << "upd(FUSED_AB)=" << upd_FUSED_AB << ", "
238 << "upd(FUSED_BC)=" << upd_FUSED_BC << ", "
239 << "upd(SCALED_A)=" << upd_SCALED_A << ", "
240 << "upd(MIXED)=" << upd_MIXED << ", "
241 << "upd(REDUCED)=" << upd_REDUCED << " | "
242 << "upd(DP1,DP2,DP3,DP4)=" << upd_DP1 << "," << upd_DP2 << "," << upd_DP3 << "," << upd_DP4
243 << " | time=" << naive_ms << " ms\n";
244
245 auto check_point = [&](int idx){
246 double a = SRC_A->retrieve(LhaID(idx))->get_val();
247 double b = SRC_B->retrieve(LhaID(idx))->get_val();
248 double c = SRC_C->retrieve(LhaID(idx))->get_val();
249 double z = FUSED_AB->retrieve(LhaID(idx))->get_val();
250 double y = FUSED_BC->retrieve(LhaID(idx))->get_val();
251 double sf = SCALE->retrieve(LhaID(1))->get_val();
252 double s = SCALED_A->retrieve(LhaID(idx))->get_val();
253 double m = MIXED->retrieve(LhaID(idx))->get_val();
254 assert(std::abs(z - (a + 2.0*b)) < 1e-9);
255 assert(std::abs(y - (b - c)) < 1e-9);
256 assert(std::abs(s - (sf * a)) < 1e-9);
257 assert(std::abs(m - (0.5*a + z + 0.1*y)) < 1e-9);
258 };
259 check_point(3); check_point(17); check_point(123);
260
261 upd_FUSED_AB = upd_FUSED_BC = upd_SCALED_A = upd_MIXED = upd_REDUCED = 0;
262 upd_DP1 = upd_DP2 = upd_DP3 = upd_DP4 = 0;
263
264 auto t2 = clock::now();
265 for (int k = 0; k < OPS; ++k) {
266 int i = k % N;
267 opt.set_value("SRC_A", LhaID(i), 0.1*k + std::sin(0.001*k));
268 opt.set_value("SRC_B", LhaID(i), 0.2*k + std::cos(0.001*k));
269 opt.set_value("SRC_C", LhaID(i), 0.05*k);
270 if ((k % 1000) == 0) opt.set_value("SCALE", LhaID(1), 1.0 + 0.001*k);
271 if ((k % 5000) == 0) {
272 LhaID nid(N + (k/5000));
273 opt.set_value("SRC_A", nid, 42.0 + k*1e-3);
274 }
275 }
276 opt.commit(true);
277 auto t3 = clock::now();
278 auto opt_ms = std::chrono::duration_cast<std::chrono::milliseconds>(t3 - t2).count();
279
280 std::cout << "[OPT ] "
281 << "upd(FUSED_AB)=" << upd_FUSED_AB << ", "
282 << "upd(FUSED_BC)=" << upd_FUSED_BC << ", "
283 << "upd(SCALED_A)=" << upd_SCALED_A << ", "
284 << "upd(MIXED)=" << upd_MIXED << ", "
285 << "upd(REDUCED)=" << upd_REDUCED << " | "
286 << "upd(DP1,DP2,DP3,DP4)=" << upd_DP1 << "," << upd_DP2 << "," << upd_DP3 << "," << upd_DP4
287 << " | time=" << opt_ms << " ms\n";
288
289 check_point(7); check_point(111); check_point(482);
290 {
291 double a1 = SRC_A->retrieve(LhaID(i1))->get_val();
292 double z1 = FUSED_AB->retrieve(LhaID(i1))->get_val();
293 double a2 = SRC_A->retrieve(LhaID(i2))->get_val();
294 double z2 = FUSED_AB->retrieve(LhaID(i2))->get_val();
295 double m1 = MIXED->retrieve(LhaID(j1))->get_val();
296 double s1 = SCALED_A->retrieve(LhaID(j1))->get_val();
297 double m2 = MIXED->retrieve(LhaID(j2))->get_val();
298 double s2 = SCALED_A->retrieve(LhaID(j2))->get_val();
299
300 assert(std::abs(DERIVED1->retrieve(LhaID(1001))->get_val() - (a1 + z1)) < 1e-9);
301 assert(std::abs(DERIVED1->retrieve(LhaID(1002))->get_val() - (a2 - z2)) < 1e-9);
302 assert(std::abs(DERIVED2->retrieve(LhaID(2001))->get_val() - (m1 + s1)) < 1e-9);
303 assert(std::abs(DERIVED2->retrieve(LhaID(2002))->get_val() - (m2 - 0.3*s2)) < 1e-9);
304 }
305
306 std::cout << "OK: test complexe — cohérence + bench.\n";
307 return 0;
308}
BlockAccessor.h
Alias-aware façade for accessing and manipulating multiple parameter blocks.
Block.h
Defines classes used to store parameters and to build derived/dependent parameter blocks.
DependentParameter.h
Defines parameters whose values are lazily computed from other parameters.
ParamOptimizer.h
Defines the ParamOptimizer class to apply batched updates to parameters.
BlockSrc
Lightweight view over a set of source blocks.
Definition SourcesView.h:71
ParamOptimizer
Helper class to batch parameter updates on one or several BlockAccessor scopes.
Definition ParamOptimizer.h:75
ParamSrc
Lightweight view over a set of source parameters keyed by ParamId.
Definition SourcesView.h:249
main
int main()
Definition main_optimizer_deeper.cpp:19
LhaID
Represents an identifier of a LHA element, possibly containing several sub-ids.
Definition LhaID.h:56
ParamId
Composite identifier for a single parameter.
Definition ParamID.h:57