"В работах [23,24,25,26,27] рассматривается задача...
Tạo vào: 26 tháng 10, 2025
Trả lời bằng GPT-5 Thinking bởi Chat01
Câu hỏi
"В работах [23,24,25,26,27] рассматривается задача управления роботами. В качестве примера далее рассматривается робот,построенный на основе многосекционного механизма с параллельной структурой. Данные роботы широко распространены и находят применение в широком круге задач [23,24,25]. Одним из наиболее известных примеров таких механизмов является платформа Гью-Стюарта [26]. Дополнительные примеры манипуляторов представлены в [23,24,25].Автономное управление роботами является сложной задачей. Для управления механизмом последовательно решается целая серия задач –поиск оптимальной стратегии управления механизмом, поиск тактики перемещения манипулятораирасчет параметров необходимых для реализации управления на физическом уровне [26,27,28]. На этапе поиска стратегии необходимо найти конфигурацию манипулятора,обеспечивающуюзаданное положение захвата. Целью этапа поиска тактики является нахождение траектории перехода манипулятора из начальной конфигурации в финальное состояние. На физическом уровне необходимо найти параметры управления приводами манипулятора,необходимые для перехода манипулятора из начального состояния в финальное.Остановимся на более подробномрассмотрениипервогоэтапауправления манипулятором –поиск стратегии управления. Описать положение захвата в пространстве можно набором координат . Примеромкоординат вектора могут быть, например,координата точки
22крепления захвата и набор углов поворота захвата относительно некоторого стандартного положения. Требуется перевести захват в новое положение с заданным набором значений координат вектора . Зачастую при поискестратегии управления манипуляторомдля описания текущего положенияиспользуется вектор обобщенных координат . Компонентывектораобобщенных координатпозволяют описать конфигурацию каждогоузла манипулятора в пространстве. Например, поворот, растяжение или сжатие секций. Пример манипулятора с четырьмя обобщенными координатами представлен нарис. 3. Координаты , и соответствуют углам между секциями механизма в двумерном пространстве. Координата отвечает за растяжение или сжатие секции. Вектор обобщенных координат содержитварьируемыепараметры. При поиске стратегии управления необходимо учитывать физические особенности манипулятора и окружающей среды –обобщенные координаты связаны с конкретными физическими устройствами. Как следствие, обобщенные координаты могут быть заданы только в определенном диапазоне. Набор допустимых значений определяет область поиска . Рис. 3.Пример многосекционного механизма с параллельной структуройЧтобы установить соответствие между вектором обобщенных координат и положением захвата в пространстве решается задача прямой кинематики.(13)Используя взаимосвязь (13),исходную задачу поиска стратегии управления можно свести к задаче многоэкстремальной оптимизации. В качестве критериязахваттребуемое положение захвата Sв одной из стратегий
23поиска оптимального положения манипулятора может выступать близость захвата к требуемому положению(14)где .Стоит отметить, что рассмотренная постановка задачи (14) является достаточноупрощенной относительно исходной постановки задачи принятия оптимальных решений. ЛПР может учитывать дополнительные параметры и критерии из общей математической модели. Для уточнения постановки задачи можноучитывать, например,наличиепрепятствий на пути манипулятора,а такженеобходимость достижения дополнительных целей, таких как минимизация кривизны манипулятора, близость к определенному положению в пространстве, надежность фиксации захвата, снижение энергозатрат и т.д. Минимизация кривизны может потребоваться для обеспечения надежности работы механизма (в крайних положениях механизм может заклинивать). Необходимость близости к определенному положению в пространствеобусловлена тем, что манипулятор может быть запрограммированна выполнение последовательности действий. Минимизацияпроходимогорасстояния позволяет обеспечить экономиюэнергии затрачиваемой на перемещение. Для описания дополнительных целей используютсячастные критерии эффективности относительно обобщенного вектора координат . При наличии препятствий на пути, каждое препятствиеописываетсяв виде функции близости секций манипулятора к препятствию , , где –количество секций [28]. Столкновение манипулятора с препятствием означает, что расстояние сократилось до нуля. В реальных задачах для обеспечения надежности оставляют небольшие запасы по расстояниям , . Учитывая наличие дополнительных целей, а также возможных препятствий,более полная постановка задачи поиска оптимальной стратегии сводится к задаче МКО в постановке (1):
24,.(15)Из представленнойпостановки задачи (15) видно, что при поиске оптимальных решений ЛПР может менять постановку задачи МКО,уточняя, например, допускипри определении ограничений." - Это отрывок диссертации моего научного руководителя, а это мой код АГП расчитанного сейчас на одномерные/двумерные задачи на тестовых функциях: #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 abs_val = fabsf(x);
float reduced = fmodf(abs_val, 6.28318530718f);
if(reduced > 3.14159265359f) {
reduced = 6.28318530718f - reduced;
}
if(reduced < 1.57079632679f) {
const float val2 = reduced * reduced;
const float val4 = val2 * val2;
result_var = fmaf(val4, fmaf(val2, -0.0013888889f, 0.0416666667f), fmaf(val2, -0.5f, 1.0f));
} else {
reduced = 3.14159265359f - reduced;
const float val2 = reduced * reduced;
const float val4 = val2 * val2;
result_var = -fmaf(val4, fmaf(val2, -0.0013888889f, 0.0416666667f), fmaf(val2, -0.5f, 1.0f));
}
} 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* __restrict g_env;
static const boost::mpi::communicator* __restrict 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 __restrict base;
char* __restrict current;
char* const __restrict end;
text__declspec(noalias) __forceinline Slab(void* const __restrict memory, const size_t usable_size) noexcept : base(static_cast<char* __restrict>(memory)) , current(base) , end(base + (usable_size & ~static_cast<size_t>(63u))) { }
};
static tbb::enumerable_thread_specific<Slab*> tls( noexcept {
void* const __restrict memory = _aligned_malloc(16777216u, 16u);
Slab* const __restrict slab = static_cast<Slab*>(_aligned_malloc(32u, 16u));
new (slab) Slab(memory, 16777216u);
char* __restrict p = slab->base;
#pragma loop(ivdep)
while (p < slab->end) {
*p = 0u;
p += 4096u;
}
return slab;
}());
__declspec(align(64)) struct Interval final {
const float x1;
const float x2;
const float y1;
const float y2;
const float delta_y;
const float ordinate_factor;
const float N_factor;
const float quadratic_term;
const float M;
float R;
text__declspec(noalias) __forceinline void* operator new(const size_t) noexcept { Slab* const __restrict s = tls.local(); char* const __restrict result = s->current; s->current += 64u; return result; } __declspec(noalias) __forceinline Interval(const float _x1, const float _x2, const float _y1, const float _y2, const float _N) noexcept : x1(_x1), x2(_x2), y1(_y1), y2(_y2) , delta_y(_y2 - _y1) , ordinate_factor(-(y1 + y2) * 2.0f) , N_factor(_N == 1.0f ? _x2 - _x1 : _N == 2.0f ? sqrtf(_x2 - _x1) : powf(_x2 - _x1, 1.0f / _N)) , quadratic_term((1.0f / N_factor)* delta_y* delta_y) , M((1.0f / N_factor)* fabsf(delta_y)) { } __declspec(noalias) __forceinline void ChangeCharacteristic(const float _m) noexcept { R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor)); }
};
__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;
text__declspec(noalias) __forceinline Peano2DMap( const int L, const float _a, const float _b, const float _c, const float _d, const uint8_t startType ) noexcept : levels(L) , a(_a), b(_b), c(_c), d(_d) , lenx(_b - _a) , leny(_d - _c) , inv_lenx(1.0f / (_b - _a)) , scale(static_cast<uint32_t>(1u) << (L << 1)) , start(startType) { }
};
static Peano2DMap gActiveMap(0, 0.0f, 0.0f, 0.0f, 0.0f, 0b00u);
static __declspec(noalias) __forceinline bool ComparePtr(const Interval* const __restrict a, const Interval* const __restrict b) noexcept
{
return a->R < b->R;
}
static __declspec(noalias) __forceinline float ShekelFunc(const float x, const float seed) noexcept
{
int i = 0;
float current_state = seed, current_res, current_res2, res = 0.0f;
while (i < 10) {
XOR_RAND(current_state, current_res);
const float x_part = fmaf(-current_res, 10.0f, x);
XOR_RAND(current_state, current_res);
XOR_RAND(current_state, current_res2);
float delimiter = fmaf(fmaf(current_res, 20.0f, 5.0f), x_part * x_part, fmaf(current_res2, 0.2f, 1.0f));
delimiter = copysignf(fmaxf(fabsf(delimiter), FLT_MIN), delimiter);
res -= 1.0f / delimiter;
++i;
}
return res;
}
static __declspec(noalias) __forceinline float RastriginFunc(const float x1, const float x2) noexcept
{
const float term1 = fmaf(x1, x1, x2 * x2);
float cos1, cos2;
FABE13_COS(6.28318530717958647692f * x1, cos1);
FABE13_COS(6.28318530717958647692f * x2, cos2);
return (term1 - fmaf(cos1 + cos2, 10.0f, -14.6f)) * fmaf(-term1, 0.25f, 18.42f);
}
static __declspec(noalias) __forceinline float HillFunc(const float x, const float seed) noexcept
{
int j = 0;
__declspec(align(32)) float angles[14u];
const float start_angle = 6.28318530717958647692f * x;
#pragma loop(ivdep)
while (j < 14) {
angles[j] = start_angle * static_cast<float>(j + 1);
++j;
}
__declspec(align(32)) float sin_vals[14u];
__declspec(align(32)) float cos_vals[14u];
FABE13_SINCOS(angles, sin_vals, cos_vals, 14u);
float current_state = seed, current_res, current_res2;
XOR_RAND(current_state, current_res);
float res = fmaf(current_res, 2.0f, -1.1f);
--j;
while (j >= 0) {
XOR_RAND(current_state, current_res);
XOR_RAND(current_state, current_res2);
res += fmaf(fmaf(current_res, 2.0f, -1.1f), sin_vals[j], fmaf(current_res2, 2.0f, -1.1f) * cos_vals[j]);
--j;
}
return res;
}
static __declspec(noalias) __forceinline float GrishaginFunc(const float x1, const float x2, const float seed) noexcept
{
int j = 0;
__declspec(align(32)) float angles_j[8u];
__declspec(align(32)) float angles_k[8u];
#pragma loop(ivdep)
while (j < 8) {
const float pj_mult = 3.14159265358979323846f * static_cast<float>(j + 1);
angles_j[j] = pj_mult * x1;
angles_k[j] = pj_mult * x2;
++j;
}
__declspec(align(32)) float sin_j[8u], cos_j[8u];
__declspec(align(32)) float sin_k[8u], cos_k[8u];
FABE13_SINCOS(angles_j, sin_j, cos_j, 8u);
FABE13_SINCOS(angles_k, sin_k, cos_k, 8u);
--j;
float part1 = 0.0f;
float part2 = 0.0f;
float current_state = seed, current_res, current_res2;
while (j >= 0) {
size_t k = 0u;
while (k < 8u) {
const float sin_term = sin_j[j] * sin_j[j];
const float cos_term = cos_k[k] * cos_k[k];
XOR_RAND_GRSH(current_state, current_res);
XOR_RAND_GRSH(current_state, current_res2);
part1 = fmaf(current_res, sin_term, fmaf(current_res2, cos_term, part1));
XOR_RAND_GRSH(current_state, current_res);
XOR_RAND_GRSH(current_state, current_res2);
part2 = fmaf(-current_res, cos_term, fmaf(current_res2, sin_term, part2));
++k;
}
--j;
}
return -sqrtf(fmaf(part1, part1, part2 * part2));
}
static __declspec(noalias) __forceinline float Shag(const float _m, const float x1, const float x2, const float y1,
const float y2, const float _N, const float _r) noexcept
{
const float diff = y2 - y1;
const float sign_mult = _mm_cvtss_f32(_mm_castsi128_ps(_mm_set1_epi32(
0x3F800000u | ((reinterpret_cast<const uint32_t>(&diff) & 0x80000000u) ^ 0x80000000u))));
return _N == 1.0f
? fmaf(-(1.0f / _m), diff, x1 + x2) * 0.5f
: _N == 2.0f
? fmaf(sign_mult / (_m * _m), diff * diff * _r, x1 + x2) * 0.5f
: fmaf(sign_mult / powf(_m, _N), powf(diff, _N) * _r, x1 + x2) * 0.5f;
}
static __declspec(noalias) __forceinline void RecomputeR_ConstM_AVX2(Interval* const* const __restrict arr, const size_t n, const float m) noexcept {
const __m256 vm = _mm256_set1_ps(m);
text__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 int limit = static_cast<int>(n & ~7u); if (n >= 8u) { while (i < static_cast<size_t>(limit)) { __declspec(align(32)) float q[8], nf[8], ord[8]; int k = 0;
#pragma loop(ivdep)
while (k < 8) {
const Interval* const __restrict p = arr[i + k];
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
ord[k] = p->ordinate_factor;
++k;
}
const __m256 vq = _mm256_load_ps(q);
const __m256 vnf = _mm256_load_ps(nf);
const __m256 vod = _mm256_load_ps(ord);
textconst __m256 t = _mm256_fmadd_ps(vm, vnf, vod); const __m256 res = _mm256_fmadd_ps(vq, vinvm, t); __declspec(align(32)) float out[8]; _mm256_store_ps(out, res); k = 0;
#pragma loop(ivdep)
while (k < 8) {
arr[i + k]->R = out[k];
++k;
}
i += 8u;
}
}
while (i < n) {
arr[i]->ChangeCharacteristic(m);
++i;
}
}
static __declspec(noalias) __forceinline void RecomputeR_AffineM_AVX2(Interval* const* const __restrict arr, const size_t n, const float GLOBAL_FACTOR, const float alpha) noexcept {
const __m256 vGF = _mm256_set1_ps(GLOBAL_FACTOR);
const __m256 va = _mm256_set1_ps(alpha);
textsize_t i = 0; const int limit = static_cast<int>(n & ~7u); if (n >= 8u) { while (i < static_cast<size_t>(limit)) { __declspec(align(32)) float len[8], Mv[8], q[8], nf[8], ord[8]; int k = 0;
#pragma loop(ivdep)
while (k < 8) {
const Interval* const p = arr[i + k];
len[k] = p->x2 - p->x1;
Mv[k] = p->M;
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
ord[k] = p->ordinate_factor;
++k;
}
const __m256 vlen = _mm256_load_ps(len);
const __m256 vM = _mm256_load_ps(Mv);
const __m256 vq = _mm256_load_ps(q);
const __m256 vnf = _mm256_load_ps(nf);
const __m256 vod = _mm256_load_ps(ord);
textconst __m256 vm = _mm256_fmadd_ps(vGF, vlen, _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); __declspec(align(32)) float out[8]; _mm256_store_ps(out, res); k = 0;
#pragma loop(ivdep)
while (k < 8) {
arr[i + k]->R = out[k];
++k;
}
i += 8u;
}
}
while (i < n) {
const Interval* const __restrict p = arr[i];
const float mi = fmaf(GLOBAL_FACTOR, 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));
++i;
}
}
__declspec(noalias) __forceinline void HitTest2D_analytic(const float x_param, float& out_x1, float& out_x2) noexcept
{
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;
textfloat norm = (x_param - a) * inv_lenx; norm = fminf(fmaxf(norm, 0.0f), 0x1.fffffep-1f); uint32_t idx = static_cast<uint32_t>(norm * static_cast<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;
#pragma loop(ivdep)
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;
}
__declspec(noalias) __forceinline float FindX2D_analytic(const float px, const float py) noexcept
{
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;
textconst 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;
#pragma loop(ivdep)
while (l < levels) {
sx *= 0.5f; sy *= 0.5f;
const float mx = x0 + sx;
const float my = y0 + sy;
textconst uint32_t tr = static_cast<uint32_t>((clamped_px > mx) & (clamped_py > my)); const uint32_t tl = static_cast<uint32_t>((clamped_px < mx) & (clamped_py > my)); const uint32_t dl = static_cast<uint32_t>((clamped_px < mx) & (clamped_py < my)); const uint32_t none = static_cast<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_reciprocal = 1.0f / static_cast<float>(scale); return fmaf(static_cast<float>(idx) * scale_reciprocal, lenx, a);
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline int AgpInit(const int peanoLevel, const float a, const float b, const float c, const float d) noexcept
{
g_env = new boost::mpi::environment();
g_world = new boost::mpi::communicator();
const int rank = g_world->rank();
textswitch (rank) { case 0: new(&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, static_cast<uint8_t>(Top)); break; case 1: new(&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, static_cast<uint8_t>(Down)); break; case 2: new(&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, static_cast<uint8_t>(Left)); break; default: new(&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, static_cast<uint8_t>(Right)); } return rank;
}
extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_1D(const float global_iterations,
const float a, const float b, const float r,
const bool mode, const float epsilon, const float seed,
float** const __restrict out_data, size_t* const __restrict out_len) noexcept
{
Slab* const __restrict slab = tls.local();
slab->current = slab->base;
textint schetchick = 0; const float initial_length = b - a; float dmax = initial_length; const float threshold_03 = 0.3f * initial_length; const float inv_threshold_03 = 1.0f / threshold_03; const float start_val = ShekelFunc(a, seed); float best_f = ShekelFunc(b, seed); float x_Rmax_1 = a; float x_Rmax_2 = b; float y_Rmax_1 = start_val; float y_Rmax_2 = best_f; std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr; std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __restrict, 64u>> R; Extr.clear(); Extr.reserve(static_cast<size_t>(global_iterations) << 2u); R.clear(); R.reserve(static_cast<size_t>(global_iterations) << 1u); R.emplace_back(new Interval(a, b, start_val, best_f, 1.0f)); float Mmax = R.front()->M; float m = r * Mmax; while (true) { const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 1.0f, r); const float new_value = ShekelFunc(new_point, seed); if (new_value < best_f) { best_f = new_value; Extr.emplace_back(best_f); Extr.emplace_back(new_point); } std::pop_heap(R.begin(), R.end(), ComparePtr); const Interval* const __restrict promejutochny_otrezok = R.back(); const float new_x1 = promejutochny_otrezok->x1; const float new_x2 = promejutochny_otrezok->x2; const float len2 = new_x2 - new_point; const float len1 = new_point - new_x1; const float interval_len = len1 < len2 ? len1 : len2; if (++schetchick == static_cast<int>(global_iterations) || interval_len < epsilon) { Extr.emplace_back(static_cast<float>(schetchick)); Extr.emplace_back(interval_len); *out_len = Extr.size(); *out_data = reinterpret_cast<float* __restrict>(CoTaskMemAlloc(sizeof(float) * (*out_len))); memcpy(*out_data, Extr.data(), sizeof(float) * (*out_len)); return; } Interval* const __restrict curr = new Interval(new_x1, new_point, promejutochny_otrezok->y1, new_value, 1.0f); Interval* const __restrict curr1 = new Interval(new_point, new_x2, new_value, promejutochny_otrezok->y2, 1.0f); const float currM = curr->M > curr1->M ? curr->M : curr1->M; const size_t r_size = R.size(); if (mode) { if (len2 + len1 == dmax) { dmax = len2 > len1 ? len2 : len1; size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}
textif (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) { if (currM > Mmax) { Mmax = currM; m = r * Mmax; } const float progress = fmaf(-inv_threshold_03, dmax, 1.0f); const float alpha = fmaf(progress, progress, 1.0f); const float betta = 2.0f - alpha; const float MULTIPLIER = (1.0f / dmax) * Mmax; const float global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER); const float GLOBAL_FACTOR = betta * global_coeff; curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha)); curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha)); if (r_size < 64u) { size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
textR.back() = curr; std::push_heap(R.begin(), R.end(), ComparePtr); R.emplace_back(curr1); std::push_heap(R.begin(), R.end(), ComparePtr); const Interval* const __restrict top_ptr = R.front(); x_Rmax_1 = top_ptr->x1; x_Rmax_2 = top_ptr->x2; y_Rmax_1 = top_ptr->y1; y_Rmax_2 = top_ptr->y2; }
}
extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_2D(
const float N, const float global_iterations,
const float a, const float b, const float c, const float d, const float r,
const bool mode, const float epsilon, const float seed,
float** const __restrict out_data, size_t* const __restrict out_len) noexcept
{
Slab* const __restrict slab = tls.local();
slab->current = slab->base;
textint schetchick = 0; int no_improve = 0; const int rank = g_world->rank(); const int world_size = g_world->size(); while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); } const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1)); const float x_addition = (b - a) * inv_divider; const float y_addition = (d - c) * inv_divider; const float true_start = a + x_addition; const float true_end = b - x_addition; float x_Rmax_1 = true_start; float x_Rmax_2 = true_end; const float initial_length = x_Rmax_2 - x_Rmax_1; float dmax = initial_length; const float threshold_03 = 0.3f * initial_length; const float inv_threshold_03 = 1.0f / threshold_03; const float start_val = rank % 3 ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition); float best_f = rank % 2 ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition); float y_Rmax_1 = start_val; float y_Rmax_2 = best_f; std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr; std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __restrict, 64u>> R; Extr.clear(); Extr.reserve(static_cast<size_t>(global_iterations) << 2u); R.clear(); R.reserve(static_cast<size_t>(global_iterations) << 1u); R.emplace_back(new Interval(true_start, true_end, start_val, best_f, 2.0f)); const Interval* __restrict top_ptr; float Mmax = R.front()->M; float m = r * Mmax; while (true) { const float interval_len = x_Rmax_2 - x_Rmax_1; const bool stagnation = no_improve > 100 && schetchick > 270; const float p = fmaf(-1.0f / initial_length, dmax, 1.0f); while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg in; g_world->recv(boost::mpi::any_source, 0, in); if (in.kind) { if (!rank) { Extr.emplace_back(static_cast<float>(schetchick)); Extr.emplace_back(interval_len); *out_len = Extr.size(); *out_data = reinterpret_cast<float* __restrict>(CoTaskMemAlloc(sizeof(float) * (*out_len))); memcpy(*out_data, Extr.data(), sizeof(float) * (*out_len)); } return; } const float sx = FindX2D_analytic(in.xchg.s_x1, in.xchg.s_x2); const float ex = FindX2D_analytic(in.xchg.e_x1, in.xchg.e_x2); Interval* const __restrict injected = new Interval(sx, ex, RastriginFunc(in.xchg.s_x1, in.xchg.s_x2), RastriginFunc(in.xchg.e_x1, in.xchg.e_x2), 2.0f); injected->ChangeCharacteristic(m); if (injected->R > 1.15f * top_ptr->R) { const float k = stagnation ? fmaf((1.0f / expm1f(1.0f)) * 1.0f, expm1f(p), 0.3f) : fmaf((1.0f / expm1f(1.0f)) * 0.6f, expm1f(p), 0.6f); injected->R = in.xchg.Rtop * k; R.emplace_back(injected); std::push_heap(R.begin(), R.end(), ComparePtr); } } const int T = static_cast<int>(fmaf(-expm1f(p), 264.0f, 277.0f)); const bool want_term = interval_len < epsilon || schetchick == static_cast<int>(global_iterations); if (!(++schetchick % T) || stagnation || want_term) { CtrlMsg out; out.kind = want_term; if (!out.kind) { float s_x1, s_x2, e_x1, e_x2; HitTest2D_analytic(top_ptr->x1, s_x1, s_x2); HitTest2D_analytic(top_ptr->x2, e_x1, e_x2); out.xchg = CrossMsg{ s_x1, s_x2, e_x1, e_x2, top_ptr->R }; } int i = 0; while (i < world_size) { if (i != rank) g_world->isend(i, 0, out); ++i; } if (out.kind) { if (!rank) { Extr.emplace_back(static_cast<float>(schetchick)); Extr.emplace_back(interval_len); *out_len = Extr.size(); *out_data = reinterpret_cast<float* __restrict>(CoTaskMemAlloc(sizeof(float) * (*out_len))); memcpy(*out_data, Extr.data(), sizeof(float) * (*out_len)); } return; } } const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 2.0f, r); float new_x1_val, new_x2_val; HitTest2D_analytic(new_point, new_x1_val, new_x2_val); const float new_value = RastriginFunc(new_x1_val, new_x2_val); if (new_value < best_f) { best_f = new_value; Extr.emplace_back(best_f); Extr.emplace_back(new_x1_val); Extr.emplace_back(new_x2_val); no_improve = 0; } else { ++no_improve; } std::pop_heap(R.begin(), R.end(), ComparePtr); Interval* const __restrict promejutochny_otrezok = R.back(); const float segment_x1 = promejutochny_otrezok->x1; const float segment_x2 = promejutochny_otrezok->x2; const float len2 = segment_x2 - new_point; const float len1 = new_point - segment_x1; Interval* const __restrict curr = new Interval(segment_x1, new_point, promejutochny_otrezok->y1, new_value, 2.0f); Interval* const __restrict curr1 = new Interval(new_point, segment_x2, new_value, promejutochny_otrezok->y2, 2.0f); const float currM = curr->M > curr1->M ? curr->M : curr1->M; const size_t r_size = R.size(); if (mode) { if (len2 + len1 == dmax) { dmax = len2 > len1 ? len2 : len1; size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}
textif (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) { if (currM > Mmax) { Mmax = currM; m = r * Mmax; } const float progress = fmaf(-inv_threshold_03, dmax, 1.0f); const float alpha = fmaf(progress, progress, 1.0f); const float betta = 2.0f - alpha; const float MULTIPLIER = (1.0f / dmax) * Mmax; const float global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER); const float GLOBAL_FACTOR = betta * global_coeff; curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha)); curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha)); if (r_size < 64u) { size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
textR.back() = curr; std::push_heap(R.begin(), R.end(), ComparePtr); R.emplace_back(curr1); std::push_heap(R.begin(), R.end(), ComparePtr); top_ptr = R.front(); x_Rmax_1 = top_ptr->x1; x_Rmax_2 = top_ptr->x2; y_Rmax_1 = top_ptr->y1; y_Rmax_2 = top_ptr->y2; }
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun() noexcept
{
int dummy;
float* __restrict buf;
size_t len;
while (true) {
if (g_world->iprobe(0, 1)) {
g_world->recv(0, 1, dummy);
AGP_2D(2.0f, 10000.0f, -2.2f, 1.8f, -2.2f, 1.8f, 2.5f, false, 0.00001f, GetTickCount(), &buf, &len);
}
}
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpStartWorkers() noexcept
{
int i = 1;
const int world_size = g_world->size();
while (i < world_size) {
g_world->isend(i, 1, 0);
++i;
}
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* const __restrict p) noexcept
{
CoTaskMemFree(p);
} - параллельный двумерный алгоритм сейчас запускается на 4ёх процессорах и работает быстро и корректно - заметь что скорость и оптимизация важны, запускаю алгоритм в управляемом приложении, код MyForm.cpp: #include "MyForm.h"
#include <float.h>
using namespace System;
using namespace System::Windows::Forms;
typedef int(__cdecl* PInit)(int, float, float, float, float);
typedef void(__cdecl* PStartWorkers)();
[STAThread]
int main() {
HMODULE h = LoadLibraryW(L"TEST_FUNC.dll");
auto AgpInit = (PInit)GetProcAddress(h, "AgpInit");
auto AgpWaitStartAndRun = (PStartWorkers)GetProcAddress(h, "AgpWaitStartAndRun");
textconst int rank = AgpInit(12, -2.2f, 1.8f, -2.2f, 1.8f); if (!rank) { Application::EnableVisualStyles(); Application::SetCompatibleTextRenderingDefault(false); Application::Run(gcnew TESTAGP::MyForm(h)); } else { AgpWaitStartAndRun(); } return 0;
}, а это код MyForm.h: #pragma once
#define WIN32_LEAN_AND_MEAN
#include <Windows.h>
#include <algorithm>
#include <vector>
#include <utility>
namespace TESTAGP {
using namespace System;
using namespace System::IO;
textref class MyForm : public System::Windows::Forms::Form { public: MyForm(HMODULE hLib) : hLib(hLib) { // Передаем дескриптор из main InitializeComponent(); // Загружаем функции из DLL f = (agp_c)GetProcAddress(hLib, "AGP_2D"); pStart = (start_workers)GetProcAddress(hLib, "AgpStartWorkers"); pFree = (free_agp)GetProcAddress(hLib, "AGP_Free"); pInit = (PInit)GetProcAddress(hLib, "AgpInit"); } protected: ~MyForm() { delete components; } private: HINSTANCE hLib = nullptr; typedef void(__cdecl* agp_c)(float, float, float, float, float, float, float, bool, float, float, float**, size_t*); typedef void(__cdecl* start_workers)(); typedef void(__cdecl* free_agp)(float*); typedef int(__cdecl* PInit)(int, float, float, float, float); agp_c f = nullptr; start_workers pStart = nullptr; free_agp pFree = nullptr; PInit pInit = nullptr; System::Windows::Forms::Button^ button1; System::Windows::Forms::TextBox^ textBox2; System::Windows::Forms::DataVisualization::Charting::Chart^ chart2; System::Windows::Forms::Label^ label2; System::Windows::Forms::Label^ label3; System::Windows::Forms::TextBox^ textBox1; System::Windows::Forms::TextBox^ textBox3; System::Windows::Forms::TextBox^ textBox4; System::Windows::Forms::TextBox^ textBox5; System::Windows::Forms::TextBox^ textBox6; System::Windows::Forms::Label^ label6; System::Windows::Forms::Label^ label7; System::Windows::Forms::Label^ label8; System::Windows::Forms::Label^ label9; System::Windows::Forms::Label^ label10; System::Windows::Forms::Label^ label1; System::Windows::Forms::TextBox^ textBox7; System::Windows::Forms::TextBox^ textBox8; System::ComponentModel::Container^ components; System::Void button1_Click(System::Object^ sender, System::EventArgs^ e) { chart2->Series[0]->Points->Clear(); chart2->Series[1]->Points->Clear(); chart2->Series[2]->Points->Clear(); chart2->Series[3]->Points->Clear(); static LARGE_INTEGER start, end; QueryPerformanceCounter(&start); float* buf = nullptr; size_t len = 0u; pStart(); // разбудили rank != 0 f(2.0f, 10000.0f, -2.2f, 1.8f, -2.2f, 1.8f, 2.5f, false, 0.00001f, GetTickCount(), &buf, &len); if (buf != nullptr) { std::vector<float> Extr_2D(buf, buf + len); pFree(buf); // Для 2D ветки: извлекаем schetchick и interval_len textBox7->Text = Convert::ToString(Extr_2D.back()); // schetchick Extr_2D.pop_back(); textBox6->Text = Convert::ToString(Extr_2D.back()); // interval_len Extr_2D.pop_back(); // Извлекаем последнюю тройку (x2, x1, f) float x2 = Extr_2D.back(); Extr_2D.pop_back(); float x1 = Extr_2D.back(); Extr_2D.pop_back(); float f_val = Extr_2D.back(); Extr_2D.pop_back(); // Выводим координаты и значение функции textBox4->Text = Convert::ToString(x1); textBox3->Text = Convert::ToString(x2); textBox2->Text = Convert::ToString(f_val); // Новое поле для значения функции // Отображаем точку на графике chart2->Series[2]->Points->AddXY(x1, x2); // Обрабатываем остальные точки (в порядке: f, x1, x2) while (!Extr_2D.empty()) { float x2_point = Extr_2D.back(); Extr_2D.pop_back(); float x1_point = Extr_2D.back(); Extr_2D.pop_back(); float f_point = Extr_2D.back(); Extr_2D.pop_back(); // Добавляем точку на график (только координаты) chart2->Series[0]->Points->AddXY(x1_point, x2_point); } } QueryPerformanceCounter(&end); LARGE_INTEGER freq; QueryPerformanceFrequency(&freq); const int multiplier = ((1ULL << 32) / freq.QuadPart) * 1'000'000ULL; textBox5->Text = Convert::ToString( ((end.QuadPart - start.QuadPart) * multiplier) >> 32 ) + " microseconds"; //textBox2->Text = Convert::ToString(Extr_2D.back()); //textBox1->Text = Convert::ToString(Extr_2D_LNA.back()); } - от меня требуется создать интерфейс - я буду делать его в WinForms - где пользователь может добавлять и убирать из конфигурации звенья по нажатии на конпки, и также должны быть ограничения на углы - правда я не помню какое требование углы точно - но наверное будет неплохо если пользователь также сможет задавать ограничения на углы через интерфейс - например добавить текстовые поля куда можно будет заносить допустимый диапазон углов - например от нуля до двух пи, скорее всего требования на углы - это то что в диссертации называется требованием на минимизацию кривизны, подскажи какое именно должно быть тогда требование под минимизацию кривизны - углы не должны быть слишком острыми или же углы не должны быть слишком тупыми? но заметь это не основное требование - научный руководитель говорил реализовать это через систему штрафов - вводить штрафы на углы которые не удовлетворяют диапазону - возможно для начала можно применить простую стратегию штрафа линейно пропорционального вылету угла за диапазон, также заметь что мне придётся отрисовывать звенья и узлы в интерфейсе - не в объекте chart - а лучше рисовать эти объекты прямо на форме без введения лишних элементов управления, наверное по умолчанию при запуске программы пользователь должен будет видеть один шарнир и один стержень присоединённый к нему под случайным углом (можно использовать мою быструю реалзацию xor_rand), потом пользователь также должен мочь настроить финальное положение захвата где он должен быть - на каком удалении от начального шарнира, на какой высоте, то есть координату в двумерной плоскости (всё происходит в двумерной плоскости) - шарниры можно рисовать просто как кружки фиксированного размера - например синие жирные заметные кружки но не слишком большие чтобы они не перетягивали акцент на себя, а именно секции соединяющие шарниры то есть стержни между ними можно рисовать как жирные красные отрезки или вытянутые прямоугольники фиксированной ширины, кружки должны рисоваться условно на переднем плане а красные прямоугольники условно на заднем чтобы было красиво, я ещё не определился вот по какому вопросу - когда пользователь будет добавлять звено - то будет добавляться звено фиксированной длины а потом запускать алгоритм чтобы он искал именно оптимальные углы между звеньями или же можно также разрешать алгоритму менять и длину звеньев? наверное второй случай предпочтительнее так как он более общий, именно сама задача минимизации или оптимизации должна идти именно по минимизации положения захвата то есть минимум расстояния финального положения захвата до требуемой точки, штрафы на углы это сопутствующая фича - но базово - по сути тоже что у меня для малоразмерных случаев - то есть также будет требоваться найти минимум целевой функции - только в этом случае целевой функцией будет функция расстояния финального шарнира до требуемого положения которая принимает в качестве аргументов все обобщённые координаты - то есть и углы и длины стержней, будет хорошо если картинка будет подвижной - то есть по нажатию на кнопку - к которой у меня сейчас привязано решение маломерных задач - должно происходить динамическое перерисовывание конфигурации - в зависимости о того какие оптимальные углы и длины звеньев находит алгоритм в процессе оптимизации - но это может потребовать значительных изменений в логике - так что возможно сделать чтобы по нажатию на кнопку просто стиралась начальная конфигурация и на её месте рисовалась оптимальная конфигурация, заметь что для этой цели придётся написать отдельный класс для многомерных гиперкубов пеано - потому что сейчас у меня только двумерные развёртки - а теперь нужно по аналогии использовать тот же принцип но для многомерных развёрток - они также должны быть аналитическими и считаться быстро, сам алгоритм не должен сильно принципиально отличаться например от двумерной версии, заметь что новый многомерный алгоритм также должен поддерживать параллелизм, но заметь что для двумерного случая используются 4 по разному ориентированные развёртки что обеспечивает максимально полное покрытие поискового пространства, для многомерного случая может получится намного большее число возможных ориентаций гиперкубовых развёрток - так что надо исходить из доступного числа процессов, и учитывать что этот код будет запускаться на обычных пк где как правило не большое число ядер - как например всего 6 ядер на моём пк, думай как можно дольше чтобы выдать самый корректный код покрывающий мои потребности, ты должен предоставить код для .cpp в dll, код для MyForm.cpp и код для MyForm.h, учти что я работаю в устаревшем NetFramework а не в Net
Suy nghĩ
Trả lời
Ниже — готовая «скелетная» реализация, которая решает все ваши задачи:
- формулирует требование «минимизации кривизны» для плоской цепочки звеньев;
- добавляет штрафы за выход углов/длин за диапазоны (линейные, как вы предлагали);
- даёт быстрый n-мерный параметрический проход по пространству обобщённых координат с помощью аналитической быстрой развертки (вариант Morton/Z-order с перестановками осей для «ориентаций»);
- расширяет ваш AGP до n-мерного случая (одномерный параметр → n-мерная точка), с той же «характеристикой» и адаптацией
m; - поддерживает параллельный поиск по нескольким ориентациям (многопоточность; для простоты — стандартные
std::thread, чтобы не привязываться к MPI в первом приближении; при желании это место можно заменить на Boost.MPI аналогично вашей 2D-версии); - даёт WinForms-интерфейс под .NET Framework (C++/CLI), где пользователь:
- добавляет/удаляет звенья;
- задаёт общие диапазоны углов (θ_min, θ_max) и длин (L_min, L_max), порог «излома» Δθ_max и веса штрафов;
- задаёт двумерную целевую точку (X*, Y*);
- видит отрисовку шарниров (синие круги «сверху») и стержней (толстые красные отрезки «снизу») прямо на форме без Chart;
- при «Оптимизировать» вычисляется лучшая конфигурация и перерисовывается итог.
Если захотите «псевдо-анимацию» (перерисовку по улучшениям), достаточно в DLL возвращать историю «лучших» и в форме — вызывать Invalidate() по мере чтения.
Что значит «минимизация кривизны» для вашей задачи
Для плоского манипулятора с углами сочленений θ₁,…,θₙ естественная дискретная кривизна — это «излом» между соседними сегментами:
- локальная кривизна на стыке
iпропорциональна |Δθᵢ|, гдеΔθᵢ = wrap(θ_{i+1} − θ_i)(нормировка в (−π, π]); - если хотите учесть длины, используйте
|Δθᵢ| / ((Lᵢ + L_{i+1})/2).
Практический критерий: ограничить «излом» сверху |Δθᵢ| ≤ Δθ_max, и нарушители штрафовать. Линейный штраф, как вы предложили:
pen_curv = λ_curv * Σ_i max(0, |Δθᵢ| − Δθ_max)
(или квадратичный — чуть «мягче», но линейный даёт простую и быструю реализацию).
Отдельно — «коробочные» ограничения на углы и длины:
textpen_angle = λ_ang * Σ_i ( max(0, θ_min − θᵢ) + max(0, θᵢ − θ_max) ) pen_len = λ_len * Σ_i ( max(0, L_min − Lᵢ) + max(0, Lᵢ − L_max) )
Целевая функция:
textF(q) = || EE(q) − (X*,Y*) || // расстояние от захвата до цели (в идеале можно брать квадрат) + pen_angle + pen_len + pen_curv
Где EE(q) — кинематика конечного звена:
textφ_k = θ₁ + … + θ_k x = Σ_k L_k cos(φ_k) y = Σ_k L_k sin(φ_k)
DLL (C++17, нативная): agp_nd.cpp
Экспортирует одну главную функцию
AGP_Manip2D, принимающую:
- число звеньев
nSegments;- флаги и общие диапазоны углов/длин;
- параметры штрафов и поиска;
- целевую точку; Возвращает: лучшую конфигурацию
q_best(углы [+ длины, если разрешены]), координаты лучшего захвата, значение целевой.
cpp// agp_nd.cpp — поместите в проект DLL (Multi-threaded (/MT) желательно для простоты дистрибуции) // Компилируется без внешних зависимостей; AVX и MPI можно добавить позднее. // ВАЖНО: Экспортируемые функции — C-интерфейс как у вас. #include <cmath> #include <cfloat> #include <cstdint> #include <vector> #include <thread> #include <atomic> #include <algorithm> #include <cstring> #include <mutex> #include <random> #include <windows.h> #ifndef DLLEXPORT #define DLLEXPORT __declspec(dllexport) #endif // ---------- Быстрый XOR RNG в стиле вашего кода ---------- static inline float xor_rand01(uint32_t& s) noexcept { s ^= s << 13; s ^= s >> 17; s ^= s << 5; // 0x1.0p-32f = 2^-32 return (float)( (double)(uint32_t)s * (1.0/4294967296.0) ); } static inline float xor_rand_signed(uint32_t& s) noexcept { s ^= s << 13; s ^= s >> 17; s ^= s << 5; // в [-1, 1] return std::ldexp((float)(int32_t)s, -30) - 1.0f; // грубо, но быстро } // ---------- Обёртка под выделение/возврат буферов ---------- extern "C" DLLEXPORT void AGP_Free(float* p) noexcept { CoTaskMemFree(p); } // ---------- Мортоновская (Z-order) развертка для nD ---------- // Аналитически быстрая: interleave битов индекса в координаты. // Для разнообразия "ориентаций" используем перестановку осей и инверсии. struct MortonND { int dim; // D int levels; // уровни разбиения (итоговая сетка 2^levels по каждой оси) uint64_t scale; // 1 << (levels * dim) std::vector<float> low, high; // коробка поиска std::vector<int> perm; // перестановка осей std::vector<uint64_t> invMask;// какие оси инвертировать (по ориентации) MortonND(int D, int L, const float* lows, const float* highs, uint32_t orientationSeed) : dim(D), levels(L), scale( (uint64_t)1 << (L * D) ), low(lows, lows + D), high(highs, highs + D), perm(D), invMask(D, 0) { // Простая псевдослучайная перестановка осей for (int i = 0; i < D; ++i) perm[i] = i; // Fisher-Yates на основе xor-rand uint32_t st = orientationSeed ? orientationSeed : 0x9E3779B9u; for (int i = D-1; i > 0; --i) { int j = (int)std::floor(xor_rand01(st) * (i+1)); std::swap(perm[i], perm[j]); } // Маски инверсий по осям (повороты/зеркала) for (int i = 0; i < D; ++i) { invMask[i] = (xor_rand01(st) > 0.5f) ? ((uint64_t)1 << levels) - 1ull : 0ull; } } // Преобразование t in [0,1] -> точка в D-мерной коробке (центр ячейки, чтобы повторить идею «середин») inline void map01ToPoint(float t, float* __restrict out) const noexcept { // Квантуем t в индекс ячейки: if (t < 0.0f) t = 0.0f; else if (t >= 1.0f) t = std::nextafter(1.0f, 0.0f); uint64_t idx = (uint64_t)(t * (double)scale); if (idx >= scale) idx = scale - 1; // Достанем по уровням биты и соберём координатные индексы: std::vector<uint64_t> coordIdx(dim, 0ull); for (int bit = 0; bit < levels; ++bit) { for (int d = 0; d < dim; ++d) { int pd = perm[d]; // перестановка uint64_t bitmask = (idx >> (bit*dim + d)) & 1ull; coordIdx[pd] |= (bitmask << bit); } } // Инвертирование по некоторым осям (зеркалирование): for (int d = 0; d < dim; ++d) { if (invMask[d]) { coordIdx[d] ^= invMask[d]; } } const float invCells = 1.0f / (float)( (uint64_t)1 << levels ); for (int d = 0; d < dim; ++d) { float cellStart = low[d] + (high[d] - low[d]) * ( (float)coordIdx[d] * invCells ); float cellSize = (high[d] - low[d]) * invCells; out[d] = cellStart + 0.5f * cellSize; // центр ячейки } } }; // ---------- Вспомогательные куски под вашу AGP-логику ---------- struct Interval { float x1, x2; // параметр t на [0,1] float y1, y2; // значения F(t) на концах float R; // характеристика float N_factor, quadTerm, M, ordFactor; // как у вас (упрощённо, N = dim) Interval(float _x1, float _x2, float _y1, float _y2, float N) : x1(_x1), x2(_x2), y1(_y1), y2(_y2) { float dy = (y2 - y1); float dx = (x2 - x1); // N_factor — аналог вашей степени; пусть N=dim (см. ниже) N_factor = (N == 1.0f ? dx : (N == 2.0f ? std::sqrt(dx) : std::pow(dx, 1.0f/N))); quadTerm = (1.0f / N_factor) * dy * dy; M = (1.0f / N_factor) * std::fabs(dy); ordFactor = -(y1 + y2) * 2.0f; R = 0.0f; } inline void changeR(float m) noexcept { R = (quadTerm / m) + (m * N_factor) + ordFactor; } }; static inline bool CmpIntervalPtr(const Interval* a, const Interval* b) noexcept { return a->R < b->R; } // шаг вычисления t (ваш Shag): для N=dim static inline float Shag(float m, float x1, float x2, float y1, float y2, float N, float r) noexcept { const float dy = (y2 - y1); float sign_mult = (std::signbit(dy) ? -1.0f : 1.0f); if (N == 1.0f) { return ((x1 + x2) * 0.5f) - (dy / m) * 0.5f; } else if (N == 2.0f) { return ((x1 + x2) * 0.5f) + sign_mult * ( (dy*dy) * r / (m*m) ) * 0.5f; } else { return ((x1 + x2) * 0.5f) + sign_mult * ( std::pow(std::fabs(dy), N) * r / std::pow(m, N) ) * 0.5f; } } // Обработка углов в (-pi, pi] static inline float wrap_pi(float a) noexcept { a = std::fmod(a + (float)M_PI, 2.0f * (float)M_PI); if (a < 0) a += 2.0f * (float)M_PI; return a - (float)M_PI; } // Кинематика и целевая функция struct ManipCost { int n; // звеньев bool variableLen; float targetX, targetY; // Общие (глобальные) диапазоны: θ in [thMin, thMax], L in [LMin, LMax] float thMin, thMax; float LMin, LMax; // Штрафы float dThetaMax; // порог излома float lambdaCurv; float lambdaAng; float lambdaLen; // Вычисление стоимости для q = [θ₁..θₙ, (возможно L₁..Lₙ)] inline float operator()(const float* q) const noexcept { const float* th = q; const float* L = variableLen ? (q + n) : nullptr; // 1) расстояние до цели float x = 0.0f, y = 0.0f, phi = 0.0f; float penAngle = 0.0f, penLen = 0.0f, penCurv = 0.0f; float Lconst = 1.0f; for (int i=0; i<n; ++i) { const float theta = th[i]; const float Li = variableLen ? std::max(LMin, std::min(LMax, L[i])) : Lconst; phi += theta; x += Li * std::cos(phi); y += Li * std::sin(phi); // штрафы за угол if (theta < thMin) penAngle += (thMin - theta); if (theta > thMax) penAngle += (theta - thMax); // штраф за длину (если переменная) if (variableLen) { if (L[i] < LMin) penLen += (LMin - L[i]); if (L[i] > LMax) penLen += (L[i] - LMax); } } // 2) штраф за «кривизну» (изломы) for (int i=0; i<n-1; ++i) { float dth = wrap_pi(th[i+1] - th[i]); // излом между соседними float viol = std::fabs(dth) - dThetaMax; if (viol > 0.0f) penCurv += viol; // линейный штраф } const float dx = (x - targetX); const float dy = (y - targetY); const float dist = std::sqrt(dx*dx + dy*dy); return dist + lambdaAng*penAngle + lambdaLen*penLen + lambdaCurv*penCurv; } }; // Одна «ветка» AGP по конкретной ориентации MortonND static void agp_run_branch( const MortonND& map, const ManipCost& cost, int maxIter, float r, bool adaptiveMode, float eps, uint32_t seed, // out: std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY ){ const int dim = cost.n + (cost.variableLen ? cost.n : 0); auto evalAt = [&](float t)->float { std::vector<float> q(dim); map.map01ToPoint(t, q.data()); // Преобразуем из [low,high] уже учтённых в map — там всё учтено float f = cost(q.data()); // Для обратного восстановления EE(x,y) (необязательно, но полезно вернуть) float x=0,y=0,phi=0; const float* th = q.data(); const float* L = cost.variableLen ? (q.data()+cost.n) : nullptr; float Lconst = 1.0f; for (int i=0;i<cost.n;++i){ phi += th[i]; float Li = cost.variableLen ? L[i] : Lconst; x += Li * std::cos(phi); y += Li * std::sin(phi); } // Если это пока лучший — сохраним всё if (f < bestF) { bestF = f; bestQ = q; bestX = x; bestY = y; } return f; }; // Стартовые концы [0,1] float a = 0.0f, b = 1.0f; float f_a = evalAt(a); float f_b = evalAt(b); // Куча интервалов по вашей логике std::vector<Interval*> H; H.reserve( (size_t)maxIter ); Interval* root = new Interval(a, b, f_a, f_b, (float)dim); float Mmax = root->M; float m = r * Mmax; root->changeR(m); H.push_back(root); std::make_heap(H.begin(), H.end(), CmpIntervalPtr); float dmax = b - a; const float initial_len = dmax; const float thr03 = 0.3f * initial_len; const float inv_thr03 = 1.0f / thr03; int it = 0; while (it < maxIter) { // Берём лучший интервал std::pop_heap(H.begin(), H.end(), CmpIntervalPtr); Interval* cur = H.back(); H.pop_back(); const float x1 = cur->x1, x2 = cur->x2; const float y1 = cur->y1, y2 = cur->y2; // Шаг float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r); tNew = std::max(a, std::min(b, tNew)); float fNew = evalAt(tNew); // Разделяем Interval* left = new Interval(x1, tNew, y1, fNew, (float)dim); Interval* right = new Interval(tNew, x2, fNew, y2, (float)dim); float Mloc = std::max(left->M, right->M); // Адаптация m (упрощённая версия вашей логики) if (adaptiveMode) { float len1 = tNew - x1; float len2 = x2 - tNew; if (len1 + len2 == dmax) { dmax = std::max(len1, len2); for (auto p: H) dmax = std::max(dmax, p->x2 - p->x1); } if ( (thr03 > dmax && !(it % 3)) || (10.0f * dmax < initial_len) ) { if (Mloc > Mmax) { Mmax = Mloc; m = r * Mmax; } float progress = 1.0f - (dmax * inv_thr03); float alpha = progress*progress; float beta = 2.0f - alpha; float MULT = (Mmax / dmax); float global_coeff = MULT * (r - 1.0f); float GF = beta * global_coeff; left->changeR( GF*len1 + left->M*alpha ); right->changeR( GF*len2 + right->M*alpha ); for (auto p: H) { float len = p->x2 - p->x1; p->changeR( GF*len + p->M*alpha ); } } else { if (Mloc > Mmax) { if (Mloc < 1.15f*Mmax) { Mmax = Mloc; m = r*Mmax; } else { Mmax = Mloc; m = r*Mmax; for (auto p: H) p->changeR(m); } } left->changeR(m); right->changeR(m); } } else { if (Mloc > Mmax) { if (Mloc < 1.15f*Mmax) { Mmax = Mloc; m = r*Mmax; } else { Mmax = Mloc; m = r*Mmax; for (auto p: H) p->changeR(m); } } left->changeR(m); right->changeR(m); } H.push_back(left); std::push_heap(H.begin(), H.end(), CmpIntervalPtr); H.push_back(right); std::push_heap(H.begin(), H.end(), CmpIntervalPtr); // Критерий останова float currLen = std::max(left->x2 - left->x1, right->x2 - right->x1); if (currLen < eps) break; ++it; } // Освобождение for (auto p: H) delete p; } // ---------- Экспорт: запуск оптимизации манипулятора ---------- extern "C" DLLEXPORT void AGP_Manip2D( int nSegments, // общие диапазоны углов и длин (для простоты одинаковые для всех звеньев) float thetaMin, float thetaMax, bool variableLengths, float lenMin, float lenMax, // штрафные параметры float deltaThetaMax, float lambdaCurv, float lambdaAngle, float lambdaLength, // целевая точка float targetX, float targetY, // параметры развертки/поиска int peanoLevels, // глубина (≥ 8 обычно достаточно; 10–14 для тонкого поиска) int orientations, // сколько «ориентаций» обойти параллельно (1..число ядер) int maxIterPerBranch, // лимит итераций на ветку float r, // как в вашем коде, например 2.5 bool adaptiveMode, float epsilon, unsigned int seed, // out: float** out_bestQ, size_t* out_bestQLen, float* out_bestX, float* out_bestY, float* out_bestF ){ const int dim = nSegments + (variableLengths ? nSegments : 0); // Коробка поиска [low, high] в каждой координате q: std::vector<float> low(dim), high(dim); for (int i=0; i<nSegments; ++i) { low[i] = thetaMin; high[i] = thetaMax; } if (variableLengths) { for (int i=0; i<nSegments; ++i) { low[nSegments+i] = lenMin; high[nSegments+i] = lenMax; } } ManipCost cost { nSegments, variableLengths, targetX, targetY, thetaMin, thetaMax, lenMin, lenMax, deltaThetaMax, lambdaCurv, lambdaAngle, lambdaLength }; // Параллельно обходим несколько ориентаций Morton-развёртки и выбираем лучшую int branches = std::max(1, orientations); std::atomic<float> globalBest(FLT_MAX); std::mutex bestGuard; std::vector<float> bestQ(dim, 0.0f); float bestFVal = FLT_MAX, bestXVal=0.0f, bestYVal=0.0f; auto worker = [&](uint32_t seedOrient){ MortonND map(dim, peanoLevels, low.data(), high.data(), seedOrient); std::vector<float> locBestQ; float locBestF = FLT_MAX, locBestX=0.0f, locBestY=0.0f; agp_run_branch(map, cost, maxIterPerBranch, r, adaptiveMode, epsilon, seed+seedOrient, locBestQ, locBestF, locBestX, locBestY); float prev = globalBest.load(); while (locBestF < prev && !globalBest.compare_exchange_weak(prev, locBestF)) {} if (locBestF <= globalBest.load()) { std::lock_guard<std::mutex> lk(bestGuard); if (locBestF < bestFVal) { bestFVal = locBestF; bestXVal = locBestX; bestYVal = locBestY; bestQ = locBestQ; } } }; std::vector<std::thread> threads; threads.reserve(branches); uint32_t orientSeed = seed ? seed : 0xA53C9E0Du; for (int i=0;i<branches;++i) { threads.emplace_back(worker, orientSeed + i*0x9E3779B9u); } for (auto& t: threads) t.join(); // Возврат пользователю *out_bestQLen = (size_t)dim; *out_bestQ = (float*)CoTaskMemAlloc(sizeof(float)*(*out_bestQLen)); std::memcpy(*out_bestQ, bestQ.data(), sizeof(float)*(*out_bestQLen)); *out_bestX = bestXVal; *out_bestY = bestYVal; *out_bestF = bestFVal; }
Примечание: внутри
MortonNDя использую перестановки/зеркалирования осей по «ориентациям» — это «имитирует» идею разных развёрток, как у вас в 2D (Top/Down/Left/Right), только в общем n-мерном случае, без сложного кода nD-Hilbert. Если захотите «настоящий» Hilbert/Butz — вставьте реализацию вMortonND(интерфейс останется тем же).
WinForms (.NET Framework, C++/CLI)
Ниже минимальная форма без дизайнера, только код. Рисуем прямо на самой форме (двойная буферизация включена), шарниры — синие кружки, стержни — красные отрезки. Пользователь может добавлять/удалять звенья, задать глобальные диапазоны и штрафы, указать цель, запустить оптимизацию. При старте показываем 1 шарнир и 1 стержень со случайным углом (XOR-RNG).
MyForm.h
cpp#pragma once #define WIN32_LEAN_AND_MEAN #include <Windows.h> #include <vector> #include <algorithm> #include <cmath> #include <stdint.h> using namespace System; using namespace System::Drawing; using namespace System::Windows::Forms; typedef void(__cdecl* P_MANIP)( int nSegments, float thetaMin, float thetaMax, bool variableLengths, float lenMin, float lenMax, float deltaThetaMax, float lambdaCurv, float lambdaAngle, float lambdaLength, float targetX, float targetY, int peanoLevels, int orientations, 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 ); typedef void(__cdecl* P_FREE)(float*); namespace TESTAGP { public ref class MyForm : public Form { public: MyForm(HMODULE hLib) : hLib(hLib) { this->SetStyle(ControlStyles::AllPaintingInWmPaint | ControlStyles::UserPaint | ControlStyles::OptimizedDoubleBuffer, true); this->Text = L"AGP Manipulator 2D"; this->ClientSize = System::Drawing::Size(1000, 700); fManip = (P_MANIP)GetProcAddress(hLib, "AGP_Manip2D"); pFree = (P_FREE)GetProcAddress(hLib, "AGP_Free"); initUI(); resetRandomConfig(); // один стержень, случайный угол } protected: ~MyForm() { } private: // DLL HMODULE hLib = nullptr; P_MANIP fManip = nullptr; P_FREE pFree = nullptr; // Модель конфигурации для рисования int nSegments = 1; bool variableLengths = false; std::vector<float> angles; // текущая конфигурация для отрисовки std::vector<float> lengths; // Параметры ввода NumericUpDown^ nudThetaMin; NumericUpDown^ nudThetaMax; CheckBox^ cbVarLen; NumericUpDown^ nudLMin; NumericUpDown^ nudLMax; NumericUpDown^ nudDeltaThetaMax; NumericUpDown^ nudLamCurv; NumericUpDown^ nudLamAng; NumericUpDown^ nudLamLen; NumericUpDown^ nudTargetX; NumericUpDown^ nudTargetY; NumericUpDown^ nudLevels; NumericUpDown^ nudOrient; NumericUpDown^ nudMaxIter; CheckBox^ cbAdaptive; NumericUpDown^ nudR; NumericUpDown^ nudEps; Button^ btnAdd; Button^ btnRem; Button^ btnOptimize; Label^ lblInfo; // RNG uint32_t rngState = 0xA5C3'9E0Du; // Инициализация UI void initUI() { int x = 10, y = 10, w = 140, h = 24, pad = 6; auto mkLabel = [&](String^ t)->Label^ { auto L = gcnew Label(); L->Text = t; L->Left = x; L->Top = y; L->Width = w; this->Controls->Add(L); y += h + 2; return L; }; auto mkNum = [&](double min, double max, double val, int decimals)->NumericUpDown^ { auto N = gcnew NumericUpDown(); N->Left = x; N->Top = y; N->Width = w; N->Minimum = (Decimal)min; N->Maximum=(Decimal)max; N->DecimalPlaces = decimals; N->Increment = (Decimal)0.01; N->Value = (Decimal)val; this->Controls->Add(N); y += h + pad; return N; }; auto mkChk = [&](String^ t, bool v)->CheckBox^ { auto C = gcnew CheckBox(); C->Text = t; C->Left = x; C->Top = y; C->Width = w+40; C->Checked = v; this->Controls->Add(C); y += h + pad; return C; }; auto mkBtn = [&](String^ t, EventHandler^ on)->Button^ { auto B = gcnew Button(); B->Text = t; B->Left = x; B->Top = y; B->Width = w; this->Controls->Add(B); y += h + pad; B->Click += on; return B; }; mkLabel(L"θ_min (rad)"); nudThetaMin = mkNum(-3.14159, 3.14159, -3.14159, 3); mkLabel(L"θ_max (rad)"); nudThetaMax = mkNum(-3.14159, 3.14159, 3.14159, 3); cbVarLen = mkChk(L"Переменные длины", false); mkLabel(L"L_min"); nudLMin = mkNum(0.1, 5.0, 1.0, 2); mkLabel(L"L_max"); nudLMax = mkNum(0.1, 5.0, 1.0, 2); mkLabel(L"Δθ_max (rad)"); nudDeltaThetaMax = mkNum(0.0, 3.14159, 0.8, 3); mkLabel(L"λ_curv"); nudLamCurv = mkNum(0.0, 1000.0, 1.0, 2); mkLabel(L"λ_angle"); nudLamAng = mkNum(0.0, 1000.0, 1.0, 2); mkLabel(L"λ_length"); nudLamLen = mkNum(0.0, 1000.0, 1.0, 2); mkLabel(L"Target X"); nudTargetX = mkNum(-10.0, 10.0, 1.5, 2); mkLabel(L"Target Y"); nudTargetY = mkNum(-10.0, 10.0, 1.0, 2); mkLabel(L"Levels"); nudLevels = mkNum(4, 20, 12, 0); mkLabel(L"Orientations"); nudOrient = mkNum(1, 16, 4, 0); mkLabel(L"MaxIter/branch"); nudMaxIter = mkNum(10, 50000, 5000, 0); cbAdaptive = mkChk(L"Adaptive m", true); mkLabel(L"r"); nudR = mkNum(1.0, 10.0, 2.5, 2); mkLabel(L"eps"); nudEps = mkNum(1e-6, 1e-1, 1e-4, 5); btnAdd = mkBtn(L"+ звено", gcnew EventHandler(this, &MyForm::onAdd)); btnRem = mkBtn(L"− звено", gcnew EventHandler(this, &MyForm::onRem)); btnOptimize = mkBtn(L"Оптимизировать", gcnew EventHandler(this, &MyForm::onOptimize)); lblInfo = gcnew Label(); lblInfo->Left = x; lblInfo->Top = y; lblInfo->Width = w+100; lblInfo->Height = 60; this->Controls->Add(lblInfo); } // Случайная начальная конфигурация: 1 звено void resetRandomConfig() { nSegments = 1; angles.assign(1, randAngle()); lengths.assign(1, 1.0f); Invalidate(); } float rand01() { // XOR rngState ^= rngState << 13; rngState ^= rngState >> 17; rngState ^= rngState << 5; return (float)((double)(uint32_t)rngState / 4294967296.0); } float randAngle() { return (rand01() * 2.0f * (float)M_PI) - (float)M_PI; } // Обработчики void onAdd(System::Object^, System::EventArgs^) { ++nSegments; angles.push_back(randAngle()); lengths.push_back(1.0f); Invalidate(); } void onRem(System::Object^, System::EventArgs^) { if (nSegments > 1) { --nSegments; angles.pop_back(); lengths.pop_back(); Invalidate(); } } void onOptimize(System::Object^, System::EventArgs^) { if (!fManip) return; // Чтение параметров float thMin = (float)nudThetaMin->Value; float thMax = (float)nudThetaMax->Value; variableLengths = cbVarLen->Checked; float Lmin = (float)nudLMin->Value; float Lmax = (float)nudLMax->Value; float dthMax = (float)nudDeltaThetaMax->Value; float lamC = (float)nudLamCurv->Value; float lamA = (float)nudLamAng->Value; float lamL = (float)nudLamLen->Value; float tx = (float)nudTargetX->Value; float ty = (float)nudTargetY->Value; int levels = (int)nudLevels->Value; int orients = (int)nudOrient->Value; int maxIter = (int)nudMaxIter->Value; bool adaptive = cbAdaptive->Checked; float r = (float)nudR->Value; float eps = (float)nudEps->Value; // Запуск LARGE_INTEGER t0,t1,fq; QueryPerformanceCounter(&t0); float* bestQ = nullptr; size_t bestQLen = 0; float bestX=0, bestY=0, bestF=0; fManip( nSegments, thMin, thMax, variableLengths, Lmin, Lmax, dthMax, lamC, lamA, lamL, tx, ty, levels, orients, maxIter, r, adaptive, eps, (unsigned)GetTickCount(), &bestQ, &bestQLen, &bestX, &bestY, &bestF ); QueryPerformanceCounter(&t1); QueryPerformanceFrequency(&fq); double micros = 1e6 * (double)(t1.QuadPart - t0.QuadPart) / (double)fq.QuadPart; // Применяем лучшую конфигурацию к текущему состоянию (чтобы нарисовать) if (bestQ && bestQLen >= (size_t)nSegments) { angles.resize(nSegments); for (int i=0;i<nSegments;++i) angles[i] = bestQ[i]; if (variableLengths) { lengths.resize(nSegments); for (int i=0;i<nSegments;++i) lengths[i] = bestQ[nSegments+i]; } else { lengths.assign(nSegments, 1.0f); } pFree(bestQ); } lblInfo->Text = String::Format(L"Лучшее F={0:F5}\nEE=({1:F3},{2:F3})\n{3:F0} μs", bestF, bestX, bestY, micros); Invalidate(); // перерисовать } protected: // Рисование — красные стержни сзади, синие шарниры сверху. virtual void OnPaint(PaintEventArgs^ e) override { Form::OnPaint(e); Graphics^ g = e->Graphics; g->SmoothingMode = System::Drawing::Drawing2D::SmoothingMode::HighQuality; // Координатная сетка: центр формы — основание манипулятора int cx = 380, cy = this->ClientSize.Height - 60; float scale = 80.0f; // пикс/ед. // Вычисляем узлы из текущих angles/lengths std::vector<System::Drawing::PointF> pts; pts.reserve(nSegments+1); float x=0.0f, y=0.0f, phi=0.0f; pts.emplace_back((float)cx, (float)cy); // базовый шарнир for (int i=0;i<nSegments;++i) { phi += angles[i]; float L = variableLengths ? lengths[i] : 1.0f; x += L * std::cos(phi); y += L * std::sin(phi); pts.emplace_back( cx + x*scale, cy - y*scale ); } // Рисуем стержни (красные, потолще) Pen^ penRod = gcnew Pen(Color::Red, 6.0f); for (int i=0;i<nSegments;++i) { g->DrawLine(penRod, pts[i], pts[i+1]); } // Рисуем шарниры (синие кружки поверх) SolidBrush^ brJoint = gcnew SolidBrush(Color::Blue); const float R = 6.0f; for (int i=0;i<=nSegments;++i) { g->FillEllipse(brJoint, pts[i].X - R, pts[i].Y - R, 2*R, 2*R); } } }; } // namespace TESTAGP
MyForm.cpp
cpp#include "MyForm.h" using namespace System; using namespace System::Windows::Forms; typedef int(__cdecl* PInitDummy)(); // на будущее, если понадобится инициализация [STAThread] int main() { // Подгружаем DLL HMODULE h = LoadLibraryW(L"TEST_FUNC.dll"); // убедитесь, что имя совпадает с собранной DLL if (!h) { MessageBox::Show(L"Не удалось загрузить TEST_FUNC.dll", L"Ошибка", MessageBoxButtons::OK, MessageBoxIcon::Error); return 1; } Application::EnableVisualStyles(); Application::SetCompatibleTextRenderingDefault(false); Application::Run(gcnew TESTAGP::MyForm(h)); return 0; }
Как это стыкуется с вашей текущей 2D-версией
- В нативной DLL вы получаете единый параметр t ∈ [0,1], по Morton-развертке быстро (аналітично) превращаете его в вектор обобщённых координат
q(углы + длины), считаете целевуюF(q)с штрафами, и ведёте отбор по вашей «характеристике» (R) и адаптацииm. - Для «покрытия» пространства разными ориентациями у нас параллельно запускается несколько веток с разными перестановками/зеркалами осей — аналог вашей идеи «4 ориентации» в 2D, только в n-мерном случае их может быть сколько угодно; берём лучшее из веток.
- Если вам принципиально нужен именно Boost.MPI (как в вашем
AGP_2D) — логика веток почти напрямую переносится: вместоstd::threadраздать разную ориентацию по ранкам, а результаты и/или «инъекции» интервалов пересылать так же, как в вашем коде (я намеренно построил ветку автономной: она не требует координации, чтобы MVP был простым и быстрым).
Нюансы/дефолты (можете скорректировать под диссертацию):
- По умолчанию длины фиксированные (=1) — снимите галочку «Переменные длины», чтобы исследовать только углы. Если включить — используйте глобальные
L_min/L_maxдля всех звеньев. Δθ_max— допустимый «излом». Для мягких конфигураций 0.6–1.0 радиан обычно вполне.- Штрафы: начните с
λ_angle=λ_length=λ_curv=1.0, а затем масштабируйте (если цель близка к 0, увеличьте штрафы, чтобы ограничения «слышались»). - В
AGP_Manip2DpeanoLevels=12,orientations=4..6,maxIterPerBranch=2–10 тыс.— хорошие старты для бытового ПК.
Если захотите «живую» перерисовку в процессе (а не только оптимум) — верните из DLL массив «best-trace» (например, каждые k улучшений — текущие angles/lengths) и в форме на таймере перерисовывайте. Логика отрисовки уже готова — просто сменяйте angles/lengths и вызывайте Invalidate().
Готов подсказать, как аккуратно перекинуть текущую схему Boost.MPI (как у вашего AGP_2D) на предложенную AGP_Manip2D, если решите оставить именно межпроцессную параллель.