#include <stdio.h>#include <math.h>#include <stdlib.h>#include <string.h>#include "clustalw.h"/* *   Prototypes *//* *   Global variables */extern sint max_aa,gap_pos1,gap_pos2;void calc_prf1(sint **profile, char **alignment, sint *gaps,  sint matrix[NUMRES][NUMRES],  sint *seq_weight, sint prf_length, sint first_seq, sint last_seq){  sint **weighting, sum2, d, i, res;   sint numseq;  sint r, pos;  int f;  float scale;  weighting = (sint **) ckalloc( (NUMRES+2) * sizeof (sint *) );  for (i=0;i<NUMRES+2;i++)    weighting[i] = (sint *) ckalloc( (prf_length+2) * sizeof (sint) );  numseq = last_seq-first_seq;  sum2 = 0;  for (i=first_seq; i<last_seq; i++)       sum2 += seq_weight[i];  for (r=0; r<prf_length; r++)   {      for (d=0; d<=max_aa; d++)        {            weighting[d][r] = 0;            for (i=first_seq; i<last_seq; i++)               if (d == alignment[i][r]) weighting[d][r] += seq_weight[i];        }      weighting[gap_pos1][r] = 0;      for (i=first_seq; i<last_seq; i++)         if (gap_pos1 == alignment[i][r]) weighting[gap_pos1][r] += seq_weight[i];      weighting[gap_pos2][r] = 0;      for (i=first_seq; i<last_seq; i++)         if (gap_pos2 == alignment[i][r]) weighting[gap_pos2][r] += seq_weight[i];   }  for (pos=0; pos< prf_length; pos++)    {      if (gaps[pos] == numseq)        {           for (res=0; res<=max_aa; res++)             {                profile[pos+1][res] = matrix[res][gap_pos1];             }           profile[pos+1][gap_pos1] = matrix[gap_pos1][gap_pos1];           profile[pos+1][gap_pos2] = matrix[gap_pos2][gap_pos1];        }      else        {           scale = (float)(numseq-gaps[pos]) / (float)numseq;           for (res=0; res<=max_aa; res++)             {                f = 0;                for (d=0; d<=max_aa; d++)                     f += (weighting[d][pos] * matrix[d][res]);                f += (weighting[gap_pos1][pos] * matrix[gap_pos1][res]);                f += (weighting[gap_pos2][pos] * matrix[gap_pos2][res]);                profile[pos+1][res] = (sint  )(((float)f / (float)sum2)*scale);             }           f = 0;           for (d=0; d<=max_aa; d++)                f += (weighting[d][pos] * matrix[d][gap_pos1]);           f += (weighting[gap_pos1][pos] * matrix[gap_pos1][gap_pos1]);           f += (weighting[gap_pos2][pos] * matrix[gap_pos2][gap_pos1]);           profile[pos+1][gap_pos1] = (sint )(((float)f / (float)sum2)*scale);           f = 0;           for (d=0; d<=max_aa; d++)                f += (weighting[d][pos] * matrix[d][gap_pos2]);           f += (weighting[gap_pos1][pos] * matrix[gap_pos1][gap_pos2]);           f += (weighting[gap_pos2][pos] * matrix[gap_pos2][gap_pos2]);           profile[pos+1][gap_pos2] = (sint )(((float)f / (float)sum2)*scale);        }    }  for (i=0;i<NUMRES+2;i++)    weighting[i]=ckfree((void *)weighting[i]);  weighting=ckfree((void *)weighting);}