in core/indigo-core/layout/src/molecule_layout_macrocycles.cpp [691:1251]
float MoleculeLayoutMacrocycles::depictionMacrocycleGreed(bool /* profi */)
{
{
profTimerStart(tt, "answerField");
// MoleculeLayoutMacrocyclesLattice molecule_layout_lattice(length);
}
if (length >= max_size)
return 1e9;
unsigned short(&minRotates)[max_size][max_size][2][max_size][max_size] = data.minRotates;
// first : number of edge
// second : summary angle of rotation (in M_PI/3 times)
// third : last rotation is contraclockwise
// fourth : x-coordinate
// fifth : y-coordinate
// value : minimum number of vertexes sticked out + count of CIS-configurations
unsigned short infinity = 60000;
int shift = -1;
int max_value = -infinity;
for (int i = 0; i < length; i++)
{
if (_edge_stereo[i] != 2)
{
int value = _edge_stereo[i] + _vertex_weight[i] + _vertex_weight[(i + 1) % length] - _vertex_weight[(i + length - 1) % length] / 2 -
_vertex_weight[(i + 2) % length] / 2;
if (shift == -1 || value > max_value)
{
shift = i;
max_value = value;
}
}
}
if (shift == -1)
return 1e9;
rotate(_vertex_weight.ptr(), length, shift);
rotate(_vertex_stereo.ptr(), length, shift);
rotate(_edge_stereo.ptr(), length, shift);
rotate(_angle_importance.ptr(), length, shift);
int dx[6];
int dy[6];
dx[0] = 1;
dy[0] = 0;
dx[1] = 0;
dy[1] = 1;
dx[2] = -1;
dy[2] = 1;
dx[3] = -1;
dy[3] = 0;
dx[4] = 0;
dy[4] = -1;
dx[5] = 1;
dy[5] = -1;
int x_left = max(init_x - length, 1);
int x_right = min(init_x + length, max_size - 2);
int y_left = max(init_y - length, 1);
int y_right = min(init_y + length, max_size - 2);
int rot_left = max(init_rot - length, 1);
int rot_right = min(init_rot + length, max_size - 2);
for (int i = 0; i <= length; i++)
for (int j = rot_left - 1; j <= rot_right + 1; j++)
for (int p = 0; p < 2; p++)
for (int k = x_left - 1; k <= x_right + 1; k++)
for (int t = y_left - 1; t <= y_right + 1; t++)
minRotates[i][j][p][k][t] = infinity;
minRotates[0][init_rot][1][init_x][init_y] = 0;
minRotates[1][init_rot][1][init_x + 1][init_y] = 0;
// int max_dist = length;
{
profTimerStart(tt, "main.for");
for (int k = 1; k < length; k++)
{
// printf("Step number %d\n", k);
int not_this_p = -1;
if (k == length - 1)
{
if (_edge_stereo[length - 1] == MoleculeCisTrans::CIS)
not_this_p = 0;
if (_edge_stereo[length - 1] == MoleculeCisTrans::TRANS)
not_this_p = 1;
}
for (int rot = rot_left; rot <= rot_right; rot++)
{
if (!_vertex_stereo[k])
{
int xchenge = dx[rot % 6];
int ychenge = dy[rot % 6];
for (int p = 0; p < 2; p++)
if (p != not_this_p)
{
for (int x = x_left; x <= x_right; x++)
{
unsigned short* ar1 = minRotates[k + 1][rot][p][x + xchenge] + ychenge;
unsigned short* ar2 = minRotates[k][rot][p][x];
for (int y = y_left; y <= y_right; y++)
{
if (ar1[y] > ar2[y])
{
ar1[y] = ar2[y];
}
}
}
}
}
else
{
for (int p = 0; p < 2; p++)
{
// trying to rotate like CIS
if (_edge_stereo[k - 1] != MoleculeCisTrans::TRANS)
if (p != not_this_p)
{
int nextRot = rot;
if (p)
nextRot++;
else
nextRot--;
int xchenge = dx[nextRot % 6];
int ychenge = dy[nextRot % 6];
int add = 1;
if (abs(_vertex_weight[k]) > WEIGHT_FACTOR)
{
if (!p && _vertex_weight[k] > 0)
add += _vertex_weight[k];
if (p && _vertex_weight[k] < 0)
add -= _vertex_weight[k];
}
for (int x = x_left; x <= x_right; x++)
{
unsigned short* ar1 = minRotates[k + 1][nextRot][p][x + xchenge] + ychenge;
unsigned short* ar2 = minRotates[k][rot][p][x];
for (int y = y_left; y <= y_right; y++)
{
if (ar1[y] > ar2[y] + add)
{
ar1[y] = static_cast<unsigned short>(ar2[y] + add);
}
}
}
}
// trying to rotate like TRANS
if (_edge_stereo[k - 1] != MoleculeCisTrans::CIS)
if ((p ^ 1) != not_this_p)
{
int nextRot = rot;
if (p)
nextRot--;
else
nextRot++;
int add = 0;
if (abs(_vertex_weight[k]) > WEIGHT_FACTOR)
{
if (p && _vertex_weight[k] > 0)
add += _vertex_weight[k];
if (!p && _vertex_weight[k] < 0)
add -= _vertex_weight[k];
}
int xchenge = dx[nextRot % 6];
int ychenge = dy[nextRot % 6];
for (int x = x_left; x <= x_right; x++)
{
unsigned short* ar1 = minRotates[k + 1][nextRot][p ^ 1][x + xchenge] + ychenge;
unsigned short* ar2 = minRotates[k][rot][p][x];
for (int y = y_left; y <= y_right; y++)
{
if (ar1[y] > ar2[y] + add)
{
ar1[y] = static_cast<unsigned short>(ar2[y] + add);
}
}
}
}
}
}
}
}
}
/* for (int rot = rot_left; rot <= rot_right; rot++)
for (int x = x_left; x <= x_right; x++)
for (int y = y_left; y <= y_right; y++)
minRotates[length][rot][1][x][y]++;*/
struct point
{
int diff;
int x;
int y;
int p;
int rot;
int type;
point(int _d, int _x, int _y, int _p, int _r, int _t)
{
diff = _d;
x = _x;
y = _y;
p = _p;
rot = _r;
type = _t;
}
bool operator<(const point& pt) const
{
return diff - pt.diff;
}
};
QS_DEF(ObjArray<point>, points);
points.clear();
int critical_diff_grid = infinity - 1;
std::multiset<int> diff_set;
int var_count = 100;
enum
{
CIRCLE_TYPE,
GRID_TYPE
};
for (int rot = init_rot; rot <= rot_right; rot++)
for (int p = 0; p < 2; p++)
for (int x = x_left; x <= x_right; x++)
for (int y = y_left; y <= y_right; y++)
if (minRotates[length][rot][p][x][y] <= infinity)
{
int curr_dif = get_diff_grid(x, y, rot, minRotates[length][rot][p][x][y]);
if (curr_dif <= critical_diff_grid)
{
diff_set.insert(curr_dif);
if (static_cast<long>(diff_set.size()) > var_count)
diff_set.erase(*diff_set.rbegin());
critical_diff_grid = *diff_set.rbegin();
}
}
for (int rot = init_rot; rot <= rot_right; rot++)
{
for (int p = 0; p < 2; p++)
{
for (int x = x_left; x <= x_right; x++)
{
for (int y = y_left; y <= y_right; y++)
{
if (minRotates[length][rot][p][x][y] < infinity)
{
int curdiff = get_diff_grid(x, y, rot, minRotates[length][rot][p][x][y]);
if (curdiff <= critical_diff_grid)
{
point curr_point(curdiff, x, y, p, rot, GRID_TYPE);
points.push(curr_point);
}
}
}
}
}
}
diff_set.clear();
int critical_diff_circle = infinity - 1;
for (int rot = rot_left; rot <= rot_right; rot++)
for (int p = 0; p < 2; p++)
for (int x = x_left; x <= x_right; x++)
for (int y = y_left; y <= y_right; y++)
if (minRotates[length][rot][p][x][y] <= infinity)
{
int curr_dif = get_diff_circle(x - init_x - length / 2, y - length / 2, rot, minRotates[length][rot][p][x][y]);
if (curr_dif <= critical_diff_circle)
{
diff_set.insert(curr_dif);
if (static_cast<long>(diff_set.size()) > var_count)
diff_set.erase(*diff_set.rbegin());
critical_diff_circle = *diff_set.rbegin();
}
}
for (int rot = rot_left; rot <= rot_right; rot++)
{
for (int p = 0; p < 2; p++)
{
for (int x = x_left; x <= x_right; x++)
{
for (int y = y_left; y <= y_right; y++)
{
if (minRotates[length][rot][p][x][y] < infinity)
{
int curdiff = get_diff_circle(x - init_x - length / 2, y - length / 2, rot, minRotates[length][rot][p][x][y]);
if (curdiff <= critical_diff_circle)
{
point curr_point(curdiff, x, y, p, rot, CIRCLE_TYPE);
points.push(curr_point);
}
}
}
}
}
}
/* for (int rot = rot_left; rot <= rot_right; rot++)
for (int x = x_left; x <= x_right; x++)
for (int y = y_left; y <= y_right; y++)
minRotates[length][rot][1][x][y]--;*/
int x_result[max_size + 1];
int y_result[max_size + 1];
int rot_result[max_size + 1];
int p_result[max_size + 1];
float bestBadness = 1e30f;
int bestIndex = 0;
for (int index = 0; index < points.size(); index++)
{
profTimerStart(tt, "selector.for");
x_result[length] = points[index].x;
y_result[length] = points[index].y;
rot_result[length] = points[index].rot;
p_result[length] = points[index].p;
for (int k = length - 1; k > 0; k--)
{
int xchenge = dx[rot_result[k + 1] % 6];
int ychenge = dy[rot_result[k + 1] % 6];
x_result[k] = x_result[k + 1] - xchenge;
y_result[k] = y_result[k + 1] - ychenge;
if (!_vertex_stereo[k])
{
p_result[k] = p_result[k + 1];
rot_result[k] = rot_result[k + 1];
}
else
{
if (p_result[k + 1])
rot_result[k] = rot_result[k + 1] - 1;
else
rot_result[k] = rot_result[k + 1] + 1;
float l = _2FLOAT(k * (sqrt(3.0) + 1.5) * M_PI / 12.);
Vec2f vec(_2FLOAT(y_result[k] - init_y), 0.f);
vec.rotate(_2FLOAT(M_PI / 3.));
vec += Vec2f(_2FLOAT(x_result[k] - init_x), 0.f);
float x = vec.length();
const float eps = 1e-3f;
float alpha = 0;
if (x > eps)
{
float L = eps;
float R = _2FLOAT(2. * M_PI - eps);
while (R - L > eps)
{
float M = (L + R) / 2;
if (M * x / (2 * sin(M / 2)) > l)
R = M;
else
L = M;
}
alpha = vec.tiltAngle2() + R / 2;
}
p_result[k] = 2;
int is_cis_better = (alpha < M_PI / 3 * (rot_result[k] - init_rot) + M_PI / length) ^ (!p_result[k + 1]);
int add = 0;
if (abs(_vertex_weight[k]) > WEIGHT_FACTOR)
{
if (!p_result[k + 1] && _vertex_weight[k] > 0)
add = _vertex_weight[k];
if (p_result[k + 1] && _vertex_weight[k] < 0)
add = -_vertex_weight[k];
}
if (!is_cis_better)
{
if (_edge_stereo[k - 1] != MoleculeCisTrans::TRANS)
{
// try CIS
if (minRotates[k][rot_result[k]][p_result[k + 1]][x_result[k]][y_result[k]] + add + 1 ==
minRotates[k + 1][rot_result[k + 1]][p_result[k + 1]][x_result[k + 1]][y_result[k + 1]])
{
p_result[k] = p_result[k + 1];
}
}
}
if (_edge_stereo[k - 1] != MoleculeCisTrans::CIS)
{
// try TRANS
if (minRotates[k][rot_result[k]][p_result[k + 1] ^ 1][x_result[k]][y_result[k]] + add ==
minRotates[k + 1][rot_result[k + 1]][p_result[k + 1]][x_result[k + 1]][y_result[k + 1]])
{
p_result[k] = p_result[k + 1] ^ 1;
}
}
if (is_cis_better)
{
if (_edge_stereo[k - 1] != MoleculeCisTrans::TRANS)
{
// try CIS
if (minRotates[k][rot_result[k]][p_result[k + 1]][x_result[k]][y_result[k]] + add + 1 ==
minRotates[k + 1][rot_result[k + 1]][p_result[k + 1]][x_result[k + 1]][y_result[k + 1]])
{
p_result[k] = p_result[k + 1];
}
}
}
if (p_result[k] == 2)
{
throw Error("Path not find (%d): %d.", length, minRotates[k + 1][rot_result[k + 1]][p_result[k + 1]][x_result[k + 1]][y_result[k + 1]]);
}
}
}
x_result[0] = init_x;
y_result[0] = init_y;
int rotateAngle[max_size];
int edgeLenght[max_size];
int vertexNumber[max_size];
int ind = 0;
for (int i = 0; i < length; i++)
{
if (_vertex_stereo[i])
vertexNumber[ind++] = i;
}
if (ind < 3)
{
ind = 0;
for (int i = 0; i < length; i++)
vertexNumber[ind++] = i;
}
vertexNumber[ind] = length;
for (int i = 0; i < ind - 1; i++)
edgeLenght[i] = vertexNumber[i + 1] - vertexNumber[i];
edgeLenght[ind - 1] = vertexNumber[0] - vertexNumber[ind - 1] + length;
for (int i = 0; i < ind; i++)
if (vertexNumber[i] != 0)
{
rotateAngle[i] = rot_result[vertexNumber[i] + 1] > rot_result[vertexNumber[i]] ? 1
: rot_result[vertexNumber[i] + 1] == rot_result[vertexNumber[i]] ? 0
: -1;
}
if (vertexNumber[0] == 0)
{
rotateAngle[0] = 1;
}
// rotateAngle[0] = -1;
Vec2f p[max_size];
for (int i = 0; i <= ind; i++)
{
p[i] = Vec2f(_2FLOAT(y_result[vertexNumber[i]]), 0.f);
p[i].rotate(_2FLOAT(M_PI / 3.));
p[i] += Vec2f(_2FLOAT(x_result[vertexNumber[i]]), 0.f);
}
for (int i = ind; i >= 0; i--)
p[i] -= p[0];
if (points[index].type == CIRCLE_TYPE)
{
if (p[ind].lengthSqr() < EPSILON)
continue;
Vec2f rotate_vector(p[ind] / p[ind].length());
for (int i = 0; i <= ind; i++)
{
p[i].rotateL(rotate_vector);
}
float max_x = p[ind].x;
if (max_x <= 0)
max_x = 1;
float radii = _2FLOAT(max_x / (2. * M_PI));
for (int i = 0; i <= ind; i++)
{
float angle = _2FLOAT(2. * M_PI * (p[i].x) / max_x);
p[i].set(radii - p[i].y / 2.f, 0.f);
p[i].rotate(angle);
}
}
/* Vec2f last_vector(p[0] - p[ind]);
for (int i = 0; i <= ind; i++) p[i] += (last_vector * vertexNumber[i]) / length;*/
// smoothing(ind, length, rotateAngle, edgeLenght, vertexNumber, p, profi);
// for (int i = 0; i <= ind; i++ ) printf("%5.5f %5.5f\n", p[i].x, p[i].y);
smoothing2(ind, length, rotateAngle, edgeLenght, vertexNumber, p);
int diff = minRotates[length][points[index].rot][points[index].p][points[index].x][points[index].y];
for (int i = 0; i < length - 1; i++)
{
if (rotateAngle[i] == 1 && rotateAngle[i + 1] == 1)
diff--;
if (rotateAngle[i] == -1 && rotateAngle[i + 1] == -1)
diff--;
}
float newBadness = badness(ind, length, rotateAngle, edgeLenght, vertexNumber, p, diff);
if (newBadness < bestBadness)
{
bestBadness = newBadness;
bestIndex = index;
for (int i = 0; i < ind; i++)
{
int nexti = (i + 1) % ind;
for (int j = vertexNumber[i], t = 0; j != vertexNumber[nexti]; j = (j + 1) % length, t++)
{
_positions[j] = (p[i] * _2FLOAT(edgeLenght[i] - t) + p[nexti] * _2FLOAT(t)) / _2FLOAT(edgeLenght[i]);
}
}
}
}
Vec2f shifted_positons[max_size];
for (int i = 0; i < length; i++)
shifted_positons[(i + shift) % length] = _positions[i];
for (int i = 0; i < length; i++)
_positions[i] = shifted_positons[i];
rotate(_vertex_weight.ptr(), length, -shift);
rotate(_vertex_stereo.ptr(), length, -shift);
rotate(_edge_stereo.ptr(), length, -shift);
rotate(_angle_importance.ptr(), length, -shift);
return bestBadness;
}