/*
    Theseus - maximum likelihood superpositioning of macromolecular structures

    Copyright (C) 2004-2014 Douglas L. Theobald

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the:

    Free Software Foundation, Inc.,
    59 Temple Place, Suite 330,
    Boston, MA  02111-1307  USA

    -/_|:|_|_\-
*/

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <time.h>
#include <limits.h>
#include "Error.h"
#include "pdbMalloc.h"
#include "Cds.h"
#include "PDBCds.h"
#include "pdbStats.h"
#include "distfit.h"
#include "DLTmath.h"
#include "pdbSSM.h"


/* SSM is Edgar-speak for Structure Similarity Matrix file type, quote:

For a file format, I suggest the following.

The file type is called SSM (structure similarity matrices).

The file contains text.

The first line is formatted as follows:

#SSM<n>

where <n>=N is the number of structures.

Following the first line are N sequences in FASTA format. Presumably you will
use the standard amino acid alphabet here, but you can use any symbols you like
-- e.g. an X at each position just as a placeholder. MUSCLE will use a sequence
for use in the output alignment and to get the length so that the similarity
matrix dimension is known in advance.

Following the FASTA data are N(N-1)/2 similarity matrices. A similarity matrix
is formatted as follows:

#MATRIX<i>,<j>
<matrix data>
#ENDMATRIX

where <i> and <j> are "structure indexes". A structure index=1, 2 ... N as
defined by the order they appear in the FASTA data section of the file.

The matrices must appear in order of increasing i, then within a given i in
order of increasing j, including only matrices in which i < j (of course, the
matrix for i,j contains the same data as the matrix for j,i). So, for example,
if N=4 then the order is:

#MATRIX1,2
...
#MATRIX1,3
...
#MATRIX1,4
...
#MATRIX2,3
...
#MATRIX2,4
...
#MATRIX3,4
...

The matrix data is formatted as follows.

The data contains L_i lines, where L_i is the length of the i'th structure.

One line contains L_j floating-point values, where L_j is the length of the j'th
structure. A floating-point value may be formatted in any way readable by the C
language atof() function.

The first value starts in the first column of the line, values are separated by
exactly one space character.
*/


SSM
*SSMInit(void)
{
    SSM            *ssm = NULL;

    ssm = (SSM *) malloc(sizeof(SSM));
    if (ssm == NULL)
    {
        perror("\n  ERROR");
        fprintf(stderr,
                "\n  ERROR198: could not allocate memory in function SSMInit(). \n");
        PrintTheseusTag();
        exit(EXIT_FAILURE);
    }

    ssm->L = NULL;
    ssm->mat = NULL;

    return(ssm);
}


void
SSMAlloc(SSM *ssm, CdsArray *cdsA)
{
    int         i, j, k, cnum;

    ssm->n = cnum = cdsA->cnum;

    ssm->L = (int *) malloc(cnum * sizeof(int));

    for (i = 0; i < cnum; ++i)
    {
        ssm->L[i] = cdsA->cds[i]->aalen;
        //printf("\nlength: %d", ssm->L[i]);
    }

    ssm->mat = (double ***) malloc(0.5 * cnum * (cnum - 1) * sizeof(double **));

    k = 0;
    for (i = 0; i < cnum; ++i)
    {
        for (j = i+1; j < cnum; ++j)
        {
            ssm->mat[k] = MatAlloc(ssm->L[i], ssm->L[j]);
            ++k;
        }
    }

//    printf("[[%d][%d]]\n", k, ssm->n);fflush(NULL);
}


void
SSMDestroy(SSM **ssm_ptr)
{
    SSM *ssm = *ssm_ptr;
    int         i;

    for (i = 0; i < ssm->n * (ssm->n - 1) / 2; ++i)
        MatDestroy(&(ssm->mat[i]));

    free(ssm->mat);
    ssm->mat = NULL;

    free(ssm->L);
    ssm->L = NULL;

    free(ssm);
    *ssm_ptr = NULL;
}


void
SSMCalc(SSM *ssm, CdsArray *cdsA)
{
    int         i, j, k, m, n, p, q;
    int         cnum = cdsA->cnum;
    double     *lnvar = malloc(cdsA->vlen * sizeof(double));
    double     *invvar = malloc(cdsA->vlen * sizeof(double));

    for (i = 0; i < cdsA->vlen; ++i)
        lnvar[i] = log(cdsA->var[i]);

    for (i = 0; i < cdsA->vlen; ++i)
        invvar[i] = 1.0 / cdsA->var[i];

    k = 0;
    for (i = 0; i < cnum; ++i)
    {
        for (j = i+1; j < cnum; ++j)
        {
            for (m = p = 0; m < cdsA->vlen; ++m)
            {
                if (cdsA->cds[i]->nu[m])
                {
                    for (n = q = 0; n < cdsA->vlen; ++n)
                    {
                        if (cdsA->cds[j]->nu[n])
                        {
/*                           printf("[%d][%d][%d]:%e %e %e %e %e %e\n", k, p, q, */
/*                                  sqrt(SqrCdsDist(cdsA->cds[i], m, cdsA->avecds, m)), */
/*                                  sqrt(SqrCdsDist(cdsA->cds[i], m, cdsA->avecds, n)), */
/*                                  sqrt(SqrCdsDist(cdsA->cds[j], n, cdsA->avecds, m)), */
/*                                  sqrt(SqrCdsDist(cdsA->cds[j], n, cdsA->avecds, n)), */
/*                                  lnvar[i], lnvar[j]); */
/*                           fflush(NULL); */

                            ssm->mat[k][p][q] =
                            (invvar[m] * SqrCdsDist(cdsA->cds[i], m, cdsA->avecds, m) +
                             invvar[n] * SqrCdsDist(cdsA->cds[i], m, cdsA->avecds, n) +
                             invvar[m] * SqrCdsDist(cdsA->cds[j], n, cdsA->avecds, m) +
                             invvar[n] * SqrCdsDist(cdsA->cds[j], n, cdsA->avecds, n) +
                             2.0 * (lnvar[i] + lnvar[j])) * -0.25;

//                           printf("[%d][%d][%d]:%e\n", k, p, q, ssm->mat[k][p][q]);
//                           fflush(NULL);

                            q++;
                        }
                    }
                    //printf("\nq:%d", q);

                    p++;
                }
            }
            //printf("\np:%d", p);
            k++;
        }
    }

    free(lnvar);
    free(invvar);
}


void
WriteSSM(SSM *ssm)
{
    int         i, j, k, m, n;
    FILE       *fp = NULL;

    fp = fopen("ssm.txt", "w");

    fprintf(fp, "#SSM%d\n", ssm->n);

    k = 0;
    for (i = 0; i < ssm->n; ++i)
    {
        for (j = i+1; j < ssm->n; ++j)
        {
            fprintf(fp, "#MATRIX%d,%d\n", i+1, j+1);

            for (m = 0; m < ssm->L[i]; ++m)
            {
                for (n = 0; n < ssm->L[j]; ++n)
                {
                    fprintf(fp, "% 14.2f ", ssm->mat[k][m][n]);
                }

                fprintf(fp, "\n");
            }

            fprintf(fp, "ENDMATRIX\n");
            k++;
        }
    }

    fprintf(fp, "\n");

    fclose(fp);
}