fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5.  
  6. #define MAXYR  200
  7. #define MAXEST 100
  8. #define MSIZE  (MAXEST * (MAXEST + 1) / 2)
  9.  
  10. #define MAXSIM  600000
  11. #define MAXSTP  500
  12. #define CLCPAR  "CLC.PAR"
  13. #define CLCDAT  "CLC.DAT"
  14. #define CLCOUT  "CLC.OUT"
  15.  
  16. void RDPARS(void);
  17. void PHOUT(double *CL, double RAWCL, int ILAST, int IY);
  18. double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR,
  19.               double *ZMULT,double *POP, double *G);
  20. double CONTRL(double D, double Y, double PLEVEL, double PSLOPE);
  21. void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT, double *FMATRX, 
  22.             int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G);
  23. double STKSIM(double RK, double R, double *POP, double *CATCH);
  24. void SORT(double *ARRAY, double *ARRAY2, int N);
  25. int  _getline(char *buff, size_t size, FILE *fp);
  26.  
  27. /* from Fortran COMMON VARIABLE */
  28. struct MANPAR {
  29. 	double PPROB;  /* Probability level */
  30. 	double PYMAX;  /* Maximum value of the productivity parameter (Y) for integration */
  31. 	double PNYSTP; /* Number of step in integration over Y */
  32. 	double PKSTEP; /* Maximum relative step size in integration over K */
  33. 	double PDSTEP; /* Target step size for integration over depletion */
  34. 	double PNBSTP; /* Number of steps for integration over the bias parameter */
  35. 	double PBMIN;  /* Minimum multiplicative bias eg 0.5 means 50% downward bias */
  36. 	double PBMAX;  /* Maximum multiplicative bias eg 1.5 means 50% upward bias */
  37. 	double PSCALE; /* Raw deviance scaling factor S = 1/PSCALE**2 */
  38. 	double PHASET; /* Number of years without surveys before phaseout invoked */
  39. 	double PHASEP; /* Phaseout: anual reduction proportion */
  40. 	double PCYCLE; /* Maximum number of years before a new CLC = number of years a CLC is valid */
  41. 	double PLEVEL; /* Internal protection level */
  42. 	double PSLOPE; /* 'Slope' of catch control law */
  43. } mp;
  44.  
  45. /* All years are scaled so 0 = a year prior to the 1st input data */
  46. struct MANDAT {
  47. 	int ISTART;    /* Year of first input data (catch or abundance data) */
  48. 	int IYEAR;     /* Year for which to set first catch limit */
  49. 	int NS;        /* Number of non-zero estimates */
  50. 	int NZ;        /* Number of zero estimates */
  51. 	double RKLO;   /* Lower bound used in integration over K */
  52. 	double RKHI;   /* Upper bound used in integration over K */
  53. } md;
  54.  
  55. int main(void)
  56. {
  57. /* ***** LOCAL VARIABLES */
  58. /*
  59. 	CHARACTER STOCK*30, FORMT*50
  60. 	INTEGER IY, INITYR, I, IS, IYRCL, ILAST, N, IN, IOUT
  61. */
  62. 	double	CATCH[MAXYR + 1],  /* Catch array, indexed by year */
  63. 			SIGHT[MAXEST + 1], /* Abundance estimate array, indexed by estimate number */
  64. 			FMATRX[MSIZE + 1], /* Information matrix (H) of the log sightings estimates */
  65. 							 /* (excludeing zero estimates, stored as a lower triangle */
  66. 			ZMULT[MAXEST + 1], /* Poisson multiplier for Nth zero estimate */
  67. 			POP[MAXYR + 1],    /* Modelled population size in year I (set in STKSIM) */
  68. 			G[MAXYR + 1],      /* set & used in DEVIAN */
  69. 			RAWCL,           /* Nominal catch limit i.e. output from the Catch Limit Algorithm */
  70. 			CL;              /* The catch limit for the area considerd */
  71. 	int 	ISYR[MAXEST + 1],  /* Year of Nth abundance estimate SIGHT(N). N=0, NS-1 */
  72. 			IZYR[MAXEST + 1],  /* Year of Nth zero estimate. N=0, NZ-1 */
  73. 			POUT;            /* Parameter determing whether phaseout may be applied if */
  74. 							 /* necessary. Test for phaseout if POUT=1. (Phaseout is not */
  75. 							 /* applied to Medium of Large area nominal catch limits) */
  76. 	int 	INITYR,			 /* Year in which first premanagement catch taken */
  77. 			ILAST,			 /* Year of the most recent abundance estimate */
  78. 			IY, I, IS, IYRCL, N, J, N1;
  79. 	double	C, TOTALC;
  80. 	int		numlen, len, ret = 0;
  81. 	char /* STOCK[30], FORMT[50],*/ buff[300];
  82.  
  83. /* ***** Read in the data required for application of the CLA */
  84.  
  85. 	/* DATA IN/2/, IOUT/3/ */
  86. 	FILE *IN = NULL, *IOUT = NULL;
  87. 	/* OPEN (IN, FILE = 'CLC.DAT') */
  88. 	if ((IN = fopen(CLCDAT, "rt")) == NULL) {
  89. 		fprintf(stderr, "Can't open " CLCDAT "\n");
  90. 		ret = 1;
  91. 		goto exit;
  92. 	}
  93. 	/* OPEN (IOUT, FILE = 'CLC.OUT') */
  94. 	if ((IOUT = fopen(CLCOUT, "w")) == NULL) {
  95. 		fprintf(stderr, "Can't open " CLCOUT "\n");
  96. 		ret = 1;
  97. 		goto exit;
  98. 	}
  99.  
  100. 	RDPARS();
  101.  
  102. 	/* READ (IN, '(A30/)') STOCK */
  103. 	_getline(buff, sizeof(buff), IN);
  104. 	/* WRITE (IOUT, '(A30/)') STOCK */
  105. 	fprintf(IOUT, "%s by C\n\n", buff);
  106. 	fprintf(stdout, "%s by C\n", buff);
  107. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  108.  
  109. /* 	Read the year of the first catch, the first year for which catch */
  110. /* limits are to be set & the phaseout option */
  111.  
  112. 	/* READ (IN, '(T30,I4)') INITYR */
  113. 	_getline(buff, sizeof(buff), IN);
  114. 	numlen = 4; len = strlen(buff);
  115. 	INITYR = atoi(&buff[len-numlen]);
  116.  
  117. 	/* WRITE (IOUT, '(A,T30,I4)') 'Year of first input data', INITYR */
  118. 	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, INITYR);
  119. 	fprintf(stdout, "%.*s %4d\n", len - numlen, buff, INITYR);
  120.  
  121. 	/* READ (IN, '(T30,I4)') IYRCL */
  122. 	_getline(buff, sizeof(buff), IN);
  123. 	numlen = 4; len = strlen(buff);
  124. 	IYRCL = atoi(&buff[len-numlen]);
  125.  
  126. 	/* WRITE (IOUT, '(A, T30, I4)') 'Year of first catch limit', IYRCL */
  127. 	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, IYRCL);
  128. 	fprintf(stdout, "%.*s %4d\n", len - numlen, buff, IYRCL);
  129.  
  130. 	/* READ (IN, '(T30, I4)') POUT */
  131. 	_getline(buff, sizeof(buff), IN);
  132. 	numlen = 4; len = strlen(buff);
  133. 	POUT = atoi(&buff[len-numlen]);
  134. 	if (POUT == 1) {
  135. 		/* WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout if necessary', 'YES' */
  136. 		fprintf(IOUT,   "Apply phaseout if necessary    YES\n");
  137. 		fprintf(stdout, "Apply phaseout if necessary    YES\n");
  138. 	} else {
  139. 		/* WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout', 'No' */
  140. 		fprintf(IOUT,   "Apply phaseout                  No\n");
  141. 		fprintf(stdout, "Apply phaseout                  No\n");
  142. 	}
  143.  
  144. /* Re-scale IYRCL such that 0 is the year prior to the first input data */
  145.  
  146. 	md.ISTART = 0;
  147. 	md.IYEAR = IYRCL - INITYR;
  148. 	/* IF (IYEAR .LE. 0 .OR. IYEAR .GT. MAXYR) STOP 'INVALID YEAR' */
  149. 	if ((md.IYEAR <= 0) || (md.IYEAR > MAXYR)) {
  150. 		fprintf(stderr, "INVALID YEAR");
  151. 		ret = 1;
  152. 		goto exit;
  153. 	}
  154.  
  155. /* Initialize the catch array */
  156.  
  157. 	/* DO 10 I = 0, MAXYR */
  158. 	for (I = 0; I <= MAXYR; I++) {
  159. 		CATCH[I] = 0.0;
  160. 	} /* 10   CONTINUE */
  161.  
  162. /* Read in the catch data, scaling each year to the initial year */
  163. 	/* READ (IN, '(// A)') FORMT (I4,F8.0) */
  164. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  165. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  166. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  167.  
  168. 	/* WRITE (IOUT, '(/A)') 'Historic catches:' */
  169. 	fprintf(IOUT, "\nHistoric catches:\n");
  170. 	fprintf(stdout, "\nHistoric catches:\n");
  171. 	TOTALC = 0;
  172. 	/* DO 20 I = 0, MAXYR */
  173. 	for (I = 0; I <= MAXYR; I++) {
  174. 		/* READ (IN, FORMT) IY, C */ 
  175. 		_getline(buff, sizeof(buff), IN);
  176. 		sscanf(buff,"%4d%7lf", &IY, &C);
  177. 		/* IF (IY .LT. 0) GO TO 25 */
  178. 		if (IY < 0) break; /* EOD */
  179.  
  180. 		/* WRITE (IOUT, FORMT) IY, C */
  181. 		fprintf(stdout, "%4d%6.f\n", IY, C);
  182. 		fprintf(IOUT, "%4d%6.f\n", IY, C);
  183. 		IY = IY - INITYR;
  184. 		/* IF (IY .LT. 0 .OR. IY .GE. IYEAR) STOP */
  185. 		if (IY < 0 || IY > md.IYEAR) {
  186. 			fprintf(stderr, " ERROR: IY is out of range(%d)\n", IY);
  187. 			ret = 1;
  188. 			goto exit;
  189. 		}
  190. 		CATCH[IY] = C;
  191. 		TOTALC = TOTALC + C;
  192. 	} /* 20   CONTINUE */
  193.  
  194. /* 25   IF (TOTALC .LT. 0) STOP ' ERROR: No historic catch input' */
  195. 	if (TOTALC < 0) {
  196. 		fprintf(stderr, " ERROR: No historic catch input\n");
  197. 		ret = 1;
  198. 		goto exit;
  199. 	}
  200.  
  201. /* Read in the non-zero sightings estimates and information matrix */
  202.  
  203. 	/* READ (IN, '(//T30, I4 / A)') NS, FORMT (I4,F10.0,10F8.0) */
  204. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  205. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  206. 	_getline(buff, sizeof(buff), IN);
  207. 	numlen = 4; len = strlen(buff);
  208. 	md.NS = atoi(&buff[len-numlen]);
  209.  
  210. 	/* WRITE (IOUT, '(/A)') 'Abundance estimates:' */
  211. 	fprintf(IOUT,   "\nAbundance estimates:\n");
  212. 	fprintf(stdout, "\nAbundance estimates:\n");
  213. 	/* IF (NS .GT. MAXEST) STOP 'ERROR: Abundance year out of range' */
  214. 	if (md.NS > MAXEST) {
  215. 		fprintf(stderr, "ERROR: Abundance year out of range\n");
  216. 		ret = 1;
  217. 		goto exit;
  218. 	}
  219. 	_getline(buff, sizeof(buff), IN);
  220. 	/* DO 30 N=0, NS-1 */
  221. 	for (N = 0; N <= md.NS - 1; N++) {
  222. 		N1 = (N * (N + 1)) / 2;
  223. 		/* READ (IN, FORMT) ISYR(N), SIGHT(N), (FMATRX(N1 + J), J = 0, N) */
  224. 		_getline(buff, sizeof(buff), IN);
  225. 		sscanf(buff, "%4d%10lf", &ISYR[N], &SIGHT[N]);
  226. 		fprintf(IOUT,   "%4d %8.f", ISYR[N], SIGHT[N]);
  227. 		fprintf(stdout, "%4d %8.f", ISYR[N], SIGHT[N]);
  228. 		len = strlen(buff);
  229. 		for (J = 0; J < N + 1; J++) {
  230. 			if (len > J * 7 + 14) {
  231. 				sscanf(&buff[0] + J * 7 + 14,"%7lf", &FMATRX[N1 + J]);
  232. 			} else {
  233. 				FMATRX[N1 + J] = 0.0;
  234. 			}
  235. 			fprintf(IOUT,   "%7.f", FMATRX[N1 + J]);
  236. 			fprintf(stdout, "%7.f", FMATRX[N1 + J]);
  237. 		}
  238. 		fprintf(IOUT,   "\n");
  239. 		fprintf(stdout, "\n");
  240. 		/* WRITE (IOUT, FORMT) ISYR(N), SIGHT(N), (FMATRX (N1 + J), J = 0, N) */
  241. 		ISYR[N] = ISYR[N] - INITYR;
  242. 		if (ISYR[N] < 0 || ISYR[N] >= md.IYEAR) {
  243. 			/*STOP 'ERROR: Sight year out of range'*/
  244. 			fprintf(stderr, "ERROR: Sight year out of range\n");
  245. 			ret = 1;
  246. 			goto exit;
  247. 		}
  248. 		if (SIGHT[N] <= 0.0) {
  249. 			/* STOP ' ERROR: Estimate not positive' */
  250. 			fprintf(stderr, "ERROR: Estimate not positive\n");
  251. 			ret = 1;
  252. 			goto exit;
  253. 		}
  254. 	} /* 30   CONTINUE */
  255.  
  256. /* Read in the Poisson multipliers for any zero sightings estimates */
  257.  
  258. 	/* READ (IN, '(/T30, I3, /A)') NZ, FORMAT */
  259. 	_getline(buff, sizeof(buff), IN);
  260. 	_getline(buff, sizeof(buff), IN);
  261. 	numlen = 3; len = strlen(buff);
  262. 	md.NZ = atoi(&buff[len-numlen]);
  263. 	if (md.NZ > 0) {
  264. 		/* WRITE (IOUT, '(/A)') 'Zero abundance estimates:' */
  265. 		fprintf(IOUT,   "\nZero abundance estimates:\n");
  266. 		fprintf(stdout, "\nZero abundance estimates:\n");
  267. 	}
  268. 	if (md.NZ > MAXEST) {
  269. 		/* STOP ' ERROR: Zero estimate array too small' */
  270. 		fprintf(stderr, "ERROR: Zero estimate array too small\n");
  271. 		ret = 1;
  272. 		goto exit;
  273. 	}
  274. 	_getline(buff, sizeof(buff), IN);
  275. 	/* DO 40 N = 0, NZ - 1 */
  276. 	for (N = 0; N <= md.NZ - 1; N++) {
  277. 		N1 = (N * (N + 1)) / 2;
  278. 		/*READ (IN, FORMT) IZYR(N), ZMULT(N)*/
  279. 		_getline(buff, sizeof(buff), IN);
  280. 		sscanf(buff, "%4d%10lf", &IZYR[N], &ZMULT[N]);
  281. 		/* WRITE (IOUT, FORMT) IZYR(N), ZMULT(N) */
  282. 		fprintf(IOUT,   "%4d%7.f\n", IZYR[N], ZMULT[N]);
  283. 		fprintf(stdout, "%4d%7.f\n", IZYR[N], ZMULT[N]);
  284. 		IZYR[N] = IZYR[N] - INITYR;
  285. 		if (IZYR[N] < 0 || IZYR[N] >= md.IYEAR) {
  286. 			/* STOP ' Sight year out of range' */
  287. 			fprintf(stderr, "Sight year out of range\n");
  288. 			ret = 1;
  289. 			goto exit;
  290. 		}
  291.  
  292. 		if (ZMULT[N] <= 0.0) {
  293. 			/* STOP ' ERROR: Multiplier not positive' */
  294. 			fprintf(stderr, "ERROR: Multiplier not positive\n");
  295. 			ret = 1;
  296. 			goto exit;
  297. 		}
  298. 	} /* 40   CONTINUE */
  299. 	/*WRITE (IOUT, '()')*/
  300. 	fprintf(IOUT,"\n");
  301. 	fprintf(stdout,"\n");
  302. /* Bound the range for the integration over K */
  303.  
  304. 	md.RKHI = 1.0E7;
  305. 	md.RKLO = 0.0;
  306. 	IS = 0; /* orignal not initialized bug? */
  307. /* ****** */
  308. /* ****** Run the CLA to obtain the nominal catch limit */
  309. /* ****** */
  310. 	RAWCL = CLIMIT(CATCH, SIGHT, FMATRX, ISYR, IZYR, ZMULT, POP, G);
  311. /* Set the catch limits for PCYCLE years. If the catch limit may be */
  312. /* subject to phaseout, call POUT to apply the phaseout rule. */
  313.  
  314. /* First set ILAST = year of the most recent abundance estimate */
  315.  
  316. 	if (md.NS > 0) ILAST = ISYR[md.NS - 1];
  317. 	/* IF (NZ .GT. 0) ILAST = MAX(IS, IZYR(md.NZ - 1)) */
  318. 	if (md.NZ > 0) ILAST = (IS > IZYR[md.NZ - 1])? IS: IZYR[md.NZ - 1];
  319.  
  320. 	/* DO 100 IY = IYEAR, IYEAR + mp.PCYCLE - f1 */
  321. 	for (IY = md.IYEAR; IY <= md.IYEAR + mp.PCYCLE - 1; IY++) {
  322. 		if (POUT == 1) {
  323. 			PHOUT(&CL, RAWCL, ILAST, IY);
  324. 		} else {
  325. 			CL = RAWCL;
  326. 		}
  327. 		/* WRITE (IOUT, '(A6, I5, A17, I6)') 'Year:', IY + INITYR,'Catch limit:', NINT(CL) */
  328. 		fprintf(IOUT,   "Year: %4d     Catch limit:%6d\n", IY + INITYR, (int)(CL + 0.5));
  329. 		fprintf(stdout, "Year: %4d     Catch limit:%6d\n", IY + INITYR, (int)(CL + 0.5));
  330.  
  331. 	} /* 100  CONTINUE */
  332.  
  333. exit:
  334. 	if (IN != NULL)   fclose(IN);
  335. 	if (IOUT != NULL) fclose(IOUT);
  336. 	return ret;
  337. 	/*	STOP */
  338. 	/*	END  */
  339. }
  340.  
  341. /*
  342. C **************************************************************************************
  343. C
  344. C	CLC version 6
  345. C
  346. C **************************************************************************************
  347. C
  348. C	31 January 1994
  349. C
  350. C **************************************************************************************
  351. */
  352. /*
  353. 	SUBROUTINE RDPARS
  354. 	Read the file of input parameters
  355. */
  356. void RDPARS(void)
  357. {
  358. 	double *p[] = {
  359. 		&mp.PPROB, &mp.PYMAX, &mp.PNYSTP,&mp.PKSTEP, &mp.PDSTEP, &mp.PBMIN, &mp.PBMAX,
  360. 		&mp.PNBSTP, &mp.PSCALE, &mp.PHASET, &mp.PHASEP, &mp.PCYCLE, &mp.PLEVEL, &mp.PSLOPE
  361. 	};
  362. 	FILE *IFILE;
  363. 	char buff[80];
  364. 	int count;
  365.  
  366. 	/* Open the input file (only read once) */
  367. 	IFILE = fopen(CLCPAR, "rt");
  368. 	if (IFILE == NULL) {
  369. 		fprintf(stderr,"Can't open " CLCPAR "\n");
  370. 		exit(1);
  371. 	}
  372. 	count = 0;
  373. 	while (_getline(buff, sizeof(buff), IFILE) != 0) {
  374. 		/* READ (IFILE, '(T30, F10.0)') PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PBMIN, PBMAX,
  375. 								 PNBSTP, PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE */
  376. 		sscanf(buff,"%*28c%lf", p[count]);
  377. 		count++;
  378. 	}
  379. 	/* CLOSE (IFILE) */
  380. 	fclose(IFILE);
  381. 	return;
  382. 	/* RETURN
  383. 	   END */
  384. }
  385.  
  386. /*
  387. 	SUBROUTINE PHOUT (CL, RAWCL, ILAST, IY)
  388. */
  389. void PHOUT(double *CL, double RAWCL, int ILAST, int IY)
  390. {
  391. /*
  392. 	PHASET Number of years without surveys before phaseout invoked
  393. 	PHASEP Phaseout annual reduction proportion
  394. 	Phaseout: Reduce catch limit if there is no survey data in the
  395. 	last PHASET years
  396. */
  397. 	if (IY >= ILAST + mp.PHASET) {
  398. 		*CL = RAWCL * (1.0 - mp.PHASEP * (IY - ILAST - mp.PHASET));
  399. 		if (*CL < 0.0) *CL = 0.0;
  400. 	} else {
  401. 		*CL = RAWCL;
  402. 	}
  403. 	/* RETURN
  404. 	   END */
  405. }
  406.  
  407. /*
  408. 	REAL FUNCTION CLIMIT (CATCH,SIGHT,FMATRX,ISYR,IZYR,ZMULT,POP,G)
  409. 	Run the CLA to obtain to nominal catch limit
  410. */
  411. double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR, double *ZMULT,double *POP, double *G)
  412. {
  413. /*
  414. 	Local variables: 
  415.  
  416. 	REAL PRES(0:MAXSIM), QRES(0:MAXSIM), SS0, SS1, SS2, SS3
  417. 	REAL SF, Y(0:MAXSTP), B(0:MAXSTP), RLGB(0:MAXSTP)
  418. 	INTEGER NB, NR
  419. 	NB	Number of bias steps
  420. 	NR	Number of steps for Y (the productivity parameter)
  421. 	SF	Deviance scale factor (SF = .5 / PSCALE**2)
  422. */
  423. 	double SF, SS0, SS1, SS2, SS3;
  424. 	double Y[MAXSTP + 1], B[MAXSTP + 1], RLGB[MAXSTP + 1];
  425. 	double *PRES = NULL, *QRES = NULL;
  426.  
  427. 	int NB, NR;
  428. 	double BINC, YINC, DK, R, RK, D, DD, P, Q, PTOT;
  429. 	int I, J, N, N2;
  430.  
  431. 	double _CLIMIT = 0.0;
  432. 	if (md.NS <= 0) return _CLIMIT;
  433.  
  434. 	PRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
  435. 	if (PRES == NULL) {
  436. 		fprintf(stderr, "Can't Allocate memory (PRES)\n");
  437. 		goto err_exit;
  438. 	}
  439. 	memset(PRES, 0, sizeof(double) * (MAXSTP + 1));
  440.  
  441. 	QRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
  442. 	if (QRES == NULL) {
  443. 		fprintf(stderr, "Can't Allocate memory (QRES)\n");
  444. 		goto err_exit;
  445. 	}
  446. 	memset(QRES, 0, sizeof(double) * (MAXSTP + 1));
  447.  
  448. 	/* Set deviance scale factor S = 1/PSCALE**2 */
  449. 	SF = 0.5 / (mp.PSCALE * mp.PSCALE);
  450.  
  451. 	/* Check the sizes of the Y and B arrays are large enough */
  452. 	if (mp.PNBSTP > MAXSTP || mp.PNYSTP > MAXSTP) {
  453. 		fprintf(stderr, "Y &/or b array sizes too small\n");
  454. 		goto err_exit;
  455. 	}
  456. 	/* Set sightings bias step sizes & their log values. BINC = Bias increment */
  457. 	NB = (int)mp.PNBSTP;
  458. 	BINC = (mp.PBMAX - mp.PBMIN) / NB;
  459. 	/* DO 10 I = 0, NB - 1 */
  460. 	for (I = 0; I <= NB - 1; I++) {
  461.         *(B + I) = mp.PBMIN + (I + 0.5) * BINC;
  462.         *(RLGB + I) = -log(*(B + I));
  463. 	} /* 10   CONTINUE */
  464.  
  465. 	/* Set productivity parameter step sizes (midpoints) */
  466. 	NR = (int)mp.PNYSTP;
  467. 	YINC = mp.PYMAX / NR;
  468. 	/* DO 20 I = 0, NR - 1 */
  469. 	for (I = 0; I <= NR - 1; I++) {
  470. 		*(Y + I) = (I + 0.5) * YINC;
  471. 	}/* 20   CONTINUE */
  472.  
  473. 	PTOT = 0;
  474. 	N = 0;
  475. 	/* DO 50 I = 0, NR - 1 */
  476. 	for (I = 0; I <= NR - 1; I++) {
  477. 		/* Set R from the productivity parameter, Y */
  478. 		R = 1.4184 * *(Y + I);
  479.  
  480. 		/* Step size for K */
  481.         DK = mp.PKSTEP;
  482.         D = 1.0;
  483.         RK = md.RKHI;
  484.  
  485. 		/* Use function STKSIM to set up the Nth population trajectory */
  486. 		/* i.e. set up the pop array */
  487.  
  488. 		while (1) {
  489. 		/* 30     IF (RK .LE. RKLO .OR. STKSIM (RK, R, POP, CATCH) .LE. 0.) GOTO 40 */
  490. 			if ((RK <= md.RKLO) || (STKSIM(RK, R, POP, CATCH) <= 0.0)) break;
  491. 			/* IF (N .GE. MAXSIM) STOP 'ERROR: TOO MANY SIMULATIONS' */
  492. 			if (N >= MAXSIM) {
  493. 				fprintf(stderr, "ERROR: TOO MANY SIMULATIONS");
  494. 				goto err_exit;
  495. 			}
  496. 			/* How much depletion covered? */
  497. 			DD = D - POP[md.IYEAR] / RK;
  498. 			D = POP[md.IYEAR] / RK;
  499. 			P = 0.0;
  500.  
  501. 			if (DD > 0.0) {
  502. 				/* Compute the internal catch limit corresponding to D and Y(I) */
  503. 				QRES[N] = CONTRL(D, Y[I], mp.PLEVEL, mp.PSLOPE) * POP[md.IYEAR];
  504.  
  505. 				/* Calculate deviance */
  506.             	DEVIAN(&SS0, &SS1, &SS2, &SS3, SIGHT, FMATRX, ISYR,IZYR, ZMULT, POP, G);
  507.  
  508. 				/* Scale likelihood and integrate over values for the bias parameter */
  509. 				/* DO 35 J = 0, NB - 1 */
  510. 				for (J = 0; J <= NB - 1; J++) {
  511. 					P = P + exp(-SF * (SS0 + *(RLGB + J) * (SS1 + *(RLGB + J) * SS2) + SS3 * *(B + J)));
  512. 				} /* 35         CONTINUE */
  513.  
  514. 				/* Calculate the weight for this point (& total area under likelihood) */
  515. 				PRES[N] = P * DD;
  516. 				PTOT = PTOT + PRES[N];
  517.  
  518. 				/* Update counter */
  519. 				N = N + 1;
  520.  
  521. 				/* Find the next K */
  522. 				DK = DK * mp.PDSTEP / DD;
  523. 				/* IF (DK .GT. PKSTEP) DK = PKSTEP */
  524. 				if (DK > mp.PKSTEP) DK = mp.PKSTEP;
  525.  
  526. 			} else {
  527. 				/* IF DD = 0 change the step size only */
  528. 				DK = mp.PKSTEP;
  529. 			}
  530. 			/* Set the new value of K */
  531. 			RK = RK / (1.0 + DK);
  532. 		} /* GOTO 30 */ /* 40     CONTINUE */
  533. 	} /* 50   CONTINUE */
  534. 	/* IF (PTOT .LE. 0.) STOP 'ERROR: PROB INTEGRATES TO ZERO' */
  535. 	if (PTOT < 0.0) {
  536. 		fprintf(stderr,"ERROR: PROB INTEGRATES TO ZERO");
  537. 		goto err_exit;
  538. 	}
  539.  
  540. 	/* Sort the QRES and PRES arrays in ascending order of QRES */
  541. 	N2 = N;
  542. 	SORT(QRES, PRES, N2);
  543.  
  544. 	/* Normalize the relative likelihoods */
  545. 	/* DO 60 I = 0, N - 1 */
  546. 	for (I = 0; I <= N - 1; I++) {
  547. 		PRES[I] = PRES[I] / PTOT;
  548. 	} /* 60   CONTINUE */
  549.  
  550. 	/* Extract the desired probability level: the nominal catch limit (NCL) */
  551. 	/* is the lower PROB% of the distribution. */
  552. 	/* First calculate PRES(I), the probability that the NCL is between */
  553. 	/* QRES(I) & QRES(I + 1). */
  554.  
  555. 	P = 0;
  556. 	/* DO 70 I = 0, N - 1 */
  557. 	for (I = 0; I <= N - 1; I++) {
  558. 		P = P + PRES[I];
  559. 		/* IF (P .GT. PPROB) GOTO 80 */
  560. 		if (P > mp.PPROB) break;
  561. 	} /* 70   CONTINUE */
  562.  
  563. 	/* Interpolate to set the nominal catch limit */
  564. 	/* 80   IF (I .GE. N - 1) THEN */
  565. 	if  (I >= N - 1) {
  566.         Q = QRES[N - 1];
  567. 	} else {
  568.         Q = (QRES[I + 1] * (mp.PPROB - P + PRES[I]) + QRES[I] * (P - mp.PPROB)) / PRES[I];
  569. 	}
  570.  
  571. 	_CLIMIT = Q;
  572.  
  573. 	if (QRES != NULL) free(QRES);
  574. 	if (PRES != NULL) free(PRES);
  575. /*	if (Y != NULL) free(Y);
  576. 	if (B != NULL) free(B);
  577. 	if (RLGB != NULL) free(RLGB);*/
  578. 	return _CLIMIT;
  579.  
  580. err_exit:
  581. 	if (QRES != NULL) free(QRES);
  582. 	if (PRES != NULL) free(PRES);
  583. /*	if (Y != NULL) free(Y);
  584. 	if (B != NULL) free(B);
  585. 	if (RLGB != NULL) free(RLGB);*/
  586. 	exit(1);
  587. 	return 0.0; /* not reached */
  588. 	/* RETURN */
  589. 	/* END */
  590. }
  591.  
  592. /*
  593. 	FUNCTION CONTRL (D, Y, PLEVEL, PSLOPE)
  594. 	Catch control law
  595. */
  596. double CONTRL(double D, double Y, double PLEVEL, double PSLOPE)
  597. {
  598. 	double ret;
  599. 	if (D < PLEVEL) {
  600. 		ret = 0.0;
  601. 	} else if (D < 1.0) {
  602. 		ret = PSLOPE * Y * (D - PLEVEL);
  603. 	} else {
  604. 		ret = PSLOPE * Y * (1.0 - PLEVEL);
  605. 	}
  606. 	return ret;
  607. 	/* END */
  608. }
  609.  
  610. /*
  611. 	SUBROUTINE DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR,
  612. 	Calculate deviance (-2 log likelihood) in terms of coefficients for 
  613. 	the bias and log bias parameters
  614. */
  615. void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT,
  616.             double *FMATRX, int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G)
  617. {
  618. 	/* REAL SS0,SS1,SS2,SIGHT(0:*),FMATRX(0:*),ZMULT(0:*),POP(0:*),G(0:*) */
  619. 	/* INTEGER ISYR(0:*), IZYR(0:*) */
  620. 	int N, J, K;
  621.  
  622. 	*SS0 = 0.0;
  623. 	*SS1 = 0.0;
  624. 	*SS2 = 0.0;
  625. 	*SS3 = 0.0;
  626.  
  627. 	/* DO 100 N = 0, NS - 1 */
  628. 	for (N = 0; N <= md.NS - 1; N++) { 
  629.  
  630. 		*(G + N) = log(*(SIGHT + N) / *(POP+*(ISYR + N)));
  631. 		K = N * (N + 1) / 2;
  632.  
  633. 		/* DO 10 J = 0, N - 1 */
  634. 		for (J = 0; J <= N - 1; J++) {
  635. 			/*1st add non diagonal contributions (which are doubled up)*/
  636. 			*SS0 = *SS0 + 2.0 * *(G + J) * *(G + N) * *(FMATRX + K);
  637. 			*SS1 = *SS1 + 2.0 * (*(G + J) + *(G + N)) * *(FMATRX + K);
  638. 			*SS2 = *SS2 + *(FMATRX + K) + *(FMATRX + K);
  639. 			K = K + 1;
  640. 		} /* 10     CONTINUE*/
  641.  
  642. 		/* Now add diagnoal contribution */
  643.         *SS0 = *SS0 + *(G + N) * *(G + N) * *(FMATRX + K);
  644.         *SS1 = *SS1 + 2.0 * *(G + N) * *(FMATRX + K);
  645.         *SS2 = *SS2 + *(FMATRX + K);
  646.  
  647. 	} /* 100  CONTINUE */
  648.  
  649. 	/* Now do the zero estimates */
  650. 	/* DO 200 N = 0, NZ - 1 */
  651. 	for (N =0; N <= md.NZ - 1; N++) {
  652. 		*SS3 = *SS3 + 2.0 * *(POP + *(IZYR + N)) / *(ZMULT + N);
  653.  	} /* 200  CONTINUE */
  654.  
  655. 	/* RETURN */
  656. 	/* END    */
  657. 	return;
  658. }
  659. /*
  660. 	FUNCTION STKSIM (RK, R, POP, CATCH)
  661. 	Calculate the stock trajectory with parameter RK and R: return
  662. 	the current stock size
  663. 	RK = notional carrying capacity; R = productivity parameter * 1.4184 
  664. */
  665. double STKSIM(double RK, double R, double *POP, double *CATCH)
  666. {
  667. 	/*INTEGER ISTART, IYEAR, NS, NZ, RKLO, RKHI */
  668. 	/* REAL CATCH(0:*), POP(0:*) */
  669. 	int I;
  670. 	double D;
  671. 	*(POP+md.ISTART) = RK;
  672. 	/*DO 10 I = md.ISTART + 1, md.IYEAR*/
  673. 	for (I = md.ISTART + 1; I <= md.IYEAR; I++) {
  674. 		D = *(POP + I - 1) / RK;
  675. 		*(POP + I) = *(POP + I - 1) * (1.0 + R * (1.0 - D * D)) - *(CATCH + I -1);
  676. 		if (*(POP + I) <= 0) {	/*GOTO 20*/
  677. 			return 0.0;
  678. 			/* 20 STKSIM = 0. */
  679. 		}
  680. 	} /* 10   CONTINUE */
  681. 	/* STKSIM = POP(md.IYEAR) */
  682. 	return *(POP + md.IYEAR);
  683. 	/* RETURN */
  684.     /* END */
  685. }
  686. /* 
  687. 	SUBROUTINE SORT (ARRAY, ARRAY2, N)
  688. 	SORT sorts a pair of arrays in ascending order of the first array
  689.    (using the heapsort algorithm)
  690. */
  691. void SORT(double *ARRAY, double *ARRAY2, int N)
  692. {
  693. 	double /*ARRAY(0:*), ARRAY2(0:*),*/ TEMP, TEMP2;
  694. 	int  /*N,*/ K, IR, I, J;
  695. 	if (N < 2) return;
  696. 	K = N / 2;
  697. 	IR = N - 1;
  698. 	while (1) { /*10*/
  699. 		if (K != 0) {
  700. 			K = K - 1;
  701. 			TEMP = *(ARRAY + K);
  702. 			TEMP2 = *(ARRAY2 + K);
  703. 		} else {
  704. 			TEMP = *(ARRAY + IR);
  705. 			TEMP2 = *(ARRAY2 + IR);
  706. 			*(ARRAY + IR) = *ARRAY;
  707. 			*(ARRAY2 + IR) = *ARRAY2;
  708. 			IR = IR - 1;
  709. 			if (IR == 0) {
  710. 				*ARRAY = TEMP;
  711. 				*ARRAY2 = TEMP2;
  712. 				return;
  713. 			}
  714. 		}
  715. 		I = K;
  716. 		J = K + K + 1;
  717. 		while (1) {/*20*/
  718. 			if (J > IR) break;
  719. 			if ((J < IR) && (*(ARRAY + J) < *(ARRAY+ J + 1))) J = J + 1;
  720. 			if (TEMP < *(ARRAY + J)) {
  721. 				*(ARRAY + I) = *(ARRAY + J);
  722. 				*(ARRAY2 + I) = *(ARRAY2 + J);
  723. 				I = J;
  724. 				J = J + I + 1;
  725. 			} else {
  726. 				J = IR + 1;
  727. 			}
  728. 		} /*GOTO 20*/
  729. 		*(ARRAY + I) = TEMP;
  730. 		*(ARRAY2 + I) = TEMP2;
  731. 	}/*GOTO 10*/
  732. 	/*END*/
  733. }
  734. int _getline(char *buff, size_t size, FILE *fp){
  735. 	char *ret;
  736. 	ret = fgets(buff, size, fp);
  737. 	if (ret != NULL) {
  738. 		if (buff[strlen(buff) - 1] == '\n')
  739. 			buff[strlen(buff) - 1] = 0;
  740. 		if (buff[strlen(buff) - 1] == '\r')
  741. 			buff[strlen(buff) - 1] = 0;
  742. 	}
  743. 	return (ret != NULL);
  744. }
  745.  
Runtime error #stdin #stdout #stderr 0s 9432KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
Can't open CLC.DAT
https://ideone.com/L5uPBe
language:
C (gcc 8.3)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!