logoChat01.ai

"#include "pch.h" #define XOR_RAND(state, result_...

Created on: November 5, 2025

Answered using GPT-5 Thinking by Chat01

Question

"#include "pch.h"

#define XOR_RAND(state, result_var)
do {
int s = (state);
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
state = s;
result_var = state * 0x1.0p-32f;
} while (0)

#define XOR_RAND_GRSH(state, result_var)
do {
int s = (state);
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
state = s;
result_var = fmaf(state, 0x1.0p-31f, -1.0f);
} while (0)

#define FABE13_COS(x, result_var)
do {
const float ax = fabsf(x);
float r = fmodf(ax, 6.28318530718f);
if (r > 3.14159265359f)
r = 6.28318530718f - r;
if (r < 1.57079632679f) {
const float t2 = r * r;
const float t4 = t2 * t2;
result_var = fmaf(t4, fmaf(t2, -0.0013888889f, 0.0416666667f), fmaf(t2, -0.5f, 1.0f));
} else {
r = 3.14159265359f - r;
const float t2 = r * r;
const float t4 = t2 * t2;
result_var = -fmaf(t4, fmaf(t2, -0.0013888889f, 0.0416666667f), fmaf(t2, -0.5f, 1.0f));
}
} while (0)

#define FABE13_SIN(x, result_var)
do {
const float x = (x);
const float ax = fabsf(x);
float r = fmodf(ax, 6.28318530718f);
bool sfl = r > 3.14159265359f;
if (sfl)
r = 6.28318530718f - r;
bool cfl = r > 1.57079632679f;
if (cfl)
r = 3.14159265359f - r;
const float t2 = r * r;
float _s = fmaf(t2, fmaf(t2, fmaf(t2, -0.0001984127f, 0.0083333333f), -0.16666666f), 1.0f) * r;
result_var = ((x < 0.0f) ^ sfl) ? -_s : _s;
} while (0)

#define FABE13_SINCOS(in, sin_out, cos_out, n)
do {
int i = 0;
const int limit = (n) & ~7;
if ((n) >= 8) {
static __declspec(align(32)) const __m256 VEC_TWOPI = _mm256_set1_ps(6.28318530718f);
static __declspec(align(32)) const __m256 VEC_PI = _mm256_set1_ps(3.14159265359f);
static __declspec(align(32)) const __m256 VEC_PI_2 = _mm256_set1_ps(1.57079632679f);
static __declspec(align(32)) const __m256 INV_TWOPI = _mm256_set1_ps(0.15915494309189535f);
static __declspec(align(32)) const __m256 BIAS = _mm256_set1_ps(12582912.0f);
static __declspec(align(32)) const __m256 VEC_COS_P5 = _mm256_set1_ps(-0.0013888889f);
static __declspec(align(32)) const __m256 VEC_COS_P3 = _mm256_set1_ps(0.0416666667f);
static __declspec(align(32)) const __m256 VEC_COS_P1 = _mm256_set1_ps(-0.5f);
static __declspec(align(32)) const __m256 VEC_COS_P0 = _mm256_set1_ps(1.0f);
static __declspec(align(32)) const __m256 VEC_SIN_P5 = _mm256_set1_ps(-0.0001984127f);
static __declspec(align(32)) const __m256 VEC_SIN_P3 = _mm256_set1_ps(0.0083333333f);
static __declspec(align(32)) const __m256 VEC_SIN_P1 = _mm256_set1_ps(-0.16666666f);
static __declspec(align(32)) const __m256 VEC_SIN_P0 = _mm256_set1_ps(1.0f);
static __declspec(align(32)) const __m256 VEC_ZERO = _mm256_setzero_ps();
while (i < limit) {
const __m256 vx = _mm256_load_ps(&(in)[i]);
const __m256 vax = _mm256_andnot_ps(_mm256_set1_ps(-0.0f), vx);
__m256 q = _mm256_fmadd_ps(vax, INV_TWOPI, BIAS);
q = _mm256_sub_ps(q, BIAS);
const __m256 r = _mm256_fnmadd_ps(VEC_TWOPI, q, vax);
const __m256 r1 = _mm256_min_ps(r, _mm256_sub_ps(VEC_TWOPI, r));
const __m256 r2 = _mm256_min_ps(r1, _mm256_sub_ps(VEC_PI, r1));
const __m256 t2 = _mm256_mul_ps(r2, r2);
const __m256 cosv = _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_COS_P5, VEC_COS_P3), VEC_COS_P1), VEC_COS_P0);
const __m256 sinv = _mm256_mul_ps(_mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_SIN_P5, VEC_SIN_P3), VEC_SIN_P1), VEC_SIN_P0), r2);
const __m256 cflip = _mm256_cmp_ps(r1, VEC_PI_2, _CMP_GT_OQ);
const __m256 sflip = _mm256_xor_ps(_mm256_cmp_ps(vx, VEC_ZERO, _CMP_LT_OQ), _mm256_cmp_ps(r, VEC_PI, _CMP_GT_OQ));
_mm256_store_ps(&(cos_out)[i], _mm256_blendv_ps(cosv, _mm256_sub_ps(VEC_ZERO, cosv), cflip));
_mm256_store_ps(&(sin_out)[i], _mm256_blendv_ps(sinv, _mm256_sub_ps(VEC_ZERO, sinv), sflip));
i += 8;
}
}
while (i < (n)) {
const float x = (in)[i];
const float ax = fabsf(x);
float q = fmaf(ax, 0.15915494309189535f, 12582912.0f);
q -= 12582912.0f;
float r = fmaf(-6.28318530718f, q, ax);
const bool sflip = r > 3.14159265359f;
if (sflip)
r = 6.28318530718f - r;
const bool cflip = r > 1.57079632679f;
if (cflip)
r = 3.14159265359f - r;
const float t2 = r * r;
const float c = fmaf(t2, fmaf(t2, fmaf(t2, -0.0013888889f, 0.0416666667f), -0.5f), 1.0f);
const float s = fmaf(t2, fmaf(t2, fmaf(t2, -0.0001984127f, 0.0083333333f), -0.16666666f), 1.0f) * r;
(cos_out)[i] = cflip ? -c : c;
(sin_out)[i] = ((x < 0.0f) ^ sflip) ? -s : s;
++i;
}
} while (0)

enum List : uint8_t {
Top = 0b00u,
Down = 0b01u,
Left = 0b10u,
Right = 0b11u
};

__declspec(align(4)) struct Step final {
const uint8_t next;
const uint8_t dx;
const uint8_t dy;
};

__declspec(align(4)) struct InvStep final {
const uint8_t q;
const uint8_t next;
};

__declspec(align(64)) static const Step g_step_tbl[4][4] = {
{ { Right, 0u, 0u }, { Top, 0u, 1u }, { Top, 1u, 1u }, { Left, 1u, 0u } },
{ { Left, 1u, 1u }, { Down, 1u, 0u }, { Down, 0u, 0u }, { Right, 0u, 1u } },
{ { Down, 1u, 1u }, { Left, 0u, 1u }, { Left, 0u, 0u }, { Top, 1u, 0u } },
{ { Top, 0u, 0u }, { Right, 1u, 0u }, { Right, 1u, 1u }, { Down, 0u, 1u } }
};

__declspec(align(64)) static const InvStep g_inv_tbl[4][4] = {
{ { 0u, Right }, { 1u, Top }, { 3u, Left }, { 2u, Top } },
{ { 2u, Down }, { 3u, Right }, { 1u, Down }, { 0u, Left } },
{ { 2u, Left }, { 1u, Left }, { 3u, Top }, { 0u, Down } },
{ { 0u, Top }, { 3u, Down }, { 1u, Right }, { 2u, Right } }
};

static const boost::mpi::environment* g_env;
static const boost::mpi::communicator* g_world;

__declspec(align(16)) struct CrossMsg final {
float s_x1, s_x2;
float e_x1, e_x2;
float Rtop;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& s_x1& s_x2& e_x1& e_x2& Rtop;
}
};

__declspec(align(16)) struct CtrlMsg final {
bool kind;
CrossMsg xchg;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& kind& xchg;
}
};

__declspec(align(16)) struct Slab final {
char* const base;
char* current;
char* const end;
__forceinline Slab(void* const memory, const size_t usable) : base((char*)memory), current(base), end(base + (usable & ~(size_t)63u)) {}
};

static tbb::enumerable_thread_specific<Slab*> tls( noexcept {
void* memory = _aligned_malloc(16777216u, 64u);
Slab* slab = (Slab*)_aligned_malloc(32u, 64u);
new (slab) Slab(memory, 16777216u);
char* p = slab->base;
while (p < slab->end) {
*p = 0;
p += 4096u;
}
return slab;
});

__declspec(align(16)) struct Peano2DMap final {
const int levels;
const float a, b, c, d;
const float lenx, leny;
const float inv_lenx;
const uint32_t scale;
const uint8_t start;
__forceinline Peano2DMap(int L, float _a, float _b, float _c, float _d, uint8_t st)
: levels(L), a(_a), b(_b), c(_c), d(_d), lenx(_b - _a), leny(_d - _c), inv_lenx(1.0f / (_b - _a)), scale((uint32_t)1u << (L << 1)), start(st) {
}
};

static Peano2DMap gActiveMap(0, 0, 0, 0, 0, 0);

__forceinline bool ComparePtr(const struct Interval* a, const struct Interval* b) noexcept;

__declspec(align(64)) struct Interval final {
const float x1, x2, y1, y2, delta_y, ordinate_factor, N_factor, quadratic_term, M;
float R;
__forceinline void* operator new(size_t) noexcept {
Slab* s = tls.local();
char* r = s->current;
s->current += 64u;
return r;
}
__forceinline Interval(float _x1, float _x2, float _y1, float _y2, float _N) noexcept
: x1(_x1), x2(_x2), y1(_y1), y2(_y2), delta_y(_y2 - _y1), ordinate_factor(-(y1 + y2) * 2.0f),
N_factor(_N == 1.0f ? _x2 - _x1 : _N == 2.0f ? sqrtf(_x2 - _x1) : powf(_x2 - _x1, 1.0f / _N)),
quadratic_term((1.0f / N_factor)* delta_y* delta_y), M((1.0f / N_factor)* fabsf(delta_y)) {
}
__forceinline void ChangeCharacteristic(float _m) noexcept {
R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor));
}
};

__forceinline bool ComparePtr(const Interval* a, const Interval* b) noexcept {
return a->R < b->R;
}

__forceinline void RecomputeR_ConstM_AVX2(Interval* const* arr, size_t n, float m) {
const __m256 vm = _mm256_set1_ps(m);
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));
size_t i = 0;
const size_t limit = n & ~7ull;
alignas(32) float q[8], nf[8], od[8], out[8];
for (; i < limit; i += 8) {
for (int k = 0; k < 8; ++k) {
const Interval* p = arr[i + k];
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
od[k] = p->ordinate_factor;
}
const __m256 vq = _mm256_load_ps(q), vnf = _mm256_load_ps(nf), vod = _mm256_load_ps(od);
const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);
_mm256_store_ps(out, res);
for (int k = 0; k < 8; ++k)
arr[i + k]->R = out[k];
}
for (; i < n; ++i)
arr[i]->ChangeCharacteristic(m);
}

__forceinline void RecomputeR_AffineM_AVX2(Interval* const* arr, size_t n, float GF, float alpha) {
const __m256 vGF = _mm256_set1_ps(GF), va = _mm256_set1_ps(alpha);
size_t i = 0;
const size_t limit = n & ~7ull;
alignas(32) float ln[8], Mv[8], q[8], nf[8], od[8], out[8];
for (; i < limit; i += 8) {
for (int k = 0; k < 8; ++k) {
const Interval* p = arr[i + k];
ln[k] = p->x2 - p->x1;
Mv[k] = p->M;
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
od[k] = p->ordinate_factor;
}
const __m256 vln = _mm256_load_ps(ln), vM = _mm256_load_ps(Mv), vq = _mm256_load_ps(q), vnf = _mm256_load_ps(nf), vod = _mm256_load_ps(od);
const __m256 vm = _mm256_fmadd_ps(vGF, vln, _mm256_mul_ps(va, vM));
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));
const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);
_mm256_store_ps(out, res);
for (int k = 0; k < 8; ++k)
arr[i + k]->R = out[k];
}
for (; i < n; ++i) {
const Interval* p = arr[i];
const float mi = fmaf(GF, (p->x2 - p->x1), p->M * alpha);
arr[i]->R = fmaf(1.0f / mi, p->quadratic_term, fmaf(mi, p->N_factor, p->ordinate_factor));
}
}

__forceinline float fast_pow_int(float v, int n) {
float r;
switch (n) {
case 3: {
float v2 = v * v;
r = v2 * v;
} break;
case 4: {
float v2 = v * v;
r = v2 * v2;
} break;
case 5: {
float v2 = v * v;
r = v2 * v2 * v;
} break;
case 6: {
float v2 = v * v;
float v4 = v2 * v2;
r = v4 * v2;
} break;
case 7: {
float v2 = v * v;
float v4 = v2 * v2;
r = v4 * v2 * v;
} break;
case 8: {
float v2 = v * v;
float v4 = v2 * v2;
r = v4 * v4;
} break;
case 9: {
float v3 = v * v * v;
float v6 = v3 * v3;
r = v6 * v3;
} break;
case 10: {
float v2 = v * v;
float v4 = v2 * v2;
float v8 = v4 * v4;
r = v8 * v2;
} break;
case 11: {
float v2 = v * v;
float v4 = v2 * v2;
float v8 = v4 * v4;
r = v8 * v2 * v;
} break;
case 12: {
float v3 = v * v * v;
float v6 = v3 * v3;
r = v6 * v6;
} break;
case 13: {
float v3 = v * v * v;
float v6 = v3 * v3;
r = v6 * v6 * v;
} break;
case 14: {
float v7 = v * v * v * v * v * v * v;
r = v7 * v7;
} break;
case 15: {
float v7 = v * v * v * v * v * v * v;
r = v7 * v7 * v;
} break;
default: {
float v2 = v * v;
float v4 = v2 * v2;
float v8 = v4 * v4;
r = v8 * v8;
}
}
return r;
}

__forceinline float Shag(float _m, float x1, float x2, float y1, float y2, float _N, float _r) {
const float diff = y2 - y1;
const float sign_mult = _mm_cvtss_f32(_mm_castsi128_ps(_mm_set1_epi32(0x3F800000u | (((((uint32_t)&diff)) & 0x80000000u) ^ 0x80000000u))));
if (_N == 1)
return fmaf(-(1.0f / _m), diff, x1 + x2) * 0.5f;
if (_N == 2)
return fmaf(sign_mult / (_m * _m), diff * diff * _r, x1 + x2) * 0.5f;
const float invmN = 1.0f / fast_pow_int(_m, _N);
const float dN = fast_pow_int(fabsf(diff), _N);
return fmaf(sign_mult * invmN, dN * _r, x1 + x2) * 0.5f;
}

struct MortonCachePerRank {
std::vector<int> permCache;
std::vector<uint64_t> invMaskCache;
uint32_t baseSeed;
};

static MortonCachePerRank g_mc;

struct MortonND final {
int dim, levels;
uint64_t scale;
std::vector<float> low, high, step, invStep, baseOff;
std::vector<int> perm;
std::vector<uint64_t> invMask, pextMask;
float invScaleLevel;
__forceinline MortonND(int D, int L, const float* lows, const float* highs, const MortonCachePerRank& mc)
: dim(D), levels(L), scale(1ull << ((uint64_t)D * L)), low(lows, lows + D), high(highs, highs + D), step(D, 0.0f), invStep(D, 0.0f), baseOff(D, 0.0f),
perm(mc.permCache.begin(), mc.permCache.begin() + D), invMask(mc.invMaskCache.begin(), mc.invMaskCache.begin() + D), pextMask(D, 0ull),
invScaleLevel(1.0f / (float)((uint64_t)1 << L)) {
uint64_t level_mask = (1ull << L) - 1ull;
for (int k = 0; k < D; ++k)
invMask[k] &= level_mask;
for (int d = 0; d < dim; ++d) {
float rng = high[d] - low[d];
float st = rng * invScaleLevel;
step[d] = st;
invStep[d] = 1.0f / st;
baseOff[d] = fmaf(0.5f, st, low[d]);
}
for (int d = 0; d < dim; ++d) {
uint64_t m = 0ull, bitpos = (uint64_t)d;
for (int b = 0; b < levels; ++b) {
m |= 1ull << bitpos;
bitpos += (uint64_t)dim;
}
pextMask[d] = m;
}
}
__forceinline void map01ToPoint(float t, float* __restrict out) const noexcept {
if (t <= 0.0f)
t = 0.0f;
else if (t >= 1.0f)
t = 0x1.fffffep-1f;
uint64_t idx = (uint64_t)((double)t * (double)scale);
if (idx >= scale)
idx = scale - 1ull;
for (int d = 0; d < dim; ++d) {
int pd = perm[d];
uint64_t bits = _pext_u64(idx, pextMask[d]);
if (invMask[pd])
bits ^= invMask[pd];
out[pd] = fmaf(step[pd], (float)bits, baseOff[pd]);
}
}
__forceinline uint64_t clamp_bits(int64_t v, int bits) const noexcept {
if (v < 0)
v = 0;
int64_t maxv = ((int64_t)1 << bits) - 1;
if (v > maxv)
v = maxv;
return (uint64_t)v;
}
__forceinline float pointToT(const float* __restrict q) const noexcept {
const int bits = levels;
uint64_t idx = 0ull;
for (int d = 0; d < dim; ++d) {
int pd = perm[d];
float v = (q[pd] - baseOff[pd]) * invStep[pd];
int64_t cell = (int64_t)_mm_cvt_ss2si(_mm_round_ss(_mm_setzero_ps(), _mm_set_ss(v), _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC));
uint64_t b = clamp_bits(cell, bits);
if (invMask[pd])
b ^= invMask[pd];
idx |= _pdep_u64(b, pextMask[d]);
}
return ((float)idx + 0.5f) / (float)scale;
}
};

struct ManipCost final {
int n;
bool variableLen;
float targetX, targetY;
float minTheta;
__forceinline ManipCost(int _n, bool _variableLen, float _targetX, float _targetY, float _minTheta)
: n(_n), variableLen(_variableLen), targetX(_targetX), targetY(_targetY), minTheta(_minTheta) {
}
__forceinline float operator()(const float* __restrict q, float& out_x, float& out_y) const noexcept {
const float* th = q;
const float* L = variableLen ? (q + n) : nullptr;
__declspec(align(64)) float phi[16], s_arr[16], c_arr[16];
float x = 0.0f, y = 0.0f, phi_acc = 0.0f, penC = 0.0f;
for (int i = 0; i < n; ++i) {
phi_acc += th[i];
phi[i] = phi_acc;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);
const float Lc = 1.0f;
if (variableLen) {
for (int i = 0; i < n; ++i) {
float Li = L[i];
x = fmaf(Li, c_arr[i], x);
y = fmaf(Li, s_arr[i], y);
}
}
else {
for (int i = 0; i < n; ++i) {
x = fmaf(Lc, c_arr[i], x);
y = fmaf(Lc, s_arr[i], y);
}
}
for (int i = 0; i < n; ++i) {
float ai = fabsf(th[i]);
float v = minTheta - ai;
float u = v > 0.0f ? v : 0.0f;
penC = fmaf(u, u, penC);
}
float dx = x - targetX;
float dy = y - targetY;
float dist = sqrtf(fmaf(dx, dx, dy * dy));
out_x = x;
out_y = y;
return dist + penC;
}
};

__forceinline void HitTest2D_analytic(float x_param, float& out_x1, float& out_x2) {
const float a = gActiveMap.a;
const float inv_lenx = gActiveMap.inv_lenx;
const uint32_t scale = gActiveMap.scale;
const uint32_t scale_minus_1 = scale - 1u;
const float lenx = gActiveMap.lenx;
const float leny = gActiveMap.leny;
const float c = gActiveMap.c;
const uint8_t start = gActiveMap.start;
const int levels = gActiveMap.levels;
float norm = (x_param - a) * inv_lenx;
norm = fminf(fmaxf(norm, 0.0f), 0x1.fffffep-1f);
uint32_t idx = (uint32_t)(norm * (float)scale);
idx = idx > scale_minus_1 ? scale_minus_1 : idx;
float sx = lenx, sy = leny;
float x1 = a, x2 = c;
uint8_t type = start;
int l = levels - 1;
while (l >= 0) {
const uint32_t q = (idx >> (l * 2)) & 3u;
const Step s = g_step_tbl[type][q];
type = s.next;
sx *= 0.5f;
sy *= 0.5f;
x1 += s.dx ? sx : 0.0f;
x2 += s.dy ? sy : 0.0f;
--l;
}
out_x1 = x1 + sx * 0.5f;
out_x2 = x2 + sy * 0.5f;
}

__forceinline float FindX2D_analytic(float px, float py) {
const float a = gActiveMap.a;
const float b = gActiveMap.b;
const float c = gActiveMap.c;
const float d = gActiveMap.d;
const float lenx = gActiveMap.lenx;
const float leny = gActiveMap.leny;
const uint32_t scale = gActiveMap.scale;
const uint8_t start = gActiveMap.start;
const int levels = gActiveMap.levels;
const float clamped_px = fminf(fmaxf(px, a), b);
const float clamped_py = fminf(fmaxf(py, c), d);
float sx = lenx, sy = leny;
float x0 = a, y0 = c;
uint32_t idx = 0u;
uint8_t type = start;
int l = 0;
while (l < levels) {
sx *= 0.5f;
sy *= 0.5f;
const float mx = x0 + sx;
const float my = y0 + sy;
const uint32_t tr = (uint32_t)((clamped_px > mx) & (clamped_py > my));
const uint32_t tl = (uint32_t)((clamped_px < mx) & (clamped_py > my));
const uint32_t dl = (uint32_t)((clamped_px < mx) & (clamped_py < my));
const uint32_t none = (uint32_t)(1u ^ (tr | tl | dl));
const uint32_t dd = (tr << 1) | tr | tl | (none << 1);
const InvStep inv = g_inv_tbl[type][dd];
type = inv.next;
idx = (idx << 2) | inv.q;
const uint32_t dx = dd >> 1;
const uint32_t dy = dd & 1u;
x0 += dx ? sx : 0.0f;
y0 += dy ? sy : 0.0f;
++l;
}
const float scale_recip = 1.0f / (float)scale;
return fmaf((float)idx * scale_recip, lenx, a);
}

__declspec(align(16)) struct MultiCrossMsg final {
float intervals[15];
uint8_t count;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& intervals& count;
}
};

__declspec(align(16)) struct BestSolutionMsg final {
float bestF;
float bestX, bestY;
float bestQ[32];
uint8_t dim;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& bestF& bestX& bestY& bestQ& dim;
}
};

__declspec(align(16)) struct CtrlMsgND final {
uint8_t kind;
CrossMsg xchg;
MultiCrossMsg multiXchg;
BestSolutionMsg bestSol;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& kind;
if (kind == 1)
ar& xchg;
else if (kind == 2)
ar& multiXchg;
else if (kind == 3)
ar& bestSol;
}
};

static __forceinline int generate_heuristic_seeds(const ManipCost& cost, int dim, float* __restrict S, int stride) {
const int n = cost.n;
const bool VL = cost.variableLen;
const float tx = cost.targetX, ty = cost.targetY;
const float phi = atan2f(ty, tx);
const float rho = sqrtf(fmaf(tx, tx, ty * ty));
float len = fminf(fmaxf(rho / (float)n, 0.1f), 2.0f);
float* s0 = S + 0 * stride;
for (int i = 0; i < n; ++i)
s0[i] = phi / (float)n;
if (VL)
for (int i = 0; i < n; ++i)
s0[n + i] = len;
float* s1 = S + 1 * stride;
for (int i = 0; i < n; ++i)
s1[i] = 0.5f * phi * ((i & 1) ? -1.0f : 1.0f);
if (VL) {
for (int i = 0; i < n; ++i)
s1[n + i] = len * (0.8f + 0.4f * (float)i / (float)n);
}
float* s2 = S + 2 * stride;
const float inv = (n > 1) ? 1.0f / (float)(n - 1) : 0.0f;
for (int i = 0; i < n; ++i) {
float pr = (float)i * inv;
s2[i] = phi * (1.0f - 0.3f * pr);
}
if (VL) {
for (int i = 0; i < n; ++i) {
float si = sinf(1.5f * (float)i);
s2[n + i] = len * (1.0f + 0.2f * si);
}
}
return 3;
}

static __forceinline void agp_run_branch_mpi(const MortonND& map, const ManipCost& cost, int maxIter, float r, bool adaptive, float eps, unsigned seed,
std::vector<Interval*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY,
float M_prior = 1e-3f) {
const int n = cost.n, dim = n + (cost.variableLen ? n : 0);
alignas(64) float q_local[32], phi[32], s_arr[32], c_arr[32], sum_s[32], sum_c[32], q_try[32];
bestQ.reserve(dim);
float x = 0.0f, y = 0.0f;
int no_improve = 0;
auto evalAt = [&](float t) -> float {
map.map01ToPoint(t, q_local);
float f = cost(q_local, x, y);
if (f < bestF * 1.05f) {
float acc = 0.0f;
for (int i = 0; i < n; ++i) {
acc += q_local[i];
phi[i] = acc;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);
float as = 0.0f, ac = 0.0f;
for (int k = n - 1; k >= 0; --k) {
const float Lk = cost.variableLen ? q_local[n + k] : 1.0f;
as += Lk * s_arr[k];
ac += Lk * c_arr[k];
sum_s[k] = as;
sum_c[k] = ac;
}
const float dx = x - cost.targetX, dy = y - cost.targetY;
float dist = sqrtf(fmaf(dx, dx, dy * dy)) + 1e-8f;
float eta = 0.1f;
for (int step = 0; step < 2; ++step) {
for (int i = 0; i < n; ++i) {
float gpen = 0.0f;
float ai = fabsf(q_local[i]);
float v = cost.minTheta - ai;
if (v > 0.0f)
gpen = -2.0f * copysignf(v, q_local[i]);
float grad = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen;
q_try[i] = q_local[i] - eta * grad;
if (q_try[i] < -3.14159265359f)
q_try[i] = -3.14159265359f;
else if (q_try[i] > 3.14159265359f)
q_try[i] = 3.14159265359f;
}
if (cost.variableLen) {
for (int i = 0; i < n; ++i) {
float grad = (dx * c_arr[i] + dy * s_arr[i]) / dist;
float v = q_local[n + i] - eta * grad;
if (v < 0.1f)
v = 0.1f;
else if (v > 2.0f)
v = 2.0f;
q_try[n + i] = v;
}
}
float x2, y2;
float f2 = cost(q_try, x2, y2);
if (f2 < f) {
memcpy(q_local, q_try, dim * sizeof(float));
f = f2;
x = x2;
y = y2;
break;
}
eta = 0.5f;
}
}
if (f < bestF) {
bestF = f;
bestQ.assign(q_local, q_local + dim);
bestX = x;
bestY = y;
no_improve = 0;
}
else
++no_improve;
return f;
};
float a = 0.0f, b = 1.0f;
float f_a = evalAt(a), f_b = evalAt(b);
const int K = (std::min)((std::max)(2 * dim, 8), 128);
H.clear();
H.reserve((size_t)maxIter + K + 16);
float Mmax = M_prior;
const int rank = g_world->rank();
const int world = g_world->size();
alignas(64) float seeds[3 * 32];
const int seedCnt = generate_heuristic_seeds(cost, dim, seeds, 32);
if (seedCnt > 0) {
const float s0 = seeds;
float t_seed = map.pointToT(s0);
float t1 = fmaxf(a, t_seed - 0.002f), t2 = fminf(b, t_seed + 0.002f);
if (t2 > t1) {
alignas(64) float q1[32], q2[32];
float x1, y1, x2, y2;
map.map01ToPoint(t1, q1);
float f1 = cost(q1, x1, y1);
map.map01ToPoint(t2, q2);
float f2 = cost(q2, x2, y2);
Interval* I = new Interval(t1, t2, f1, f2, (float)dim);
Mmax = (std::max)(Mmax, I->M);
I->ChangeCharacteristic(r * Mmax);
I->R = 0.85f;
H.emplace_back(I);
if (f1 < bestF) {
bestF = f1;
bestQ.assign(q1, q1 + dim);
bestX = x1;
bestY = y1;
}
if (f2 < bestF) {
bestF = f2;
bestQ.assign(q2, q2 + dim);
bestX = x2;
bestY = y2;
}
}
}
float prev_t = a, prev_f = f_a;
for (int k = 1; k <= K; ++k) {
float t = a + (b - a) * ((float)k / (K + 1)) + (float)rank / (float)(world * (K + 1));
float f = evalAt(t);
Interval I = new Interval(prev_t, t, prev_f, f, (float)dim);
Mmax = (std::max)(Mmax, I->M);
I->ChangeCharacteristic(r * Mmax);
H.emplace_back(I);
std::push_heap(H.begin(), H.end(), ComparePtr);
prev_t = t;
prev_f = f;
}
Interval* tail = new Interval(prev_t, b, prev_f, f_b, (float)dim);
Mmax = (std::max)(Mmax, tail->M);
tail->ChangeCharacteristic(r * Mmax);
H.emplace_back(tail);
std::push_heap(H.begin(), H.end(), ComparePtr);
if (seedCnt > 1) {
for (int i = 1; i < seedCnt; ++i) {
const float* s = seeds + i * 32;
float t_seed = map.pointToT(s);
float t1 = fmaxf(a, t_seed - 0.001f), t2 = fminf(b, t_seed + 0.001f);
if (t2 > t1) {
float f1 = evalAt(t1), f2 = evalAt(t2);
Interval* J = new Interval(t1, t2, f1, f2, (float)dim);
Mmax = (std::max)(Mmax, J->M);
J->ChangeCharacteristic(r * Mmax);
J->R = 0.90f;
H.emplace_back(J);
std::push_heap(H.begin(), H.end(), ComparePtr);
}
}
}
float dmax = b - a, initial_len = dmax, thr03 = 0.3f * initial_len, inv_thr03 = 1.0f / thr03;
int it = 0;
while (it < maxIter) {
if ((it % 100) == 0 && no_improve > 50 && !bestQ.empty()) {
float t_best = map.pointToT(bestQ.data());
for (int i = 0; i < 2; ++i) {
float off = (i == 0) ? 0.01f : -0.01f;
float t_seed = fminf(b, fmaxf(a, t_best + off));
float f_seed = evalAt(t_seed);
Interval J = new Interval(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed, (float)dim);
J->ChangeCharacteristic(r * Mmax);
J->R = 0.9f;
H.emplace_back(J);
std::push_heap(H.begin(), H.end(), ComparePtr);
}
no_improve = 0;
}
const float p = fmaf(-1.0f / initial_len, dmax, 1.0f);
bool stagnation = (no_improve > 100) && (it > 270);
float r_eff = fmaxf(1.0f, r * (0.7f + 0.3f * (1.0f - p)));
std::pop_heap(H.begin(), H.end(), ComparePtr);
Interval cur = H.back();
H.pop_back();
const float x1 = cur->x1, x2 = cur->x2, y1 = cur->y1, y2 = cur->y2;
float m = r_eff * Mmax;
float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r);
tNew = fminf(fmaxf(tNew, a), b);
float fNew = evalAt(tNew);
Interval* L = new Interval(x1, tNew, y1, fNew, (float)dim);
Interval* Rv = new Interval(tNew, x2, fNew, y2, (float)dim);
float Mloc = (std::max)(L->M, Rv->M);
if (adaptive) {
float len1 = tNew - x1, len2 = x2 - tNew;
if (len1 + len2 == dmax) {
dmax = (std::max)(len1, len2);
for (auto pI : H) {
float Ls = pI->x2 - pI->x1;
if (Ls > dmax)
dmax = Ls;
}
}
if ((thr03 > dmax && !(it % 3)) || (10.0f * dmax < initial_len)) {
if (Mloc > Mmax) {
Mmax = Mloc;
m = r_eff * Mmax;
}
const float progress = fmaf(-dmax, inv_thr03, 1.0f), alpha = progress * progress, beta = fmaf(-alpha, 1.0f, 2.0f);
const float MULT = (1.0f / dmax) * Mmax, global_coeff = fmaf(MULT, r_eff, -MULT), GF = fmaf(beta, global_coeff, 0.0f);
L->ChangeCharacteristic(fmaf(GF, len1, L->M * alpha));
Rv->ChangeCharacteristic(fmaf(GF, len2, Rv->M * alpha));
size_t sz = H.size();
RecomputeR_AffineM_AVX2(H.data(), sz, GF, alpha);
std::make_heap(H.begin(), H.end(), ComparePtr);
}
else {
if (Mloc > Mmax) {
Mmax = Mloc;
m = r_eff * Mmax;
L->ChangeCharacteristic(m);
Rv->ChangeCharacteristic(m);
if (Mloc > 1.15f * Mmax) {
size_t sz = H.size();
RecomputeR_ConstM_AVX2(H.data(), sz, m);
std::make_heap(H.begin(), H.end(), ComparePtr);
}
}
else {
L->ChangeCharacteristic(m);
Rv->ChangeCharacteristic(m);
}
}
}
else {
if (Mloc > Mmax) {
Mmax = Mloc;
float m2 = r_eff * Mmax;
L->ChangeCharacteristic(m2);
Rv->ChangeCharacteristic(m2);
if (Mloc > 1.15f * Mmax) {
size_t sz = H.size();
RecomputeR_ConstM_AVX2(H.data(), sz, m2);
std::make_heap(H.begin(), H.end(), ComparePtr);
}
}
else {
L->ChangeCharacteristic(m);
Rv->ChangeCharacteristic(m);
}
}
H.push_back(L);
std::push_heap(H.begin(), H.end(), ComparePtr);
H.push_back(Rv);
std::push_heap(H.begin(), H.end(), ComparePtr);
if (!H.empty()) {
_mm_prefetch((const char*)H[0], _MM_HINT_T0);
if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0);
}
Interval* top = H.front();
float interval_len = top->x2 - top->x1;
const int T = (int)fmaf(-expm1f(p), 264.0f, 277.0f);
bool want_term = (exp2f(log2f(interval_len) / (float)dim) < eps) || (it == maxIter - 1);
if (!(it % T) || want_term) {
CtrlMsgND out;
out.kind = want_term ? 0 : 2;
if (!want_term) {
uint8_t cnt = (uint8_t)((H.size() >= 3) ? 3 : H.size());
out.multiXchg.count = cnt;
float* dest = out.multiXchg.intervals;
Interval* t1 = H[0];
Interval* t2 = (H.size() > 1 ? H[1] : H[0]);
Interval* t3 = (H.size() > 2 ? H[2] : H[H.size() - 1]);
Interval* tops[3] = { t1, t2, t3 };
for (uint8_t i2 = 0; i2 < cnt; ++i2) {
Interval* Tt = tops[i2];
dest[0] = Tt->x1;
dest[1] = 0.0f;
dest[2] = Tt->x2;
dest[3] = 0.0f;
dest[4] = Tt->R;
dest += 5;
}
}
for (int i2 = 0; i2 < world; ++i2)
if (i2 != rank)
g_world->isend(i2, 0, out);
if (want_term)
break;
}
if (!(it % 500) && !bestQ.empty()) {
CtrlMsgND out;
out.kind = 3;
out.bestSol.bestF = bestF;
out.bestSol.bestX = bestX;
out.bestSol.bestY = bestY;
out.bestSol.dim = (uint8_t)bestQ.size();
memcpy(out.bestSol.bestQ, bestQ.data(), bestQ.size() * sizeof(float));
for (int i2 = 0; i2 < world; ++i2)
if (i2 != rank)
g_world->isend(i2, 0, out);
}
while (g_world->iprobe(boost::mpi::any_source, 0)) {
CtrlMsgND in;
g_world->recv(boost::mpi::any_source, 0, in);
if (in.kind == 0) {
if (!rank)
break;
else
return;
}
else if (in.kind == 1) {
float sx = in.xchg.s_x1, ex = in.xchg.e_x1;
if (sx < 0.0f)
sx = 0.0f;
if (ex > 1.0f)
ex = 1.0f;
if (ex > sx) {
alignas(64) float tmp[32];
float tx, ty;
map.map01ToPoint(sx, tmp);
float y1 = cost(tmp, tx, ty);
map.map01ToPoint(ex, tmp);
float y2 = cost(tmp, tx, ty);
Interval* inj = new Interval(sx, ex, y1, y2, (float)dim);
inj->ChangeCharacteristic(r * Mmax);
if (!H.empty()) {
_mm_prefetch((const char*)H[0], _MM_HINT_T0);
if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0);
}
Interval* topH = H.front();
if (inj->R > 1.15f * topH->R) {
float p = fmaf(-1.0f / initial_len, dmax, 1.0f);
float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p), 0.3f) : fmaf(0.3491860241215959f, expm1f(p), 0.6f);
inj->R = in.xchg.Rtop * k;
H.emplace_back(inj);
std::push_heap(H.begin(), H.end(), ComparePtr);
}
}
}
else if (in.kind == 2) {
const MultiCrossMsg& mX = in.multiXchg;
for (uint8_t ii = 0; ii < mX.count; ++ii) {
const float* d = &mX.intervals[ii * 5];
float sx = d[0], ex = d[2];
if (sx < 0.0f)
sx = 0.0f;
if (ex > 1.0f)
ex = 1.0f;
if (ex > sx) {
alignas(64) float tmp[32];
float tx, ty;
map.map01ToPoint(sx, tmp);
float y1 = cost(tmp, tx, ty);
map.map01ToPoint(ex, tmp);
float y2 = cost(tmp, tx, ty);
Interval* inj = new Interval(sx, ex, y1, y2, (float)dim);
inj->ChangeCharacteristic(r * Mmax);
if (!H.empty()) {
_mm_prefetch((const char*)H[0], _MM_HINT_T0);
if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0);
}
Interval* topH = H.front();
if (inj->R > 1.15f * topH->R) {
float p = fmaf(-1.0f / initial_len, dmax, 1.0f);
float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p), 0.3f) : fmaf(0.3491860241215959f, expm1f(p), 0.6f);
inj->R = d[4] * k;
H.emplace_back(inj);
std::push_heap(H.begin(), H.end(), ComparePtr);
}
}
}
}
else if (in.kind == 3) {
const BestSolutionMsg& bm = in.bestSol;
if (bm.bestF < bestF * 1.15f) {
alignas(64) float tmp_q[32];
memcpy(tmp_q, bm.bestQ, bm.dim * sizeof(float));
float t_best = map.pointToT(tmp_q);
float t1 = fmaxf(a, t_best - 0.001f), t2 = fminf(b, t_best + 0.001f);
if (t2 > t1) {
alignas(64) float tq1[32], tq2[32];
float xx1, yy1, xx2, yy2;
map.map01ToPoint(t1, tq1);
float f1 = cost(tq1, xx1, yy1);
map.map01ToPoint(t2, tq2);
float f2 = cost(tq2, xx2, yy2);
Interval* I = new Interval(t1, t2, f1, f2, (float)dim);
I->ChangeCharacteristic(r * Mmax);
I->R *= 0.90f;
H.emplace_back(I);
std::push_heap(H.begin(), H.end(), ComparePtr);
}
if (bm.bestF < bestF) {
bestF = bm.bestF;
bestX = bm.bestX;
bestY = bm.bestY;
bestQ.assign(bm.bestQ, bm.bestQ + bm.dim);
}
}
}
}
++it;
}
}

struct BestPacket {
float bestF;
int dim;
float bestX;
float bestY;
template<typename Archive>
void serialize(Archive& ar, const unsigned int) {
ar& bestF& dim& bestX& bestY;
}
};

extern "C" __declspec(dllexport) __declspec(noalias) void AGP_Manip2D(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon,
unsigned int seed, float** out_bestQ, size_t* out_bestQLen, float* out_bestX,
float* out_bestY, float* out_bestF) {
const int dim = nSegments + (variableLengths ? nSegments : 0);
g_mc.permCache.resize(dim);
for (int i = 0; i < dim; ++i)
g_mc.permCache[i] = i;
uint32_t s = g_mc.baseSeed;
for (int i = dim - 1; i > 0; --i) {
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
uint32_t j = s % (uint32_t)(i + 1);
std::swap(g_mc.permCache[i], g_mc.permCache[j]);
}
g_mc.invMaskCache.resize(dim);
for (int k = 0; k < dim; ++k) {
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
g_mc.invMaskCache[k] = (uint64_t)s;
}
const float thetaMin = -3.14159265359f, thetaMax = 3.14159265359f, lenMin = 0.1f, lenMax = 2.0f;
std::vector<float> low;
low.reserve(dim);
std::vector<float> high;
high.reserve(dim);
for (int i = 0; i < nSegments; ++i) {
low.push_back(thetaMin);
high.push_back(thetaMax);
}
if (variableLengths) {
for (int i = 0; i < nSegments; ++i) {
low.push_back(lenMin);
high.push_back(lenMax);
}
}
ManipCost cost(nSegments, variableLengths, targetX, targetY, minTheta);
const int rank = g_world->rank(), world = g_world->size();
std::vector<float> bestQ;
float bestF = FLT_MAX, bx = 0.0f, by = 0.0f;
const int levels0 = (std::min)(peanoLevels, 8), maxIter0 = (int)(maxIterPerBranch * 0.2f);
MortonND map0(dim, levels0, low.data(), high.data(), g_mc);
std::vector<Interval*> H_coarse;
std::vector<float> bestQ_coarse;
float bestF_coarse = FLT_MAX, bx_coarse = 0.0f, by_coarse = 0.0f;
float M_prior = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << levels0));
if (variableLengths)
M_prior += 1.41421356237f * (1.0f / (float)(1u << levels0));
agp_run_branch_mpi(map0, cost, maxIter0, r, adaptiveMode, epsilon, seed, H_coarse, bestQ_coarse, bestF_coarse, bx_coarse, by_coarse, M_prior);
if (bestF_coarse < bestF) {
bestF = bestF_coarse;
bestQ = bestQ_coarse;
bx = bx_coarse;
by = by_coarse;
}
if (levels0 < peanoLevels) {
MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc);
std::vector<Interval*> H_fine;
std::vector<float> bestQ_fine = bestQ;
float bestF_fine = bestF, bx_fine = bx, by_fine = by;
float M_prior_fine = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << peanoLevels));
if (variableLengths)
M_prior_fine += 1.41421356237f * (1.0f / (float)(1u << peanoLevels));
if (!H_coarse.empty()) {
std::sort(H_coarse.begin(), H_coarse.end(), [](const Interval* a, const Interval* b) { return a->R < b->R; });
const size_t topCount = (size_t)(H_coarse.size() * 0.3f);
for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) {
const Interval* C = H_coarse[i];
alignas(64) float q1[32], q2[32];
float x1, y1, x2, y2;
map1.map01ToPoint(C->x1, q1);
float f1 = cost(q1, x1, y1);
map1.map01ToPoint(C->x2, q2);
float f2 = cost(q2, x2, y2);
H_fine.push_back(new Interval(C->x1, C->x2, f1, f2, (float)dim));
if (f1 < bestF_fine) {
bestF_fine = f1;
bestQ_fine.assign(q1, q1 + dim);
bx_fine = x1;
by_fine = y1;
}
if (f2 < bestF_fine) {
bestF_fine = f2;
bestQ_fine.assign(q2, q2 + dim);
bx_fine = x2;
by_fine = y2;
}
}
std::make_heap(H_fine.begin(), H_fine.end(), ComparePtr);
}
agp_run_branch_mpi(map1, cost, maxIterPerBranch - maxIter0, r, adaptiveMode, epsilon, seed, H_fine, bestQ_fine, bestF_fine, bx_fine, by_fine,
M_prior_fine);
if (bestF_fine < bestF) {
bestF = bestF_fine;
bestQ = bestQ_fine;
bx = bx_fine;
by = by_fine;
}
}
BestPacket me{ bestF, dim, bx, by };
if (!rank) {
std::vector<float> winnerQ = bestQ;
float winF = bestF, wx = bx, wy = by;
for (int i = 1; i < world; ++i) {
BestPacket bp;
g_world->recv(i, 2, bp);
std::vector<float> qin;
g_world->recv(i, 3, qin);
if (bp.bestF < winF) {
winF = bp.bestF;
wx = bp.bestX;
wy = bp.bestY;
winnerQ = qin;
}
}
*out_bestQLen = winnerQ.size();
out_bestQ = (float)CoTaskMemAlloc(sizeof(float) * (*out_bestQLen));
memcpy(*out_bestQ, winnerQ.data(), sizeof(float) * (*out_bestQLen));
*out_bestX = wx;
*out_bestY = wy;
*out_bestF = winF;
}
else {
g_world->send(0, 2, me);
g_world->send(0, 3, bestQ);
}
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline int AgpInit(int peanoLevel, float a, float b, float c, float d) {
g_env = new boost::mpi::environment();
g_world = new boost::mpi::communicator();
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
const int rank = g_world->rank();
const int world_size = g_world->size();
if (world_size == 4) {
new (&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, rank & 3);
}
g_mc.baseSeed = fmaf(0x9E3779B9u, rank, 0x9E3779B9u);
return rank;
}

__forceinline float ShekelFunc(float x, float seed) {
int i = 0;
float st = seed, r1, r2, res = 0.0f;
while (i < 10) {
XOR_RAND(st, r1);
float xp = fmaf(-r1, 10.0f, x);
XOR_RAND(st, r1);
XOR_RAND(st, r2);
float d = fmaf(fmaf(r1, 20.0f, 5.0f), xp * xp, fmaf(r2, 0.2f, 1.0f));
d = copysignf(fmaxf(fabsf(d), FLT_MIN), d);
res -= 1.0f / d;
++i;
}
return res;
}

__forceinline float RastriginFunc(float x1, float x2) {
const float t = fmaf(x1, x1, x2 * x2);
float c1, c2;
FABE13_COS(6.28318530717958647692f * x1, c1);
FABE13_COS(6.28318530717958647692f * x2, c2);
return (t - fmaf(c1 + c2, 10.0f, -14.6f)) * fmaf(-t, 0.25f, 18.42f);
}

__forceinline float HillFunc(float x, float seed) {
int j = 0;
__declspec(align(32)) float ang[14u];
float st = 6.28318530717958647692f * x;
while (j < 14) {
ang[j] = st * (float)(j + 1);
++j;
}
__declspec(align(32)) float sv[14u], cv[14u];
FABE13_SINCOS(ang, sv, cv, 14u);
float state = seed, r1, r2;
XOR_RAND(state, r1);
float res = fmaf(r1, 2.0f, -1.1f);
--j;
while (j >= 0) {
XOR_RAND(state, r1);
XOR_RAND(state, r2);
res += fmaf(fmaf(r1, 2.0f, -1.1f), sv[j], fmaf(r2, 2.0f, -1.1f) * cv[j]);
--j;
}
return res;
}

__forceinline float GrishaginFunc(float x1, float x2, float seed) {
int j = 0;
__declspec(align(32)) float aj[8u], ak[8u];
while (j < 8) {
float pj = 3.14159265358979323846f * (float)(j + 1);
aj[j] = pj * x1;
ak[j] = pj * x2;
++j;
}
__declspec(align(32)) float sj[8u], cj[8u], sk[8u], ck[8u];
FABE13_SINCOS(aj, sj, cj, 8u);
FABE13_SINCOS(ak, sk, ck, 8u);
--j;
float p1 = 0.0f, p2 = 0.0f;
float st = seed, r1, r2;
while (j >= 0) {
size_t k = 0u;
while (k < 8u) {
float s = sj[j] * sj[j];
float c = ck[k] * ck[k];
XOR_RAND_GRSH(st, r1);
XOR_RAND_GRSH(st, r2);
p1 = fmaf(r1, s, fmaf(r2, c, p1));
XOR_RAND_GRSH(st, r1);
XOR_RAND_GRSH(st, r2);
p2 = fmaf(-r1, c, fmaf(r2, s, p2));
++k;
}
--j;
}
return -sqrtf(fmaf(p1, p1, p2 * p2));
}

extern "C" __declspec(dllexport) __declspec(noalias) void AGP_1D(float global_iterations, float a, float b, float r, bool mode, float epsilon, float seed,
float** out_data, size_t* out_len) {
Slab* slab = tls.local();
slab->current = slab->base;
int schetchick = 0;
const float initial_length = b - a;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float inv_threshold_03 = 1.0f / threshold_03;
const float start_val = ShekelFunc(a, seed);
float best_f = ShekelFunc(b, seed);
float x_Rmax_1 = a;
float x_Rmax_2 = b;
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval*, boost::alignment::aligned_allocator<Interval*, 64u>> R;
Extr.reserve((size_t)global_iterations << 2u);
R.reserve((size_t)global_iterations << 1u);
R.emplace_back(new Interval(a, b, start_val, best_f, 1.0f));
float Mmax = R.front()->M;
float m = r * Mmax;
while (true) {
const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 1.0f, r);
const float new_value = ShekelFunc(new_point, seed);
if (new_value < best_f) {
best_f = new_value;
Extr.emplace_back(best_f);
Extr.emplace_back(new_point);
}
std::pop_heap(R.begin(), R.end(), ComparePtr);
const Interval* pro = R.back();
const float new_x1 = pro->x1;
const float new_x2 = pro->x2;
const float len2 = new_x2 - new_point;
const float len1 = new_point - new_x1;
const float interval_len = len1 < len2 ? len1 : len2;
if (++schetchick == (int)global_iterations || interval_len < epsilon) {
Extr.emplace_back((float)schetchick);
Extr.emplace_back(interval_len);
*out_len = Extr.size();
out_data = (float)CoTaskMemAlloc(sizeof(float) * (out_len));
memcpy(out_data, Extr.data(), sizeof(float) * (out_len));
return;
}
Interval curr = new Interval(new_x1, new_point, pro->y1, new_value, 1.0f);
Interval curr1 = new Interval(new_point, new_x2, new_value, pro->y2, 1.0f);
const float currM = curr->M > curr1->M ? curr->M : curr1->M;
const size_t r_size = R.size();
if (mode) {
if (len2 + len1 == dmax) {
dmax = len2 > len1 ? len2 : len1;
for (auto p : R) {
float L = p->x2 - p->x1;
if (L > dmax)
dmax = L;
}
}
if (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) {
if (currM > Mmax) {
Mmax = currM;
m = r * Mmax;
}
const float progress = fmaf(-inv_threshold_03, dmax, 1.0f);
const float alpha = progress * progress;
const float betta = 2.0f - alpha;
const float MULT = (1.0f / dmax) * Mmax;
const float global_coeff = fmaf(MULT, r, -MULT);
const float GF = betta * global_coeff;
curr->ChangeCharacteristic(fmaf(GF, len1, curr->M * alpha));
curr1->ChangeCharacteristic(fmaf(GF, len2, curr1->M * alpha));
RecomputeR_AffineM_AVX2(R.data(), r_size, GF, alpha);
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
R.back() = curr;
std::push_heap(R.begin(), R.end(), ComparePtr);
R.emplace_back(curr1);
std::push_heap(R.begin(), R.end(), ComparePtr);
const Interval top = R.front();
x_Rmax_1 = top->x1;
x_Rmax_2 = top->x2;
y_Rmax_1 = top->y1;
y_Rmax_2 = top->y2;
}
}

extern "C" __declspec(dllexport) __declspec(noalias) void AGP_2D(const float N, const float global_iterations, const float a, const float b, const float c,
const float d, const float r, const bool mode, const float epsilon, const float seed,
float** const __restrict out_data, size_t* const __restrict out_len) noexcept {
Slab* const __restrict slab = tls.local();
slab->current = slab->base;
int schetchick = 0;
int no_improve = 0;
const int rank = g_world->rank();
const int world_size = g_world->size();
while (g_world->iprobe(boost::mpi::any_source, 0)) {
CtrlMsg dummy;
g_world->recv(boost::mpi::any_source, 0, dummy);
}
const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1));
const float x_addition = (b - a) * inv_divider;
const float y_addition = (d - c) * inv_divider;
const float true_start = a + x_addition;
const float true_end = b - x_addition;
float x_Rmax_1 = true_start;
float x_Rmax_2 = true_end;
const float initial_length = x_Rmax_2 - x_Rmax_1;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float inv_threshold_03 = 1.0f / threshold_03;
const float start_val = rank % 3 ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition);
float best_f = rank % 2 ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition);
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __restrict, 64u>> R;
Extr.clear();
Extr.reserve(static_cast<size_t>(global_iterations) << 2u);
R.clear();
R.reserve(static_cast<size_t>(global_iterations) << 1u);
R.emplace_back(new Interval(true_start, true_end, start_val, best_f, 2.0f));
const Interval* __restrict top_ptr;
float Mmax = R.front()->M;
float m = r * Mmax;
while (true) {
const float interval_len = x_Rmax_2 - x_Rmax_1;
const bool stagnation = no_improve > 100 && schetchick > 270;
const float p = fmaf(-1.0f / initial_length, dmax, 1.0f);
while (g_world->iprobe(boost::mpi::any_source, 0)) {
CtrlMsg in;
g_world->recv(boost::mpi::any_source, 0, in);
if (in.kind) {
if (!rank) {
Extr.emplace_back(static_cast<float>(schetchick));
Extr.emplace_back(interval_len);
*out_len = Extr.size();
out_data = reinterpret_cast<float __restrict>(CoTaskMemAlloc(sizeof(float) * (*out_len)));
memcpy(*out_data, Extr.data(), sizeof(float) * (out_len));
}
return;
}
const float sx = FindX2D_analytic(in.xchg.s_x1, in.xchg.s_x2);
const float ex = FindX2D_analytic(in.xchg.e_x1, in.xchg.e_x2);
Interval const __restrict injected = new Interval(sx, ex, RastriginFunc(in.xchg.s_x1, in.xchg.s_x2), RastriginFunc(in.xchg.e_x1, in.xchg.e_x2), 2.0f);
injected->ChangeCharacteristic(m);
if (injected->R > 1.15f * top_ptr->R) {
const float k = stagnation ? fmaf(0.5819767068693265f, expm1f(p), 0.3f) : fmaf(0.3491860241215959f, expm1f(p), 0.6f);
injected->R = in.xchg.Rtop * k;
R.emplace_back(injected);
std::push_heap(R.begin(), R.end(), ComparePtr);
}
}
const int T = static_cast<int>(fmaf(-expm1f(p), 264.0f, 277.0f));
const bool want_term = interval_len < epsilon || schetchick == static_cast<int>(global_iterations);
if (!(++schetchick % T) || stagnation || want_term) {
CtrlMsg out;
out.kind = want_term;
if (!out.kind) {
float s_x1, s_x2, e_x1, e_x2;
HitTest2D_analytic(top_ptr->x1, s_x1, s_x2);
HitTest2D_analytic(top_ptr->x2, e_x1, e_x2);
out.xchg = CrossMsg{ s_x1, s_x2, e_x1, e_x2, top_ptr->R };
}
for (int i = 0; i < world_size; ++i) {
if (i != rank)
g_world->isend(i, 0, out);
}
if (out.kind) {
if (!rank) {
Extr.emplace_back(static_cast<float>(schetchick));
Extr.emplace_back(interval_len);
*out_len = Extr.size();
out_data = reinterpret_cast<float __restrict>(CoTaskMemAlloc(sizeof(float) * (out_len)));
memcpy(out_data, Extr.data(), sizeof(float) * (out_len));
}
return;
}
}
const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 2.0f, r);
float new_x1_val, new_x2_val;
HitTest2D_analytic(new_point, new_x1_val, new_x2_val);
const float new_value = RastriginFunc(new_x1_val, new_x2_val);
if (new_value < best_f) {
best_f = new_value;
Extr.emplace_back(best_f);
Extr.emplace_back(new_x1_val);
Extr.emplace_back(new_x2_val);
no_improve = 0;
}
else {
++no_improve;
}
std::pop_heap(R.begin(), R.end(), ComparePtr);
Interval const __restrict promejutochny_otrezok = R.back();
const float segment_x1 = promejutochny_otrezok->x1;
const float segment_x2 = promejutochny_otrezok->x2;
const float len2 = segment_x2 - new_point;
const float len1 = new_point - segment_x1;
Interval const __restrict curr = new Interval(segment_x1, new_point, promejutochny_otrezok->y1, new_value, 2.0f);
Interval const __restrict curr1 = new Interval(new_point, segment_x2, new_value, promejutochny_otrezok->y2, 2.0f);
const float currM = (std::max)(curr->M, curr1->M);
const size_t r_size = R.size();
if (mode) {
if (len2 + len1 == dmax) {
dmax = (std::max)(len1, len2);
for (auto pI : R) {
float L = pI->x2 - pI->x1;
if (L > dmax)
dmax = L;
}
}
if (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) {
if (currM > Mmax) {
Mmax = currM;
m = r * Mmax;
}
const float progress = fmaf(-inv_threshold_03, dmax, 1.0f);
const float alpha = progress * progress;
const float betta = 2.0f - alpha;
const float MULTIPLIER = (1.0f / dmax) * Mmax;
const float global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER);
const float GLOBAL_FACTOR = betta * global_coeff;
curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha));
curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha));
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
R.back() = curr;
std::push_heap(R.begin(), R.end(), ComparePtr);
R.emplace_back(curr1);
std::push_heap(R.begin(), R.end(), ComparePtr);
top_ptr = R.front();
x_Rmax_1 = top_ptr->x1;
x_Rmax_2 = top_ptr->x2;
y_Rmax_1 = top_ptr->y1;
y_Rmax_2 = top_ptr->y2;
}
}

struct RunParams {
int nSegments;
uint8_t varLen;
float minTheta;
float tx, ty;
int levels, maxIter;
float r;
uint8_t adaptive;
float eps;
unsigned int seed;
template<typename Archive>
void serialize(Archive& ar, const unsigned int) {
ar& nSegments& varLen& minTheta& tx& ty& levels& maxIter& r& adaptive& eps& seed;
}
};

extern "C" __declspec(dllexport) __declspec(noalias) void AgpStartManipND(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon,
unsigned int seed) {
RunParams p;
p.nSegments = nSegments;
p.varLen = (uint8_t)variableLengths;
p.minTheta = minTheta;
p.tx = targetX;
p.ty = targetY;
p.levels = peanoLevels;
p.maxIter = maxIterPerBranch;
p.r = r;
p.adaptive = (uint8_t)adaptiveMode;
p.eps = epsilon;
p.seed = seed;
const int rank = g_world->rank();
const int world = g_world->size();
for (int i = 1; i < world; ++i) {
if (i != rank)
g_world->isend(i, 1, p);
}
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun() {
while (true) {
if (g_world->iprobe(0, 1)) {
RunParams p;
g_world->recv(0, 1, p);
float* q = nullptr;
size_t qlen = 0;
float bx = 0, by = 0, bf = 0;
AGP_Manip2D(p.nSegments, (bool)p.varLen, p.minTheta, p.tx, p.ty, p.levels, p.maxIter, p.r, (bool)p.adaptive, p.eps, p.seed, &q, &qlen, &bx, &by, &bf);
if (q)
CoTaskMemFree(q);
}
Sleep(0);
}
}

extern "C" __declspec(dllexport) __declspec(noalias) void AgpStartWorkers() {
int i = 1;
const int world = g_world->size();
while (i < world) {
g_world->isend(i, 1, 0);
++i;
}
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* p) {
CoTaskMemFree(p);
}" - это мой алгоритм глобального поиска - рассмотри именно его многомерную реализацию и глубоко проанализируй саму вычислительную схему - верна ли она? что можно предпринять чтобы улучшить сходимость? укажи на слабые места многомерной схемы которые могут препятствовать сходимости, дай мне конкретные рецепты и внедрения в код - что нужно опробовать в первую очередь, также заметь что одномерные и двумерные алгоритмы уже практически вылизаны до идеала - ты можешь ориентироваться на них, также deepseek дал мне следующий ответ на тот же вопрос: "Глубокий анализ вашего многомерного алгоритма глобального поиска показывает несколько фундаментальных проблем, препятствующих сходимости:

Критические проблемы сходимости

  1. Некорректное использование кривых Пеано
    cpp
    __forceinline void map01ToPoint(float t, float* __restrict out) const noexcept
    Проблема: Кривые Пеано сохраняют лишь частичную близость в многомерном пространстве. При отображении 1D→ND близкие точки в параметре t могут оказаться далеко в исходном пространстве, что нарушает основное предположение метода Липшица.

Решение: Используйте отображения с лучшими свойствами сохранения близости (кривые Гильберта) или перейдите к настоящему многомерному поиску.

  1. Некорректная характеристика R в многомерном случае
    cpp
    R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor));
    Проблема: Характеристика вычисляется для 1D интервалов, но используется для многомерной функции. Это математически некорректно - константа Липшица в ND пространстве не сводится к 1D оценке.

  2. Неадекватная оценка константы Липшица
    cpp
    M = (1.0f / N_factor)* fabsf(delta_y)
    Проблема: Для многомерных функций оценка М должна учитывать градиент по всем направлениям, а не только разность значений.

Конкретные улучшения

  1. Переход к настоящему многомерному поиску
    cpp
    struct NDInterval {
    std::vector<float> min, max; // Границы гиперпрямоугольника
    float f_min, f_max; // Лучшее и худшее значение в области
    float diameter; // Диаметр области
    // Многомерная характеристика
    double compute_characteristic(float L) const {
    return f_min - L * diameter; // Для минимизации
    }
    };

  2. Корректная оценка константы Липшица
    cpp
    float estimate_lipschitz_nd(const std::vector<NDInterval>& intervals) {
    float max_gradient = 0.0f;
    for (const auto& interval : intervals) {
    // Оцениваем градиент по разностям в соседних точках
    float grad_est = (interval.f_max - interval.f_min) / interval.diameter;
    max_gradient = std::max(max_gradient, grad_est);
    }
    return max_gradient * safety_factor;
    }

  3. Улучшенная стратегия деления
    cpp
    void split_nd_interval(const NDInterval& original,
    std::vector<NDInterval>& result) {
    // Делим по самой длинной координате
    int split_dim = 0;
    float max_length = 0.0f;
    for (int i = 0; i < original.min.size(); ++i) {
    float length = original.max[i] - original.min[i];
    if (length > max_length) {
    max_length = length;
    split_dim = i;
    }
    }

    // Создаем два под-интервала
    NDInterval left = original, right = original;
    float mid = (original.min[split_dim] + original.max[split_dim]) * 0.5f;

    left.max[split_dim] = mid;
    right.min[split_dim] = mid;

    // Переоцениваем функции в новых интервалах
    reevaluate_interval(left);
    reevaluate_interval(right);

    result.push_back(left);
    result.push_back(right);
    }

  4. Локальная оптимизация с учетом многомерности
    cpp
    void local_refinement_nd(const std::vector<float>& point,
    ManipCost& cost,
    std::vector<float>& improved) {
    // Используйте квази-Ньютоновские методы (BFGS)
    // или trust-region методы для локального улучшения
    const float epsilon = 1e-6f;
    std::vector<float> gradient = compute_numerical_gradient(cost, point, epsilon);

    // Простой градиентный спуск с ограничениями
    improved = point;
    float step_size = 0.01f;
    for (int iter = 0; iter < 100; ++iter) {
    std::vector<float> new_point = improved;
    for (int i = 0; i < new_point.size(); ++i) {
    new_point[i] -= step_size * gradient[i];
    // Проекция на допустимую область
    new_point[i] = std::max(lower_bounds[i],
    std::min(upper_bounds[i], new_point[i]));
    }

    text
     float new_value = cost(new_point.data(), ...);
 if (new_value < cost(improved.data(), ...)) {
     improved = new_point;
     gradient = compute_numerical_gradient(cost, improved, epsilon);
 } else {
     step_size *= 0.5f;
 }

    }
    }
    Рекомендуемая архитектура

  5. Многомерное дерево поиска
    cpp
    class NDGlobalSearch {
    private:
    std::priority_queue<NDInterval> active_intervals_;
    std::vector<float> best_point_;
    float best_value_;

public:
void search(int max_iterations) {
while (!active_intervals_.empty() && iterations < max_iterations) {
NDInterval current = active_intervals_.top();
active_intervals_.pop();

text
        if (current.promising(best_value_)) {
            auto children = split_interval(current);
            for (auto& child : children) {
                if (child.f_min < best_value_) {
                    best_value_ = child.f_min;
                    best_point_ = child.best_point;
                }
                if (child.promising(best_value_)) {
                    active_intervals_.push(child);
                }
            }
        }
    }
}

};
2. Баланс глобального и локального поиска
Используйте критерий current.diameter > epsilon для глобального поиска

При малых диаметрах переключайтесь на локальную оптимизацию

Комбинируйте с мета-эвристиками (мультистарт, байесовская оптимизация)

  1. Параллельная стратегия для MPI
    cpp
    // Вместо обмена 1D интервалами
    void exchange_promising_regions() {
    // Обмен гиперпрямоугольниками с лучшими характеристиками
    // Координационная информация о многообещающих регионах
    }
    Вывод
    Текущая схема фундаментально ограничена использованием 1D методов для многомерных задач. Для достижения сходимости необходимо перейти к настоящим многомерным методам ветвей и границ с корректной оценкой характеристик и стратегией деления многомерных областей." и потом ещё "Понимаю! Вы используете АГП Стронгина с редукцией размерности через кривые Мортона. Это действительно классический подход, но даже в этой схеме есть несколько критических мест, влияющих на сходимость.

Анализ текущей реализации АГП с редукцией

  1. Проблемы с отображением Мортона
    cpp
    __forceinline void map01ToPoint(float t, float* __restrict out) const noexcept
    Проблема: Кривая Мортона имеет худшие свойства сохранения близости, чем кривая Гильберта. Близкие точки в 1D могут соответствовать далеким точкам в ND пространстве.

Решение:

cpp
// Замена Мортона на Гильберта
class HilbertND {
// Реализация кривой Гильберта с лучшими свойствами непрерывности
void map01ToPoint(float t, float* out) const;
};
2. Некорректная адаптивная оценка константы Липшица
cpp
__forceinline void RecomputeR_AffineM_AVX2(Interval* const* arr, size_t n, float GF, float alpha)
Проблема: Адаптивная оценка M не учитывает геометрию многомерного пространства. В ND пространстве изменение функции зависит от направления.

Улучшение:

cpp
float estimate_effective_lipschitz(const Interval* interval,
const MortonND& map,
const ManipCost& cost) {
// Оцениваем градиент в нескольких точках вдоль кривой
float t_mid = 0.5f * (interval->x1 + interval->x2);
float t_delta = 0.01f * (interval->x2 - interval->x1);

text
alignas(32) float q1[32], q2[32];
float x1, y1, x2, y2;

map.map01ToPoint(t_mid - t_delta, q1);
float f1 = cost(q1, x1, y1);
map.map01ToPoint(t_mid + t_delta, q2);
float f2 = cost(q2, x2, y2);

// Учитываем также длину пути в исходном пространстве
float path_length = compute_nd_distance(q1, q2, map.dim);
return fabsf(f2 - f1) / (path_length + 1e-8f);

}
3. Неоптимальная стратегия деления
Проблема: Деление всегда происходит по параметру t, но не учитывается геометрия целевой функции.

Улучшение - адаптивное направление деления:

cpp
struct EnhancedInterval : public Interval {
std::vector<float> sensitivity; // Чувствительность по разным измерениям
int best_split_dim; // Направление для деления
};

void compute_sensitivity(EnhancedInterval* interval,
const MortonND& map,
const ManipCost& cost) {
// Оцениваем, вдоль каких измерений функция меняется сильнее
float t_mid = 0.5f * (interval->x1 + interval->x2);
alignas(32) float q_base[32], q_test[32];
map.map01ToPoint(t_mid, q_base);

text
for (int dim = 0; dim < map.dim; ++dim) {
    memcpy(q_test, q_base, map.dim * sizeof(float));
    float delta = 0.01f * (map.high[dim] - map.low[dim]);
    q_test[dim] += delta;
    
    float fx, fy, fdx, fdy;
    cost(q_base, fx, fy);
    cost(q_test, fdx, fdy);
    
    interval->sensitivity[dim] = fabsf(fdx - fx) / delta;
}

}
4. Слабая локальная оптимизация
cpp
// Текущая реализация слишком простая
float eta = 0.1f;
for (int step = 0; step < 2; ++step) {
// ...
}
Улучшение - более агрессивная локальная оптимизация:

cpp
void enhanced_local_optimization(float* q, int dim, ManipCost& cost,
int max_iter = 50) {
float current_f, x, y;
float best_f = cost(q, x, y);
float best_q[32];
memcpy(best_q, q, dim * sizeof(float));

text
// Используем адаптивный шаг
float step_size = 0.1f;
for (int iter = 0; iter < max_iter; ++iter) {
    // Численный градиент
    alignas(32) float grad[32] = {0};
    compute_numerical_gradient(cost, q, dim, grad);
    
    // Проекция градиента с учетом ограничений
    project_gradient(q, grad, dim, cost.minTheta);
    
    // Поиск с регулировкой шага
    bool improved = false;
    for (int ls = 0; ls < 10; ++ls) {
        float new_q[32];
        for (int i = 0; i < dim; ++i) {
            new_q[i] = q[i] - step_size * grad[i];
        }
        project_to_feasible(new_q, dim, cost);
        
        float new_f, new_x, new_y;
        new_f = cost(new_q, new_x, new_y);
        
        if (new_f < best_f - 1e-6f) {
            best_f = new_f;
            memcpy(best_q, new_q, dim * sizeof(float));
            memcpy(q, new_q, dim * sizeof(float));
            step_size *= 1.2f; // Увеличиваем шаг
            improved = true;
            break;
        }
    }
    
    if (!improved) {
        step_size *= 0.5f;
        if (step_size < 1e-8f) break;
    }
}

memcpy(q, best_q, dim * sizeof(float));

}
5. Неэффективный обмен информацией в MPI
cpp
CtrlMsgND out;
out.kind = want_term ? 0 : 2;
Проблема: Обмен только лучшими интервалами, но не информацией о многообещающих регионах в исходном пространстве.

Улучшение:

cpp
struct RegionInfo {
float center[32]; // Центр региона в исходном пространстве
float radius; // Радиус региона
float f_value; // Значение функции
int dim; // Размерность
};

// Обмен информацией о перспективных регионах
void exchange_promising_regions(const std::vector<Interval*>& H,
const MortonND& map) {
std::vector<RegionInfo> promising_regions;

text
// Выбираем топ-3 интервала и преобразуем в исходное пространство
for (int i = 0; i < std::min(3, (int)H.size()); ++i) {
    RegionInfo region;
    float t_center = 0.5f * (H[i]->x1 + H[i]->x2);
    map.map01ToPoint(t_center, region.center);
    region.radius = 0.5f * (H[i]->x2 - H[i]->x1);
    region.f_value = 0.5f * (H[i]->y1 + H[i]->y2);
    region.dim = map.dim;
    promising_regions.push_back(region);
}

// Рассылаем другим процессам
// ...

}
6. Улучшенная инициализация
cpp
// Текущая эвристика может быть улучшена
__forceinline int generate_heuristic_seeds(const ManipCost& cost, int dim,
float* S, int stride)
Улучшение - использование Latin Hypercube Sampling:

cpp
void latin_hypercube_sampling(int n_samples, int dim,
const float* low, const float* high,
float* samples) {
// Реализация LHS для лучшего покрытия пространства
for (int d = 0; d < dim; ++d) {
std::vector<float> positions(n_samples);
for (int i = 0; i < n_samples; ++i) {
positions[i] = (i + 0.5f) / n_samples;
}
std::shuffle(positions.begin(), positions.end(), std::mt19937{});
for (int i = 0; i < n_samples; ++i) {
samples[i * dim + d] = low[d] + positions[i] * (high[d] - low[d]);
}
}
}
Ключевые рекомендации:
Замените Мортона на Гильберта - это даст немедленное улучшение сходимости

Улучшите локальную оптимизацию - текущая слишком слабая для многомерных задач

Введите адаптивное направление деления - учитывайте геометрию целевой функции

Улучшите обмен информацией в MPI - передавайте информацию о регионах в исходном пространстве

Используйте LHS для инициализации - лучшее покрытие начального пространства

Эти изменения должны значительно улучшить сходимость без изменения фундаментальной схемы АГП с редукцией размерности." какие из этих советов действительно правда и к каким можно прислушаться? я не хотел бы отказываться от кривых мортона так как и кривые пеано тоже не сохраняют близость точек в многомере - так что этот совет кажется мне бессмысленным - особенно учитывая то что с кривыми мортона я удобно работаю через побитовые сдвиги, рассмотри многомерную схему в целом как бы издалека - я не могу изменить базовый принцип алгоритма - многомерный агп стронгина с редукцией размерности, я должен по сути использовать именно одномерный поиск всегда с редукцией размерности, но в остальном полная творческая свобода - информацию о градиенте можно использовать - подготовь мне конкретные рецепты которые стоит опробовать и подробно объясни почему - думай долго для более корректного ответа, дольше чем ты обычно думаешь - для меня это очень важно, жду конкретные рекомендации по внедрению в код + подробные объяснения почему это стоит внедрить, мой приоритет - скорость работы алгоритма + гарантия сходимости

Thought

Exploring the N-factor calculation

Okay, the N_factor fudge uses different exponent formulas depending on the dimension, like len^(1/N). They calculate 'R' based on that and a quadratic term. This seems to follow the principles from Hölder's extension. The use of N in calculations matches the generalization of reducing ND to 1D.

Clarifying the mapping and consistency

The mapping being used is Morton instead of strictly Peano, but it's still space-filling, meaning the overall approach holds. The Hölder exponent remains 1/N. To estimate the constant H, they calculate M using |Δf| and N_factor, which is consistent with the scheme, so the issues with R being incorrect might not be valid.

Addressing potential weaknesses of the mapping

Okay, Morton is being used here instead of a continuous Peano curve, which causes some discontinuity at cell boundaries. The mapping from t to ND can have jumps because of how it's discretized. However, Hölder continuity remains consistent if the mapping is handled as a continuous process despite these discontinuities.

Examining the impact of Morton versus continuous curves

Morton introduces discontinuity in mapping, unlike continuous Peano curves like Hilbert or Lebesgue. This leads to large variations in f(t) and potentially impacts Lipschitz/Hölder assumptions. Though the code works with cell centers in Morton order, the decomposition may slow convergence due to less local proximity compared to alternatives.

Refining the Morton approach

The user prefers using Morton, so Gray-code ordering could reduce jumps—though it doesn't guarantee Euclidean adjacency. For estimating Lipschitz constants, the code leverages local gradients and discrete steps, but Morton’s non-continuous nature complicates things. There’s potential to improve by using ND diameter instead of 1D length for bounding, enhancing robustness while keeping Morton’s bitwise benefits intact. This could involve modifying N_factor.

Optimizing with Morton

To keep Morton, we'll first compute a cheap upper bound for the ND diameter of a given interval. By using bitwise operations and precomputed step arrays, we can estimate the bounding box. This allows us to treat the ND diameter as a key metric for Lipschitz bounds, thereby improving our interval estimates without changing the core algorithm.

Improving algorithm efficiency

To improve interval computation and optimization, I’ll propose adjustments like using Morton indices with Gray coding and optimizing M estimation with robust quantile smoothing. I’ll also enhance seeding with Latin hypercube sampling for better initial points and consider local optimizer improvements. For MPI messaging, I'll suggest sending bounding boxes with diameters for better remote computation without changing too much.

Tuning parameters and optimizations

For better efficiency, I’m proposing adjustments like adapting the r parameter schedule to shrink as necessary and leveraging 64-bit scaling for large dimensions in Morton encoding. I’ll also focus on interval diameter computation, gradient-based updates, and optimization techniques like local refinement and robust seeding. I’ll integrate Mawby's ideas with necessary code tweaks and boundary constraints.

Optimizing algorithm

I’m accepting changes that improve local optimization, seeding, and MPI region exchange, along with gradient-driven M updates. Specifically, I’ll focus on incorporating Gray-coded Morton mapping, calculating diameter-aware R, and refining the local optimizer. I'll provide code patches that include updates for MortonND, diameter calculations, and adjusted Interval struct. Additionally, I’ll prioritize optimizing M and the gradient-based updates while ensuring low computational cost per interval.

Optimizing Shag and Interval functions

I’m considering modifying the Shag function to include diameter-based characteristics. We'll compute _N based on interval.N_factor. Changing the Interval struct and Gray-encoding Morton indices seems useful for better locality. I’ll adjust mappings and diameter calculations; ensure compatibility without breaking theoretical results. The changes should focus on simplifying code with toggles for "diameter mode" and Gray-coded indexing, optimizing performance and spatial coherence.

Thought for 2m 29s

Answer

Коротко: ваша многомерная схема — это ровно «АГП Стронгина по Пеано/Мортону» с гёльдеровской степенью α=1/dim\alpha = 1/\text{dim} (у вас это через N_factor = (x2-x1)^{1/N}) и с характеристикой вида
R=(Δf)2m dα+m dα−2(f1+f2)R = \frac{(\Delta f)^2}{m\,d^\alpha} + m\,d^\alpha - 2(f_1+f_2).
Это математически корректная редукция если характеристика использует «правильную» размерность интервала. Проблема вашей реализации не в «самой идее редукции», и не в Мортоне как таковом, а в том, что в ND варианте вы оцениваете «размер» области только через 1D-длину по tt. Для Пеано/Хильберта это терпимо (они непрерывны), а для Моронта — скачки соседних tt часто означают далекие точки в Rd\mathbb{R}^d. В результате RR получается слишком оптимистичным для «плохих» отрезков и вы тратите итерации на неэффективные деления.

Ниже — конкретные изменения, которые сохраняют базовый принцип (всегда 1D-поиск по tt), не требуют отказа от Моронта, и дают на практике ощутимый прирост сходимости/скорости без утраты гарантий.

Что из советов deepseek имеет смысл

  • ❌ «Переходите на настоящий ND-поиск, откажитесь от Пеано/Моронта» — не подходит вашим требованиям.
  • ⚠️ «Смените Моронта на Хильберта» — не обязательно. Можно получить почти тот же выигрыш, оставшись на Моронте, если исправить метрику в RR и/или применить Gray-кодирование индекса (см. ниже).
  • ✅ «Используйте информацию о градиенте для оценки константы Липшица» — да, и это у вас почти «на ладони»: вы уже считаете производные для локального шага.
  • ✅ «Улучшите локальную догладку» — да, но удерживаем её очень лёгкой (5–10 итераций с шагом Барзилай–Борвейна и проекцией на коробку).
  • ✅ «Лучшие сиды (LHS)» — да, 4–8 дешёвых сидов на грубом уровне дают частые выигрыши.

Рецепты (в порядке приоритета) + точечные патчи

1) ND-осведомлённая характеристика RR: заменяем 1D-длину на ND-диаметр блока Моронта

Идея. Интервал по tt соответствует блоку индексoв Моронта. Для этого блока можно быстро (чисто побитово) верхне ограничить пер-координатный размах и, значит, диаметр гиперпрямоугольника в исходном пространстве. Подставим этот диаметр DD вместо (x2−x1)1/N(x_2-x_1)^{1/N} в формулу RR и в оценку MM.

Быстрый подсчёт диаметра. Пусть для подотрезка индексы i1,i2i_1, i_2 (в Моронте). Для каждого измерения dd возьмём ваш pextMask[d] и посчитаем, сколько разрядов этой координаты различаются между i1i_1 и i2i_2:

cpp
// внутри MortonND
__forceinline float block_diameter(uint64_t i1, uint64_t i2) const noexcept {
    if (i1 > i2) std::swap(i1, i2);
    double s2 = 0.0;
    for (int d = 0; d < dim; ++d) {
        const int pd = perm[d];
        const uint64_t varying = (i1 ^ i2) & pextMask[d];
        const int nfree = _mm_popcnt_u64(varying);              // число "свободных" уровней по измерению d
        const double range = (double)step[pd] * ((nfree ? ((1ull << nfree) - 1ull) : 0ull));
        s2 += range * range;
    }
    return (float)std::sqrt(s2);
}

Это верхняя оценка реального размаха координат среди всех ячеек, которые может покрыть отрезок индексов. Дёшево (несколько popcnt) и стабильно.

Минимальные правки в структуре интервала. Чтобы не ломать остальной код, ослабьте const у трёх полей и добавьте «метрику»:

cpp
// было: const float ... N_factor, quadratic_term, M;
struct Interval final {
    float x1, x2, y1, y2, delta_y, ordinate_factor;
    float N_factor, quadratic_term, M;   // ← снимаем const
    float R;
    uint64_t i1, i2;
    float diam;
    // ...
    __forceinline Interval(float _x1, float _x2, float _y1, float _y2) noexcept
        : x1(_x1), x2(_x2), y1(_y1), y2(_y2), delta_y(_y2 - _y1),
          ordinate_factor(-(y1 + y2) * 2.0f), N_factor(0), quadratic_term(0), M(0), R(0), i1(0), i2(0), diam(0) {}
    __forceinline void set_metric(float d_alpha) noexcept {
        // d_alpha = D^alpha, для альфа=1 просто D
        N_factor = d_alpha;
        quadratic_term = (1.0f / N_factor) * delta_y * delta_y;
        M = (1.0f / N_factor) * fabsf(delta_y);
    }
    __forceinline void ChangeCharacteristic(float m) noexcept {
        R = fmaf(1.0f / m, quadratic_term, fmaf(m, N_factor, ordinate_factor));
    }
};

Где заполнять. В agp_run_branch_mpi(...) сразу после создания интервала:

cpp
auto t_to_idx = [&](float t) -> uint64_t {
    float tt = t <= 0.0f ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t);
    uint64_t idx = (uint64_t)((double)tt * (double)map.scale);
    if (idx >= map.scale) idx = map.scale - 1ull;
#if AGP_USE_GRAY_ORDER
    idx = gray_encode(idx);
#endif
    return idx;
};

// пример для I на [t1, t2]
Interval* I = new Interval(t1, t2, f1, f2);
I->i1 = t_to_idx(t1);
I->i2 = t_to_idx(t2);
I->diam = map.block_diameter(I->i1, I->i2);
I->set_metric(/* d_alpha = */ I->diam /* α = 1 */);
Mmax = std::max(Mmax, I->M);
I->ChangeCharacteristic(r_eff * Mmax);

Почему это важно. Для любой Липшицевой ff в Rd\mathbb{R}^d: ∣f(x)−f(y)∣≤L∥x−y∥|f(x)-f(y)| \le L \|x-y\|. Отсюда для вашего 1D отрезка по tt, соответствующего ND блоку, гарантированная оценка идёт через ND-диаметр блока. Даже с «рваными» соседями Моронта ваша характеристика становится корректной и стабильной, потому что «размер области» учтён по геометрии задачи, а не по формальной длине tt.

Гарантии. Метод Стронгина с m≥Lm \ge L и такой характеристикой сохраняет нижние оценки и те же условия сходимости. Вы лишь меняете «мерку» dαd^\alpha на DD, что для Липшица — естественно.

2) Gray-кодирование индекса Моронта (не меняя Моронта)

Это дёшево и сильно сглаживает последовательность по tt (соседние индексы отличаются в 1 бите ⇒ меньше «прыжков» по координатам).

cpp
static __forceinline uint64_t gray_encode(uint64_t x) noexcept { return x ^ (x >> 1); }
static __forceinline uint64_t gray_decode(uint64_t g) noexcept {
    // стандартная инверсия Грея
    g ^= g >> 32; g ^= g >> 16; g ^= g >> 8; g ^= g >> 4; g ^= g >> 2; g ^= g >> 1;
    return g;
}

// В MortonND:
bool use_gray;
__forceinline void map01ToPoint(float t, float* __restrict out) const noexcept {
    // ...
    uint64_t idx = (uint64_t)((double)t * (double)scale);
    if (idx >= scale) idx = scale - 1ull;
    if (use_gray) idx = gray_encode(idx);
    // дальше как было: _pext_u64(idx, pextMask[d]) ...
}

__forceinline float pointToT(const float* __restrict q) const noexcept {
    // как у вас, собираем idx через _pdep_u64(b, pextMask[d])
    uint64_t idx = ...;
    if (use_gray) idx = gray_decode(idx);
    return ((float)idx + 0.5f) / (float)scale;
}

В конструкторе: use_gray = true; (или флагом).

Почему это стоит делать. Это почти даром, а «локальность по tt» заметно улучшается даже для Моронта. В сочетании с пунктом 1) — особенно эффективно.

3) Градиент-драйв обновления MM (почти бесплатно)

Вы уже считаете суммарные синусы/косинусы для манипулятора и делаете 1–2 локальных шага. Из этого блока можно почти даром получить верхнюю оценку локальной константы Липшица и агрессивнее (но безопасно) поднимать Mmax⁡M_{\max}.

Добавьте внутри evalAt после блока, где вы считаете sum_s/sum_c, следующее:

cpp
// оценка ||∇f|| в физических единицах (углы в рад, длины в тех же ед.)
float grad2 = 0.0f;

// по углам
for (int i = 0; i < n; ++i) {
    float gpen = 0.0f;
    float ai = fabsf(q_local[i]);
    float v = cost.minTheta - ai;
    if (v > 0.0f) gpen = -2.0f * copysignf(v, q_local[i]); // как у вас
    float gtheta = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen;
    grad2 = fmaf(gtheta, gtheta, grad2);
}

// по длинам (если включены)
if (cost.variableLen) {
    for (int i = 0; i < n; ++i) {
        float gL = (dx * c_arr[i] + dy * s_arr[i]) / dist;
        grad2 = fmaf(gL, gL, grad2);
    }
}

float Lcand = sqrtf(grad2) + 1e-8f;

// Плавное, но неубывающее обновление
Mmax = (std::max)(Mmax, 0.9f * Mmax + 0.1f * Lcand);

Почему это корректно. Для Липшицевой ff L≥∥∇f(x)∥L \ge \|\nabla f(x)\| почти всюду (гладкость у вас есть). Значит Lcand — валидная «локальная» нижняя оценка для глобального LL. Включая её в Mmax, вы ускоряете зажатие без потери гарантии (мы берём максимум).

Синхронизация с диаметровой метрикой. Эта оценка живёт в тех же «физических» единицах, что и диаметр DD, так что они совместимы.

4) Лёгкая локальная догладка (BB-шаг, 5–10 итераций)

Сейчас у вас два полушага с убыванием eta. Этого мало в многомере. Добавьте компактную «BB-догладку» только когда точка «обещающая» (например, f < 1.1*bestF) и/или когда диаметр выделенного интервала мал.

Пример замены вашего мини-блока:

cpp
auto project_box = [&](float* q){
    for (int i = 0; i < n; ++i) q[i] = fminf( 3.14159265359f, fmaxf(-3.14159265359f, q[i]));
    if (cost.variableLen) {
        for (int i = 0; i < n; ++i) q[n+i] = fminf(2.0f, fmaxf(0.1f, q[n+i]));
    }
};

auto grad_at = [&](const float* qsrc, float* __restrict g)->void{
    // используем уже посчитанные phi, s_arr, c_arr как у вас;
    // аккуратно пересчитать для qsrc при необходимости
    float xg, yg;
    float ftmp = cost(qsrc, xg, yg);
    // заполняем g аналогично коду выше (см. пункт 3)
    // ...
};

if (f < bestF * 1.1f) {
    alignas(64) float gk[32], qk[32], qprev[32], gprev[32];
    memcpy(qk, q_local, dim * sizeof(float));
    grad_at(qk, gk);
    float alpha = 0.1f;

    for (int k = 0; k < 8; ++k) {                // 8 итераций хватает
        memcpy(qprev, qk, dim * sizeof(float));
        memcpy(gprev, gk, dim * sizeof(float));

        // шаг
        for (int i = 0; i < dim; ++i) qk[i] = qk[i] - alpha * gk[i];
        project_box(qk);

        float x2, y2;
        float f2 = cost(qk, x2, y2);
        if (f2 < f) {
            // BB шаг
            float sTy = 0.0f, sTs = 0.0f;
            for (int i = 0; i < dim; ++i) {
                float si = qk[i] - qprev[i];
                sTs = fmaf(si, si, sTs);
            }
            grad_at(qk, gk);
            for (int i = 0; i < dim; ++i) {
                float yi = gk[i] - gprev[i];
                sTy = fmaf(qk[i] - qprev[i], yi, sTy);
            }
            if (sTy > 1e-12f) alpha = fminf(0.5f, fmaxf(1e-3f, sTs / sTy));
            f = f2; x = x2; y = y2;
        } else {
            alpha *= 0.5f;                       // откат шага
            if (alpha < 1e-6f) break;
        }
    }
    memcpy(q_local, qk, dim * sizeof(float));
}

Почему это быстро. Мы уже в «горячей» точке, вычисляемые массивы phi/s/c у вас есть; 8 итераций почти не влияют на перф, но часто дают релевантное сужение и новое bestF.

5) LHS-сидирование на грубом уровне (2–8 сидов)

Добавьте к generate_heuristic_seeds лёгкий Latin Hypercube Sampling (в нормированных координатах [low,high][low, high]) и сделайте вокруг каждого сида крошечный интервал по tt (±0.001…0.003) как вы уже делаете для ваших сидов.

Скетч (в духе вашего кода, без STL-тяжёлостей):

cpp
static __forceinline int generate_lhs_seeds_lite(const MortonND& map, int dim, float* __restrict S, int stride, uint32_t seed) {
    const int ns = (std::min)(dim + 2, 8); // 2..8 сидов
    uint32_t st = seed;
    for (int s = 0; s < ns; ++s) {
        for (int d = 0; d < dim; ++d) {
            // позиция (i + U)/ns
            st ^= st << 13; st ^= st >> 17; st ^= st << 5;
            float u = (st & 0xFFFFFF) * 5.9604645e-8f; // ~U[0,1)
            float pos = ( (float)s + u ) / (float)ns;
            int pd = map.perm[d];
            float lo = map.low[pd], hi = map.high[pd];
            S[s*stride + d] = fmaf(pos, (hi-lo), lo);
        }
    }
    return ns;
}

Дальше — как с вашими «эвристическими» сидами: перевод map.pointToT, инициализация коротких интервалов в H.

6) Небольшой «пересмотр» mm: пер-масштаб по «числу свободных бит»

Вместо одного глобального Mmax добавьте 8–12 «карманов» M_by_span[k], где k=\sum_d \text{popcnt}((i1^i2)&pextMask[d]) — «уровень» блога. При апдейте Mmax = max(Mmax, M_by_span[k]). Это недорого (один popcnt на измерение), но часто даёт более «честный» m для мелких/крупных блоков.

Что менять минимально в существующей логике

  • Где вы создаёте Interval (все места): теперь вызывайте конструктор без NN и подставляйте метрику через set_metric(diam) (см. п.1).

  • В формуле Shag(...) логика под NN у вас «универсальная». Если хотите полной консистентности с диаметровой метрикой, можно вызывать её с _N=1.0f (Липшицев случай). Это одна строка в местах вызова:

    cpp
    float tNew = Shag(m, x1, x2, y1, y2, /*N=*/1.0f, r);

    Если не хотите трогать — оставьте как есть: характеристика уже станет «ND-осведомлённой», этого обычно достаточно.

  • RecomputeR_* работать не мешаем: они читают N_factor/quadratic_term, мы их уже переопределили на основе диаметра — всё совместимо.

  • MPI-обмен оставляем прежним. На стороне получателя вы тоже создаёте Interval — не забудьте повторить заполнялку i1/i2/diam/set_metric(...), чтобы R считался в той же метрике.

Почему это даст выигрыш (интуитивно)

  • Главное узкое место многомерной редукции — то, что «короткий» по tt интервал может покрывать большой по ∥x∥\|x\| регион при Моронте. Ваша текущая RR этого не видит ⇒ переоценивает привлекательность «рваных» кусочков, ⇒ лишние деления и поздние улучшения bestF.
    Лечение: метрика через ND-диаметр (п.1) делает характеристику адекватной геометрии. Это самый «лёгкий» и самый «бьющий» патч.

  • Gray-код уменьшает число «рваных» соседей по tt ⇒ меньше экстремальных Δf на близких tt ⇒ стабильнее оценка MM и RR.

  • Градиент в MM даёт честный «локальный Липшиц» в единицах вашей физики. Это резко снижает риск «недооценки» mm и раздувания «ложных» RR.

  • Лёгкая догладка вынимает «дешёвые» улучшения рядом с найденными кандидатами (где у вас уже всё вычислено), не превращаясь в тяжёлый локальный оптимизатор.

Что попробовать в первую очередь (план эксперимента)

  1. Включить ND-диаметр в RR (п.1) без других изменений.
    Ожидаемо: быстрее растёт bestF, меньше пустых делений, равномернее распределение MM.

  2. Включить Gray-код (п.2).
    Ожидаемо: меньшая дисперсия Δf на соседних отрезках ⇒ стабильнее сортировка по RR.

  3. Включить градиент-апдейт MM (п.3).
    Ожидаемо: заметное сокращение стадии «разгона» на старте, особенно в больших dim.

  4. Добавить BB-догладку 8 итераций (п.4).
    Ожидаемо: частые «мелкие» улучшения, заметное снижение числа крупных делений при тех же итерациях.

  5. LHS-сидирование 4–8 точек (п.5) — только на грубом уровне.
    Ожидаемо: более стабильный старт на сложных целевых.

  6. (Опционально) M-карманы по «числу свободных бит» (п.6).

Все пункты совместимы и включаются независимо — удобно «аблируются» флагами.

Небольшие корректности/мелочи

  • В pointToT вы округляете к центру ячейки (+0.5f). Это хорошо. Убедитесь, что при Gray-коде применяется gray_decode перед делением на scale (как в патче).

  • В map01ToPoint вы «клипуете» t в [0, 1-ε]. Для консистентности с целочисленным индексом лучше всегда делать:
    idx = min( (uint64_t)((double)t * scale), scale-1 ).

  • В agp_run_branch_mpi вы иногда формируете интервалы «почти нулевой длины». Если диаметр по п.1 даст D=0, то set_metric обеспечит корректность (M=∞ не будет, потому что мы клипуем и nfree==0 ⇒ range=0). На всякий случай можно поставить float D = fmaxf(map.block_diameter(...), 1e-12f);.

Если хотите, я могу собрать эти кусочки в единый «патч» под условные #define AGP_USE_DIAMETER и #define AGP_USE_GRAY_ORDER, чтобы вы могли по одному включать/выключать и смотреть профили. Но даже без глобальной переделки кода эти вставки интегрируются буквально в 3–5 мест и обычно дают самый большой прирост по «скорость+сходимость» в ND-редукции по Моронту.

Share this Q&A