/*-------------------------------------------------------------

Copyright (C) 2000 Rolf Backofen, Peter Clote. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/



/* Include section */
#include "general.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>


/* function declarations */
/* calc square distance-matrix from a vector of polar requirements */
void calc_matrix(float PC[NUM_ACIDS], float matrix[NUM_ACIDS][NUM_ACIDS]);
void print_pc_code(float PC[NUM_ACIDS],gcode code);
int printblockstruct( gcode block_code, int *fixed );

/* reads a vector    */
int read_non_fixed_vector(char *filename, float PC[NUM_ACIDS], 
                          int fixed[NUM_ACIDS], float V[5]);

/* sets the value of matrix nij (= digulio in Peters version)
   calculates and returns MS (including sumnij). */
Fivetuple setnij_and_ms(gcode code,
                        float matrix[NUM_ACIDS][NUM_ACIDS],
                        int nij[NUM_ACIDS][NUM_ACIDS]);


/********************************************************************** 

            General Routines for measuring Fault Tolerance 

***********************************************************************/


/* sets the value of matrix nij calculates and
   returns MS (including sumnij).

   after return:

   nij gives how often we go from block i to block j by one mutation of
   the triplet code.
   
   MS.sumnij is the sum over all nij values

   MS.whole is mean square distance of polar requirements of
   blocknumbers, where the triplet code changes by one base

   MS.one/MS.two/MS.three is mean square distance ..., where
   the triplet changes in first/second/third base.
*/
Fivetuple setnij_and_ms(gcode code,
     float matrix[NUM_ACIDS][NUM_ACIDS],
     int nij[NUM_ACIDS][NUM_ACIDS])
{
  int i, j, k, m;  /* loop variables for genetic code*/
  int sumnij = 0; /* error tolerance of the given code */
  int num1 = 0, num2 = 0, num3 = 0;
  Fivetuple MS;
       
  MS.whole = MS.one = MS.two = MS.three = 0.0;

  /* initialize matrix for DiGiulio coefficients */
  for (i=0;i<NUM_ACIDS;i++)
    for (j=0;j<NUM_ACIDS;j++)
       nij[i][j] = 0;

  /* compute nij coefficients */
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      for(k = 0; k < 4; k++)
        {
          if (code[i][j][k] != STOP )
            for (m = 0; m < 4; m++)
              {
                if  ( i != m ) 
                  {
                    if  (code[m][j][k] != STOP ) {
                      MS.one += matrix[code[i][j][k]][code[m][j][k]]; 
                      num1++;
                      nij[code[i][j][k]][code[m][j][k]]++; 
                      sumnij++;
                    }
                  }
                if  ( j != m )
                  { 
                    if  (code[i][m][k] != STOP ) {
                     MS.two += matrix[code[i][j][k]][code[i][m][k]]; 
                     num2++;
                     nij[code[i][j][k]][code[i][m][k]]++; 
                     sumnij++;
                    }
                  }
                if   ( k != m ) 
                  {
                    if  (code[i][j][m] != STOP) {
                      MS.three += matrix[code[i][j][k]][code[i][j][m]]; 
                      num3++;
                      nij[code[i][j][k]][code[i][j][m]]++; 
                      sumnij++;
                    }
                  }
              }
        }
  
  MS.whole = (MS.one + MS.two + MS.three)/(num1 + num2 + num3); 
  MS.one = MS.one/num1;
  MS.two = MS.two/num2;
  MS.three = MS.three/num3;

  MS.sumnij = sumnij;

  return MS;
}

/**************
  main function
***************/

int main (int argc, char* argv[])
{
  int i,j,k,num; /* loop variables */

  /* nij-matrix */
  int nij[NUM_ACIDS][NUM_ACIDS]; /* computes N_{i,j} */
  int s, sumnij; 
  
  /* polar requirements vector for blocks,
     corresponding square distance-matrix  */
  float PC[NUM_ACIDS];
  float matrix[NUM_ACIDS][NUM_ACIDS];

  /* mean standard deviation */
  Fivetuple MS;

  /* matrix and vector for lagrange system */
  /* this is for every i: sum_{j} N_{i,j} for i,j not fixed. */
  int nijSums[NUM_ACIDS]; 

  float A[16][16];  /* 16 x 16 linear equation matrix */
  float B[16];    /* 16 x 1 vector */

  /* this relates line in matrix A to block numbers,
     the values in vector fixed do not appear in the
     ascending values of a_blocknums */
  int a_blocknums[15] =
  { 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 14, 15, 16, 17, 18};

  /* fixed blocks and their values  (currently genetic code)*/
  float V[5]       = { 5.0, 6.6,  10.0, 10.1,  7.9 };
  int fixednums[5] = {   0,   6,    12,   13,   19 };
  /* will be calculated: 1 if block number is fixed */
  int fixed[NUM_ACIDS];
  int nextfixed;

  /**********************************************
   generated fixed vector

   the array fixednums gives the indices of entries in fixed 
   that are set to 1 by following code
  ************************************************/
  nextfixed = 0;
  for (i=0;i<NUM_ACIDS;i++) 
    {
      if (i == fixednums[nextfixed])
        {
          fixed[i] = 1;
          nextfixed++;
        }
      else
        fixed[i] = 0;
    }
  
  printf("\n    Block Structure   \n");
  printblockstruct( block_code,fixed );

  /**********************************************
  error checking routines.
  Read in file for similarity matrix: eg. wac.mat.
  ************************************************/

  if ((argc != 2) && (argc != 7))
    {
      fprintf(stderr, "%s: Insufficient arguments!\n", argv[0]);
      fprintf(stderr, "usage: %s <polarity_vector_file> or %s <polarity_vector_file> v1 v2 v3 v4 v5 \n", argv[0], argv[0]);
      exit(1);
    }
  else if (argc == 7)
    for (i=0;i<5;i++)
      V[i] = strtod(argv[i+2],NULL);

  
  /* read polarrequirements vector  */
  read_non_fixed_vector(argv[1],PC,fixed,V);

  /* compute square distance-matrix             */
  calc_matrix(PC, matrix);

 /* print polar-requirement vector including fixed values  */
  printf("\n\n PC Vector  PC genetic Code\n");
  for (i=0;i<NUM_ACIDS;i++)
    {
      printf("%8.2f      %8.2f\n",PC[i],alpha_PC[block_alpha[i]]);
    }

  printf("\n\n POLAR REQUIREMENTS CODE\n");
  print_pc_code(PC,block_code);

  MS = setnij_and_ms(block_code, matrix, nij);

  printf("MS0/FT0: %f\n",MS.whole); 
  printf("MS1/FT1: %f\n",MS.one); 
  printf("MS2/FT2: %f\n",MS.two); 
  printf("MS3/FT3: %f\n",MS.three); 

  printf("\n\n\n");

  /* Compute sum_j N_{i,j} for each i*/
  for (i=0;i<NUM_ACIDS;i++) 
    {
      nijSums[i] = 0; 
      for (j=0;j<NUM_ACIDS;j++)
        nijSums[i] += nij[i][j]; 
    }

  /* computes fixednums */
  printf("DiGiulio SumNIJ is: %d\n\n",MS.sumnij);
  printf("DiGiulio NUM_ACIDS x NUM_ACIDS n_ij matrix\n\n\n    ");
  for (i=0;i<NUM_ACIDS;i++) printf("%2d ",i); 
  printf("\n\n");
  for (i=0;i<NUM_ACIDS;i++) {
    printf("%2d: ",i);
    for (j=0;j<NUM_ACIDS;j++)
      printf("%2d ",nij[i][j]);
    printf("\n");
  }

  
  /* calculate A for 0..14 (represents the variables different from 
     lambda)
  */
  for (i=0;i<15;i++)
    for (j=0;j<15;j++)
      if (i==j)
        A[i][i] =  (4.0 / (float) MS.sumnij)*nijSums[a_blocknums[i]];
      else
        A[i][j] = (-4.0 / (float) MS.sumnij)
	  * nij[a_blocknums[i]][a_blocknums[j]];

  /* calculate A for line 15   */
  for (j=0;j<15;j++)
    A[15][j] = 1.0;
  /* calculate A for column 15   */
  for (i=0;i<15;i++)
    A[i][15] = 1.0;
  /* calculate A(15,15)   */
  A[15][15] = 0.0;

  /* print out A matrix */
  printf("\n\n Final A matrix \n");
  for (i=0;i<16;i++) {
    for (j=0;j<16;j++)
      printf("%10.5f",A[i][j]);
    printf("\n");
    }
  
  /* print out B matrix */
  printf("\n\n Final B matrix \n");
  for (i=0;i<15;i++) {
    B[i] = 0;    
    for (j=0;j<5;j++)
        B[i] +=  nij[a_blocknums[i]][fixednums[j]]*V[j];
    B[i] = (4.0 / (float) MS.sumnij)*B[i];
    printf("%10.5f\n",B[i]);
    }
  B[15] = 148.4 - (V[0] +  V[1] + V[2] + V[3] + V[4]);
  printf("%10.5f\n",B[15]);
   
  return 0;
}


/* This function reads a polar requirement vector from a file
   <filename>. The array PC is initialized with the values from the
   file for nonfixed blocks and values from array V for fixed blocks. 
   
   returns 0 in case of success, otherwise 1
*/
int read_non_fixed_vector(char *filename, float PC[NUM_ACIDS], 
                          int fixed[NUM_ACIDS], float V[5])
{
  FILE *fp = fopen(filename, "r");
  int i; /* counters for loops */
  int f=0; /* counts the number of times, a fixed value was used. */
  
  if (fp == NULL)
  {
    fprintf(stderr, "Could not polar requirements file!\n");
    return 1;
  }

  for (i = 0; i < NUM_ACIDS; i++)
    {
      if (fixed[i])
        {
          PC[i] = V[f];
          f++;
        }
      else
        {
          if (fscanf(fp, "%f ", &PC[i]) == EOF)
            {
              fprintf(stderr, "File contains no polar req. vector!\n");
              fclose(fp);
              return 1;
            }
        }
    }
  
  return 0;
}

/* calculates square distance-matrix from PC-vector
 */
void calc_matrix(float PC[NUM_ACIDS], float matrix[NUM_ACIDS][NUM_ACIDS])
{
  int i,j;

  for (i=0;i<NUM_ACIDS;i++)
    for (j=0;j<NUM_ACIDS;j++)
      matrix[i][j]=(PC[i]-PC[j])*(PC[i]-PC[j]);
}

/* print code table with polar requirement from array <PC> for every
   triplet code (translated by <code> to blocknumber)*/
void print_pc_code(float PC[NUM_ACIDS],gcode code)
{
  int i,j,k; 
  for (i=0;i<4;i++)
    {
      for (j=0;j<4;j++)
        {
          for (k=0;k<4;k++)
            if (code[i][j][k] == STOP)
              printf("    STOP"); 
            else 
              printf("%8.4f",PC[code[i][j][k]]);
          printf("\n");
        }
      printf("\n\n");
    }
}



/* prints a blockstructure */
int printblockstruct( gcode block_code, int *fixed )
{
  int i,j,k; 
  for (i=0;i<4;i++)
    {
    for (j=0;j<4;j++)
      {
      for (k=0;k<4;k++) 
        {
          if (block_code[i][j][k] == STOP)
            printf("%s\t","STOP"); 
          else {
            if (fixed[block_code[i][j][k]] == 1) 
              printf("%s\t","FIX");
            else
              printf("%s\t","VAR"); 
          }
        }
      printf("\n");
      }
    printf("\n\n");
    }
  return 0;
}