461 lines
12 KiB
C
461 lines
12 KiB
C
#include <math.h>
|
|
#ifdef BENCHMARK
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#endif
|
|
|
|
#include "rng.h"
|
|
#include "rbtree.h"
|
|
#include "metric.h"
|
|
#include "utils.h"
|
|
#include "cryst.h"
|
|
#ifdef BENCHMARK
|
|
#include "bench_cryst.h"
|
|
#endif
|
|
#include "int_cryst.h"
|
|
|
|
#ifdef BENCHMARK
|
|
struct sweep {
|
|
float pos;
|
|
unsigned idx;
|
|
int end;
|
|
};
|
|
|
|
struct cpair {
|
|
unsigned idx[2];
|
|
};
|
|
|
|
struct bump_t {
|
|
struct rbtree active;
|
|
struct sweep *poss;
|
|
struct cpair *bumps[3];
|
|
unsigned cnt[3], n;
|
|
};
|
|
|
|
struct bump_t *bump_mk (unsigned axes, unsigned n) {
|
|
struct bump_t *state;
|
|
unsigned i, j;
|
|
|
|
if (!(state = malloc (sizeof (struct bump_t)))) goto state_err;
|
|
for (i = 0; i < 3; ++i) state->bumps[i] = NULL;
|
|
for (i = 0; i < 3; ++i) if (axes >> i & 0x1 && !(
|
|
state->bumps[i] = calloc (n * n, sizeof (struct cpair))
|
|
)) goto poss_err;
|
|
if (!(state->poss = malloc (2 * n * sizeof (struct sweep)))) goto poss_err;
|
|
|
|
for (j = i = 0; i < n; ++i) {
|
|
state->poss[j++] = (struct sweep) { .idx = i, .end = 0 };
|
|
state->poss[j++] = (struct sweep) { .idx = i, .end = 1 };
|
|
}
|
|
rbtree_mk (&state->active);
|
|
state->n = n;
|
|
return state;
|
|
|
|
free (state->poss);
|
|
poss_err:
|
|
for (i = 0; i < 3; ++i) if (state->bumps[i]) free (state->bumps[i]);
|
|
free (state);
|
|
state_err:
|
|
return NULL;
|
|
}
|
|
|
|
void bump_fin (struct bump_t *state) {
|
|
unsigned i;
|
|
free (state->poss);
|
|
for (i = 0; i < 3; ++i) if (state->bumps[i]) free (state->bumps[i]);
|
|
free (state);
|
|
}
|
|
|
|
int bump_eval_orig (struct crystal const *cr, narrow_f *narrow, signed *ret_) {
|
|
struct bump_t const *state = cr->bbump;
|
|
signed ret = 0, score;
|
|
for (unsigned i = 0; i < state->n; ++i) {
|
|
for (unsigned j = i + 1; j < state->n; ++j) {
|
|
if ((*narrow) (cr, i, j, &score)) ret += score;
|
|
}
|
|
}
|
|
*ret_ = ret;
|
|
return 1;
|
|
}
|
|
|
|
int bump_eval_bc (struct crystal const *cr, narrow_f *narrow, signed *ret_) {
|
|
struct bump_t const *state = cr->bbump;
|
|
struct cr_ball const *balls = cr->atoms->scats->balls;
|
|
signed ret = 0, score;
|
|
for (unsigned i = 0; i < state->n; ++i) if (balls[i].flag & 0x2) {
|
|
for (unsigned j = 0; j < state->n; ++j) {
|
|
if (i == j || (balls[j].flag & 0x2 && i > j)) continue;
|
|
if ((*narrow) (cr, i, j, &score)) ret += score;
|
|
}
|
|
}
|
|
*ret_ = ret;
|
|
return 1;
|
|
}
|
|
|
|
static int sweep_cmp (void const *p1, void const *p2) {
|
|
struct sweep const *s1 = p1, *s2 = p2;
|
|
float const x = s1->pos - s2->pos;
|
|
return x != 0.0 ? x > 0.0 : s2->end - s1->end;
|
|
}
|
|
|
|
static int cpair_cmp (void const *p1, void const *p2) {
|
|
struct cpair const *c1 = p1, *c2 = p2;
|
|
return c1->idx[0] != c2->idx[0] ?
|
|
c1->idx[0] - c2->idx[0] : c1->idx[1] - c2->idx[1];
|
|
}
|
|
|
|
static int bump_axis_bb (
|
|
struct bump_t *state, struct cr_ball const balls[], unsigned axis
|
|
) {
|
|
struct rbtree *active = &state->active;
|
|
struct rbnode *node, *tmp;
|
|
struct cpair *bump = state->bumps[axis];
|
|
unsigned nn = 2 * state->n, cnt, i, j;
|
|
|
|
for (j = i = 0; i < state->n; ++i) {
|
|
state->poss[j++] = (struct sweep)
|
|
{ .idx = i, .end = 0, .pos = balls[i].bound[axis][0] };
|
|
state->poss[j++] = (struct sweep)
|
|
{ .idx = i, .end = 1, .pos = balls[i].bound[axis][1] };
|
|
}
|
|
qsort (state->poss, nn, sizeof (struct sweep), &sweep_cmp);
|
|
cnt = 0;
|
|
|
|
for (i = 0; i < nn; ++i) {
|
|
j = state->poss[i].idx;
|
|
if (state->poss[i].end) {
|
|
if ((tmp = rbtree_get (active, j))) rbtree_rm (active, tmp);
|
|
} else {
|
|
if (!(tmp = malloc (sizeof (struct rbnode)))) return 0;
|
|
for (
|
|
node = rbtree_first (active); node;
|
|
node = rbtree_next (active, node)
|
|
) if ((balls[j].flag | balls[node->peer].flag) & 0x2) {
|
|
bump[cnt++] = j < node->peer ?
|
|
(struct cpair) {{ j, node->peer }} :
|
|
(struct cpair) {{ node->peer, j }};
|
|
}
|
|
tmp->peer = j;
|
|
rbtree_add (active, tmp);
|
|
}
|
|
}
|
|
|
|
// Here duplicates might be introduced.
|
|
for (i = 0; active->root != &active->nil; ++i) {
|
|
j = state->poss[i].idx;
|
|
if (state->poss[i].end) {
|
|
if ((tmp = rbtree_get (active, j))) rbtree_rm (active, tmp);
|
|
} else if (!rbtree_get (active, j)) for (
|
|
node = rbtree_first (active); node;
|
|
node = rbtree_next (active, node)
|
|
) if ((balls[j].flag | balls[node->peer].flag) & 0x2) {
|
|
bump[cnt++] = j < node->peer ?
|
|
(struct cpair) {{ j, node->peer }} :
|
|
(struct cpair) {{ node->peer, j }};
|
|
}
|
|
}
|
|
|
|
rbtree_clear (active);
|
|
qsort (bump, cnt, sizeof (struct cpair), &cpair_cmp);
|
|
state->cnt[axis] = cnt;
|
|
return 1;
|
|
}
|
|
|
|
int bump_eval_bb (struct crystal const *cr, narrow_f *narrow, signed *ret_) {
|
|
struct bump_t *state = cr->bbump;
|
|
struct cr_ball const *balls = cr->atoms->scats->balls;
|
|
struct cpair *cur[3], *end[3], cp = {{ 0, 0 }};
|
|
unsigned ptr[4], max = 0, i, j;
|
|
signed ret = 0, score;
|
|
|
|
for (i = 0; i < 3; ++i) {
|
|
if (state->bumps[i]) {
|
|
if (!bump_axis_bb (state, balls, i)) return 0;
|
|
if (!state->cnt[i]) goto retn;
|
|
end[max] = state->bumps[i] + state->cnt[i];
|
|
cur[max++] = state->bumps[i];
|
|
}
|
|
ptr[i] = i;
|
|
}
|
|
for (i = 1; i < max; ++i) for (
|
|
j = i; j && cpair_cmp (cur[ptr[j]], cur[ptr[j - 1]]) < 0; --j
|
|
) {
|
|
ptr[3] = ptr[j - 1]; ptr[j - 1] = ptr[j]; ptr[j] = ptr[3];
|
|
}
|
|
|
|
for (j = 0;;) {
|
|
struct cpair tmp = *cur[ptr[0]];
|
|
if (!cpair_cmp (&cp, &tmp)) ++j;
|
|
else {
|
|
cp = tmp;
|
|
j = 1;
|
|
}
|
|
// `ret' is only used in narrow_count(), so add 1 instead of `score'.
|
|
if (j == max && (*narrow) (cr, cp.idx[0], cp.idx[1], &score)) ++ret;
|
|
|
|
// Also handle duplicates.
|
|
do if (++cur[ptr[0]] == end[ptr[0]]) {
|
|
// See above for why `++ret' is used.
|
|
for (i = 1; i < max && !cpair_cmp (&cp, cur[ptr[i]]); ++i) if (
|
|
++j == max && (*narrow) (cr, cp.idx[0], cp.idx[1], &score)
|
|
) ++ret;
|
|
goto retn;
|
|
} while (!cpair_cmp (&tmp, cur[ptr[0]]));
|
|
|
|
for (
|
|
i = 1;
|
|
i < max && cpair_cmp (cur[ptr[i]], cur[ptr[i - 1]]) < 0;
|
|
++i
|
|
) {
|
|
ptr[3] = ptr[i - 1]; ptr[i - 1] = ptr[i]; ptr[i] = ptr[3];
|
|
}
|
|
}
|
|
|
|
retn:
|
|
*ret_ = ret;
|
|
return 1;
|
|
}
|
|
#endif
|
|
|
|
void cryst_step (struct crystal *cr, float len) {
|
|
struct cr_pos *pos = cr->poss + cr->perm[cr->nPos[1]];
|
|
pos->pos[0] += len * pos->bound[2];
|
|
if (pos->pos[0] > pos->bound[1]) pos->pos[0] -= pos->bound[2];
|
|
else if (pos->pos[0] < pos->bound[0]) pos->pos[0] += pos->bound[2];
|
|
pos->flag = 1;
|
|
|
|
if (++cr->nPos[1] == cr->nPos[0]) {
|
|
rng_shuf (cr->rng, cr->perm, cr->nPos[0]);
|
|
cr->nPos[1] = 0;
|
|
}
|
|
}
|
|
|
|
void cryst_ack (struct crystal *cr, int keep) {
|
|
struct cr_pos *poss = cr->poss;
|
|
for (unsigned i = 0; i < cr->nPos[0]; ++i) if (poss[i].flag) {
|
|
if (keep) {
|
|
poss[i].pos[1] = poss[i].pos[0];
|
|
poss[i].flag = 0;
|
|
} else poss[i].pos[0] = poss[i].pos[1];
|
|
}
|
|
}
|
|
|
|
static void spec_eval (struct crystal *cr) {
|
|
for (unsigned i = 0; i < cr->nHill[0]; ++i) {
|
|
struct cr_hill *hill = cr->hills + i;
|
|
hill->val = 0.0;
|
|
for (unsigned j = 0; j < hill->hill.n; ++j) {
|
|
float sc[2] = {
|
|
(float) hill->peaks[j].sc[0] / cr->scales[0],
|
|
(float) hill->peaks[j].sc[1] / cr->scales[0]
|
|
};
|
|
hill->val += hill->peaks[j].peak.weight *
|
|
(sc[0] * sc[0] + sc[1] * sc[1]);
|
|
}
|
|
}
|
|
}
|
|
|
|
#ifdef BENCHMARK
|
|
static void spec_eval_bc (struct crystal *cr) {
|
|
for (unsigned i = 0; i < cr->nHill[0]; ++i) {
|
|
struct cr_hill *hill = cr->hills + i;
|
|
hill->val = 0.0;
|
|
for (unsigned j = 0; j < hill->hill.n; ++j) {
|
|
float *sc = hill->peaks[j].bsc;
|
|
hill->val += hill->peaks[j].peak.weight *
|
|
(sc[0] * sc[0] + sc[1] * sc[1]);
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
|
|
static float tv_eval (struct crystal *cr) {
|
|
float ret = 0.0, sum = 0.0;
|
|
unsigned i;
|
|
for (i = 0; i < cr->nHill[0]; ++i) sum += cr->hills[i].val;
|
|
for (i = 0; i < cr->nHill[0]; ++i) {
|
|
cr->hills[i].val /= sum;
|
|
ret += fabs (cr->hills[i].val - cr->hills[i].hill.val0);
|
|
}
|
|
return 0.5 * ret;
|
|
}
|
|
|
|
static float metric_dist (
|
|
metric const *met, float const a[3], float const b[3]
|
|
) {
|
|
float dif[3];
|
|
for (unsigned i = 0; i < 3; ++i) dif[i] = b[i] - a[i];
|
|
return sqrt (metric_get (met, dif));
|
|
}
|
|
|
|
#ifdef BENCHMARK
|
|
int narrow_orig (
|
|
struct crystal const *cr, unsigned i, unsigned j, signed *ret_
|
|
) {
|
|
struct cr_ball const *balls = cr->atoms->scats->balls;
|
|
float ret = metric_dist (cr->met, balls[i].xyz, balls[j].xyz) /
|
|
balls[i].dists[balls[j].elem];
|
|
*ret_ = cr->scales[1] * (ret > cr->ctl[1] ? 0.0 : (
|
|
ret < cr->ctl[0] ? 1.0 : cr->ctl[2] * (cr->ctl[1] - ret)
|
|
));
|
|
return ret < cr->ctl[1];
|
|
}
|
|
|
|
#else
|
|
static
|
|
#endif
|
|
int narrow_cryst (
|
|
struct crystal const *cr, unsigned i, unsigned j, signed *ret_
|
|
) {
|
|
struct cr_ball const *balls = cr->atoms->scats->balls;
|
|
float ret = metric_dist (cr->met, balls[i].xyz, balls[j].xyz) /
|
|
balls[i].dists[balls[j].elem];
|
|
*ret_ = cr->scales[1] * (balls[i].n + balls[j].n) *
|
|
(ret > cr->ctl[1] ? 0.0 : (ret < cr->ctl[0] ? 1.0 :
|
|
cr->ctl[2] * (cr->ctl[1] - ret)
|
|
));
|
|
return ret < cr->ctl[1];
|
|
}
|
|
|
|
static int bump_eval (struct crystal *cr, narrow_f *narrow) {
|
|
struct rbforest *bump = &cr->bump;
|
|
struct cr_ball *balls = cr->atoms->scats->balls;
|
|
signed score;
|
|
unsigned i, j;
|
|
for (i = 0; i < bump->n; ++i) {
|
|
if (balls[i].flag & 0x1) rbforest_rm (bump, i);
|
|
}
|
|
for (i = 0; i < bump->n; ++i) {
|
|
unsigned flag0 = balls[i].flag;
|
|
if (!(flag0 & 0x2)) continue;
|
|
for (j = 0; j < bump->n; ++j) {
|
|
unsigned flag1 = balls[j].flag;
|
|
if (i == j || (flag1 & 0x2 && i > j) || !(
|
|
(flag0 | flag1) & 0x1 && (*narrow) (cr, i, j, &score)
|
|
)) continue;
|
|
if (!rbforest_add (bump, i, j, score)) return 0;
|
|
}
|
|
}
|
|
return 1;
|
|
}
|
|
|
|
static void score_eval (struct crystal *cr) {
|
|
cr->scores[1] = tv_eval (cr);
|
|
if ((cr->scores[2] /= 2 * cr->nAtom[1]) > 1.0) cr->scores[2] = 1.0;
|
|
cr->scores[0] = cr->scales[2] * cr->scores[2] +
|
|
(1 - cr->scales[2]) * cr->scores[1];
|
|
}
|
|
|
|
int cryst_eval (struct crystal *cr) {
|
|
struct cr_pos *pos = cr->poss;
|
|
struct cr_atom *atom = cr->atoms;
|
|
struct cr_ball *ball = atom->scats->balls;
|
|
unsigned i;
|
|
for (i = 0; i < cr->nPos[0]; ++i, ++pos) if (pos->flag) {
|
|
pos->atom->xyz[pos->axis] = torus_pos (pos->pos[0]);
|
|
pos->atom->flag = 1;
|
|
}
|
|
for (i = 0; i < cr->nAtom[0]; ++i, ++atom) if (atom->flag) {
|
|
atom_eval (cr, atom);
|
|
atom->flag = 0;
|
|
}
|
|
spec_eval (cr);
|
|
if (cr->scales[2] != 0.0) {
|
|
if (!bump_eval (cr, &narrow_cryst)) return 0;
|
|
for (i = 0; i < cr->nAtom[1]; ++i, ++ball) ball->flag &= 0x2;
|
|
}
|
|
cr->scores[2] = cr->bump.score / cr->scales[1];
|
|
score_eval (cr);
|
|
return 1;
|
|
}
|
|
|
|
#ifdef BENCHMARK
|
|
static void pre_eval (struct crystal *cr, int sc) {
|
|
struct cr_peak *peak = cr->hills->peaks;
|
|
struct cr_pos *pos = cr->poss;
|
|
unsigned i;
|
|
for (i = 0; i < cr->nPos[0]; ++i, ++pos) {
|
|
pos->atom->xyz[pos->axis] = pos->pos[0];
|
|
}
|
|
for (i = 0; i < cr->nHill[1]; ++i, ++peak) memset
|
|
(peak->bsc, 0, 2 * sizeof (float));
|
|
for (i = 0; i < cr->nAtom[0]; ++i) atom_eval_bb (cr, cr->atoms + i, sc);
|
|
}
|
|
|
|
int cryst_eval_bb (struct crystal *cr, broad_f *broad, narrow_f *narrow) {
|
|
pre_eval (cr, 0);
|
|
if (broad && narrow) {
|
|
signed score;
|
|
if (!(*broad) (cr, narrow, &score)) return 0;
|
|
cr->scores[0] = score;
|
|
}
|
|
return 1;
|
|
}
|
|
|
|
int cryst_eval_bc (struct crystal *cr) {
|
|
signed score;
|
|
pre_eval (cr, 1);
|
|
spec_eval_bc (cr);
|
|
if (!bump_eval_bc (cr, &narrow_cryst, &score)) return 0;
|
|
cr->scores[2] = score / cr->scales[1];
|
|
score_eval (cr);
|
|
return 1;
|
|
}
|
|
#endif
|
|
|
|
unsigned cryst_dof (struct crystal const *cr) {
|
|
return cr->nPos[0];
|
|
}
|
|
|
|
void cryst_dump (struct crystal const *cr, uint32_t buf[]) {
|
|
for (unsigned i = 0; i < cr->nPos[0]; ++i) {
|
|
buf[i] = hftonl (cr->poss[i].pos[0]);
|
|
}
|
|
}
|
|
|
|
int cryst_load (struct crystal *cr, uint32_t const buf[]) {
|
|
struct cr_pos *poss = cr->poss;
|
|
for (unsigned i = 0; i < cr->nPos[0]; ++i) {
|
|
poss[i].pos[0] = nltohf (buf[i]);
|
|
poss[i].flag = 1;
|
|
if (
|
|
!isfinite (poss[i].pos[0]) ||
|
|
poss[i].pos[0] < poss[i].bound[0] ||
|
|
poss[i].pos[0] > poss[i].bound[1]
|
|
) return 0;
|
|
}
|
|
return 1;
|
|
}
|
|
|
|
float const *cryst_scores (struct crystal const *cr) {
|
|
return cr->scores;
|
|
}
|
|
|
|
#ifdef BENCHMARK
|
|
unsigned cryst_nball (struct crystal const *cr) {
|
|
return cr->nAtom[1];
|
|
}
|
|
|
|
void cryst_set (struct crystal *cr, float const xyzs[][3]) {
|
|
struct cr_pos *poss = cr->poss;
|
|
for (unsigned i = 0; i < cr->nPos[0]; ++i) {
|
|
poss[i].pos[0] = xyzs[poss[i].atom - cr->atoms][poss[i].axis];
|
|
poss[i].flag = 1;
|
|
}
|
|
}
|
|
|
|
void cryst_write (struct crystal const *cr) {
|
|
printf ("%g %g\n", cr->scores[1], cr->scores[2]);
|
|
for (unsigned i = 0;;) {
|
|
float *xyz = cr->atoms[i].scats->balls->xyz;
|
|
printf ("%g,%g,%g", xyz[0], xyz[1], xyz[2]);
|
|
if (++i < cr->nAtom[0]) printf (" ");
|
|
else {
|
|
printf ("\n");
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
|