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

#define MAXYR  200
#define MAXEST 100
#define MSIZE  (MAXEST * (MAXEST + 1) / 2)

#define MAXSIM  600000
#define MAXSTP  500

#define IN stdin
#define IOUT stdout

void PHOUT(double *CL, double RAWCL, int ILAST, int IY);
double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR,
              double *ZMULT,double *POP, double *G);
double CONTRL(double D, double Y, double PLEVEL, double PSLOPE);
void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT, double *FMATRX, 
            int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G);
double STKSIM(double RK, double R, double *POP, double *CATCH);
void SORT(double *ARRAY, double *ARRAY2, int N);
int  _getline(char *buff, size_t size, FILE *fp);

/* from Fortran COMMON VARIABLE */
struct MANPAR {
	double PPROB;  /* Probability level */
	double PYMAX;  /* Maximum value of the productivity parameter (Y) for integration */
	double PNYSTP; /* Number of step in integration over Y */
	double PKSTEP; /* Maximum relative step size in integration over K */
	double PDSTEP; /* Target step size for integration over depletion */
	double PNBSTP; /* Number of steps for integration over the bias parameter */
	double PBMIN;  /* Minimum multiplicative bias eg 0.5 means 50% downward bias */
	double PBMAX;  /* Maximum multiplicative bias eg 1.5 means 50% upward bias */
	double PSCALE; /* Raw deviance scaling factor S = 1/PSCALE**2 */
	double PHASET; /* Number of years without surveys before phaseout invoked */
	double PHASEP; /* Phaseout: anual reduction proportion */
	double PCYCLE; /* Maximum number of years before a new CLC=number of years a CLC is valid */
	double PLEVEL; /* Internal protection level */
	double PSLOPE; /* 'Slope' of catch control law */
} mp = {
	0.4102, 0.05,  200.0,  0.05,  0.01,
	100.0,  0.0,   1.6667, 4.0,   8.0,
	0.2,    5.0,   0.54,   3.0
};
/* All years are scaled so 0 = a year prior to the 1st input data */
struct MANDAT {
	int ISTART;    /* Year of first input data (catch or abundance data) */
	int IYEAR;     /* Year for which to set first catch limit */
	int NS;        /* Number of non-zero estimates */
	int NZ;        /* Number of zero estimates */
	double RKLO;   /* Lower bound used in integration over K */
	double RKHI;   /* Upper bound used in integration over K */
} md;

int main(void)
{
/* ***** LOCAL VARIABLES */
/*
	CHARACTER STOCK*30, FORMT*50
	INTEGER IY, INITYR, I, IS, IYRCL, ILAST, N, IN, IOUT
*/
	double	CATCH[MAXYR + 1],  /* Catch array, indexed by year */
			SIGHT[MAXEST + 1], /* Abundance estimate array, indexed by estimate number */
			FMATRX[MSIZE + 1], /* Information matrix (H) of the log sightings estimates */
							 /* (excludeing zero estimates, stored as a lower triangle */
			ZMULT[MAXEST + 1], /* Poisson multiplier for Nth zero estimate */
			POP[MAXYR + 1],    /* Modelled population size in year I (set in STKSIM) */
			G[MAXYR + 1],      /* set & used in DEVIAN */
			RAWCL,           /* Nominal catch limit i.e. output from the Catch Limit Algorithm */
			CL;              /* The catch limit for the area considerd */
	int 	ISYR[MAXEST + 1],  /* Year of Nth abundance estimate SIGHT(N). N=0, NS-1 */
			IZYR[MAXEST + 1],  /* Year of Nth zero estimate. N=0, NZ-1 */
			POUT;            /* Parameter determing whether phaseout may be applied if */
							 /* necessary. Test for phaseout if POUT=1. (Phaseout is not */
							 /* applied to Medium of Large area nominal catch limits) */
	int 	INITYR,			 /* Year in which first premanagement catch taken */
			ILAST,			 /* Year of the most recent abundance estimate */
			IY, I, IS, IYRCL, N, J, N1;
	double	C, TOTALC;
	int		numlen, len, ret = 0;
	char /* STOCK[30], FORMT[50],*/ buff[300];

/* ***** Read in the data required for application of the CLA */

	_getline(buff, sizeof(buff), IN);
	fprintf(IOUT, "%s by C\n\n", buff);
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */

/* 	Read the year of the first catch, the first year for which catch */
/* limits are to be set & the phaseout option */

	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	INITYR = atoi(&buff[len-numlen]);

	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, INITYR);

	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	IYRCL = atoi(&buff[len-numlen]);

	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, IYRCL);

	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	POUT = atoi(&buff[len-numlen]);
	if (POUT == 1) {
		fprintf(IOUT,   "Apply phaseout if necessary    YES\n");
	} else {
		fprintf(IOUT,   "Apply phaseout                  No\n");
	}

/* Re-scale IYRCL such that 0 is the year prior to the first input data */

	md.ISTART = 0;
	md.IYEAR = IYRCL - INITYR;
	if ((md.IYEAR <= 0) || (md.IYEAR > MAXYR)) {
		fprintf(stderr, "INVALID YEAR");
		ret = 1;
		goto exit;
	}

/* Initialize the catch array */

	for (I = 0; I <= MAXYR; I++) {
		CATCH[I] = 0.0;
	}

/* Read in the catch data, scaling each year to the initial year */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */

	fprintf(IOUT, "\nHistoric catches:\n");
	TOTALC = 0;
	for (I = 0; I <= MAXYR; I++) {
		_getline(buff, sizeof(buff), IN);
		sscanf(buff,"%4d%7lf", &IY, &C);
		if (IY < 0) break; /* EOD */

		fprintf(IOUT, "%4d%6.f\n", IY, C);
		IY = IY - INITYR;
		if (IY < 0 || IY > md.IYEAR) {
			fprintf(stderr, " ERROR: IY is out of range(%d)\n", IY);
			ret = 1;
			goto exit;
		}
		CATCH[IY] = C;
		TOTALC = TOTALC + C;
	}

	if (TOTALC < 0) {
		fprintf(stderr, " ERROR: No historic catch input\n");
		ret = 1;
		goto exit;
	}

/* Read in the non-zero sightings estimates and information matrix */

	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	md.NS = atoi(&buff[len-numlen]);

	fprintf(IOUT,   "\nAbundance estimates:\n");
	if (md.NS > MAXEST) {
		fprintf(stderr, "ERROR: Abundance year out of range\n");
		ret = 1;
		goto exit;
	}
	_getline(buff, sizeof(buff), IN);
	for (N = 0; N <= md.NS - 1; N++) {
		N1 = (N * (N + 1)) / 2;
		_getline(buff, sizeof(buff), IN);
		sscanf(buff, "%4d%10lf", &ISYR[N], &SIGHT[N]);
		fprintf(IOUT,   "%4d %8.f", ISYR[N], SIGHT[N]);
		len = strlen(buff);
		for (J = 0; J < N + 1; J++) {
			if (len > J * 7 + 14) {
				sscanf(&buff[0] + J * 7 + 14,"%7lf", &FMATRX[N1 + J]);
			} else {
				FMATRX[N1 + J] = 0.0;
			}
			fprintf(IOUT,   "%7.f", FMATRX[N1 + J]);
		}
		fprintf(IOUT,   "\n");
		ISYR[N] = ISYR[N] - INITYR;
		if (ISYR[N] < 0 || ISYR[N] >= md.IYEAR) {
			fprintf(stderr, "ERROR: Sight year out of range\n");
			ret = 1;
			goto exit;
		}
		if (SIGHT[N] <= 0.0) {
			fprintf(stderr, "ERROR: Estimate not positive\n");
			ret = 1;
			goto exit;
		}
	}

/* Read in the Poisson multipliers for any zero sightings estimates */

	_getline(buff, sizeof(buff), IN);
	_getline(buff, sizeof(buff), IN);
	numlen = 3; len = strlen(buff);
	md.NZ = atoi(&buff[len-numlen]);
	if (md.NZ > 0) {
		fprintf(IOUT,   "\nZero abundance estimates:\n");
	}
	if (md.NZ > MAXEST) {
		fprintf(stderr, "ERROR: Zero estimate array too small\n");
		ret = 1;
		goto exit;
	}
	_getline(buff, sizeof(buff), IN);
	for (N = 0; N <= md.NZ - 1; N++) {
		N1 = (N * (N + 1)) / 2;
		_getline(buff, sizeof(buff), IN);
		sscanf(buff, "%4d%10lf", &IZYR[N], &ZMULT[N]);
		fprintf(IOUT,   "%4d%7.f\n", IZYR[N], ZMULT[N]);
		IZYR[N] = IZYR[N] - INITYR;
		if (IZYR[N] < 0 || IZYR[N] >= md.IYEAR) {
			fprintf(stderr, "Sight year out of range\n");
			ret = 1;
			goto exit;
		}
                                    
		if (ZMULT[N] <= 0.0) {
			fprintf(stderr, "ERROR: Multiplier not positive\n");
			ret = 1;
			goto exit;
		}
	}
	fprintf(IOUT,"\n");
/* Bound the range for the integration over K */

	md.RKHI = 1.0E7;
	md.RKLO = 0.0;
	IS = 0; /* orignal not initialized bug? */

/* ****** Run the CLA to obtain the nominal catch limit */
	RAWCL = CLIMIT(CATCH, SIGHT, FMATRX, ISYR, IZYR, ZMULT, POP, G);
/* Set the catch limits for PCYCLE years. If the catch limit may be */
/* subject to phaseout, call POUT to apply the phaseout rule. */

/* First set ILAST = year of the most recent abundance estimate */

	if (md.NS > 0) ILAST = ISYR[md.NS - 1];
	if (md.NZ > 0) ILAST = (IS > IZYR[md.NZ - 1])? IS: IZYR[md.NZ - 1];

	for (IY = md.IYEAR; IY <= md.IYEAR + mp.PCYCLE - 1; IY++) {
		if (POUT == 1) {
			PHOUT(&CL, RAWCL, ILAST, IY);
		} else {
			CL = RAWCL;
		}
		fprintf(IOUT,   "Year: %4d     Catch limit:%6d\n", IY + INITYR, (int)(CL + 0.5));

	}

exit:
	return ret;
}

/*
C **************************************************************************************
C
C	CLC version 6
C
C **************************************************************************************
C
C	31 January 1994
C
C **************************************************************************************
*/
/*
	SUBROUTINE PHOUT (CL, RAWCL, ILAST, IY)
*/
void PHOUT(double *CL, double RAWCL, int ILAST, int IY)
{
/*
	PHASET Number of years without surveys before phaseout invoked
	PHASEP Phaseout annual reduction proportion
	Phaseout: Reduce catch limit if there is no survey data in the
	last PHASET years
*/
	if (IY >= ILAST + mp.PHASET) {
		*CL = RAWCL * (1.0 - mp.PHASEP * (IY - ILAST - mp.PHASET));
		if (*CL < 0.0) *CL = 0.0;
	} else {
		*CL = RAWCL;
	}
	/* RETURN
	   END */
}

/*
	REAL FUNCTION CLIMIT (CATCH,SIGHT,FMATRX,ISYR,IZYR,ZMULT,POP,G)
	Run the CLA to obtain to nominal catch limit
*/
double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR, double *ZMULT,double *POP, double *G)
{
	double SF, SS0, SS1, SS2, SS3;
	double Y[MAXSTP + 1], B[MAXSTP + 1], RLGB[MAXSTP + 1];
	double *PRES = NULL, *QRES = NULL;

	int NB, NR;
	double BINC, YINC, DK, R, RK, D, DD, P, Q, PTOT;
	int I, J, N, N2;

	double _CLIMIT = 0.0;
	if (md.NS <= 0) return _CLIMIT;

	PRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
	if (PRES == NULL) {
		fprintf(stderr, "Can't Allocate memory (PRES)\n");
		goto err_exit;
	}
	memset(PRES, 0, sizeof(double) * (MAXSTP + 1));

	QRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
	if (QRES == NULL) {
		fprintf(stderr, "Can't Allocate memory (QRES)\n");
		goto err_exit;
	}
	memset(QRES, 0, sizeof(double) * (MAXSTP + 1));

	SF = 0.5 / (mp.PSCALE * mp.PSCALE);

	/* Check the sizes of the Y and B arrays are large enough */
	if (mp.PNBSTP > MAXSTP || mp.PNYSTP > MAXSTP) {
		fprintf(stderr, "Y &/or b array sizes too small\n");
		goto err_exit;
	}
	/* Set sightings bias step sizes & their log values. BINC = Bias increment */
	NB = (int)mp.PNBSTP;
	BINC = (mp.PBMAX - mp.PBMIN) / NB;
	for (I = 0; I <= NB - 1; I++) {
        *(B + I) = mp.PBMIN + (I + 0.5) * BINC;
        *(RLGB + I) = -log(*(B + I));
	}

	/* Set productivity parameter step sizes (midpoints) */
	NR = (int)mp.PNYSTP;
	YINC = mp.PYMAX / NR;
	for (I = 0; I <= NR - 1; I++) {
		*(Y + I) = (I + 0.5) * YINC;
	}

	PTOT = 0;
	N = 0;
	for (I = 0; I <= NR - 1; I++) {
		/* Set R from the productivity parameter, Y */
		R = 1.4184 * *(Y + I);

		/* Step size for K */
        DK = mp.PKSTEP;
        D = 1.0;
        RK = md.RKHI;

		/* Use function STKSIM to set up the Nth population trajectory */
		/* i.e. set up the pop array */

		while (1) {
			if ((RK <= md.RKLO) || (STKSIM(RK, R, POP, CATCH) <= 0.0)) break;
			if (N >= MAXSIM) {
				fprintf(stderr, "ERROR: TOO MANY SIMULATIONS");
				goto err_exit;
			}
			/* How much depletion covered? */
			DD = D - POP[md.IYEAR] / RK;
			D = POP[md.IYEAR] / RK;
			P = 0.0;

			if (DD > 0.0) {
				/* Compute the internal catch limit corresponding to D and Y(I) */
				QRES[N] = CONTRL(D, Y[I], mp.PLEVEL, mp.PSLOPE) * POP[md.IYEAR];

				/* Calculate deviance */
            	DEVIAN(&SS0, &SS1, &SS2, &SS3, SIGHT, FMATRX, ISYR,IZYR, ZMULT, POP, G);

				/* Scale likelihood and integrate over values for the bias parameter */
				for (J = 0; J <= NB - 1; J++) {
					P = P + exp(-SF * (SS0 + *(RLGB + J) * (SS1 + *(RLGB + J) * SS2) + SS3 * *(B + J)));
				}

				/* Calculate the weight for this point (& total area under likelihood) */
				PRES[N] = P * DD;
				PTOT = PTOT + PRES[N];

				/* Update counter */
				N = N + 1;

				/* Find the next K */
				DK = DK * mp.PDSTEP / DD;
				if (DK > mp.PKSTEP) DK = mp.PKSTEP;

			} else {
				DK = mp.PKSTEP;
			}
			/* Set the new value of K */
			RK = RK / (1.0 + DK);
		}
	}
	if (PTOT < 0.0) {
		fprintf(stderr,"ERROR: PROB INTEGRATES TO ZERO");
		goto err_exit;
	}

	/* Sort the QRES and PRES arrays in ascending order of QRES */
	N2 = N;
	SORT(QRES, PRES, N2);

	/* Normalize the relative likelihoods */
	for (I = 0; I <= N - 1; I++) {
		PRES[I] = PRES[I] / PTOT;
	}

	/* Extract the desired probability level: the nominal catch limit (NCL) */
	/* is the lower PROB% of the distribution. */
	/* First calculate PRES(I), the probability that the NCL is between */
	/* QRES(I) & QRES(I + 1). */

	P = 0;
	for (I = 0; I <= N - 1; I++) {
		P = P + PRES[I];
		if (P > mp.PPROB) break;
	}

	/* Interpolate to set the nominal catch limit */
	if  (I >= N - 1) {
        Q = QRES[N - 1];
	} else {
        Q = (QRES[I + 1] * (mp.PPROB - P + PRES[I]) + QRES[I] * (P - mp.PPROB)) / PRES[I];
	}

	_CLIMIT = Q;

	if (QRES != NULL) free(QRES);
	if (PRES != NULL) free(PRES);
	return _CLIMIT;

err_exit:
	if (QRES != NULL) free(QRES);
	if (PRES != NULL) free(PRES);
	exit(1);
	return 0.0; /* not reached */
}

/*
	FUNCTION CONTRL (D, Y, PLEVEL, PSLOPE)
	Catch control law
*/
double CONTRL(double D, double Y, double PLEVEL, double PSLOPE)
{
	double ret;
	if (D < PLEVEL) {
		ret = 0.0;
	} else if (D < 1.0) {
		ret = PSLOPE * Y * (D - PLEVEL);
	} else {
		ret = PSLOPE * Y * (1.0 - PLEVEL);
	}
	return ret;
	/* END */
}

/*
	SUBROUTINE DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR,
	Calculate deviance (-2 log likelihood) in terms of coefficients for 
	the bias and log bias parameters
*/
void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT,
            double *FMATRX, int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G)
{
	int N, J, K;

	*SS0 = 0.0;
	*SS1 = 0.0;
	*SS2 = 0.0;
	*SS3 = 0.0;

	for (N = 0; N <= md.NS - 1; N++) { 

		*(G + N) = log(*(SIGHT + N) / *(POP+*(ISYR + N)));
		K = N * (N + 1) / 2;

		for (J = 0; J <= N - 1; J++) {
			/*1st add non diagonal contributions (which are doubled up)*/
			*SS0 = *SS0 + 2.0 * *(G + J) * *(G + N) * *(FMATRX + K);
			*SS1 = *SS1 + 2.0 * (*(G + J) + *(G + N)) * *(FMATRX + K);
			*SS2 = *SS2 + *(FMATRX + K) + *(FMATRX + K);
			K = K + 1;
		}

		/* Now add diagnoal contribution */
        *SS0 = *SS0 + *(G + N) * *(G + N) * *(FMATRX + K);
        *SS1 = *SS1 + 2.0 * *(G + N) * *(FMATRX + K);
        *SS2 = *SS2 + *(FMATRX + K);

	}

	/* Now do the zero estimates */
	for (N =0; N <= md.NZ - 1; N++) {
		*SS3 = *SS3 + 2.0 * *(POP + *(IZYR + N)) / *(ZMULT + N);
 	}

	return;
}
/*
	FUNCTION STKSIM (RK, R, POP, CATCH)
	Calculate the stock trajectory with parameter RK and R: return
	the current stock size
	RK = notional carrying capacity; R = productivity parameter * 1.4184 
*/
double STKSIM(double RK, double R, double *POP, double *CATCH)
{
	int I;
	double D;
	*(POP+md.ISTART) = RK;
	for (I = md.ISTART + 1; I <= md.IYEAR; I++) {
		D = *(POP + I - 1) / RK;
		*(POP + I) = *(POP + I - 1) * (1.0 + R * (1.0 - D * D)) - *(CATCH + I -1);
		if (*(POP + I) <= 0) {
			return 0.0;
		}
	}

	return *(POP + md.IYEAR);
}
/* 
	SUBROUTINE SORT (ARRAY, ARRAY2, N)
	SORT sorts a pair of arrays in ascending order of the first array
   (using the heapsort algorithm)
*/
void SORT(double *ARRAY, double *ARRAY2, int N)
{
	double /*ARRAY(0:*), ARRAY2(0:*),*/ TEMP, TEMP2;
	int  /*N,*/ K, IR, I, J;
	if (N < 2) return;
	K = N / 2;
	IR = N - 1;
	while (1) { /*10*/
		if (K != 0) {
			K = K - 1;
			TEMP = *(ARRAY + K);
			TEMP2 = *(ARRAY2 + K);
		} else {
			TEMP = *(ARRAY + IR);
			TEMP2 = *(ARRAY2 + IR);
			*(ARRAY + IR) = *ARRAY;
			*(ARRAY2 + IR) = *ARRAY2;
			IR = IR - 1;
			if (IR == 0) {
				*ARRAY = TEMP;
				*ARRAY2 = TEMP2;
				return;
			}
		}
		I = K;
		J = K + K + 1;
		while (1) {/*20*/
			if (J > IR) break;
			if ((J < IR) && (*(ARRAY + J) < *(ARRAY+ J + 1))) J = J + 1;
			if (TEMP < *(ARRAY + J)) {
				*(ARRAY + I) = *(ARRAY + J);
				*(ARRAY2 + I) = *(ARRAY2 + J);
				I = J;
				J = J + I + 1;
			} else {
				J = IR + 1;
			}
		} /*GOTO 20*/
		*(ARRAY + I) = TEMP;
		*(ARRAY2 + I) = TEMP2;
	}/*GOTO 10*/
	/*END*/
}
int _getline(char *buff, size_t size, FILE *fp){
	char *ret;
	ret = fgets(buff, size, fp);
	if (ret != NULL) {
		if (buff[strlen(buff) - 1] == '\n')
			buff[strlen(buff) - 1] = 0;
		if (buff[strlen(buff) - 1] == '\r')
			buff[strlen(buff) - 1] = 0;
	}
	return (ret != NULL);
}