home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 474 / ems.arc / PC-PLAN.FOR < prev    next >
Text File  |  1986-02-21  |  12KB  |  418 lines
  1. C
  2. C        IPMAX -- LARGEST NUMBER OF THREATMENTS THAT CAN BE PERMUTED
  3. C                 IPERM IS (IPMAX,FACTORIAL(IPMAX))
  4. C        LMAX  -- LARGEST NUMBER OF TREATMENTS THAT CAN BE ASSIGNED RANDOMLY
  5. C        ISMAX -- LARGEST NUMBER OF SUBJECTS PER BLOCK
  6. C        NPERPG -- NUMBER OF SUBJECTS PER PRINTED PAGE FOR TREATMENTS
  7. C                     ASSIGNED AT RANDOM
  8. C
  9. C        TREATMENT LABELS CAN CONTAIN UP TO 15 CHARACTERS. SEE ALSO
  10. C              SUBROUTINE RPERMC
  11. C        BLOCK LABELS CAN CONTAIN UP TO 40 CHARACTERS.
  12. C
  13.       DIMENSION KTR(20), IP(400), IPERM(5,120), IPERM0(5), IPERM1(5)
  14.       CHARACTER*15 TLAB(20), SLAB(400), SLABP(20)
  15.       CHARACTER*40 FNAME, BLAB(20)
  16.       CHARACTER*1 QUERY
  17.       LOGICAL QFILE
  18.       DOUBLE PRECISION DCSEED, DTSEED
  19.       DATA IIN /0/, IOUT/0/, IWOUT/10/
  20.       DATA LMAX /20/, ISMAX /400/, IPMAX /5/, NPERPG /80/
  21. C
  22.       QFILE = .FALSE.
  23.       DO 70 I = 1, IPMAX
  24.       IPERM0(I) = I
  25.       IPERM1(I) = I
  26.    70 CONTINUE
  27.       DO 80 I = 1, LMAX
  28.       BLAB(I) = '                    '
  29.       TLAB(I) = '                    '
  30.       SLABP(I) = '                    '
  31.    80 CONTINUE
  32.       DO 90 I = 1, ISMAX
  33.    90 SLAB(I) = '                    '
  34. C
  35.       WRITE(IOUT,91)
  36.    91 FORMAT(/' Do you wish to save the results of this session?  ',
  37.      *  '(Y or N):  '$)
  38.       READ (IIN,200) QUERY
  39.       IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QFILE = .TRUE.
  40.       IF (.NOT. QFILE) GOTO 93
  41.       WRITE (IOUT,92)
  42.    92 FORMAT (' Enter filename:  '$)
  43.       READ (IIN,'(A)') FNAME
  44.       OPEN (IWOUT, FILE = FNAME, STATUS = 'NEW')
  45. C
  46. C        SET SEEDS
  47. C
  48.    93 CALL SEEDER(IX,IY,IZ,IIN,IOUT)
  49.       WRITE(IOUT,94) IX,IY,IZ
  50.       IF (QFILE) WRITE(IWOUT,94) IX,IY,IZ
  51.    94 FORMAT(/' The initial seeds for the random number generator',
  52.      *   ' are:  '/10X,3I15)
  53. C
  54.  
  55.   100 WRITE (IOUT,110) LMAX
  56.   110 FORMAT(/' Enter the number of blocks [.LE.',I2,']:  '$)
  57.       READ (IIN,*,ERR=100) NBLK
  58.       IF (NBLK.LT.1 .OR. NBLK.GT.LMAX) GOTO 100
  59.       IF (NBLK.EQ.1) GOTO 150
  60. C
  61.       WRITE (IOUT,130)
  62.   130 FORMAT (' Enter block labels, one per line:')
  63.       READ (IIN,'(A)') (BLAB(I),I=1,NBLK)
  64. C
  65.   150 WRITE (IOUT,160)
  66.   160 FORMAT (/' Enter number of subjects per block:  '$)
  67.       READ (IIN,*,ERR=150) NSUBJ
  68.       IF (NSUBJ.GT.ISMAX) STOP 'Too many subjects...recompile'
  69. C
  70.   170 WRITE (IOUT,180)
  71.   180 FORMAT (/' Do you wish to assign treatments to subjects (A),'/
  72.      *   '         or permute treatments within subject (P)?:  '$)
  73.       READ (IIN,200) QUERY
  74.   200 FORMAT(A1)
  75.       IF (QUERY.NE.'A' .AND. QUERY.NE.'a' .AND. 
  76.      *    QUERY.NE.'P' .AND. QUERY.NE.'p') GOTO 170
  77. C
  78.       IF (QUERY.EQ.'P' .OR. QUERY.EQ.'p') GOTO 1300
  79. C
  80. C        ASSIGN TREATMENTS TO SUBJECTS
  81. C
  82.   300 WRITE (IOUT,310) LMAX
  83.   310 FORMAT(/' Enter the number of treatments [.LE.',I2,']:  '$)
  84.       READ (IIN,*,ERR=300) NTR
  85.       IF (NTR.LT.2 .OR. NTR.GT.LMAX) GOTO 300
  86.       FNTR = NTR
  87. C
  88.       WRITE (IOUT,330)
  89.   330 FORMAT (' Enter treatment labels, one per line:')
  90.       READ (IIN,'(A)') (TLAB(I),I=1,NTR)
  91. C
  92.   340 WRITE (IOUT,350)
  93.   350 FORMAT (/' Are treatments:'//
  94.      * '      1 -- assigned at random without any forced balancing,'/
  95.      * '      2 -- assigned at random with an equal number per'
  96.      * ' treatment, or'/
  97.      * '      3 -- assigned at random with a prespecified number per',
  98.      * ' treatment?:  '$)
  99. C
  100.       READ (IIN,*,ERR=340) NPICK
  101.       IF (NPICK.LT.1 .OR. NPICK.GT.3) GOTO 340
  102.       IF (NPICK - 2) 400, 351, 360
  103. C
  104.   351 IF (NSUBJ.EQ.NSUBJ/NTR*NTR) GOTO 400
  105.       WRITE (IOUT,352)
  106.   352 FORMAT(/'** Numbers of subjects and treatments do not mesh'/)
  107.       GOTO 340
  108. C
  109.   360 WRITE (IOUT,370)
  110.   370 FORMAT (' Enter the number of subjects for each treatment:  '$)
  111.       READ (IIN,*,ERR=360) (KTR(I),I=1,NTR)
  112.       KTRSUM = 0
  113.       DO 380 I = 1, NTR
  114.   380 KTRSUM = KTRSUM + KTR(I)
  115.       IF (KTRSUM.EQ.NSUBJ) GOTO 400
  116.       WRITE (IOUT,390)
  117.   390 FORMAT(/'** Total number of subjects and number of subjects'/
  118.      *        '     per treatment do not mesh.'/)
  119.       GOTO 360
  120. C
  121.   400 DO 800 I = 1, NBLK
  122. C
  123.       GOTO (410,430,450), NPICK
  124. C
  125. C        TOTALLY AT RANDOM
  126. C
  127.   410 DO 420 K = 1, NSUBJ
  128.       RAND = RANDOM(IX,IY,IZ)
  129.       INDEX = 1. + FNTR * RAND
  130.       SLAB(K) = TLAB(INDEX)
  131.   420 CONTINUE
  132.       GOTO 500
  133. C
  134. C        EQUAL NUMBER PER GROUP
  135. C
  136.   430 NPERTR = NSUBJ / NTR
  137.       K = 0
  138.       DO 440 I1 = 1, NTR
  139.       DO 440 I2 = 1, NPERTR
  140.       K = K + 1
  141.       SLAB(K) = TLAB(I1)
  142.   440 CONTINUE
  143. C
  144.       CALL RPERMC(SLAB,NSUBJ,IX,IY,IZ)
  145.       GOTO 500
  146. C
  147. C        PRESPECIFIED NUMBER PER GROUP
  148. C
  149.   450 K = 0
  150.       DO 460 I1 = 1, NTR
  151.       NPERTR = KTR(I1)
  152.       DO 460 I2 = 1, NPERTR
  153.       K = K + 1
  154.       SLAB(K) = TLAB(I1)
  155.   460 CONTINUE
  156. C
  157.       CALL RPERMC(SLAB,NSUBJ,IX,IY,IZ)
  158. C
  159.   500 IF (NBLK.GT.1) WRITE(IOUT,510) BLAB(I)
  160.       IF (NBLK.GT.1 .AND. QFILE) WRITE(IWOUT,510) BLAB(I)
  161.   510 FORMAT (////5X,A40)
  162. C
  163.       NPAGE = 1 + (NSUBJ - 1) / NPERPG
  164.       DO 700 IPDEX = 1, NPAGE
  165. C
  166.       ISTART = 1 + (IPDEX - 1) * NPERPG
  167.       IEND = MIN0(NSUBJ,IPDEX*NPERPG)
  168.       NHALF = (ISTART + IEND) / 2
  169.       INCR = NHALF - ISTART + 1
  170. C
  171.       WRITE(IOUT,515)
  172.       IF (QFILE) WRITE(IWOUT,515) 
  173.   515 FORMAT (///)
  174. C
  175.       DO 600 IW = ISTART, NHALF - 1
  176.       WRITE(IOUT,520) IW, SLAB(IW), IW+INCR, SLAB(IW+INCR)
  177.       IF (QFILE)
  178.      *   WRITE(IWOUT,520) IW, SLAB(IW), IW+INCR, SLAB(IW+INCR)
  179.   520 FORMAT(8X,I3,'.',5X,A15,10X,I3,'.',5X,A15)
  180.   600 CONTINUE
  181.       IODD = MOD(IEND - ISTART + 1,2)
  182.       IF (IODD.EQ.0) WRITE(IOUT,520) NHALF, SLAB(NHALF),
  183.      *   NHALF+INCR, SLAB(NHALF+INCR)
  184.       IF (IODD.EQ.1) WRITE(IOUT,520) NHALF, SLAB(NHALF)
  185.       IF (IODD.EQ.0 .AND. QFILE) WRITE(IWOUT,520) NHALF, SLAB(NHALF),
  186.      *   NHALF+INCR, SLAB(NHALF+INCR)
  187.       IF (IODD.EQ.1 .AND. QFILE) WRITE(IWOUT,520) NHALF, SLAB(NHALF)
  188.   700 CONTINUE
  189.   800 CONTINUE
  190.       STOP ' '
  191. C
  192. C        PERMUTE WITHIN SUBJECT
  193. C
  194.  1300 WRITE (IOUT,1310) IPMAX
  195.  1310 FORMAT(/' Enter the number of treatments [.LE.',I1,']:  '$)
  196.       READ (IIN,*,ERR=1300) NTR
  197.       IF (NTR.LT.2 .OR. NTR.GT.IPMAX) GOTO 1300
  198.       FNTR = NTR
  199.       NPERM = 1
  200.       DO 1320 I = 2, NTR
  201.  1320 NPERM = NPERM * I
  202.       FNPERM = NPERM
  203. C
  204.       DO 1325 J = 1, NPERM
  205.       CALL PERMUT(IPERM1,IPERM0,NTR)
  206.       DO 1325 K = 1, NTR
  207.       IPERM(K,J)= IPERM1(K)
  208.  1325 CONTINUE
  209. C
  210.       WRITE (IOUT,1330)
  211.  1330 FORMAT (' Enter treatment labels, one per line:')
  212.       READ (IIN,'(A)') (TLAB(I),I=1,NTR)
  213. C
  214.       NPICK = 1
  215.       IF (NSUBJ.NE.NSUBJ/NPERM*NPERM) GOTO 1360
  216.  1340 WRITE (IOUT,1350)
  217.  1350 FORMAT (' Are permutations:'//
  218.      * '      1 -- assigned at random without any forced balancing'
  219.      * ' or'/
  220.      * '      2 -- assigned at random with an equal number of',
  221.      * ' each?:  '$)
  222. C
  223.       READ (IIN,*,ERR=1340) NPICK
  224.       IF (NPICK.LT.1 .OR. NPICK.GT.2) GOTO 1340
  225. C
  226.  1360 DO 1800 I = 1, NBLK
  227. C
  228.       GOTO (1410,1430), NPICK
  229. C
  230. C        TOTALLY AT RANDOM
  231. C
  232.  1410 DO 1420 K = 1, NSUBJ
  233.       RAND = RANDOM(IX,IY,IZ)
  234.  1420 IP(K) = 1. + FNPERM * RAND
  235. C
  236.       GOTO 1500
  237. C
  238. C        EQUAL NUMBER PER GROUP
  239. C
  240.  1430 NPERTR = NSUBJ / NPERM
  241.       K = 0
  242.       DO 1440 I1 = 1, NPERM
  243.       DO 1440 I2 = 1, NPERTR
  244.       K = K + 1
  245.       IP(K) = I1
  246.  1440 CONTINUE
  247. C
  248.       CALL RPERMI(IP,NSUBJ,IX,IY,IZ)
  249. C
  250. C        OUTPUT 
  251. C
  252.  1500 IF (NBLK.GT.1) WRITE(IOUT,510) BLAB(I)
  253.       IF (NBLK.GT.1 .AND. QFILE) WRITE(IWOUT,510) BLAB(I)
  254. C
  255.       DO 1600 I1 = 1, NSUBJ
  256.       DO 1510 K = 1, NTR
  257.  1510 SLABP(K) = TLAB(IPERM(K,IP(I1)))
  258. C
  259.       WRITE(IOUT,1520) I1, (SLABP(K),K=1,NTR)
  260.       IF (QFILE) WRITE(IWOUT,1520) I1, (SLABP(K),K=1,NTR)
  261.  1520 FORMAT (8X,I3,'.',5X,A15,5X,A15,5X,A15,:/
  262.      *   (25X,A15,5X,A15,5X,A15))
  263. C
  264.  1600 CONTINUE
  265.  1800 CONTINUE
  266.  
  267.       STOP ' '
  268.       END
  269. C
  270.       SUBROUTINE RPERMI(LIST,N,IX,IY,IZ)
  271. C
  272. C        A RANDOM PERMUTATION OF THE ITEMS IN LIST
  273. C
  274.       DIMENSION LIST(N)
  275. C
  276.       NM1 = N - 1
  277.       DO 10 I = 1, NM1
  278.       II = N - I + 1
  279.       J = FLOAT(II) * RANDOM(IX,IY,IZ) + 1.
  280.       IDUMMY = LIST(II)
  281.       LIST(II) = LIST(J)
  282.       LIST(J) = IDUMMY
  283.    10 CONTINUE
  284.       RETURN
  285.       END
  286. C
  287.       SUBROUTINE RPERMC(LIST,N,IX,IY,IZ)
  288. C
  289. C        A RANDOM PERMUTATION OF THE ITEMS IN LIST
  290. C
  291.       CHARACTER*15 LIST(N),IDUMMY
  292. C
  293.       NM1 = N - 1
  294.       DO 10 I = 1, NM1
  295.       II = N - I + 1
  296.       J = FLOAT(II) * RANDOM(IX,IY,IZ) + 1.
  297.       IDUMMY = LIST(II)
  298.       LIST(II) = LIST(J)
  299.       LIST(J) = IDUMMY
  300.    10 CONTINUE
  301.       RETURN
  302.       END
  303. C
  304.       SUBROUTINE PERMUT(IX,IP,N)
  305. C
  306. C        CACM ALGORITHM 86
  307. C        PERMUTE
  308. C        BY PECK, J.E.L. AND G.F. SCHRACK
  309. C
  310. C        TRANSLATED INTO FORTRAN BY
  311. C                   Gerard E. Dallal
  312. C                   USDA-Human Nutrition Research Center on Aging
  313. C                              at Tufts University
  314. C                   711 Washington Street
  315. C                   Boston, MA  02111
  316. C
  317.       DIMENSION IX(N), IP(N)
  318.       DO 100 K = 2, N
  319.       IT = IX(1)
  320.       KM = K - 1
  321.       DO 50 I = 1, KM
  322.    50 IX(I) = IX(I+1)
  323.       IX(K) = IT
  324.       IP(K) = IP(K) - 1
  325.       IF(IP(K).NE.0) RETURN
  326.       IP(K) = K
  327.   100 CONTINUE
  328.       RETURN
  329.       END
  330. C
  331.       FUNCTION RANDOM(IX,IY,IZ)
  332. C
  333. C        ALGORITHM AS 183  APPL. STATIST. (1982) VOL.31, NO.2
  334. C
  335. C        RETURNS A PSEUDO-RANDOM NUMBER RECTANGULARLY DISTRIBUTED
  336. C        BETWEEN 0 AND 1.
  337. C
  338. C        IX, IY AND IZ SHOULD BE SET TO INTEGER VALUES BETWEEN
  339. C        1 AND 30000 BEFORE FIRST ENTRY
  340. C
  341.       IX = 171 * MOD(IX, 177) - 2  * (IX / 177)
  342.       IY = 172 * MOD(IY, 176) - 35 * (IY / 176)
  343.       IZ = 170 * MOD(IZ, 178) - 63 * (IZ / 178)
  344. C
  345.       IF (IX .LT. 0)  IX = IX + 30269
  346.       IF (IY .LT. 0)  IY = IY + 30307
  347.       IF (IZ .LT. 0)  IZ = IZ + 30323
  348. C
  349. C        IF INTEGER ARITHMETIC UP TO 5212632 IS AVAILABLE,
  350. C        THE PRECEDING 6 STATEMENTS MAY BE REPLACED BY
  351. C
  352. C        IX = MOD(171 * IX, 30269)
  353. C        IY = MOD(172 * IY, 30307)
  354. C        IZ = MOD(170 * IZ, 30323)
  355. C
  356. C        ON SOME MACHINES, THIS MAY SLIGHTLY INCREASE
  357. C        THE SPEED.  THE RESULTS WILL BE IDENTICAL.
  358. C
  359.       RANDOM = AMOD(FLOAT(IX) / 30269.0 + FLOAT(IY) / 30307.0 +
  360.      *              FLOAT(IZ) / 30323.0, 1.0)
  361.       RETURN
  362.       END
  363.  
  364. C  DATTIM.FOR program - To access the date and time:
  365.  
  366.       INTERFACE TO SUBROUTINE TIME (N,STR)
  367.       CHARACTER*10 STR [NEAR,REFERENCE]
  368.       INTEGER*2 N [VALUE]
  369.       END
  370. C
  371. C
  372.       SUBROUTINE SEEDER(IX,IY,IZ,IIN,IOUT)
  373. C
  374. C        SEEDS MUST BE BETWEEN 1 AND 30000
  375. C
  376.       CHARACTER*1 QUERY
  377.       CHARACTER*10 TSTR
  378.       DIMENSION ITIME(8)
  379. C
  380.       WRITE(IOUT,50)
  381.    50 FORMAT(/' Do you wish to enter the random number generator''s',
  382.      *   ' seeds yourself? (Y or N):  '$)
  383.       READ (IIN,60) QUERY
  384.    60 FORMAT(A1)
  385.       IF (QUERY.NE.'Y' .AND. QUERY.NE.'y') GOTO 200
  386. C
  387. C        ENTER SEEDS
  388. C
  389.    70 WRITE(IOUT,80)
  390.    80 FORMAT(/' Enter three seeds between 1 and 30000 separated',
  391.      *   ' by commas:')
  392.       READ(IIN,*,ERR=70) IX,IY,IZ
  393.       IF (IX.LE.1 .OR. IX.GE.30000 .OR. IY.LE.1 .OR. IY.GE.30000 
  394.      *   .OR. IZ.LE.1 .OR. IZ.GE.30000) GOTO 70
  395.       RETURN
  396. C
  397. C         GET SEEDS FROM CLOCK
  398. C
  399.   200 CALL TIME (10,TSTR)
  400.       DO 210 I = 1,8
  401.   210 ITIME(I) = ICHAR(TSTR(I:I)) - 48
  402. C
  403.       IX = ITIME(8) + 10 * (ITIME(7) + 10 * (ITIME(5) + 10 * (ITIME(4) +
  404.      *   10 * (ITIME(2) + 10 * ITIME(1)))))
  405.  
  406.       IY = ITIME(2) + 10 * (ITIME(1) + 10 * (ITIME(8) + 10 * (ITIME(7) +
  407.      *   10 * (ITIME(5) + 10 * ITIME(4)))))
  408.  
  409.       IZ = ITIME(5) + 10 * (ITIME(4) + 10 * (ITIME(2) + 10 * (ITIME(1) +
  410.      *   10 * (ITIME(8) + 10 * ITIME(7)))))
  411.  
  412.       IY = 1 + MOD(IY,29999)
  413.       IZ = 1 + MOD(IZ,29999)
  414.  
  415.       RETURN
  416.       END
  417.  
  418.