/*
 * 20 May 2002
 * $Id: geo_gap.c,v 1.1 2007/09/28 16:57:05 mmundry Exp $
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "coord.h"
#include "geo_gap.h"
#include "mprintf.h"

/* ---------------- Constants --------------------------------
 */
static const float R_CN = 1.32;  /* Ideal bond distance */
static const float TRIG = 2.0;   /* Minimum dist to worry about. Below
                             this, there is no gap */
/* ---------------- coord_geo_gap -----------------------------
 * Implement gap penalties for a sequence to structure alignment,
 * based on geometry of model.
 * c is the coordinate set we are looking at.
 * scale is a coefficient which we multiply our final result by.
 * max is the maximum distance. This is a cap. If the calculated
 * distance is more than max, treat it as max.
 * In the code (hardwired), if the C..N distance is less than
 * TRIG (2.0 Angstrom), do not treat it as a gap.
 * num_gap is a pointer into which we return the number of gaps.
 * We want to give the caller several numbers back
 *   gap penalty using quadratic penalties
 *               using linear penalties
 *               using crazy logistic function (rising part)
 *   the number of gaps
 * To see what our logistic function looks like, feed this to gnuplot:
     c1 = 2;
     c2 = 4;
     f(x) = 1. / (1. + c2 * exp (c1 - ( (x < 1.8) ? 0 : x)));
     set xrange [0:15]
     plot f(x);
 * More notes..
 * Some example gap sizes from real structures
 *  gap in residues / gap in Angstrom
 *  1  / 3 to 3.5
 *  2  / 4.2 to about 7
 *  3  / Anywhere from 4 to 9
 */
int
coord_geo_gap (struct coord *c,
               float *quadratic, float *linear, float *logistic,
               unsigned int *num_gap, const float scale, const float max)
{
    const char *this_sub = "coord_geo_gap";
    struct RPoint *rn = c->rp_n + 1;
    struct RPoint *rc = c->rp_c;
    struct RPoint *rclast;
    unsigned ngap = 0;
    float quad    = 0.0;
    float lin     = 0.0;
    float logi    = 0.0;
    const float TRIG2 = TRIG * TRIG;  /* Don't worry. This is optimised away */
    const float max2 = max * max;
    const float k1 = 0.0,
                k2 = 1;
    if (c->size < 2) {
        err_printf (this_sub, "Coordinates too small ( < 2 residues)\n");
        return EXIT_FAILURE;
    }
    rclast = rc + c->size - 1;
    for ( ; rc < rclast; rc++, rn++) {
        float x = rc->x - rn->x;
        float y = rc->y - rn->y;
        float z = rc->z - rn->z;
        float d2 = (x * x) + (y * y) + (z * z);
        float d;
        if (d2 < TRIG2)       /* If dist is less than our threshold */
            continue;         /* there is no penalty */

        if (d2 > max2) {      /* If distance is too big, then just use */
            d2 = max2;        /* the maximum we were given. */
            d  = max;
        } else {
            d = sqrt (d2);    /* Real distance between atoms */
        }
        ngap++;             /* Whatever happens, we track the number of gaps */
        d -= R_CN;                                  /* This is our violation */
        logi += 1 / (1 + k2 * exp (k1 - d/2 ));   /* logistic function */
        lin  += d;
        quad += (d * d);
    }

    *num_gap   = ngap;            /* These are really the return values */
    *quadratic = quad * scale;
    *linear    = lin  * scale;
    *logistic  = logi * scale;
    return EXIT_SUCCESS;
}