/*
    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 <stdio.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "CovMat.h"
#include "DLTutils.h"
#include "DLTmath.h"
#include "Error.h"
#include "PCAstats.h"
#include "pdbIO.h"
#include "pdbMalloc.h"
#include "Cds.h"
#include "PDBCds.h"
#include "pdbStats.h"
#include "pdbUtils.h"


void
CalcPCA(CdsArray *cdsA)
{
    int            i, j;
    int            vlen = (int) cdsA->vlen;
    int            pcanum;
    double       **CovMat = cdsA->CovMat;
    double         sum, runsum;
    PDBCds        *pdbave = NULL;
    char           pcafile_name[256];
    FILE          *pcavecs_fp = NULL, *pcastats_fp = NULL;
    double         biggest, bstick;
    char           aster;
    char          *cov_name = NULL, *cor_name = NULL, *pcvecs_name = NULL, *pcstats_name = NULL;


    pdbave = cdsA->pdbA->avecds;
    pcanum = algo->pca;

    if (pcanum > cdsA->vlen)
        pcanum = algo->pca = cdsA->vlen;

    if (pcanum > cdsA->cnum - 1)
        pcanum = algo->pca = cdsA->cnum - 1;

    cov_name = mystrcat(algo->rootname, "_cov.mat");
    PrintCovMatGnuPlot((const double **) CovMat, vlen, cov_name);

    /* convert it to a correlation matrix */
    if (algo->cormat)
    {
        CovMat2CorMat(CovMat, vlen);
        cor_name = mystrcat(algo->rootname, "_cor.mat");
        PrintCovMatGnuPlot((const double **) CovMat, vlen, cor_name);
    }

//    write_C_mat((const double **) CovMat, vlen, 8, 0);

    /* find the total variance */
    sum = 0.0;
    for (i = 0; i < vlen; ++i)
       sum += CovMat[i][i];

    cdsA->pcamat = MatAlloc(vlen, vlen);
    cdsA->pcavals = malloc(vlen * sizeof(double));

    EigenGSLDest(CovMat, vlen, cdsA->pcavals, cdsA->pcamat, 1);
    MatTransIp(cdsA->pcamat, vlen);
    //MatPrint(cdsA->pcamat, vlen);

    /* LAPACK DSYEVR() computes selected eigenvalues, and optionally, eigenvectors of a
       real symmetric matrix.  Find all eigenvalues (w[]) and eigenvectors (mat[][]).
       The pcanum eigenvalues are in the first pcanum elements of the w[] vector,
       ordered smallest to biggest.  Weird and horrible, but true (and in fact
       makes sense if you think about it). */
    //dsyevr_dest(CovMat, vlen, lower, upper, cdsA->pcavals, cdsA->pcamat, 1e-8);

    //MatPrint(cdsA->pcamat, vlen);

    pcvecs_name = mystrcat(algo->rootname, "_pcvecs.txt");
    pcstats_name = mystrcat(algo->rootname, "_pcstats.txt");
    pcavecs_fp = myfopen(pcvecs_name, "w");
    pcastats_fp = myfopen(pcstats_name, "w");

    if (pcavecs_fp == NULL || pcastats_fp == NULL)
    {
        fprintf(stderr, "\n  ERROR: Could not open PCA files \n");
        PrintTheseusTag();
        exit(EXIT_FAILURE);
    }

    runsum = 0.0;
    fprintf(pcastats_fp, "eigenv   bstick     raw    raw_%%  cumul_%% \n");

    for (i = 0; i < pcanum; ++i)
    {
        bstick = 0.0;
        for (j = i+1; j <= vlen; ++j)
            bstick += (1.0 / (double) j);

        if (bstick < cdsA->pcavals[i] && cdsA->pcavals[i] > 1.0)
            aster = '*';
        else
            aster = ' ';

        runsum += cdsA->pcavals[i];
        fprintf(pcastats_fp, "%-6i %8.3f %8.3f %8.3f %8.3f %c\n",
                i+1, bstick, cdsA->pcavals[i],
                cdsA->pcavals[i] * 100.0 / sum, runsum * 100.0 / sum, aster);
    }
    fputc('\n', pcastats_fp);

    fprintf(pcavecs_fp, "atom ");
    for (j = 0; j < pcanum; ++j)
        fprintf(pcavecs_fp, "     %3d ", j+1);
    fputc('\n', pcavecs_fp);

    for (i = 0; i < vlen; ++i)
    {
        fprintf(pcavecs_fp, "%-4d ", i+1);

        for (j = 0; j < pcanum; ++j)
            fprintf(pcavecs_fp, "%8.3f ", sqrt(cdsA->pcavals[j]) * cdsA->pcamat[j][i]);

        fputc('\n', pcavecs_fp);
    }

    CopyCds2PDB(pdbave, cdsA->avecds);

    for (i = 0; i < pcanum; ++i)
    {
        //printf("\neigenval %d: %g\n", i, cdsA->pcavals[i]);
        /* find largest absolute value in the eigenvector PCA */
        biggest = -DBL_MAX;
        for (j = 0; j < vlen; ++j)
            if (biggest < fabs(cdsA->pcamat[i][j]))
                biggest = fabs(cdsA->pcamat[i][j]);

        /* rescale (for rasmol really) so that the largest eigenvalue component
           is = 99.99, i.e. the largest value allowable in the b-value column
           of a PDB file */
        /* biggest = 1.0; */
        for (j = 0; j < vlen; ++j)
            cdsA->pdbA->avecds->tempFactor[j] = cdsA->pcamat[i][j] * (99.99 / biggest);

        sprintf(pcafile_name, "%s_pc%d_ave.pdb", algo->rootname, i+1);

        //strcpy(pcafile_name, mystrcat(algo->rootname, "_pc"));
        //pcafile_name[11] = '\0';
        //strncat(pcafile_name, itoa(i+1, &numstring[0], 10), 5);
        //strncat(pcafile_name, "_ave.pdb", 8);

        WriteAvePDBCdsFile(cdsA->pdbA, pcafile_name);
    }

    if (pcanum == vlen)
    {
        char *pcvecs_mat_name = mystrcat(algo->rootname, "_pcvecs.mat");

        for (i = 0; i < vlen; ++i)
            for (j = 0; j < pcanum; ++j)
                cdsA->pcamat[j][i] *= sqrt(cdsA->pcavals[j]);

        PrintCovMatGnuPlot((const double **) cdsA->pcamat, vlen, pcvecs_mat_name);
        free(pcvecs_mat_name);
    }

    fclose(pcastats_fp);
    fclose(pcavecs_fp);

    if (cov_name != NULL)
        free(cov_name);

    if (cor_name != NULL)
        free(cor_name);

    if (pcvecs_name != NULL)
        free(pcvecs_name);

    if (pcstats_name != NULL)
        free(pcstats_name);
}


void
Calc3NPCA(CdsArray *cdsA)
{
    int                      i, j;
    int                      vlen = (int) 3 * cdsA->vlen;
    double                 **mat = NULL;
    int                      pcanum;
    double                 **evecs = NULL, *evals = NULL;
    double                   sum, runsum;
    PDBCds               *pdbave = NULL;
    char                     pcafile_name[256];
    FILE                    *pcavecs_fp = NULL, *pcastats_fp = NULL;
    double                   biggest, bstick;
    char                     aster;

    pdbave = cdsA->pdbA->avecds;
    mat = MatAlloc(vlen, vlen);

    if (algo->pca > cdsA->cnum - 1)
        pcanum = algo->pca = cdsA->cnum - 1;
    else
        pcanum = algo->pca;

    /* copy over the covariance matrix */
    memcpy(mat[0], cdsA->FullCovMat[0], vlen * vlen * sizeof(double));
    /* MatPrint(cdsA->FullCovMat, vlen); */
    /* fflush(NULL); */

    /* convert it to a correlation matrix */
    if (algo->cormat)
        CovMat2CorMat(mat, vlen);

    /* find the total variance */
    sum = 0.0;
    for (i = 0; i < vlen; ++i)
       sum += mat[i][i];

    evecs = MatAlloc(vlen, vlen);
    evals = malloc(vlen * sizeof(double));

    /* LAPACK DSYEVR() computes selected eigenvalues, and optionally, eigenvectors of a
       real symmetric matrix.  Find all eigenvalues (w[]) and eigenvectors (mat[][]).
       The pcanum eigenvalues are in the first pcanum elements of the w[] vector,
       ordered smallest to biggest.  Weird and horrible, but true (and in fact
       makes sense if you think about it). */
    //dsyevr_opt_dest(mat, vlen, lower, upper, evals, evecs, 1e-8);

    EigenGSLDest(mat, vlen, evals, evecs, 1);
    MatTransIp(cdsA->pcamat, vlen);

    PrintCovMatGnuPlot((const double **) evecs, vlen, "evecs.mat");

    pcavecs_fp = fopen("pcavecs.txt", "w");
    pcastats_fp = fopen("pcastats.txt", "w");
    if (pcavecs_fp == NULL || pcastats_fp == NULL)
    {
        fprintf(stderr, "\n  ERROR1000: Could not open PCA files \n");
        PrintTheseusTag();
        exit(EXIT_FAILURE);
    }

    runsum = 0.0;
    fprintf(pcastats_fp, "eigenv   bstick     raw    raw_%%  cumul_%% \n");

    for (i = 0; i < pcanum; ++i)
    {
        bstick = 0.0;
        for (j = i+1; j <= vlen; ++j)
            bstick += (1.0 / (double) j);

        if (bstick < evals[i] && evals[i] > 1.0)
            aster = '*';
        else
            aster = ' ';

        runsum += cdsA->pcavals[i];
        fprintf(pcastats_fp, "%-6i %8.3f %8.3f %8.3f %8.3f %c\n",
                i+1, bstick, evals[i],
                evals[i] * 100.0 / sum, runsum * 100.0 / sum, aster);
    }
    fputc('\n', pcastats_fp);

    fprintf(pcavecs_fp, "atom ");
    for (j = 0; j < pcanum; ++j)
        fprintf(pcavecs_fp, "     %3d ", j+1);
    fputc('\n', pcavecs_fp);

    for (i = 0; i < vlen; ++i)
    {
        fprintf(pcavecs_fp, "%-4d ", i+1);

        for (j = 0; j < pcanum; ++j)
            fprintf(pcavecs_fp, "%8.3f ", sqrt(evals[j]) * evecs[j][i]);

        fputc('\n', pcavecs_fp);
    }

    CopyCds2PDB(pdbave, cdsA->avecds);

    for (i = 0; i < pcanum; ++i)
    {
        /* find largest absolute value in the eigenvector PCA */
        biggest = -DBL_MAX;
        for (j = 0; j < vlen; ++j)
            if (biggest < fabs(evecs[i][j]))
                biggest = fabs(evecs[i][j]);

        /* rescale (for rasmol really) so that the largest eigenvalue component
           is = 99.99, i.e. the largest value allowable in the b-value column
           of a PDB file */
        for (j = 0; j < vlen; ++j)
            cdsA->pdbA->avecds->tempFactor[j] = evecs[i][j] * (99.99 / biggest);

        sprintf(pcafile_name, "%s_pc%d", algo->rootname, i+1);

        //strncpy(pcafile_name, mystrcat(algo->rootname, "_pc"), 11);
        //pcafile_name[11] = '\0';
        //strncat(pcafile_name, itoa(i+1, &numstring[0], 10), 5);

        WriteAvePDBCdsFile(cdsA->pdbA, pcafile_name);
    }

    cdsA->pcamat = evecs; /* DLT debug -- this should be copied or dealt with better */
    cdsA->pcavals = evals;

    fclose(pcastats_fp);
    fclose(pcavecs_fp);

    MatDestroy(&mat);
}


void
WritePCAFile(PDBCdsArray *parray, CdsArray *cdsA, const char *outfile_root)
{
    FILE           *pdbfile = NULL;
    char            pcafile_name[256];
    int             i, j, k, m;
    double          biggest = -DBL_MAX;
    const double  **mat = (const double **) cdsA->pcamat;
    const int       cvlen = cdsA->vlen;
    const int       pvlen = parray->vlen;
    double          tempFactor;
    //char            numstring[5];
    char            covcor_str[16] = "correlation";

    /* find largest absolute value in the eigenvector PCA */
    for (i = 0; i < algo->pca; ++i)
    {
        biggest = -DBL_MAX;
        for (j = 0; j < cvlen; ++j)
        {
            /* printf("\n%3d %3d % f", i, j, mat[i][j]); */
            if (biggest < fabs(mat[algo->pca - 1 - i][j]))
                biggest = fabs(mat[algo->pca - 1 - i][j]);
        }
/*     } */
        biggest = 99.99 / biggest;

/*     for (i = 0; i < algo->pca; ++i) */
/*     { */
        /* rescale (for rasmol really) so that the largest eigenvector component
           is = 99.99, i.e. the largest value allowable in the b-value column
           of a PDB file */
       if (algo->atoms == 0)
       {
            m = 0;
            for (j = 0; j < cvlen; ++j)
            {
                /* skip inital PDBCds that may have been selected out */ /* DLT debug fix */
                while (/* strncmp(cdsA->cds[0]->resName[j], parray->cds[0]->resName[m], 3) != 0 || */
                       cdsA->cds[0]->chainID[j] != parray->cds[0]->chainID[m] ||
                       cdsA->cds[0]->resSeq[j] != parray->cds[0]->resSeq[m])
                {
                    ++m;

                    if (m >= pvlen)
                        break;
                }

                /* while they match, set the B-factor according to the given PC */
                while (/* strncmp(cdsA->cds[0]->resName[j], parray->cds[0]->resName[m], 3) == 0 && */
                       cdsA->cds[0]->chainID[j] == parray->cds[0]->chainID[m] &&
                       cdsA->cds[0]->resSeq[j] == parray->cds[0]->resSeq[m])
                {
                    tempFactor = mat[algo->pca - 1 - i][j] * biggest;
                    /* printf("\n%4d %4d % f", j, m, tempFactor); */
                    for (k = 0; k < parray->cnum; ++k)
                        parray->cds[k]->tempFactor[m] = tempFactor;

                    ++m;

                    if (m >= pvlen)
                        break;
                }

                if (m >= pvlen)
                    break;
            }
        }
        else
        {
            for (j = 0; j < cvlen; ++j)
            {
                tempFactor = mat[algo->pca - 1 - i][j] * biggest;

                for (k = 0; k < parray->cnum; ++k)
                    parray->cds[k]->tempFactor[j] = tempFactor;
            }
        }

        sprintf(pcafile_name, "%s_pc%d.pdb", outfile_root, i+1);

/*         strncpy(pcafile_name, outfile_root, strlen(outfile_root)); */
/*         pcafile_name[strlen(outfile_root)] = '\0'; */
/*         strncat(pcafile_name, "_pc", 4); */
/*         strncat(pcafile_name, itoa(i+1, &numstring[0], 10), 5); */
/*         strcat(pcafile_name, ".pdb"); */

        pdbfile = fopen(pcafile_name, "w");
        if (pdbfile ==NULL)
        {
            perror("\n  ERROR");
            fprintf(stderr,
                    "\n  ERROR99: could not open file '%s' for writing. \n", pcafile_name);
            PrintTheseusTag();
            exit(EXIT_FAILURE);
        }

        if (algo->cormat == 0)
            strncpy(covcor_str, "covariance", 10);

        fprintf(pdbfile, "REMARK ===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-=\n");
        fprintf(pdbfile, "REMARK + File made by Douglas Theobald's THESEUS program\n");
        fprintf(pdbfile, "REMARK + Multiple maximum likelihood superpositioning\n");
        fprintf(pdbfile, "REMARK + Principal component %3d of %s matrix in B-factor column\n", i+1, covcor_str);
        fprintf(pdbfile, "REMARK + All B-factors scaled by %12.3f\n", biggest);
        fprintf(pdbfile, "REMARK + dtheobald@brandeis.edu\n");
        fprintf(pdbfile, "REMARK =-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===\n");

        for (j = 0; j < parray->cnum; ++j)
        {
            fprintf(pdbfile, "MODEL %8d\n", j+1);
            PrintPDBCds(pdbfile, parray->cds[j]);
            fprintf(pdbfile, "ENDMDL\n");
        }
        fprintf(pdbfile, "END\n");

        fclose(pdbfile);
    }

    PCADestroy(cdsA);
}


/* Writes out PDB format files with the PC eigenvector in the coordinate field.
   Used by CNS/XPLOR in my morph script to incrementally add the PC to the
   median PDB model. */
void
WritePCAMorphFile(PDBCdsArray *parray, CdsArray *cdsA, const char *outfile_root)
{
    FILE           *pdbfile = NULL;
    char            pcafile_name[256];
    int             i, j, m, pcanum;
    double        **vecs = cdsA->pcamat;
    double         *vals = cdsA->pcavals;
    const int       cvlen = 3 * cdsA->vlen;
    const int       pvlen = parray->vlen;
    PDBCds         *pcacds = NULL;
    char            covcor_str[16] = "correlation";

    pcacds = PDBCdsInit();
    PDBCdsAlloc(pcacds, pvlen);
    PDBCdsCopyAll(pcacds, parray->cds[0]);

    memset(pcacds->x, 0, pvlen * sizeof(double));
    memset(pcacds->y, 0, pvlen * sizeof(double));
    memset(pcacds->z, 0, pvlen * sizeof(double));

    if (algo->pca > cdsA->cnum - 1)
        pcanum = algo->pca = cdsA->cnum - 1;
    else
        pcanum = algo->pca;

CalcFullCovMat(cdsA);

    /* Multiply each PCA by the sqrt of the corresponding eigenvalue.
       If correlation matrix was used, we need to get back into std deviation space,
       so multiply by the sqrt of the corresponding variance */
    if (algo->cormat)
    {
        for (i = 0; i < pcanum; ++i)
            for (j = 0; j < cvlen; ++j)
                /* vecs[i][j] *= sqrt(cdsA->FullCovMat[j][j] * vals[i]); */
                vecs[i][j] *= sqrt(cdsA->FullCovMat[j][j]);
    }
    else if (algo->cormat == 0)
    {
        for (i = 0; i < pcanum; ++i)
            for (j = 0; j < cvlen; ++j)
                vecs[i][j] *= sqrt(vals[i]);
    }

    for (i = 0; i < pcanum; ++i)
    {
        m = 0;
        for (j = 0; j < cvlen; j += 3)
        {
            /* skip inital PDBCds that may have been selected out */
            while (strncmp(cdsA->cds[0]->resName[j/3], parray->cds[0]->resName[m], 3) != 0 ||
                   cdsA->cds[0]->chainID[j/3] != parray->cds[0]->chainID[m] ||
                   cdsA->cds[0]->resSeq[j/3] != parray->cds[0]->resSeq[m])
            {
                ++m;
            }

            /* while they match */
            while (strncmp(cdsA->cds[0]->resName[j/3], parray->cds[0]->resName[m], 3) == 0 &&
                   cdsA->cds[0]->chainID[j/3] == parray->cds[0]->chainID[m] &&
                   cdsA->cds[0]->resSeq[j/3] == parray->cds[0]->resSeq[m])
            {
                pcacds->x[m] = vecs[pcanum - 1 - i][j+0];
                pcacds->y[m] = vecs[pcanum - 1 - i][j+1];
                pcacds->z[m] = vecs[pcanum - 1 - i][j+2];

                ++m;

                if (m >= pvlen)
                    break;
            }

            if (m >= pvlen)
                break;
        }

        sprintf(pcafile_name, "%s_pca%d_morph.pdb", outfile_root, i+1);

/*         strncpy(pcafile_name, outfile_root, strlen(outfile_root)); */
/*         pcafile_name[strlen(outfile_root)] = '\0'; */
/*         strncat(pcafile_name, "_pca", 4); */
/*         strncat(pcafile_name, itoa(i+1, &numstring[0], 10), 5); */
/*         strncat(pcafile_name, "_morph", 6); */
/*         strcat(pcafile_name, ".pdb"); */

        pdbfile = fopen(pcafile_name, "w");
        if (pdbfile == NULL)
        {
            perror("\n  ERROR");
            fprintf(stderr,
                    "\n  ERROR99: could not open file '%s' for writing. \n", pcafile_name);
            PrintTheseusTag();
            exit(EXIT_FAILURE);
        }

        if (algo->cormat == 0)
            strncpy(covcor_str, "covariance", 10);

        fprintf(pdbfile, "REMARK ===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-=\n");
        fprintf(pdbfile, "REMARK + File made by Douglas Theobald's THESEUS program \n");
        fprintf(pdbfile, "REMARK + Multiple maximum likelihood superpositioning \n");
        fprintf(pdbfile, "REMARK + Principal component #%d of %s matrix, one SD in cds fields \n", i+1, covcor_str);
        fprintf(pdbfile, "REMARK + dtheobald@brandeis.edu \n");
        fprintf(pdbfile, "REMARK =-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===\n");
        PrintPDBCds(pdbfile, pcacds);
        fprintf(pdbfile, "END\n");

        fclose(pdbfile);
    }

    PDBCdsDestroy(&pcacds);
    PCADestroy(cdsA);
}


void
WritePCAProjections(PDBCdsArray *parray, CdsArray *cdsA, const char *outfile_root)
{
    FILE           *pdbfile = NULL;
    char            pcafile_name[256];
    int             i, j, m, pcanum;
    double          c;
    double        **vecs = cdsA->pcamat;
    double         *vals = cdsA->pcavals;
    const int       cvlen = 3 * cdsA->vlen;
    const int       pvlen = parray->vlen;
    PDBCds         *pcacds = NULL;
    //char            numstring[5], tmpstring[5];

    pcacds = PDBCdsInit();
    PDBCdsAlloc(pcacds, pvlen);
    PDBCdsCopyAll(pcacds, parray->cds[0]);

    if (algo->pca > cdsA->cnum - 1)
        pcanum = algo->pca = cdsA->cnum - 1;
    else
        pcanum = algo->pca;

CalcFullCovMat(cdsA);

    /* Multiply each PCA by the sqrt of the corresponding eigenvalue.
       If correlation matrix was used, we need to get back into std deviation space,
       so multiply by the sqrt of the corresponding variance */
    if (algo->cormat)
    {
        for (i = 0; i < pcanum; ++i)
            for (j = 0; j < cvlen; ++j)
                vecs[i][j] *= sqrt(cdsA->FullCovMat[j][j] * vals[i]);
    }
    else if (algo->cormat == 0)
    {
        for (i = 0; i < pcanum; ++i)
            for (j = 0; j < cvlen; ++j)
                vecs[i][j] *= sqrt(vals[i]);
    }

/*     for (j = 1; j <= cvlen; j++) */
/*         printf("\n%f", vals[cvlen - j]); */

/*     for (j = 0; j < cvlen; j++) */
/*         printf("\n%f", vecs[cvlen-1][j]); */

    for (i = 0; i < pcanum; ++i)
    {
        for (c = -3.0; c <= 3.0; c += 0.2)
        {
            m = 0;
            for (j = 0; j < cvlen; j += 3)
            {
                /* skip inital PDBCds that may have been selected out */
                while (strncmp(cdsA->cds[0]->resName[j/3], parray->cds[0]->resName[m], 3) ||
                       cdsA->cds[0]->chainID[j/3] != parray->cds[0]->chainID[m] ||
                       cdsA->cds[0]->resSeq[j/3] != parray->cds[0]->resSeq[m])
                {
                    ++m;
                }

                /* while they match */
                while (strncmp(cdsA->cds[0]->resName[j/3], parray->cds[0]->resName[m], 3) == 0 &&
                       cdsA->cds[0]->chainID[j/3] == parray->cds[0]->chainID[m] &&
                       cdsA->cds[0]->resSeq[j/3] == parray->cds[0]->resSeq[m])
                {
                    pcacds->x[m] = parray->cds[cdsA->cnum/2]->x[m] +
                                      c * vecs[pcanum - 1 - i][j+0];
                    pcacds->y[m] = parray->cds[cdsA->cnum/2]->y[m] +
                                      c * vecs[pcanum - 1 - i][j+1];
                    pcacds->z[m] = parray->cds[cdsA->cnum/2]->z[m] +
                                      c * vecs[pcanum - 1 - i][j+2];

                    ++m;

                    if (m >= pvlen)
                        break;
                }

                if (m >= pvlen)
                    break;
            }

            sprintf(pcafile_name, "%s_pca%d_%+3.1f.pdb", outfile_root, i+1, c);

/*             strncpy(pcafile_name, outfile_root, strlen(outfile_root)); */
/*             pcafile_name[strlen(outfile_root)] = '\0'; */
/*             strncat(pcafile_name, "_pca", 4); */
/*             strncat(pcafile_name, itoa(i+1, &numstring[0], 10), 5); */
/*             sprintf(tmpstring, "_%+3.1f", c); */
/*             strncat(pcafile_name, tmpstring, 5); */
/*             strcat(pcafile_name, ".pdb"); */

            pdbfile = fopen(pcafile_name, "w");
            if (pdbfile == NULL)
            {
                perror("\n  ERROR");
                fprintf(stderr,
                        "\n  ERROR99: could not open file '%s' for writing. \n", pcafile_name);
                PrintTheseusTag();
                exit(EXIT_FAILURE);
            }

            fprintf(pdbfile, "REMARK ===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-=\n");
            fprintf(pdbfile, "REMARK + File made by Douglas Theobald's THESEUS program \n");
            fprintf(pdbfile, "REMARK + Multiple maximum likelihood superpositioning \n");
            fprintf(pdbfile, "REMARK + Principal component %d of correlation matrix in B-factor column \n", i+1);
            fprintf(pdbfile, "REMARK + dtheobald@brandeis.edu \n");
            fprintf(pdbfile, "REMARK =-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===\n");
            PrintPDBCds(pdbfile, pcacds);
            fprintf(pdbfile, "END\n");

            fclose(pdbfile);
        }
    }

    PDBCdsDestroy(&pcacds);
    PCADestroy(cdsA);
}