C-------------------------------- DEMO -----------------------------------------
C
C Note: the original code source has been lost, this one was typed again by
C       JRL on 03-DEC-2006 from an hard copy of the article:
C
C       @article{LobryJR1991,
C         Author = {Lobry, J.R. and Rosso, L. and Flandrois, J.-P.}, 
C         Title = {A {FORTRAN} subroutine for the determination of parameter 
C                 confidence limits in non-linear models},
C         Journal = {Binary},
C         Volume = {3},
C         Pages = {86-93},
C         Year = {1991} }
C
C This program compiles for me without warning with g77 :
C GNU Fortran (GCC) 3.4.2
C
C VARIABLES
C----------
C
C SSR    NAME OF AN USER SUPPLIED SUBROUTINE WHICH RETURNS THE SUM OF
C        SQUARED RESIDUALS BETWEEN OBSERVED AND THEORETICAL VALUES
C
C PN     THE NUMBER OF PARAMETERS IN THE MODEL
C
C POPT   VECTOR WHICH CONTAINS THE THE PARAMETER VALUES SUCH THAT
C        THE SSR IS MINIMUM
C
C PINF   VECTOR WHICH CONTAINS THE PARAMETER LOWER BOUNDS TO
C        GENERATE STARTING POINTS
C
C PSUP   VECTOR WHICH CONTAINS THE PARAMETER UPPER BOUNDS TO
C        GENERATE STARTING POINTS
C
C THRES  THRESHOL VALUE FOR THE SSR WHICH MAY BE OBTAINED FROM:
C
C        THRES = SSR(POPT)*(1+PN*F(alpha;PN,N-PN)/(N-PN)),
C
C        WHERE N IS THE NUMBER OF POINTS IN THE DATA SET.
C        ACCORDING TO BEALE (1960):
C
C        @article{BealeEML1960,
C          Author = {Beale, E.M.L.},
C          Title = {Confidence regions in non-linear estimation},
C          Journal = {Journal of the Royal Statistical Society},
C          Year = {1960},
C          Volume = {22B},
C          Pages = {41-88} }
C
C MAXPT  MAXIMUM NUMBER OF POINTS TO BE GENERATED
C
C SEED   AN INITIAL CONDITION FOR THE PSEUDO-RANDOM NUMBERS
C        GENERATOR
C
C FILN   THE NAME OF THE OUTPUT FILE TO WRITE COORDINATES OF POINTS
C
C EPS    DETERMINES THE ACCURACY OF RESULTS DURING CONVERGENCE TO THE
C        THRESHOLD, THE STOP CRITERION IS THAT (THRES-SSR)/THRES < EPS
C
C MAXSSR A STARTING POINT IS REJECTED WHEN THE RATIO OF THE SSR AT THIS
C        POINT TO THE THRESHOLD VALUE IS GREATER THAN MAXSSR
C        (SSR/THRES > MAXSSR). THIS IS TO AVOID STARTING FROM A POINT TOO
C        FAR FROM THE OPTIMAL POINT: IN THIS CASE CONVERGENCE IS TOO
C        LONG AND SUBJECT TO FAILURE DUE TO NUMERICAL IMPRECISION
C
C MAXIT  MAXIMUM NUMBER OF ITERATIONS ALLOWED FOR CONVERGENCE TO
C        THRESHOLD
C
C-------------------------------------------------------------------------------
      PROGRAM DEMO
      EXTERNAL SSR
      REAL*8 SSR,POPT(8),PINF(8),PSUP(8),THRES,SEED,EPS,MAXSSR
      INTEGER PN,MAXPT,MAXIT
      CHARACTER*30 FILN
C
C INITIALISATION OF VARIABLES
C
      PN=4
C POPT(1): MuOpt
      PINF(1)=1.3
      POPT(1)=1.3915
      PSUP(1)=1.5
C POPT(2): Tmin
      PINF(2)=285.5
      POPT(2)=289.398
      PSUP(2)=292.0
C POPT(3): Topt
      PINF(3)=311.8
      POPT(3)=313.25
      PSUP(3)=314.5
C POPT(4): Tmax
      PINF(4)=319.9
      POPT(4)=320.234
      PSUP(4)=320.8
      
      THRES=0.0395811*(1+(4.0*3.36)/(15.0-4.0))
      MAXPT=10000
      SEED=0.1071966
      FILN="OUTDEMO"
      EPS=0.0001
      MAXSSR=1.0E+05
      MAXIT=500
C
C CALLING SUBROUTINE CRUNCH:
C
      CALL CRUNCH(SSR,PN,POPT,PINF,PSUP,THRES,
     &       MAXPT,SEED,FILN,EPS,MAXSSR,MAXIT)
      END
C----------------------------------- SSR --------------------------------------
C
C PURPOSE
C -------
C THIS FUNCTION RETURNS THE SUM OF SQUARED RESIDUALS BETWEEN OBSERVED
C AND THEORETICAL VALUES. FOR CLARITY, EXPERIMENTAL DATA HAVE BEEN
C INCLUDED DIRECTLY INSIDE THE PROGRAM. IT IS MORE USUAL, HOWEVER, TO
C READ DATA FROM A FILE AT THE LEVEL OF THE MAIN PROGRAM AND THEN TO
C PASS VALUES TO THE SUBROUTINE SSR VIA COMMON VARIABLES.
C
C Note added in Dec 2006: the included dataset of specific growth rate
C versus temperature for Escherichia coli was deduced from figure 2
C in the paper:
C
C @article{RatkowskyDA1983,
C   Author = {Ratkowsky, D.A. and Lowry, R.K. and McMeekin, T.A. and 
C     Stokes, A.N. and Chandler, R.E.},
C   Title = {Model for bacterial culture growth rate throughout the entire
C     biokinetic range},
C   Journal = {Journal of Bacteriology},
C   Volume = {154},
C   Pages = {1222-1226},
C Year = {1983} }
C
C COMPUTATION OF THE SSR IN THE CASE OF THE CARDINAL TEMPERATURE MODEL
C T1 = Muopt, T2 = Tmin, T3 = Topt, T4 = Tmax
C
C VARIABLES
C ---------
C
C T   INPUT VECTOR WHICH CONTAINS PARAMETER VALUES
C
C P   INPUT VARIABLE WHICH CONTAINS THE NUMBER OF PARAMETERS IN THE MODEL
C
C-------------------------------------------------------------------------------
      REAL*8 FUNCTION SSR(T,P)
      INTEGER P
      REAL*8 T(8)
      REAL*8 A,B,T1,T2,T3,T4
      T1=T(1)
      T2=T(2)
      T3=T(3)
      T4=T(4)
      A=T2+T4
      B=T2*T4
      SSR=
     &(0.25-T1*(1-((294-T3)**2)/((294-T3)**2+294*(A-294)-B)))**2+
     &(0.56-T1*(1-((296-T3)**2)/((296-T3)**2+296*(A-296)-B)))**2+
     &(0.61-T1*(1-((298-T3)**2)/((298-T3)**2+298*(A-298)-B)))**2+
     &(0.79-T1*(1-((300-T3)**2)/((300-T3)**2+300*(A-300)-B)))**2+
     &(0.94-T1*(1-((302-T3)**2)/((302-T3)**2+302*(A-302)-B)))**2+
     &(1.04-T1*(1-((304-T3)**2)/((304-T3)**2+304*(A-304)-B)))**2+
     &(1.16-T1*(1-((306-T3)**2)/((306-T3)**2+306*(A-306)-B)))**2+
     &(1.23-T1*(1-((308-T3)**2)/((308-T3)**2+308*(A-308)-B)))**2+
     &(1.36-T1*(1-((310-T3)**2)/((310-T3)**2+310*(A-310)-B)))**2+
     &(1.32-T1*(1-((312-T3)**2)/((312-T3)**2+312*(A-312)-B)))**2+
     &(1.36-T1*(1-((314-T3)**2)/((314-T3)**2+314*(A-314)-B)))**2+
     &(1.34-T1*(1-((316-T3)**2)/((316-T3)**2+316*(A-316)-B)))**2+
     &(0.96-T1*(1-((318-T3)**2)/((318-T3)**2+318*(A-318)-B)))**2+
     &(0.83-T1*(1-((319-T3)**2)/((319-T3)**2+319*(A-319)-B)))**2+
     &(0.16-T1*(1-((320-T3)**2)/((320-T3)**2+320*(A-320)-B)))**2
      END
C--------------------------- CRUNCH -------------------------------------------
C
C PURPOSE
C ------
C THE SUBROUTINE CRUNCH WRITES COORDINATES OF POINTS SUCH THAT
C THE SSR FUNCTION EQUALS THE THRESHOLD VALUE IN THE OUTPUT FILE
C
C------------------------------------------------------------------------------
      SUBROUTINE CRUNCH(SSR,PN,POPT,PINF,PSUP,THRES,
     &       MAXPT,SEED,FILN,EPS,MAXSSR,MAXIT)
      EXTERNAL SSR
      REAL*8 SSR,POPT(8),PINF(8),PSUP(8),THRES,SEED,EPS,MAXSSR
      INTEGER PN,MAXPT,MAXIT
      CHARACTER*30 FILN
C
C DECLARATION OF LOCAL VARIABLES
C
      REAL*8 COOR(8),TETA(8),DPINF(8),DPSUP(8),COORSV(8),TETASV(8)
      REAL*8 LOWER,UPPER,SSRLOW,SSRUP,NEW,SSRNEW
      INTEGER ITERNBR
      LOGICAL OK
C
C CHECKING POPT VALUES
C
      IF(SSR(POPT,PN).GT.THRES)THEN
        WRITE(*,*)"ERROR FROM ROUTINE CRUNCH: SSR IS GREATER"
        WRITE(*,*)"THAN THRESHOLD WITH OPTIMAL PARAMETER VALUES"
        WRITE(*,*)"THRESHOLD = ",THRES
        WRITE(*,*)"SSR = ",SSR(POPT,PN)
        STOP
      END IF
C
C OPENING OUTPUT FILE
C
      OPEN(UNIT=34,FILE=FILN,STATUS="NEW")
C
C INTERMEDIARY COMPUTATIONS
C
      DO 50,I=1,PN
        DPINF(I)=(POPT(I)-PINF(I))*SQRT(DBLE(PN))
        DPSUP(I)=(PSUP(I)-POPT(I))*SQRT(DBLE(PN))
50    CONTINUE
C
C MAIN LOOP OF THE PROGRAM.
C REPETITION OF THE SAME WITH MAXPT DIFFERENT STARTING POINTS
C 
      I=1
1     I=I+1
C
C CHOOSING A POINT AT RANDOM AT THE SURFACE OF A HYPERSHERE
C OF RADIUS 1.0 WITH ITS CENTER AT THE ORIGIN:
C
      CALL SPHERE(PN,COOR,SEED)
C
C COMPUTATION OF THE CORRESPONDING PARAMETER VALUES (SCALING)
C
      DO 110,J=1,PN
        IF(COOR(J).GT.0.0)THEN
          TETA(J)=POPT(J)+COOR(J)*DPSUP(J)
        ELSE
          TETA(J)=POPT(J)+COOR(J)*DPSUP(J)
        END IF
110   CONTINUE
      SSRUP=SSR(TETA,PN)
      OK=.TRUE.
C
C TESTING IF THE POINT IS OUTSIDE THE CONFIDENCE REGION:
C
      IF(SSRUP.LT.THRES)THEN
        WRITE(*,*)"WARNING 1 FROM SUBROUTINE CRUNCH"
        WRITE(*,*)"SSR LES THAN THRESHOLD WITH:"
        WRITE(*,*)"PARAM VALUE  % DISTANCE FROM OPTIMAL PARAM VALUE"
        DO 120,K=1,PN
          WRITE(*,1000)TETA(K),100*COOR(K)*SQRT(DBLE(PN))
120     CONTINUE
        WRITE(*,*)"-----------------------------------------"
        OK=.FALSE.
        I=I-1
      END IF
C
C TESTING IF SSR/THRES IS NOT TOO BIG:
C
      IF(SSRUP/THRES.GT.MAXSSR)THEN
        WRITE(*,*)"WARNING 2 FROM SUBROUTINE CRUNCH"
        WRITE(*,*)"SSR/THRES GREATER THAN:",MAXSSR
        WRITE(*,*)"PARAM VALUE  % DISTANCE FROM OPTIMAL PARAM VALUE"
        DO 125,K=1,PN
          WRITE(*,1000)TETA(K),100*COOR(K)*SQRT(DBLE(PN))
125     CONTINUE
        WRITE(*,*)"-----------------------------------------"
        OK=.FALSE.
        I=I-1
      END IF
C
C SAVING COORDINATES
C
      DO 115,K=1,PN
        TETASV(K)=TETA(K)
        COORSV(K)=COOR(K)
115   CONTINUE
C
C CONVERGENCE TO THRESHOLD, INITIALIZATION:
C
      IF(OK)THEN
        LOWER=0.0
        UPPER=1.0
        DO 130,J=1,PN
          TETA(J)=POPT(J)
130     CONTINUE
        SSRLOW=SSR(TETA,PN)
        ITERNBR=0
C
C CONVERGENCE TO THRESHOLD, MAIN LOOP:
C
200     NEW=LOWER+(UPPER-LOWER)*SQRT((THRES-SSRLOW)/(SSRUP-SSRLOW))
        ITERNBR=ITERNBR+1
        DO 150,J=1,PN
          IF(COOR(J).GT.0.0)THEN
            TETA(J)=POPT(J)+NEW*COOR(J)*DPSUP(J)
          ELSE
            TETA(J)=POPT(J)+NEW*COOR(J)*DPINF(J)
          END IF
150     CONTINUE
        SSRNEW=SSR(TETA,PN)
        IF(SSRNEW.GT.THRES)THEN
          UPPER=NEW
        ELSE
          LOWER=NEW
        END IF
        IF(ITERNBR.LT.MAXIT)THEN
          IF(ABS(THRES-SSRNEW)/THRES.LT.EPS)THEN
            WRITE(34,1000)(TETA(K),K=1,PN)
          ELSE
            GO TO 200
          END IF
        ELSE
          WRITE(*,*)"WARNING 3 FROM SUBROUTINE CRUNCH"
          WRITE(*,*)"NUMBER OF ITERATIONS GREATER THAN",MAXIT
          WRITE(*,*)"PARAM VALUE  % DISTANCE FROM OPTIMAL PARAM VALUE"
          DO 135,K=1,PN
            WRITE(*,1000)TETA(K),100*COOR(K)*SQRT(DBLE(PN))
135       CONTINUE
          WRITE(*,*)"-----------------------------------------"
          I=I-1
        END IF
      END IF
      IF(I.LE.MAXPT)GO TO 1
      CLOSE(34)
1000  FORMAT('',8(F12.6))
C
C Note added in Dec. 2006: I'm not sure that this output format is
C correct, I can't read well it on the hard copy of the listing.
C
      END
C------------------------------ SPHERE ----------------------------------------
C
C PURPOSE
C -------
C THE SUBROUTINE SPHERE RETURNS RANDOM COORDINATES FOR A POINT
C BELONGING TO A HYPERSPHERE OF RADIUS 1.0 WITH ITS CENTER
C AT THE ORIGIN (0,0,...,0)
C
C VARIABLES
C ---------
C N     SPACE DIMENSION: NUMBER OF COORDINATES TO GENERATE.
C       N IS IN THE RANGE [2,8]
C
C COOR  COORDINATES OF POINTS RETURNED BY SPHETEST (Note added in
C       Dec. 2006: this should be a typo for SPHERE)
C
C SEED  SEED FOR THE PSEUDO-RANDOM GENERATOR,
C       ITS VALUE HAS TO BE INITIALIZED BY THE CALLING PROGRAM
C       WITHIN [0-1] AND NOT FURTHER MODIFIED
C
C COMMENTS
C --------
C THE N COORDINATES OF A POINT ARE CHOSEN AT RANDOM WITHIN A
C HYPERCUBE. THE DISTANCE FROM THIS POINT TO THE ORIGIN (0,0,...,0) IS
C COMPUTED. IF THE DISTANCE IS LESS THAN 1.0, THEN THE POINT IS INSIDE THE
C HYPERSHERE, AND THE COORDINATES ARE NORMALIZED SO THAT THE
C DISTANCE TO THE ORIGIN EQUALS ONE. WHEN THE POINT IS OUTSIDE THE
C HYPERSPHERE, ANOTHER POINT IS CHOSEN AT RANDOM.
C THIS IS TO OBTAIN A HOMOGENEOUS DISTRIBUTION OF THE
C POINTS ON THE SURFACE OF THE HYPERSPHERE.
C
C THE OBSERVED MEAN NUMBER OF TRIALS TO OBTAIN A POINT WITHIN A
C HYPERSHERE IS GIVEN HERE AS FUNCTION OF THE SPACE DIMENSION N:
C
C N  NUMBER OF TRIALS (MEAN OVER 1000)
C 2  1.284
C 3  1.970
C 4  3.243
C 5  6.010
C 6 12.193
C 7 26.699
C 8 62.023
C
C-------------------------------------------------------------------------------
      SUBROUTINE SPHERE(N,COOR,SEED)
      INTEGER N
      REAL*8 COOR(8),SEED,SUM

5     SUM=0.0
C
C CHOICE OF N COORDINATES AT RANDOM
C
      DO 10,I=1,N
        COOR(I)=2.0*RANDOM(SEED)-1.0
        SUM=SUM+COOR(I)*COOR(I)
10    CONTINUE
C
C TESTING WHETHER THE POINT IS INSIDE THE HYPERSHERE:
C
      IF(ABS(SQRT(SUM)).LE.1.0)THEN
C
C Note added in Dec. 2006: I don't understand why the absolute value is
C taken in the line right above. Extreme defensive programing? A sum
C of squares is always positive, this shouldn't be useful. Anyway, it
C doesn't hurt, I keep it as in the original code.
C
        DO 20,I=1,N
          COOR(I)=COOR(I)/SQRT(SUM)
20      CONTINUE
      ELSE
        GO TO 5
      END IF
      END
C------------------------------ RANDOM ----------------------------------------
C
C PURPOSE
C -------
C THE FUNCTION RANDOM RETURNS A PSEUDO-RANDOM NUMBER WITHIN [0-1]
C
C VARIABLES
C ---------
C
C X  SEED FOR THE PSEUDO-RANDOM GENERATOR,
C    ITS VALUE HAS TO BE INITIALIZED BY THE CALLING PROGRAM
C    WITHIN [0-1] AND NOT FURTHER MODIFIED
C
C NOTE ADDED IN DEC 2006
C ----------------------
C
C This generator is known to behave satisfactorily for up to 100,000
C successive calls:
C
C @article{ChasseJL1974,
C   Author = {Chass{\'e}, J.-L. and Debouzie, D.},
C   Title = {Utilisation des tests de {K}iveliovitch et {V}ialar dans
C     l'{\'e}tude de quelques g{\'e}n{\'e}rateurs de nombres pseudo-
C     al{\'e}atoires},
C   Journal = {Revue de Statistique Appliqu{\'e}e},
C   Volume = {22},
C   Pages = {83-90},
C Year = {1974} }
C
C Enter ?".Random.seed" at your R prompt for references to better
C pseudo-random generators.
C
C------------------------------------------------------------------------------
      FUNCTION RANDOM(X)
      REAL*8 X,RANDOM
      RANDOM=MOD(262144*3125*X,262144.0)/262144.0
      X=RANDOM
      END

     