C*****FILE PVAX06.F C SIM APPLICATION FOR STATISTICS OF POPULATION TRAJECTORY C.....NOMINALLY COMPLETE PVA PROJECTION, SAMPLING JOINT DISTRIBUTION C OF PARAMETER UNCERTAINTY), WITH C POPULATION MODELED AS C * NO SPATIAL STRUCTURE C * NO AGE STRUCTURE C * NO EXPLICIT ENVIRONMENTAL COVARIATE C * DENSITY DEPENDENCE C ** HARD POPULATION CEILING C * DEMOGRAPHIC STOCHASTICITY C ** SPECIFIED DIRECTLY AS ADDITIVE COMPONENT OF C TOTAL VARIANCE IN FACTOR OF INCREASE (LAMBDA) C * UNATTRIBUTED ENVIRONMENTAL VARIATION C ** LOGNORMAL DISTRIBUTION OF FACTOR OF INCREASE C * PARAMETER UNCERTAINTY C ** JOINT (POSTERIOR) DISTRIBUTION FOR MEAN AND STANDARD C DEVIATION OF ENVIRONMENTAL VARIATION (LOGNORMAL) C OF FACTOR OF INCREASE C * QUASI-EXTINCTION C ** FIRST PASSAGE TIME REPORTED AND TERMINATES TRAJECTORY C * SUB TIME-STEP PROCESSES C ** NONE; FACTOR OF INCREASE IS NOT DECOMPOSED INTO C BIRTHS AND DEATHS C * NUMERICAL C ** POPULATION STATE IS REPRESENTED AS INTEGER COUNT C ** DYNAMICS ARE ROUNDED (NOT TRUNCATED) BACK TO INTEGER C ** TRUNCATE POPULATION AT 10E+6 TO AVOID INTEGER OVERFLOW C.....EACH CALL TO SIMSM SIMULATES ONE COMPLETE TRAJECTORY C.....REPORTS AFTER IT TIME STEPS C PRV(1) PROBABILITY OF (QUASI-)EXTINCTION WITHIN TIME HORIZON C (REPORTED AS MEAN IN MARGINAL SUMMARY REPORT; C IGNORE HISTOGRAM--SUPPRESS WITH 0 BINS) C PRV(2) TIME TO FIRST PASSAGE, FROM ABOVE, OF QUASI-EXTINCTION C THRESHOLD, FOR TRAJECTORIES WHICH ACHIEVE FIRST PASSAGE C (CALCULATION WILL BE SLOW IF PROBABILITY OF C (QUASI-)EXTINCTION IS LOW; SUPPRESS BY SETTING SWITCH C TO 1, OR 0 TO SUPPRESS BOTH THIS AND TERMINAL POPULATION, C AND THEN SUPPRESS HISTOGRAM BY 0 BINS) C PRV(3) TERMINAL POPULATION SIZE FOR TRAJECTORIES WHICH DID NOT C ACHIEVE FIRST PASSAGE) C (CALCULATION WILL BE SLOW IF PROBABILITY OF C (QUASI-)EXTINCTION IS HIGH; SUPPRESS BY SETTING SWITCH C TO 2, OR 0 TO SUPPRESS BOTH THIS AND TIME TO C (QUASI-)EXTINCTION AND THEN SUPPRESS HISTOGRAM BY 0 BINS) C.....RECEIVE FROM JOB FILE C SPV(1) N1 INITIAL POPULATION C SPV(2) IT TIME HORIZON C SPV(3) NQ (QUASI-)EXTINCTION THRESHOLD (MAY BE 0) C SPV(4) NKH HARD POPULATION CEILING C SPV(5) DS STD OF DEMOGRAPHIC STOCHASTICITY IN LAMBDA AT C ONE INDIVIDUAL C SPV(6) ISW CUMULATION SWITCH C 0 SUPPRESSES CUMULATION OF TIME TO (QUASI-)EXTINCTION C AND TERMINAL POPULATION SIZE C 1 SUPPRESSES CUMULATION OF TIME TO (QUASI-)EXTINCTION C BUT NOT TERMINAL POPULATION SIZE C 2 SUPPRESSES CUMULATION OF TERMINAL POPULATION SIZE C BUT NOT TIME TO (QUASI-)EXTINCTION C 3 CUMULATES BOTH TIME TO (QUASI-)EXTINCTION C AND TERMINAL POPULATION SIZE C OTHER SPECIFICATIONS C EXPECTS A FILE NAMED PVAX05.SMP WHICH WILL BE A (POSSIBLY RENAMED) C OUTPUT SMP FILE (FROM MTG RUNNING APPLICATION SPECIFIC FILE C RPVAX05) CONTAINING JOINT POSTERIOR SAMPLE OF C UXF MEAN OF LOGNORMAL DISTRIBUTION OF ENVIRONMENTAL C VARIATION IN LAMBDA C SXF SXB OF LOGNORMAL DISTRIBUTION OF ENVIRONMENTAL C VARIATION IN LAMBDA C WHERE UXF IS STORED AS PARAMETER 1 AND SXF C IS STORED AS PARAMETER 2 C THIS SMP FILE IS NOT NAMED IN THE JOB FILE, AND READING IT C IS NOT MANAGED BY THE SIM SHELL; IT IS MANAGED AND READ BY C CALLS TO SUBROUTINES SMPFI AND SMPFR IN THE CODE BELOW. C***** C VERSION OF MAR 16, 2011 C***** SUBROUTINE SIMIN(LUNS,LUNF,NDT,MDT,DTV,SPV,CFNAM,FWV,IWV) C*****JOB-SPECIFIC DEDICATED SUBROUTINE THAT INITIALIZES VALUES C FOR PROGRAM SIM C.....GIVEN C LUNS DEVICE NUMBER FOR SCREEN OUTPUT, PRESUMED OPEN C LUNF DEVICE NUMBER FOR OUTPUT FILE, PRESUMED OPEN C NDT NUMBER OF OBSERVATIONS IN DATA C MDT NUMBER OF VARIABLES PER OBSERVATION IN DATA C DTV FLOATING POINT ARRAY OF DATA ELEMENTS, IN SINGLE SUBSCRIPT C STORAGE, WHERE ELEMENT IJ=(I-1)*MDT+J IS THE J-TH VARIABLE C FROM OBSERVATION I. C SPV VECTOR OF FIXED PARAMETER VALUES READ FROM JOB FILE C.....RETURNS C CFNAM CHARACTER*8 STRING IDENTIFYING THIS SOURCE FILE, TO MATCH C AGAINST CODE NAME SPECIFIED IN JOB FILE C DTV POSSIBLY REARRANGED, TRANSFORMED, SUBSET OR AUGMENTED; WILL C BE COMMUNICATED TO SUBROUTINE SIMSM IN THIS FORM C SPV POSSIBLY REARRANGED, TRANSFORMED, SUBSET OR AUGMENTED; WILL C BE COMMUNICATED TO SUBROUTINE SIMSM IN THIS FORM C FWV VECTOR AVAILABLE FOR FLOATING POINT WORKSPACE AND COMMUNICATING C ADDITIONAL FLOATING POINT VALUES TO SUBROUTINE SIMSM C IWV VECTOR AVAILABLE FOR INTEGER WORKSPACE AND COMMUNICATING C INTEGER VALUES TO SUBROUTINE SIMSM C OPTIONALLY MAY WRITE ANY DESIRED MESSAGE TO SCREEN (UNIT LUNS) OR C THE OUTPUT FILE (UNIT LUNF) C***** CHARACTER*8 CFNAM DIMENSION DTV(1),SPV(1),FWV(1),IWV(1) C========BEGIN BLOCK OF JOB-SPECIFIC CODE=============================== C.....SELF-IDENTIFY WITH CHARACTER STRING OF UP TO 8 CHARACTERS C IN SINGLE QUOTES (MANDATORY). C DO NOT CHANGE THE NAME OF THE VARIABLE CFNAM C..... CFNAM='PVAX06 ' C.....CONVERT INITIAL POPULATION TO INTEGER, LOAD IN IWV(1) EN1=SPV(1) N1=IFIX(EN1+0.01) IWV(1)=N1 C.....CONVERT TIME HORIZON TO INTEGER, LOAD IN IWV(2) T=SPV(2) IT=IFIX(T+0.01) IWV(2)=IT C.....CONVERT THRESHOLD TO INTEGER, LOAD IN IWV(3) ENQ=SPV(3) NQ=IFIX(ENQ+0.01) IWV(3)=NQ C.....CONVERT HARD POPULATION CEILNG TO INTEGER, LOAD IN IWV(5) FNKH=SPV(4) NKH=IFIX(FNKH+0.01) IWV(5)=NKH C.....CONVERT THE DEMOGRAPHIC STOCHASTICITY STD TO VARIANCE, IN PLACE DS=SPV(5) DS2=DS*DS SPV(5)=DS2 C.....CONVERT SWITCH TO INTEGER, RECODE, AND LOAD IN IWV(4) SW=SPV(6) ISW=IFIX(SW+0.01)-1 IWV(4)=ISW C.....SET UP FOR READING THE SMP FILE LUNC=5 LUNP=8 CALL SMPFI(LUNS,LUNC,LUNF,LUNP,CFNAM,NPT,MPS,NHPS,IWV(11),NSMP, * FWV) IWV(6)=NPT IWV(7)=MPS IWV(8)=NHPS IWV(9)=NSMP C========END THIS BLOCK OF JOB SPECIFIC CODE============================ 9999 RETURN END C----- SUBROUTINE SIMSM(NDT,MDT,DTV,SPV,FWV,IWV,SMV,ISD,PRV) C*****JOB SPECIFIC DEDICATED PROGRAM TO SAMPLE PROCESS FOR PROGRAM SIM C.....GIVEN C NDT NUMBER OF OBSERVATIONS IN DATA C MDT NUMBER OF VARIABLES PER OBSERVATION IN DATA C DTV FLOATING POINT DATA AS POSSIBLY REARRANGED, TRANSFORMED, C SUBSET OR AUGMENTED BY JOB-SPECIFIC SUBROUTINE SIMIN C SPV VECTOR OF FIXED PARAMETER VALUES, AS SPECIFIED IN THE C JOB FILE AND POSSIBLY MODIFIED IN SIMIN C FWV VECTOR OF FLOATING WORKSPACE, AND POSSIBLE FLOATING POINT C VALUES COMMUNICATED FROM SUBROUTINE SIMIN C IWV VECTOR OF INTEGER WORKSPACE AND POSSIBLE INTEGER VALUES C COMMUNICATED FROM SUBROUTINE SIMIN C SMV FLOATING POINT ARRAY OF ONE NEW SAMPLE FROM THE SMP FILE C INDEXED BY THE PARAMETER DESIGNATION IN THE FILE HEADER C (I.E., THE PARAMETER DESIGNATION IN THE PROGRAM THAT GENERATED C THE SMP FILE) C ISD SEED VECTOR FOR RANDOM NUMBER GENERATOR (DO NOT MODIFY) C.....RETURNS C PRV VECTOR OF VALUES (UP TO 3 ELEMENTS) GENERATED ON THIS C REALIZATION FOR CUMULATION OF JOINT DISTRIBUTION C***** DIMENSION DTV(1),SPV(1),FWV(1),IWV(1),SMV(1),ISD(1),PRV(1) C========BEGIN BLOCK OF JOB-SPECIFIC CODE=============================== C.....UNLOAD JOB SPECIFICATIONS: C INITIAL POPULATION C..... N1=IWV(1) C.....TIME HORIZON IT=IWV(2) C.....POPULATION THRESHOLD NQ=IWV(3) C.....HARD POPULATION CEILING NKH=IWV(5) C.....DEMOGRAPHIC STOCHASTICITY VARIANCE COMPONENT DS2=SPV(5) C.....RETRIAL SWITCH ISW=IWV(4) C.....DOWNLOAD OF JOB SPECIFICATIONS COMPLETE: C UNLOAD SMP FILE CONTROL SPECIFICATIONS C..... LUNP=8 NPT=IWV(6) MPS=IWV(7) NHPS=IWV(8) NSMP=IWV(9) C.....INITIALIZE TRIAL OUTCOMES FLAG NTRJ=0 IX=0 IH=0 C.....START A TRIAL TRAJECTORY C INITIALIZE FIRST PASSAGE FLAG, AT NOT C..... 1000 IFP=0 C.....INITIALIZE TRAJECTORY N=N1 C.....PARAMETER UNCERTAINTY C SAMPLE FROM THE SMP FILE JOINT VALUES OF THE MEAN AND THE STD OF C THE (LOGNORMAL) DISTRIBUTION OF ENVIRONMENTAL VARIATION IN LAMBA, C FOR THIS TRAJECTORY C..... CALL SMPFR(LUNP,NPT,MPS,NHPS,IWV(11),NSMP,FWV,SMV) IWV(9)=NSMP UXF=SMV(1) SXF=SMV(2) C.....LOOP THROUGH TIME STEPS DO 1140 I=1,IT C.....ADD THE DEMOGRAPHIC STOCHASTICITY VARIANCE COMPONENT TO ENVIRONMENTAL FN=FLOAT(N) DV=SXF*SXF+DS2/FN SXF=SQRT(DV) C.....LOG SPACE PARAMETERS FOR THE TOTAL VARIATION IN LAMBDA CALL LNSPC(UXF,SXF,ULF,SLF) C.....SAMPLE LOGNORMAL FOR LAMBDA THIS TIME STEP CALL LNORG(ULF,SLF,FA,ISD) C.....APPLY FACTOR OF INCREASE IN FLOATING POINT TO ADVANCE POPULATION FN=FA*FN C.....ROUND TO INTEGER N=IFIX(FN+0.5) C.....ENFORCE HARD CEILING IF (N-NKH) 1100,1100,1010 1010 N=NKH C.....CHECK FOR FIRST PASSAGE 1100 IF (N-NQ) 1110,1110,1120 1110 IFP=I GO TO 1150 C.....CHECK FOR APPROACHING OVERFLOW 10,000,000 1120 IF (N-10000000) 1140,1140,1130 1130 N=10000000 GO TO 1150 C.....TIME STEP COMPLETE 1140 CONTINUE C.....TRAJECTORY COMPLETE: C.....CHECK FOR FIRST TRIAL 1150 IF (NTRJ) 9999,1160,1190 1160 NTRJ=1 IF (IFP) 9999,1170,1180 1170 PRV(1)=0.0 GO TO 1190 1180 PRV(1)=1.0 C.....CHECK RETRIAL SWITCH C RECODED TO: -1 SUPPRESSES TIME AND POPULATION SIZE C 0 SUPPRESSES TIME ONLY C 1 SUPPRESSES POPULATION SIZE ONLY C 2 NONE SUPPRESSED C..... 1190 IF (ISW) 1200,1210,1230 C.....TIME AND POPULATION SIZE SUPPRESSED 1200 PRV(2)=-9.9 PRV(3)=-9.9 GO TO 9999 C.....TIME ONLY SUPPRESSED 1210 IF (IFP) 9999,1220,1000 1220 PRV(2)=-9.9 PRV(3)=FLOAT(N) GO TO 9999 C.....POPULATION OR NEITHER SUPPRESSED 1230 IF (ISW-1) 9999,1240,1260 C.....POPULATION SUPPRESSED 1240 IF (IFP) 9999,1000,1250 1250 PRV(2)=FLOAT(IFP) PRV(3)=-9.9 GO TO 9999 C.....NEITHER SUPPRESSED 1260 IF (IX) 9999,1270,1290 1270 IF (IFP) 9999,1290,1280 1280 IX=1 PRV(2)=FLOAT(IFP) IF (IH) 9999,1300,9999 1290 IF (IH) 9999,1300,1000 1300 IF (IFP) 9999,1310,1000 1310 IH=1 PRV(3)=FLOAT(N) IF (IX) 9999,1000,9999 C========END THIS BLOCK OF JOB SPECIFIC CODE============================ 9999 RETURN END C***** C END OF FILE C*****