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

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

#define TRUE			1
#define	FALSE			0
#define	FILE_NAME_SZ		81
#define	HYDRO_CHARGE_TB		"cnnrchrg.hyd"
#define	HYDRO_CHARGE_TB_SZ	21
#define	RADIAN			(3.14159265/180)
#define	LINE_SZ			101
#define	MARGIN			6
#define	MAX_PROTEIN_CNT		1000
#define AVE_SZ			2

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

typedef struct {
    float hydrophobicity, charge;
    char residue;
} HydroCharge;

typedef struct {
    char *pProt;
    float *pAveAlpha, *pAveBeta;
} ResPos;

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

static HydroCharge hydroChargeTb[HYDRO_CHARGE_TB_SZ];
static ResPos protArray[MAX_PROTEIN_CNT];
static int  protCnt, longestProt;

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

CleanUp()
{
    int i;

    if (protCnt > 0)
    {
	for (i=0; i<protCnt; i++)
	{
	    if (protArray[i].pProt != NULL)
	    {
		free(protArray[i].pProt);
		free(protArray[i].pAveAlpha);
		free(protArray[i].pAveBeta);
	    }
	}
    }
}

GetLine(fp, pLine)
FILE *fp;
char *pLine;
{
    while (!feof(fp))
    {
	if ((*pLine = fgetc(fp)) == '\n')
	    break;
	else
	    pLine++;
    }

    *pLine = '\0';
}

int FindLine(fp, pStr)
FILE *fp;
char *pStr;
{
    char line[LINE_SZ];

    while (!feof(fp))
    {
	GetLine(fp, line);
	if (strstr(line, pStr) != NULL)
	    return TRUE;
    }

    return FALSE;
}

GetHydroChargeTable()
{
    FILE *fp;
    int i;
    float hydro, charge;
    char str[80];

    if ((fp = fopen(HYDRO_CHARGE_TB, "r")) == NULL)
    {
	printf("Could not open hydrophobicity and charge table: %s\n",
		HYDRO_CHARGE_TB);
	exit(1);
    }

    for (i=0; i<HYDRO_CHARGE_TB_SZ; i++)
    {
	fscanf(fp, "%f", &hydroChargeTb[i].hydrophobicity);
	fscanf(fp, "%f", &hydroChargeTb[i].charge);
	fscanf(fp, "%s", str);
	hydroChargeTb[i].residue = str[0];

	if (feof(fp))
	{
	    printf("There should be 21 entries in %s file\n", HYDRO_CHARGE_TB);
	    exit(1);
	}
    }

    fclose(fp);
/*
for (i=0; i<HYDRO_CHARGE_TB_SZ; i++)
printf("%1.2f, %1.2f, %c\n", hydroChargeTb[i].hydrophobicity, 
hydroChargeTb[i].charge, hydroChargeTb[i].residue);
*/
}

FILE* OpenProteinFile()
{
    FILE *fp;
    char fileName[FILE_NAME_SZ];

    printf("Enter protein file name: ");
    scanf("%s", fileName);

    if ((fp = fopen(fileName, "r")) == NULL)
	printf("Could not open %s\n", fileName);

    return fp;
}

FILE* OpenOutputFile()
{
    FILE *fp;
    char fileName[FILE_NAME_SZ];

    printf("Enter output file name: ");
    scanf("%s", fileName);

    if ((fp = fopen(fileName, "w")) == NULL)
	printf("Could not open %s\n", fileName);

    return fp;
}

float GetAngle(pStr)
char *pStr;
{
    float angle;

    printf("%s", pStr);
    scanf("%f", &angle);
    return angle;
}

int GetSlidingWindowSize()
{
    int winSz;

    for (;;)
    {
	printf("Enter window length: ");
	scanf("%i", &winSz);
	if (winSz > 0)
	    break;
    }

    return winSz;
}

int GetProtCnt(fp)
FILE *fp;
{
    int cnt;
    char line[LINE_SZ];

    cnt = 0;
    for (;;)
    {
	GetLine(fp, line);
	if (strlen(line) < 5 || strstr(line, "-- Results --") != NULL)
	    break;
	else
	    cnt++;
    }

    return cnt;
}

int GetLongestProteinLength(fp)
FILE *fp;
{
    char line[LINE_SZ];
    int cnt;

    cnt = 0;
    for (;;)
    {
	GetLine(fp, line);
	if (strstr(line, "-- Results --") != NULL)
	    break;

	cnt = cnt + ((strlen(line) - MARGIN) / 2);
	for (;;)
	{
	    GetLine(fp, line);
	    if (strstr(line, "-- Results --") != NULL)
		return cnt;

	    if (strlen(line) < 5 || feof(fp))
		break;
	}
    }

    return cnt;
}

int GetProtArray(fp)
FILE *fp;
{
    char line[LINE_SZ];
    long curPos;
    int i,j, cnt1, cnt2;

    if (!FindLine(fp, "gap penalty"))
    {
	printf("The input file is not the correct output of ree program.\n");
	return FALSE;
    }

    GetLine(fp, line);	/* Skip a blank line */
    curPos = ftell(fp);
    protCnt = GetProtCnt(fp);
    fseek(fp, curPos, 0);
    longestProt = GetLongestProteinLength(fp);
    fseek(fp, curPos, 0);

    for (i=0; i<protCnt; i++)
    {
	protArray[i].pProt = (char *) malloc(longestProt * sizeof(char) + 1);
	protArray[i].pAveAlpha = (float *) malloc(longestProt * sizeof(float));
	protArray[i].pAveBeta = (float *) malloc(longestProt * sizeof(float));
	if (protArray[i].pProt == NULL || protArray[i].pAveAlpha == NULL ||
	protArray[i].pAveBeta == NULL)
	{
	    printf("Could not allocate memory.\n");
	    CleanUp();
	    exit(1);
	}

	memset(protArray[i].pAveAlpha, 0, longestProt * sizeof(float));
	memset(protArray[i].pAveBeta, 0, longestProt * sizeof(float));
    }

    cnt1 = 0;
    for (;;)
    {
	for (i=0; i<protCnt; i++)
	{
	    GetLine(fp, line);
	    cnt2 = cnt1;
	    for (j=MARGIN; j<strlen(line); j+=2)
	    {
		*(protArray[i].pProt + cnt2) = line[j];
		cnt2++;
	    }
	}
        cnt1 = cnt2;
	GetLine(fp, line);	/* Skip blank Line */
	if (strstr(line, "-- Results --") != NULL)
	    break;
    }

    return TRUE;
}

double GetResidueHydrophobicity(residue)
char residue;
{
    int i;

    /* Search for the residue's hydrophobicity */
    for (i=0; i<HYDRO_CHARGE_TB_SZ; i++)
    {
	if (residue == hydroChargeTb[i].residue)
	    return (double) hydroChargeTb[i].hydrophobicity;
    }

    /* Search for unknown residue's value */
    for (i=0; i<HYDRO_CHARGE_TB_SZ; i++)
    {
	if (residue == 'U')
	    return (double) hydroChargeTb[i].hydrophobicity;
    }

    /* All else, return zero */
    printf("Could not find hydrophobicity value for residue %c. Zero is used.\n");
    return 0.0;
}

GetSlidingWindowData(pWin, alphaAngle, betaAngle, pAveAlpha, pAveBeta)
char *pWin;
float alphaAngle, betaAngle;
float *pAveAlpha, *pAveBeta;
{
    int i, winSz;
    double sumAlphaX, sumAlphaY, sumBetaX, sumBetaY, hValue, sumHValue, angle;
    float aveH;
    char residue;

    sumAlphaX = sumAlphaY = sumBetaX = sumBetaY = sumHValue = 0;
    winSz = strlen(pWin);
    for (i=0; i<winSz; i++)
    {
	residue = *(pWin + i);
	hValue = GetResidueHydrophobicity(residue);
	sumHValue += hValue;

	angle = alphaAngle * RADIAN * i;
	sumAlphaX += (hValue * cos(angle));
	sumAlphaY += (hValue * sin(angle));

	angle = betaAngle * RADIAN * i;
	sumBetaX += (hValue * cos(angle));
	sumBetaY += (hValue * sin(angle));
    }

    aveH = (float) (sumHValue / winSz);
    *pAveAlpha = (float) sqrt((sumAlphaX * sumAlphaX) + (sumAlphaY * sumAlphaY));
    *pAveAlpha = *pAveAlpha / winSz;
    *pAveBeta = (float) sqrt((sumBetaX * sumBetaX) + (sumBetaY * sumBetaY));
    *pAveBeta = *pAveBeta / winSz;
}

int GetNextWinCenterOffset(protArrayIndex, offset)
int protArrayIndex, offset;
{
    for (;;)
    {
	offset++;
	if (*(protArray[protArrayIndex].pProt + offset) != ' ')
	    return offset;
    }
}

int FirstResidue(protArrayIndex)
int protArrayIndex;
{
    int cnt=0;

    for (;;)
    {
	if (*(protArray[protArrayIndex].pProt + cnt) != ' ')
	    return cnt;
	cnt++;
    }
}

GetProtAve(protArrayIndex, winSz, alpha, beta)
int protArrayIndex, winSz;
float alpha, beta;
{
    char *pWin, *pTmp;
    int i, cnt, len, offset, curOffset;
    float aveAlpha, aveBeta;

    len = strlen(protArray[protArrayIndex].pProt);
    cnt = 0;
    pTmp = (char*) malloc(len + 1);
    pWin = (char*) malloc(winSz + 1);
    if (pTmp == NULL || pWin == NULL)
    {
	printf("Could not allocate memory.\n");
	CleanUp();
	exit(1);
    }

    for (i=0; i<len; i++)	/* Remove blank space */
    {
	if (*(protArray[protArrayIndex].pProt + i) != ' ')
	{
	    *(pTmp + cnt) = *(protArray[protArrayIndex].pProt + i);
	    cnt++;
	}
    }
    *(pTmp + cnt) = '\0';

    offset = FirstResidue(protArrayIndex) + (winSz/2);
    *(pWin + winSz) = '\0';
    for (i=0; i<=cnt-winSz; i++)
    {
	strncpy(pWin, pTmp + i, winSz);
	GetSlidingWindowData(pWin, alpha, beta, &aveAlpha, &aveBeta);
	*(protArray[protArrayIndex].pAveAlpha + offset) = aveAlpha;
	*(protArray[protArrayIndex].pAveBeta + offset) = aveBeta;
	offset = GetNextWinCenterOffset(protArrayIndex, offset);
    }

    free(pTmp);
    free(pWin);
}

OutputData(fp, winSz)
FILE *fp;
int winSz;
{
    int i, j, nonGapCnt;
    float sumAlpha, sumBeta, prevAveBeta, aveBeta;

    fprintf(fp, "Alignment     Average Hydrophobic Moment    Percent of\n");
    fprintf(fp, "Position         Alpha         Beta           Nongap\n");
    for (i=0; i<longestProt; i++)
    {
	nonGapCnt = 0;
	for (j=0; j<protCnt; j++)
	{
	    if (*(protArray[j].pProt + i) != ' ')
		nonGapCnt++;
	}

	fprintf(fp, "   %3i           ", i+1);
	sumAlpha = sumBeta = 0;
	for (j=0; j<protCnt; j++)
	{
	    sumAlpha += *(protArray[j].pAveAlpha + i);
	    sumBeta += *(protArray[j].pAveBeta + i);
	}

	fprintf(fp, "%1.2f          ", sumAlpha/nonGapCnt);
	aveBeta = sumBeta/nonGapCnt;
	if (i > winSz/2)
	    fprintf(fp, "%1.2f            ", (aveBeta + prevAveBeta) / 2);
	else
	    fprintf(fp, "0.0             ");
	fprintf(fp, "%2.f\n", (float) nonGapCnt*100.0/protCnt);

	prevAveBeta = aveBeta;
    }
}

main()
{
    FILE *fpProtein, *fpOutput;
    float alphaAngle, betaAngle; 
    int winSz, i;

    protCnt = 0;
    for (i=0; i<MAX_PROTEIN_CNT; i++)
	protArray[i].pProt = NULL;

    GetHydroChargeTable();

    if ((fpProtein = OpenProteinFile()) == NULL)
	exit(1);

    if ((fpOutput = OpenOutputFile()) == NULL)
	exit(1);

    alphaAngle = GetAngle("Enter Alpha Angle Value: ");
    betaAngle = GetAngle("Enter Beta Angle Value: ");
    winSz = GetSlidingWindowSize();

    if (!GetProtArray(fpProtein))
    {
	printf("Could not retrieve proteins array.\n");
	CleanUp();
    }

    for (i=0; i<protCnt; i++)
	GetProtAve(i, winSz, alphaAngle, betaAngle);

    OutputData(fpOutput, winSz);

    CleanUp();
    fclose(fpOutput);
    fclose(fpProtein);
}