///предыдущий код такой же как будет во второй вер...
생성일: 2025년 11월 7일
답변에 사용된 모델 GPT-5 Thinking by Chat01
질문
///предыдущий код такой же как будет во второй версии
__declspec(align(64)) struct Interval final {
const float x1, x2, y1, y2, delta_y, ordinate_factor;
float N_factor, quadratic_term, M; // Убрали const
float R;
uint64_t i1, i2; // Добавили индексы Мортона
float diam; // Добавили диаметр блока
text__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(0), quadratic_term(0), M(0), R(0), i1(0), i2(0), diam(0) { // Старая метрика для обратной совместимости if (_N == 1.0f) { N_factor = _x2 - _x1; } else if (_N == 2.0f) { N_factor = sqrtf(_x2 - _x1); } else { N_factor = 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 set_metric(float d_alpha) noexcept { 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)); }
};
__forceinline bool ComparePtr(const Interval* a, const Interval* b) noexcept {
return a->R < b->R;
}
///AVX2 функции и функция шага такие же как будут во второй версии кода
struct MortonCachePerRank {
std::vector<int> permCache;
std::vector<uint64_t> invMaskCache;
uint32_t baseSeed;
};
static MortonCachePerRank g_mc;
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;
}
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;
bool use_gray; // Добавили флаг Gray-кодирования
text__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)), use_gray(true) { // Включили Gray по умолчанию 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 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); } __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; // Применяем Gray-кодирование если включено if (use_gray) idx = gray_encode(idx); 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]); } // Декодируем Gray если включено if (use_gray) idx = gray_decode(idx); 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;
float Mmax = M_prior;
int no_improve = 0;
text// Функция преобразования t в индекс Мортона (без Gray для вычисления диаметра) 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; return idx; }; 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; // ДОБАВЛЕНО: Оценка градиента для обновления Mmax 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 grad = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen; grad2 = fmaf(grad, grad, grad2); } if (cost.variableLen) { for (int i = 0; i < n; ++i) { float grad = (dx * c_arr[i] + dy * s_arr[i]) / dist; grad2 = fmaf(grad, grad, grad2); } } float Lcand = sqrtf(grad2) + 1e-8f; Mmax = (std::max)(Mmax, 0.9f * Mmax + 0.1f * Lcand); // КОНЕЦ ДОБАВЛЕНИЯ 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); 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2); I->diam = map.block_diameter(I->i1, I->i2); I->set_metric(I->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t); I->diam = map.block_diameter(I->i1, I->i2); I->set_metric(I->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b); tail->diam = map.block_diameter(tail->i1, tail->i2); tail->set_metric(tail->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика J->i1 = t_to_idx(t1); J->i2 = t_to_idx(t2); J->diam = map.block_diameter(J->i1, J->i2); J->set_metric(J->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f); J->diam = map.block_diameter(J->i1, J->i2); J->set_metric(J->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика для новых интервалов L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew); L->diam = map.block_diameter(L->i1, L->i2); L->set_metric(L->diam); // α = 1 Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2); Rv->set_metric(Rv->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->set_metric(inj->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->set_metric(inj->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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); // ДОБАВЛЕНО: ND-осведомлённая характеристика I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2); I->diam = map.block_diameter(I->i1, I->i2); I->set_metric(I->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ 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);
text// Функция преобразования t в индекс Мортона для fine уровня auto t_to_idx_fine = [&](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)map1.scale); if (idx >= map1.scale) idx = map1.scale - 1ull; return idx; }; 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); Interval* I = new Interval(C->x1, C->x2, f1, f2, (float)dim); // ДОБАВЛЕНО: ND-осведомлённая характеристика для fine уровня I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2); I->diam = map1.block_diameter(I->i1, I->i2); I->set_metric(I->diam); // α = 1 // КОНЕЦ ДОБАВЛЕНИЯ H_fine.push_back(I); 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); }
}
///тестовые функции, 1D, 2D алгоритма и.т.д как во второй версии кода - код с началом внедрения твоих советов "1) ND-осведомлённая характеристика
R
R: заменяем 1D-длину на ND-диаметр блока Моронта
Идея. Интервал по
t
t соответствует блоку индексoв Моронта. Для этого блока можно быстро (чисто побитово) верхне ограничить пер-координатный размах и, значит, диаметр гиперпрямоугольника в исходном пространстве. Подставим этот диаметр
D
D вместо
(
x
2
−
x
1
)
1
/
N
(x
2
−x
1
)
1/N
в формулу
R
R и в оценку
M
M.
Быстрый подсчёт диаметра. Пусть для подотрезка индексы
i
1
,
i
2
i
1
,i
2
(в Моронте). Для каждого измерения
d
d возьмём ваш pextMask[d] и посчитаем, сколько разрядов этой координаты различаются между
i
1
i
1
и
i
2
i
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);
Почему это важно. Для любой Липшицевой
f
f в
R
d
R
d
:
∣
f
(
x
)
−
f
(
y
)
∣
≤
L
∥
x
−
y
∥
∣f(x)−f(y)∣≤L∥x−y∥. Отсюда для вашего 1D отрезка по
t
t, соответствующего ND блоку, гарантированная оценка идёт через ND-диаметр блока. Даже с «рваными» соседями Моронта ваша характеристика становится корректной и стабильной, потому что «размер области» учтён по геометрии задачи, а не по формальной длине
t
t.
Гарантии. Метод Стронгина с
m
≥
L
m≥L и такой характеристикой сохраняет нижние оценки и те же условия сходимости. Вы лишь меняете «мерку»
d
α
d
α
на
D
D, что для Липшица — естественно.
- Gray-кодирование индекса Моронта (не меняя Моронта)
Это дёшево и сильно сглаживает последовательность по
t
t (соседние индексы отличаются в 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; (или флагом).
Почему это стоит делать. Это почти даром, а «локальность по
t
t» заметно улучшается даже для Моронта. В сочетании с пунктом 1) — особенно эффективно.
- Градиент-драйв обновления
M
M (почти бесплатно)
Вы уже считаете суммарные синусы/косинусы для манипулятора и делаете 1–2 локальных шага. Из этого блока можно почти даром получить верхнюю оценку локальной константы Липшица и агрессивнее (но безопасно) поднимать
M
max
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);
Почему это корректно. Для Липшицевой
f
f
L
≥
∥
∇
f
(
x
)
∥
L≥∥∇f(x)∥ почти всюду (гладкость у вас есть). Значит Lcand — валидная «локальная» нижняя оценка для глобального
L
L. Включая её в Mmax, вы ускоряете зажатие без потери гарантии (мы берём максимум).
Синхронизация с диаметровой метрикой. Эта оценка живёт в тех же «физических» единицах, что и диаметр
D
D, так что они совместимы.
- Лёгкая локальная догладка (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;
textfor (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.
- LHS-сидирование на грубом уровне (2–8 сидов)
Добавьте к generate_heuristic_seeds лёгкий Latin Hypercube Sampling (в нормированных координатах
[
l
o
w
,
h
i
g
h
]
[low,high]) и сделайте вокруг каждого сида крошечный интервал по
t
t (±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.
- Небольшой «пересмотр»
m
m: пер-масштаб по «числу свободных бит»
Вместо одного глобального 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 для мелких/крупных блоков.", код после внедрения советов: "#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);
// ==================== СТРУКТУРА ДЛЯ 1D/2D СЛУЧАЕВ ====================
__declspec(align(64)) struct Interval1D final {
const float x1, x2, y1, y2, delta_y, ordinate_factor, N_factor, quadratic_term, M;
float R;
text__forceinline void* operator new(size_t) noexcept { Slab* s = tls.local(); char* r = s->current; s->current += 64u; return r; } __forceinline Interval1D(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 : sqrtf(_x2 - _x1)), 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)); }
};
struct MortonCachePerRank {
std::vector<int> permCache;
std::vector<uint64_t> invMaskCache;
uint32_t baseSeed;
};
static MortonCachePerRank g_mc;
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;
}
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;
bool use_gray;
text__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)), use_gray(true) { 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 float block_diameter(uint64_t i1, uint64_t i2) const noexcept { if (i1 > i2) std::swap(i1, i2); float s2 = 0.0f; 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); const float range = step[pd] * ((nfree ? ((1ull << nfree) - 1ull) : 0ull)); s2 += range * range; } return sqrtf(s2); } __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; if (use_gray) idx = gray_encode(idx); 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]); } if (use_gray) idx = gray_decode(idx); return ((float)idx + 0.5f) / (float)scale; }
};
// ==================== СТРУКТУРА ДЛЯ ND СЛУЧАЕВ ====================
__declspec(align(64)) struct IntervalND final {
const float x1, x2, y1, y2, delta_y, ordinate_factor;
float N_factor, quadratic_term, M;
float R;
uint64_t i1, i2;
float diam;
int span_level; // NEW: уровень блога
text__forceinline void* operator new(size_t) noexcept { Slab* s = tls.local(); char* r = s->current; s->current += 64u; return r; } __forceinline IntervalND(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), span_level(0) { } // NEW: вычисляем уровень блога __forceinline void compute_span_level(const MortonND& map) noexcept { span_level = 0; for (int d = 0; d < map.dim; ++d) { uint64_t varying = (i1 ^ i2) & map.pextMask[d]; span_level += _mm_popcnt_u64(varying); } span_level = (std::min)(span_level, 11); } __forceinline void set_metric(float d_alpha) noexcept { 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)); }
};
// Функции сравнения
__forceinline bool ComparePtr1D(const Interval1D* a, const Interval1D* b) noexcept {
return a->R < b->R;
}
__forceinline bool ComparePtrND(const IntervalND* a, const IntervalND* b) noexcept {
return a->R < b->R;
}
// AVX2 функции для Interval1D
__forceinline void RecomputeR_ConstM_AVX2_1D(Interval1D* 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 Interval1D* 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_1D(Interval1D* 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 Interval1D* 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 Interval1D* 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));
}
}
// AVX2 функции для IntervalND
__forceinline void RecomputeR_ConstM_AVX2_ND(IntervalND* 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 IntervalND* 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_ND(IntervalND* 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 IntervalND* 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 IntervalND* 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 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_lhs_seeds_lite(const MortonND& map, const int dim, float* __restrict S, int stride, uint32_t seed) {
int temp_dim = dim;
const int ns = --temp_dim * temp_dim;
uint32_t st = seed;
text// Генерируем случайные перестановки для каждого измерения (настоящий LHS) alignas(32) int permutations[32][16]; // Максимум 32 измерения × 8 сидов for (int d = 0; d < dim; ++d) { // Инициализация последовательности страт for (int s = 0; s < ns; ++s) { permutations[d][s] = s; } // Перемешивание Фишера-Йейтса для каждого измерения for (int s = ns - 1; s > 0; --s) { st ^= st << 13; st ^= st >> 17; st ^= st << 5; int j = st % (s + 1); std::swap(permutations[d][s], permutations[d][j]); } } // Настоящий LHS с перестановками между измерениями for (int s = 0; s < ns; ++s) { for (int d = 0; d < dim; ++d) { st ^= st << 13; st ^= st >> 17; st ^= st << 5; float u = (st & 0xFFFFFF) * 5.9604645e-8f; // Настоящий LHS: разная перестановка для каждого измерения int stratum = permutations[d][s]; float pos = ((float)stratum + 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;
}
static __forceinline int generate_heuristic_seeds(const ManipCost& cost, const MortonND& map, int dim, float* __restrict S, int stride, uint32_t seed) {
const int n = cost.n;
const bool VL = cost.variableLen;
const float tx = cost.targetX, ty = cost.targetY;
text// Существующие эвристические сиды (3 штуки) int total_seeds = 0; // Сид 0: равномерное распределение углов float* s0 = S + total_seeds * stride; float phi = atan2f(ty, tx); float rho = sqrtf(fmaf(tx, tx, ty * ty)); float len = fminf(fmaxf(rho / (float)n, 0.1f), 2.0f); 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; total_seeds++; // Сид 1: альтернирующие углы float* s1 = S + total_seeds * 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); } total_seeds++; // Сид 2: линейно убывающие углы float* s2 = S + total_seeds * 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); } } total_seeds++; // Добавляем LHS-сиды (2-8 штук) int lhs_count = generate_lhs_seeds_lite(map, dim, S + total_seeds * stride, stride, seed); total_seeds += lhs_count; return total_seeds;
}
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<IntervalND*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY,
float M_prior = 1e-3f) {
textconst int n = cost.n, dim = n + (cost.variableLen ? n : 0); // NEW: карманы для более точного вычисления Mmax alignas(32) float M_by_span[12]; for (int i = 0; i < 12; ++i) M_by_span[i] = M_prior; float Mmax = M_prior; 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 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; return idx; }; float a = 0.0f, b = 1.0f; float dmax = b - a, initial_len = dmax, thr03 = 0.3f * initial_len, inv_thr03 = 1.0f / thr03; int it = 0; auto evalAt = [&](float t) -> float { map.map01ToPoint(t, q_local); float f = cost(q_local, x, y); bool should_optimize = (f < bestF * 1.25f) || ((thr03 > dmax && !(it % 3)) || (10.0f * dmax < initial_len)); if (should_optimize) { 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) { const float* th = qsrc; const float* L = cost.variableLen ? (qsrc + n) : nullptr; __declspec(align(64)) float phi_loc[32], s_arr_loc[32], c_arr_loc[32]; float x_loc = 0.0f, y_loc = 0.0f, phi_acc = 0.0f; for (int i = 0; i < n; ++i) { phi_acc += th[i]; phi_loc[i] = phi_acc; } FABE13_SINCOS(phi_loc, s_arr_loc, c_arr_loc, n); if (cost.variableLen) { for (int i = 0; i < n; ++i) { float Li = L[i]; x_loc = fmaf(Li, c_arr_loc[i], x_loc); y_loc = fmaf(Li, s_arr_loc[i], y_loc); } } else { for (int i = 0; i < n; ++i) { x_loc = c_arr_loc[i] + x_loc; y_loc = s_arr_loc[i] + y_loc; } } alignas(64) float sum_s_loc[32], sum_c_loc[32]; float as = 0.0f, ac = 0.0f; for (int k = n - 1; k >= 0; --k) { const float Lk = cost.variableLen ? qsrc[n + k] : 1.0f; as += Lk * s_arr_loc[k]; ac += Lk * c_arr_loc[k]; sum_s_loc[k] = as; sum_c_loc[k] = ac; } float dx_loc = x_loc - cost.targetX; float dy_loc = y_loc - cost.targetY; float dist_loc = sqrtf(fmaf(dx_loc, dx_loc, dy_loc * dy_loc)); for (int i = 0; i < n; ++i) { float gpen = 0.0f; float ai = fabsf(qsrc[i]); float v = cost.minTheta - ai; if (v > 0.0f) gpen = -2.0f * copysignf(v, qsrc[i]); g[i] = (dx_loc * (-sum_s_loc[i]) + dy_loc * (sum_c_loc[i])) / dist_loc + gpen; } if (cost.variableLen) { for (int i = 0; i < n; ++i) { g[n + i] = (dx_loc * c_arr_loc[i] + dy_loc * s_arr_loc[i]) / dist_loc; } } }; alignas(64) float gk[32], qk[32], qprev[32], gprev[32]; memcpy(qk, q_local, dim * sizeof(float)); grad_at(qk, gk); float grad2 = 0.0f; for (int i = 0; i < dim; ++i) { grad2 = fmaf(gk[i], gk[i], grad2); } float Lcand = sqrtf(grad2); Mmax = (std::max)(Mmax, 0.7f * Mmax + 0.3f * Lcand); float alpha = 0.125f; for (int k = 0; k < 3; ++k) { 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) { float sTy = 0.0f, sTs = 0.0f; for (int i = 0; i < dim; ++i) { float si = qk[i] - qprev[i]; float yi = gk[i] - gprev[i]; sTs = fmaf(si, si, sTs); sTy = fmaf(si, yi, sTy); } alpha = fminf(0.3f, sTs / sTy); grad_at(qk, gk); f = f2; x = x2; y = y2; } else { alpha *= 0.5f; memcpy(qk, qprev, dim * sizeof(float)); memcpy(gk, gprev, dim * sizeof(float)); } } memcpy(q_local, qk, dim * sizeof(float)); } 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 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); const int rank = g_world->rank(); const int world = g_world->size(); alignas(64) float seeds[11 * 32]; const int seedCnt = generate_heuristic_seeds(cost, map, dim, seeds, 32, seed + rank * 7919u); // Инициализация из сидов for (int i = 0; i < seedCnt; ++i) { const float* s = seeds + i * 32; float t_seed = map.pointToT(s); float interval_size; if (i < 3) { interval_size = 0.0004f * (float)dim; } else { float start_val = 0.00031f * (float)dim; float end_val = 0.00025f * (float)dim; interval_size = start_val * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_val / start_val) * (float)(i - 3)); } float t1 = fmaxf(a, t_seed - interval_size); float t2 = fminf(b, t_seed + interval_size); 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); IntervalND* I = new IntervalND(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->compute_span_level(map); // NEW I->set_metric(I->diam); // NEW: обновляем Mmax через карманы M_by_span[I->span_level] = (std::max)(M_by_span[I->span_level], I->M); Mmax = (std::max)(Mmax, I->M); // Ключевая строка: Mmax = max(Mmax, M_by_span[k]) I->ChangeCharacteristic(r * Mmax); // Используем глобальный Mmax float multiplier; if (i < 3) { multiplier = fmaf(0.01f, (float)dim, 0.85f); } else { float start_mult = 0.214f * (float)dim; float end_mult = 0.174f * (float)dim; multiplier = start_mult * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_mult / start_mult) * (float)(i - 3)); } I->R *= multiplier; H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); 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); IntervalND* I = new IntervalND(prev_t, t, prev_f, f); I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t); I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); // NEW I->set_metric(I->diam); // NEW: обновляем Mmax через карманы M_by_span[I->span_level] = (std::max)(M_by_span[I->span_level], I->M); Mmax = (std::max)(Mmax, I->M); // Mmax = max(Mmax, M_by_span[k]) I->ChangeCharacteristic(r * Mmax); // Используем глобальный Mmax H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); prev_t = t; prev_f = f; } // Хвостовой интервал IntervalND* tail = new IntervalND(prev_t, b, prev_f, f_b); tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b); tail->diam = map.block_diameter(tail->i1, tail->i2); tail->compute_span_level(map); // NEW tail->set_metric(tail->diam); // NEW: обновляем Mmax через карманы M_by_span[tail->span_level] = (std::max)(M_by_span[tail->span_level], tail->M); Mmax = (std::max)(Mmax, tail->M); // Mmax = max(Mmax, M_by_span[k]) tail->ChangeCharacteristic(r * Mmax); // Используем глобальный Mmax H.emplace_back(tail); std::push_heap(H.begin(), H.end(), ComparePtrND); // Основной цикл 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); IntervalND* J = new IntervalND(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed); J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f); J->diam = map.block_diameter(J->i1, J->i2); J->compute_span_level(map); // NEW J->set_metric(J->diam); // NEW: обновляем Mmax через карманы M_by_span[J->span_level] = (std::max)(M_by_span[J->span_level], J->M); Mmax = (std::max)(Mmax, J->M); // Mmax = max(Mmax, M_by_span[k]) J->ChangeCharacteristic(r * Mmax); // Используем глобальный Mmax J->R *= 0.9f; H.emplace_back(J); std::push_heap(H.begin(), H.end(), ComparePtrND); } 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(), ComparePtrND); IntervalND* 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; // Используем глобальный Mmax (как и было) float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r); tNew = fminf(fmaxf(tNew, a), b); float fNew = evalAt(tNew); IntervalND* L = new IntervalND(x1, tNew, y1, fNew); IntervalND* Rv = new IntervalND(tNew, x2, fNew, y2); L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew); L->diam = map.block_diameter(L->i1, L->i2); L->compute_span_level(map); // NEW L->set_metric(L->diam); Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2); Rv->compute_span_level(map); // NEW Rv->set_metric(Rv->diam); // NEW: обновляем Mmax через карманы M_by_span[L->span_level] = (std::max)(M_by_span[L->span_level], L->M); M_by_span[Rv->span_level] = (std::max)(M_by_span[Rv->span_level], Rv->M); Mmax = (std::max)(Mmax, (std::max)(L->M, Rv->M)); // Mmax = max(Mmax, M_by_span[k]) 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_ND(H.data(), sz, GF, alpha); std::make_heap(H.begin(), H.end(), ComparePtrND); } 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_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); } } 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_ND(H.data(), sz, m2); std::make_heap(H.begin(), H.end(), ComparePtrND); } } else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); } } H.push_back(L); std::push_heap(H.begin(), H.end(), ComparePtrND); H.push_back(Rv); std::push_heap(H.begin(), H.end(), ComparePtrND); 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); } IntervalND* 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; IntervalND* t1 = H[0]; IntervalND* t2 = (H.size() > 1 ? H[1] : H[0]); IntervalND* t3 = (H.size() > 2 ? H[2] : H[H.size() - 1]); IntervalND* tops[3] = { t1, t2, t3 }; for (uint8_t i2 = 0; i2 < cnt; ++i2) { IntervalND* 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); IntervalND* inj = new IntervalND(sx, ex, y1, y2); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); // NEW inj->set_metric(inj->diam); // NEW: обновляем Mmax через карманы M_by_span[inj->span_level] = (std::max)(M_by_span[inj->span_level], inj->M); Mmax = (std::max)(Mmax, inj->M); inj->ChangeCharacteristic(r * Mmax); // Используем обновленный 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); } IntervalND* 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(), ComparePtrND); } } } 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); IntervalND* inj = new IntervalND(sx, ex, y1, y2); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); // NEW inj->set_metric(inj->diam); // NEW: обновляем Mmax через карманы M_by_span[inj->span_level] = (std::max)(M_by_span[inj->span_level], inj->M); Mmax = (std::max)(Mmax, inj->M); inj->ChangeCharacteristic(r * Mmax); // Используем обновленный 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); } IntervalND* 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(), ComparePtrND); } } } } 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); IntervalND* I = new IntervalND(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->compute_span_level(map); // NEW I->set_metric(I->diam); // NEW: обновляем Mmax через карманы M_by_span[I->span_level] = (std::max)(M_by_span[I->span_level], I->M); Mmax = (std::max)(Mmax, I->M); I->ChangeCharacteristic(r * Mmax); // Используем обновленный Mmax I->R *= 0.90f; H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); } 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) {
textconst 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, high; low.reserve(dim); 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; // Coarse level const int levels0 = (std::min)(peanoLevels, 8), maxIter0 = (int)(maxIterPerBranch * 0.2f); MortonND map0(dim, levels0, low.data(), high.data(), g_mc); std::vector<IntervalND*> 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; } // Fine level (если нужно) if (levels0 < peanoLevels) { MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc); std::vector<IntervalND*> 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 IntervalND* a, const IntervalND* b) { return a->R < b->R; }); const size_t topCount = (size_t)(H_coarse.size() * 0.3f); auto t_to_idx_fine = [&](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)map1.scale); if (idx >= map1.scale) idx = map1.scale - 1ull; return idx; }; // Переносим лучшие интервалы на fine уровень for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) { const IntervalND* 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); IntervalND* I = new IntervalND(C->x1, C->x2, f1, f2); I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2); I->diam = map1.block_diameter(I->i1, I->i2); I->set_metric(I->diam); H_fine.push_back(I); 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(), ComparePtrND); } // Запускаем fine уровень с улучшенной версией 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); }
}" - сейчас минимум ищется стабильнее чем до внедрения советов - но всё ещё присутствуёт лёгкая нестабильность - посмотри на код "издалека" - проверь правильно ли концептуально внедрены советы - действительно ли они внедрены так как задумывались? особенно внимательно проверь правильность внедрения последнего совета с "карманами", отрывок разговора с deepseek: ""Небольшой «пересмотр»
m
m: пер-масштаб по «числу свободных бит»
Вместо одного глобального 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 для мелких/крупных блоков." совет chat01ai применительно к коду: "#include "pch.h"
... далее привожу свой код" как внедрить - что это даст на практике и объяснение на пальцах суть предлагаемого изменения - внеси конкретные изменения в код для теста, это мой алгоритм глобальной оптимизации для нахождения оптимального положения захвата робота манипулятора - обрати внимание я пишу высокооптимизированный скоростной код - твои правки должны быть в том же стиле - главное скорость при математической гарантиии сходимости
Мы собираемся внедрить предложенное изменение: вместо одного глобального Mmax использовать 8-12 "карманов" M_by_span[k], где k = сумма по d popcnt((i1 XOR i2) & pextMask[d]) - уровень блога.
Идея в том, чтобы для каждого интервала вычислять "уровень" (количество свободных бит в разности индексов i1 и i2 по всем измерениям) и хранить отдельные Mmax для каждого уровня.
Это позволит более точно оценивать константу Липшица для блоков разного размера, особенно для мелких и крупных блоков.
Полная реализация для agp_run_branch_mpi
cpp
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<IntervalND*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY,
float M_prior = 1e-3f) {
textconst int n = cost.n, dim = n + (cost.variableLen ? n : 0); // NEW: карманы для Mmax по уровням (12 карманов: 0-11) alignas(32) float M_by_span[12]; for (int i = 0; i < 12; ++i) M_by_span[i] = M_prior; float Mmax = M_prior; 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 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; return idx; }; float a = 0.0f, b = 1.0f; float dmax = b - a, initial_len = dmax, thr03 = 0.3f * initial_len, inv_thr03 = 1.0f / thr03; int it = 0; // ... evalAt функция остается без изменений ... 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); const int rank = g_world->rank(); const int world = g_world->size(); alignas(64) float seeds[11 * 32]; const int seedCnt = generate_heuristic_seeds(cost, map, dim, seeds, 32, seed + rank * 7919u); // Инициализация из сидов с использованием карманов for (int i = 0; i < seedCnt; ++i) { const float* s = seeds + i * 32; float t_seed = map.pointToT(s); float interval_size; if (i < 3) { interval_size = 0.0004f * (float)dim; } else { float start_val = 0.00031f * (float)dim; float end_val = 0.00025f * (float)dim; interval_size = start_val * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_val / start_val) * (float)(i - 3)); } float t1 = fmaxf(a, t_seed - interval_size); float t2 = fminf(b, t_seed + interval_size); 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); IntervalND* I = new IntervalND(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->compute_span_level(map); // NEW I->set_metric(I->diam); // NEW: обновляем соответствующий карман M_by_span[I->span_level] = (std::max)(M_by_span[I->span_level], I->M); Mmax = (std::max)(Mmax, I->M); I->ChangeCharacteristic(r * M_by_span[I->span_level]); // NEW: используем карманный Mmax float multiplier; if (i < 3) { multiplier = fmaf(0.01f, (float)dim, 0.85f); } else { float start_mult = 0.214f * (float)dim; float end_mult = 0.174f * (float)dim; multiplier = start_mult * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_mult / start_mult) * (float)(i - 3)); } I->R *= multiplier; H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); 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); IntervalND* I = new IntervalND(prev_t, t, prev_f, f); I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t); I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); // NEW I->set_metric(I->diam); // NEW: обновляем карман M_by_span[I->span_level] = (std::max)(M_by_span[I->span_level], I->M); Mmax = (std::max)(Mmax, I->M); I->ChangeCharacteristic(r * M_by_span[I->span_level]); // NEW H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); prev_t = t; prev_f = f; } // Хвостовой интервал с использованием карманов IntervalND* tail = new IntervalND(prev_t, b, prev_f, f_b); tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b); tail->diam = map.block_diameter(tail->i1, tail->i2); tail->compute_span_level(map); // NEW tail->set_metric(tail->diam); // NEW: обновляем карман M_by_span[tail->span_level] = (std::max)(M_by_span[tail->span_level], tail->M); Mmax = (std::max)(Mmax, tail->M); tail->ChangeCharacteristic(r * M_by_span[tail->span_level]); // NEW H.emplace_back(tail); std::push_heap(H.begin(), H.end(), ComparePtrND); // Основной цикл 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); IntervalND* J = new IntervalND(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed); J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f); J->diam = map.block_diameter(J->i1, J->i2); J->compute_span_level(map); // NEW J->set_metric(J->diam); // NEW: обновляем карман M_by_span[J->span_level] = (std::max)(M_by_span[J->span_level], J->M); Mmax = (std::max)(Mmax, J->M); J->ChangeCharacteristic(r * M_by_span[J->span_level]); // NEW J->R *= 0.9f; H.emplace_back(J); std::push_heap(H.begin(), H.end(), ComparePtrND); } 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(), ComparePtrND); IntervalND* cur = H.back(); H.pop_back(); const float x1 = cur->x1, x2 = cur->x2, y1 = cur->y1, y2 = cur->y2; float m = r_eff * M_by_span[cur->span_level]; // NEW: используем карманный Mmax float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r); tNew = fminf(fmaxf(tNew, a), b); float fNew = evalAt(tNew); IntervalND* L = new IntervalND(x1, tNew, y1, fNew); IntervalND* Rv = new IntervalND(tNew, x2, fNew, y2); L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew); L->diam = map.block_diameter(L->i1, L->i2); L->compute_span_level(map); // NEW L->set_metric(L->diam); Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2); Rv->compute_span_level(map); // NEW Rv->set_metric(Rv->diam); // NEW: обновляем карманы M_by_span[L->span_level] = (std::max)(M_by_span[L->span_level], L->M); M_by_span[Rv->span_level] = (std::max)(M_by_span[Rv->span_level], Rv->M); Mmax = (std::max)(Mmax, (std::max)(L->M, Rv->M)); 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)) { // В адаптивном режиме используем индивидуальные m для каждого интервала 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); // Пересчитываем R для всех интервалов в адаптивном режиме for (auto pI : H) { float global_coeff = fmaf(MULT * M_by_span[pI->span_level], r_eff, -MULT * M_by_span[pI->span_level]); float GF = beta * global_coeff; pI->ChangeCharacteristic(fmaf(GF, (pI->x2 - pI->x1), pI->M * alpha)); } // Для новых интервалов float GF_L = beta * fmaf(MULT * M_by_span[L->span_level], r_eff, -MULT * M_by_span[L->span_level]); float GF_Rv = beta * fmaf(MULT * M_by_span[Rv->span_level], r_eff, -MULT * M_by_span[Rv->span_level]); L->ChangeCharacteristic(fmaf(GF_L, len1, L->M * alpha)); Rv->ChangeCharacteristic(fmaf(GF_Rv, len2, Rv->M * alpha)); std::make_heap(H.begin(), H.end(), ComparePtrND); } else { // Неадаптивная часть в адаптивном режиме - используем карманы L->ChangeCharacteristic(r_eff * M_by_span[L->span_level]); Rv->ChangeCharacteristic(r_eff * M_by_span[Rv->span_level]); if (Mloc > 1.15f * Mmax) { // Глобальный пересчет с использованием карманов for (auto pI : H) { pI->ChangeCharacteristic(r_eff * M_by_span[pI->span_level]); } std::make_heap(H.begin(), H.end(), ComparePtrND); } } } else { // Неадаптивный режим - только карманы L->ChangeCharacteristic(r_eff * M_by_span[L->span_level]); Rv->ChangeCharacteristic(r_eff * M_by_span[Rv->span_level]); if (Mloc > 1.15f * Mmax) { for (auto pI : H) { pI->ChangeCharacteristic(r_eff * M_by_span[pI->span_level]); } std::make_heap(H.begin(), H.end(), ComparePtrND); } } H.push_back(L); std::push_heap(H.begin(), H.end(), ComparePtrND); H.push_back(Rv); std::push_heap(H.begin(), H.end(), ComparePtrND); // ... остальная часть основного цикла без изменений ... // (обмен сообщениями, проверка условий остановки и т.д.) ++it; }
}
Модификация для AGP_Manip2D
cpp
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) {
textconst 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, high; low.reserve(dim); 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; // Coarse level const int levels0 = (std::min)(peanoLevels, 8), maxIter0 = (int)(maxIterPerBranch * 0.2f); MortonND map0(dim, levels0, low.data(), high.data(), g_mc); std::vector<IntervalND*> 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)); // NEW: передаем M_prior в функцию, которая теперь использует карманы 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; } // Fine level (если нужно) if (levels0 < peanoLevels) { MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc); std::vector<IntervalND*> 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 IntervalND* a, const IntervalND* b) { return a->R < b->R; }); const size_t topCount = (size_t)(H_coarse.size() * 0.3f); auto t_to_idx_fine = [&](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)map1.scale); if (idx >= map1.scale) idx = map1.scale - 1ull; return idx; }; // NEW: инициализируем H_fine с правильными span_level для fine уровня for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) { const IntervalND* 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); IntervalND* I = new IntervalND(C->x1, C->x2, f1, f2); I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2); I->diam = map1.block_diameter(I->i1, I->i2); I->compute_span_level(map1); // NEW: пересчитываем для fine уровня I->set_metric(I->diam); H_fine.push_back(I); 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(), ComparePtrND); } // NEW: запускаем fine уровень с карманами 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; } } // ... остальная часть функции без изменений (MPI сбор результатов) ...
}
я думаю ты неправильно понял идею chat01ai ведь там говорилось "При апдейте Mmax = max(Mmax, M_by_span[k])" - то есть Mmax имеет тот же смысл - так что нужно и дальше использовать Mmax - особая путанница в адаптивной схеме там ты делаешь полный пересчёт в местах где у меня например был пересчёт только двух новых интервалов - ты не должен менять существующую логику ты должен лишь добавить эти карманы как метод для более честного пересчёта Mmax, приведи с учётом этих правок полный код без пропусков этих двух функций
Мы должны внедрить идею "карманов" для Mmax, но при этом сохранить существующую логику, особенно в адаптивной схеме.
Идея: вместо одного Mmax, мы имеем массив M_by_span[12] (для уровней 0..11). При создании или обновлении интервала мы обновляем M_by_span[k] для соответствующего уровня k (который вычисляется как сумма popcnt по измерениям) и затем Mmax = max(Mmax, M_by_span[k]).
Таким образом, Mmax остается глобальным максимумом, но при вычислении характеристики интервала мы используем Mmax (а не M_by_span[k]). Однако, в предложении chat01ai говорится использовать M_by_span[k] для вычисления m для интервала.
Но в исходном предложении: "При апдейте Mmax = max(Mmax, M_by_span[k])", то есть Mmax обновляется как максимум по всем карманам, а затем в формуле m = r * Mmax.
Однако, в следующем предложении: "но часто даёт более «честный» m для мелких/крупных блоков" - это может означать, что для каждого интервала используется свой M_by_span[k] вместо глобального Mmax.
Уточнение: в исходном предложении не совсем ясно, что именно используется в формуле m. Но из контекста "честный m" следует, что для интервала уровня k используется m = r * M_by_span[k]." - особенно тщательно проверь метод внедрения карманов - его корректность во второй присланной версии кода, если найдёшь концептуальные несоответствия - приведи конкретные правки в код, также в целом рассмотри многомерный алгоритм для робота манипулятора и предложи конкретные алгоритмические улучшения которые могут ускорить сходимость - также приведи конкретный код, долго думай, тщательно анализируй вторую версию алгоритма изменённую по твоим советам - ищи там несоответствия с логикой советов и логические ошибки - приведи полную исправленную версию кода с улучшенной сходимостью
사고
답변
крутая работа — видно, что ты аккуратно внедрил диаметр Мортона, Gray-код и «локальный» градиент для подкачки . Лёгкая нестабильность, которую ты видишь, действительно могла появиться из-за двух вещей:
-
«карманы» в твоей второй версии стали использоваться вместо глобального (то есть ты подставлял в характеристику интервала ). По исходной идее chat01ai карманы нужны только для более честного апдейта глобального :
,
а дальше везде в формуле характеристики используется один общий . Это сохраняет гарантии Стронгина и не ломает твою адаптивную ветку. -
Во второй версии у тебя сильно выросло количество LHS-сидов (формула
--temp_dim * temp_dim) при фиксированном стекеfloat seeds[11*32]. Это легко приводит к выходу за границы массива и «дрожанию» поведения. В коде ниже я поставил безопасный запас (64 сида на стек), а правильнее — держать LHS ≤ 8 (как и предлагалось изначально).
Ниже — две функции целиком, с минимально-инвазивной интеграцией «карманов» только как механизма апдейта . Логика адаптивного/неадаптивного пересчёта оставлена как в твоём коде; я лишь исправил внутреннее условие с «1.15×» так, чтобы оно сравнивалось с предыдущим (иначе оно всегда ложное сразу после присваивания).
agp_run_branch_mpi (исправленная)
cppstatic __forceinline void agp_run_branch_mpi( const MortonND& map, const ManipCost& cost, int maxIter, float r, bool adaptive, float eps, unsigned seed, std::vector<IntervalND*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY, float M_prior = 1e-3f) { const int n = cost.n; const int dim = n + (cost.variableLen ? n : 0); // ---- NEW: "карманы" по уровню span (0..11), только для апдейта Mmax ---- alignas(32) float M_by_span[12]; for (int i = 0; i < 12; ++i) M_by_span[i] = M_prior; float Mmax = M_prior; // ------------------------------------------------------------------------ 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 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; // ВАЖНО: без Gray для оценки диаметра/уровня return idx; }; auto update_pockets_and_Mmax = [&](IntervalND* I) __forceinline { // апдейтим карман уровня и глобальный Mmax const int k = I->span_level; // 0..11 (см. compute_span_level) if (I->M > M_by_span[k]) M_by_span[k] = I->M; if (M_by_span[k] > Mmax) Mmax = M_by_span[k]; }; float a = 0.0f, b = 1.0f; auto evalAt = [&](float t) -> float { map.map01ToPoint(t, q_local); float f = cost(q_local, x, y); // локальная мини-догладка только когда «обещает» if (f < bestF * 1.25f) { // считаем суммарные sin/cos и градиент (как у тебя) 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 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 g = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen; grad2 = fmaf(g, g, grad2); } if (cost.variableLen) { for (int i = 0; i < n; ++i) { float g = (dx * c_arr[i] + dy * s_arr[i]) / dist; grad2 = fmaf(g, g, grad2); } } // подкачка глобального Mmax (как и раньше) float Lcand = sqrtf(grad2) + 1e-8f; float Mtmp = 0.7f * Mmax + 0.3f * Lcand; if (Mtmp > Mmax) Mmax = Mtmp; // пара «пробных» шагов (твой блок, оставлен по сути без изменений) float eta = 0.125f; for (int step = 0; step < 3; ++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 g = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen; q_try[i] = q_local[i] - eta * g; 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 g = (dx * c_arr[i] + dy * s_arr[i]) / dist; float v = q_local[n + i] - eta * g; 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 f_a = evalAt(a), f_b = evalAt(b); const int K = (std::min)((std::max)(2 * dim, 8), 128); // ВНИМАНИЕ: сохраняем старую логику — чистим H тут H.clear(); H.reserve((size_t)maxIter + K + 16); const int rank = g_world->rank(); const int world = g_world->size(); // безопасный стек под сиды (до 64 сидов по 32 float) alignas(64) float seeds[64 * 32]; const int seedCnt = generate_heuristic_seeds(cost, map, dim, seeds, 32, seed + rank * 7919u); // --- инициализация по сидами (как у тебя, но: m=r*Mmax и с апдейтом «карманов») --- for (int i = 0; i < seedCnt && i < 64; ++i) { const float* s = seeds + i * 32; float t_seed = map.pointToT(s); float interval_size; if (i < 3) interval_size = 0.0004f * (float)dim; else { float start_val = 0.00031f * (float)dim; float end_val = 0.00025f * (float)dim; interval_size = start_val * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_val / start_val) * (float)(i - 3)); } float t1 = fmaxf(a, t_seed - interval_size); float t2 = fminf(b, t_seed + interval_size); if (t2 <= t1) continue; 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); IntervalND* I = new IntervalND(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->compute_span_level(map); I->set_metric(I->diam); update_pockets_and_Mmax(I); I->ChangeCharacteristic(r * Mmax); // ← глобальный Mmax if (i < 3) I->R *= fmaf(0.01f, (float)dim, 0.85f); else { float start_mult = 0.214f * (float)dim; float end_mult = 0.174f * (float)dim; float mult = start_mult * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_mult / start_mult) * (float)(i - 3)); I->R *= mult; } H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); 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); IntervalND* I = new IntervalND(prev_t, t, prev_f, f); I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t); I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); I->set_metric(I->diam); update_pockets_and_Mmax(I); I->ChangeCharacteristic(r * Mmax); // ← глобальный Mmax H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); prev_t = t; prev_f = f; } // --- хвост --- IntervalND* tail = new IntervalND(prev_t, b, prev_f, f_b); tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b); tail->diam = map.block_diameter(tail->i1, tail->i2); tail->compute_span_level(map); tail->set_metric(tail->diam); update_pockets_and_Mmax(tail); tail->ChangeCharacteristic(r * Mmax); // ← глобальный Mmax H.emplace_back(tail); std::push_heap(H.begin(), H.end(), ComparePtrND); // --- рабочие переменные цикла --- 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); IntervalND* J = new IntervalND(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed); J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f); J->diam = map.block_diameter(J->i1, J->i2); J->compute_span_level(map); J->set_metric(J->diam); update_pockets_and_Mmax(J); J->ChangeCharacteristic(r * Mmax); // ← глобальный Mmax J->R *= 0.9f; H.emplace_back(J); std::push_heap(H.begin(), H.end(), ComparePtrND); } 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(), ComparePtrND); IntervalND* 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; // всегда глобальный Mmax float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r); tNew = fminf(fmaxf(tNew, a), b); float fNew = evalAt(tNew); IntervalND* L = new IntervalND(x1, tNew, y1, fNew); IntervalND* Rv = new IntervalND(tNew, x2, fNew, y2); L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew); L->diam = map.block_diameter(L->i1, L->i2); L->compute_span_level(map); L->set_metric(L->diam); Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2); Rv->compute_span_level(map); Rv->set_metric(Rv->diam); float Mloc = (std::max)(L->M, Rv->M); // апдейтим карманы и Mmax ТОЛЬКО для апдейта глобального Mmax update_pockets_and_Mmax(L); update_pockets_and_Mmax(Rv); const float prevMmax = Mmax; if (Mloc > Mmax) Mmax = Mloc; // как и было логически m = r_eff * Mmax; // m всегда через глобальный Mmax // --- твоя адаптивная логика без изменений по структуре --- 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)) { // "Аффинная" подстройка как у тебя const float progress = fmaf(-dmax, inv_thr03, 1.0f); const float alpha = progress * progress; const float beta = fmaf(-alpha, 1.0f, 2.0f); const float MULT = (1.0f / dmax) * Mmax; const float global_coeff = fmaf(MULT, r_eff, -MULT); const float 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_ND(H.data(), sz, GF, alpha); std::make_heap(H.begin(), H.end(), ComparePtrND); } else { if (Mloc > prevMmax) { // NB: проверяем с прошлым Mmax, чтобы "1.15×" работало по назначению L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); if (Mloc > 1.15f * prevMmax) { size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); } } else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); } } } else { if (Mloc > prevMmax) { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); if (Mloc > 1.15f * prevMmax) { size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); } } else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); } } H.push_back(L); std::push_heap(H.begin(), H.end(), ComparePtrND); H.push_back(Rv); std::push_heap(H.begin(), H.end(), ComparePtrND); 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); } IntervalND* 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; IntervalND* t1 = H[0]; IntervalND* t2 = (H.size() > 1 ? H[1] : H[0]); IntervalND* t3 = (H.size() > 2 ? H[2] : H[H.size() - 1]); IntervalND* tops[3] = { t1, t2, t3 }; for (uint8_t i2 = 0; i2 < cnt; ++i2) { IntervalND* 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 y1i = cost(tmp, tx, ty); map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty); IntervalND* inj = new IntervalND(sx, ex, y1i, y2i); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); inj->set_metric(inj->diam); update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax); // ← глобальный 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); } IntervalND* topH = H.front(); if (inj->R > 1.15f * topH->R) { float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f); float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f) : fmaf(0.3491860241215959f, expm1f(p2), 0.6f); inj->R = in.xchg.Rtop * k; H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND); } } } 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 y1i = cost(tmp, tx, ty); map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty); IntervalND* inj = new IntervalND(sx, ex, y1i, y2i); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); inj->set_metric(inj->diam); update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax); // ← глобальный 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); } IntervalND* topH = H.front(); if (inj->R > 1.15f * topH->R) { float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f); float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f) : fmaf(0.3491860241215959f, expm1f(p2), 0.6f); inj->R = d[4] * k; H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND); } } } } 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); IntervalND* I = new IntervalND(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->compute_span_level(map); I->set_metric(I->diam); update_pockets_and_Mmax(I); I->ChangeCharacteristic(r * Mmax); // ← глобальный Mmax I->R *= 0.90f; H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); } if (bm.bestF < bestF) { bestF = bm.bestF; bestX = bm.bestX; bestY = bm.bestY; bestQ.assign(bm.bestQ, bm.bestQ + bm.dim); } } } } ++it; } }
AGP_Manip2D (исправленная)
cppextern "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; // --- coarse level --- const int levels0 = (std::min)(peanoLevels, 8); const int maxIter0 = (int)(maxIterPerBranch * 0.2f); MortonND map0(dim, levels0, low.data(), high.data(), g_mc); std::vector<IntervalND*> 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; } // --- fine level (если нужно) --- if (levels0 < peanoLevels) { MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc); std::vector<IntervalND*> 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)); // NB: как и в твоём коде — предварительное заполнение H_fine здесь if (!H_coarse.empty()) { std::sort(H_coarse.begin(), H_coarse.end(), [](const IntervalND* a, const IntervalND* b) { return a->R < b->R; }); const size_t topCount = (size_t)(H_coarse.size() * 0.3f); auto t_to_idx_fine = [&](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)map1.scale); if (idx >= map1.scale) idx = map1.scale - 1ull; return idx; // без Gray — только для оценки диапазона }; for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) { const IntervalND* 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); IntervalND* I = new IntervalND(C->x1, C->x2, f1, f2); I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2); I->diam = map1.block_diameter(I->i1, I->i2); I->set_metric(I->diam); // span_level для fine не обязателен тут; будет вычислен в run H_fine.push_back(I); 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(), ComparePtrND); } // ВНИМАНИЕ: внутри agp_run_branch_mpi есть H.clear(); // это повторяет твою текущую логику (предварительное H_fine обнуляется). // Если хочешь реально использовать pre-seed из coarse → просто убери H.clear() в run-функции. 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); } }
Что именно поменялось и почему это стабилизирует
-
Карманы теперь используются строго по назначению: они аккумулируют локальные оценки на разных уровнях «свободных бит», но в формуле характеристики всегда стоит один глобальный . Это сохраняет гарантию корректной нижней оценки и общую структуру твоего алгоритма (включая адаптивную ветку с
RecomputeR_*). -
Сравнение с 1.15× теперь делается относительно предыдущего (
prevMmax). Раньше условие послеMmax=Mlocвсегда было ложным. -
Стек под сиды расширен до 64×32 float, чтобы исключить переполнение, если генератор сидов вдруг вернёт больше 11. (Рекомендация: держи LHS ≤ 8; если хочешь — верну прежнюю реализацию LHS на 2..8 сидов.)
-
t_to_idx для диаметра/уровня — без Gray (как ты и задумывал), а сама траектория
map01ToPoint— с Gray. Это корректно: диаметр нужен как грубая верхняя оценка геометрического блока, а Gray служит для «сглаживания» прохода по .
Опционально (одна строка, даёт быстрый профит)
Если ты действительно хотел переносить «топ-интервалы» с coarse на fine, убери одну строку в начале agp_run_branch_mpi:
cpp// H.clear(); // ← просто закомментируй
Тогда предварительно собранный H_fine реально пойдёт в работу (сейчас он обнуляется). Это обычно ускоряет «выход на цель» на fine-уровне.
Если захочешь, могу прислать аккуратный патч на generate_lhs_seeds_lite, чтобы он гарантированно возвращал ≤8 сидов и не ломал выравнивание.