#!/bin/sh
# This is a shell archive (produced by shar 3.49)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#
# made 04/21/1995 11:02 UTC by shelp@simpson
# Source directory /sunusers/stats/students/stats3/shelp/S/library/Dopt/tmp
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 8177 -rw-r--r-- Dopt.d
# 28206 -rw-r--r-- Dopt.f
# 5660 -rw-r--r-- Dopt.q
# 134 -rw-r--r-- First.lib.q
# 135 -rw-r--r-- First.lib2.q
# 189 -rw-r--r-- INDEX
# 859 -rwxr-xr-x INS
# 3054 -rw-r--r-- INSTALLATION
# 8603 -rw-r--r-- NOTES
# 70 -rw-r--r-- README
# 1138 -rw-r--r-- STATEMENT
# 385 -rwxr-xr-x TEST
# 1705 -rw-r--r-- master.out
# 1447 -rw-r--r-- test.q
#
# ============= Dopt.d ==============
if test -f 'Dopt.d' -a X"$1" != X"-c"; then
echo 'x - skipping Dopt.d (File already exists)'
else
echo 'x - extracting Dopt.d (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'Dopt.d' &&
.BG
.FN Dopt
.TL
Finding D-optimal experimental designs.
.DN
A function to generate approximately D-optimal experimental
designs using a Fortran routine of Alan Miller and Nam-Ky Nguyen
(AS295). The algorithm is a combination of the exchange
algorithm of Fedorov and a supplementary algorithm of Cook &
Nachtsheim (1989) for swapping treatments between blocks.
.CS
Dopt(treatments, blocks,
X data.treatments = sys.parent(),
X data.blocks = sys.parent(),
X initial.points = rep(0, npoints),
X fixed.points = rep(F, npoints),
X RNG = as.integer(runif(3) * 30269),
X eps = .Machine$single.eps,
X runs)
.RA
.AG treatments
A number, factor, matrix or formula specifying how the full
design matrix of candidate points is to be constructed. (If this
argument is a single number it is interpreted as a number of
simple treatments to be applied in a single factor experiment.)
.AG blocks
A number, factor, vector or formula specifying the number of runs
in the experiment and the block structure, if any. (If this
argument is a single number it specifies the number of runs in a
single block experiment.)
.OA
.AG data.treatments
Optional data frame to supply variables for the treatments
specification.
.AG data.blocks
Optional data frame to supply variables for the blocks
specification.
.AG initial.points
A vector of numbers specifying, in order, rows of the treatment
design matrix to be used in the initial trial design. A zero or
missing value denotes an unspecified run. Any unspecified runs
of the initial design are chosen at random.
.AG fixed.points
Logical vector; T if the corresponding initial point is to be
kept in the experiment, F if it may be replaced by the algorithm.
.AG RNG
Vector of three positive integers less than 30269 used to
initialize the Wichmann Hill random number generator incorporated
into the routine.
.AG eps
Tolerance for zero singular values to detect rank deficiency.
.AG runs
Alternative name for the `blocks' argument.
.RT
A data frame specifying the chosen experimental design. This
will have
X
A factor specifying the block structure, (or variables from the
blocks data frame needed to reconstruct it if a formula for
blocks is given).
X
Variables (or design submatrix) from the treatments specification
specifying the treatments to be used on each run.
X
A numeric vector labelled "Point.No" giving the row number of the
treatment design matrix used for each run of the experiment.
X
If the experiment included fixed runs, a logical vector labelled "Fixed"
specifying which runs were fixed in advance at their initial points.
X
The data frame has three additional attributes, "log(det)" giving
the log of the final determinant as returned by the algorithm,
which may be useful in comparing two runs for the same design,
"RNG", the three starting integers used for the local random
number generators for this run of the algorithm and "CRE" a
vector of `canonical relative efficiencies'.
.SE
None.
.DT
The program of Miller and Nguyen takes a treatments design
matrix, T, of "candidate points" and forms another matrix, X0,
from it by selecting rows, possibly with repitition. If the
experiment has k blocks, and X is the matrix got from X0 by
subtracting block means from each column, then X0 is chosen to
make det(X'X) is as large as possible. If this is maximal the
design is D-optimal. The program uses an exchange algorithm due
to Fedorov, with modifications for computational efficiency, and
exact optimality is not guaranteed so several runs of the
procedure may be useful.
X
Only simple block structures are allowed. Zero blocks implies no
intercept term and the matrices X0 and X above are the same. One
block implies an intercept term, so the columns of X have mean
zero.
X
The present routine is an S interface to the Miller Nguyen
program.
X
The `blocks' argument has two purposes. Through its length (or
explicitly, if it is a single number), it prescribes the number
of runs in the experiment. If it is a vector, it is coerced to
factor and used to define a block structure for the experiment.
If it is a formula it is not interpreted in the ordinary way: any
vectors present are coerced to factors and all terms are then
used to define a "subclasses" block structure. Thus ~a + b, ~a*b
and ~a/b all define the same block structure in this case. Note
that non-simple block structures such as for latin square and
row-column designs, or designs for known plot covariates are
not accommodated. If there are no blocks (see NOTE below) the
blocks argument is only used to define the number of runs.
X
The treatments argument, together with its data.treatments data
frame, if any, is used to construct the design matrix T of
candidate points. It either an S formula or effectively coerced
to such, and is interpreted in the usual way. The rows of the
design matrix, T, or `points' are numbered 1 to ncand.
X
NOTE: In Dopt the only way to specify that the design is not to
have an intercept term is to give the treatments specification as
a formula and explicitly to omit the intercept term, as in ~x +
x^2 - 1 for a quadratic regression through the origin.
X
The optional argument `initial.points' is a vector usually of
length equal to the number of runs, with integer entries
specifying rows of T to be used at that position in the initial
trial matrix X0. Any zero or missing values entries specify rows
that are to be initially chosen at random.
X
The `fixed.points' argument is a logical vector specifying whether
or not the corresponding row of X is to be held fixed at its
initially specified value. Thus the routine may be used to
augment an existing design with additional points.
X
The CRE attribute of the result gives the non-zero eigenvalues of
XX'X relative to (T'T * npoints/ncand). Under some circumstances
for an ideal design choice they should all be unity. In most
cases they will be above and below unity and their harmonic mean
gives one global measure of how well the design succeeds in
capturing potential information, with connexions to A-efficiency.
X
.SH REFERENCES
Cook, R. D. and Nachtsheim, C. J. (1989) Computer-aided blocking
of factorial and response-surface designs. Technometrics, vol.
31, 339-346.
X
Miller, A. J. and Nguyen, N-K. (1994). Algorithm AS 295: A
Fedorov exchange algorithm for D-optimal design, Applied
Statistics, vol. 43, 669-678. (Held in statlib as apstat/295).
X
Nguyen, N-K. and Miller, A. J. (1992). A review of some exchange
algorithms for constructing discrete D-optimal designs.
Comput. Statist. & Data Anal., vol.14, 489-498.
.SA
S+DOX
.EX
1. Approximate optimal design of a cubic regression with 12 points.
X
> x <- 0:100
> dsn <- Dopt(~x + x^2 + x^3, runs = 12)
> sort(dsn$x)
X [1] 0 0 0 27 28 28 72 72 73 100 100 100
X
2. A block design for 7 treatments in 7 blocks of size 3:
X
> tr <- factor(1:7)
> bl <- factor(rep(1:7, rep(3,7)))
> dsn <- Dopt(~tr, ~bl)
>
> # look at the concurrency matrix:
>
> crossprod(table(dsn$bl, dsn$tr))
X 1 2 3 4 5 6 7
1 3 1 1 1 1 1 1
2 1 3 1 1 1 1 1
3 1 1 3 1 1 1 1
4 1 1 1 3 1 1 1
5 1 1 1 1 3 1 1
6 1 1 1 1 1 3 1
7 1 1 1 1 1 1 3
> # in this case the design is a BIB and hence optimal.
X
3. An approximate main effect plan for a 2^11 experiment in 12
X runs. (In this case a Plackett Burman design could be used
X for comparison). We also ensure that the first run has all
X factors at their lowest level.
X
> l <- c("-","+")
> fac <- expand.grid(l,l,l,l,l,l,l,l,l,l,l)
> names(lis) <- LETTERS[-6][1:11] # omit F as factor name
> dsn <- Dopt(~., 12, fac, initial = 1, fixed = T)
> dsn
X A B C D E G H I J K L Point.No Fixed
X 1 - - - - - - - - - - - 1 T
X 2 - + - - + + + - - - + 1139 F
X 3 + - + - + + - - + - - 310 F
X 4 - + - + - + - - + + - 811 F
X 5 + - - + + + + + - + - 762 F
X 6 - - + + + - - - - + + 1565 F
X 7 + + - - + - - + + + + 1940 F
X 8 - + + + + - + + + - - 479 F
X 9 + + + - - - + - - + - 584 F
10 - - + - - + + + + + + 2021 F
11 + - - + - - + - + - + 1354 F
12 + + + + - + - + - - + 1200 F
>
> # in large cases several tries of the algorithm might be useful.
X
.KW dox
.WR
SHAR_EOF
chmod 0644 Dopt.d ||
echo 'restore of Dopt.d failed'
Wc_c="`wc -c < 'Dopt.d'`"
test 8177 -eq "$Wc_c" ||
echo 'Dopt.d: original size 8177, current size' "$Wc_c"
fi
# ============= Dopt.f ==============
if test -f 'Dopt.f' -a X"$1" != X"-c"; then
echo 'x - skipping Dopt.f (File already exists)'
else
echo 'x - extracting Dopt.f (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'Dopt.f' &&
C (C) COPYRIGHT, 1994. Royal Statistical Society, London.
C
C Requests to incorporate any or all of this software into any
C commercial product must be addressed to
C
C The Secretary,
C The Royal Statistical Society,
C 25 Enford St,
C LONDON, UK, W1H2BH
X
C-------------------------------------------------------------------------
X
X SUBROUTINE FEDROV(X, DIM1, NCAND, KIN, N, NBLOCK, IN, BLKSIZ, K,
X + RSTART, NRBAR, D, RBAR, PICKED, LNDET, XX, TOL, ZPZ, WK,
X + IFAULT)
X
C [Formerly named DOPT, but renamed here to avoid loader table
C conflicts. W. Venables.]
C
C
C Heuristic algorithm to pick N rows of X out of NCAND to maximize
C the determinant of X'X, using the Fedorov exchange algorithm.
C
C Code written by Alan Miller and Nam Nguyen from CSIRO, Melbourne,
C Australia and to be published in Applied Statistics.
C
C [Now appeared as AS295, Applied Statistics, vol 43, 669-677,
C (1994) W. Venables.]
C
C Auxiliary routines required: A random number generator (RAND),
C and routines CLEAR & REGCF from AS 274.
C
C [All supplied herewith. W. Venables.]
C
X INTEGER DIM1, NCAND, KIN, N, NBLOCK, IN(*), BLKSIZ(*), K, NRBAR,
X + PICKED(N), IFAULT
X LOGICAL RSTART
X REAL X(DIM1, KIN)
X DOUBLE PRECISION D(K), RBAR(NRBAR), LNDET, XX(K), TOL(K), WK(K),
X + ZPZ(NCAND, *)
C
C Local variables
C
X INTEGER I, J, NIN, POINT, CASE, NB, BLOCK, L, POS, BEST, FIRST,
X + LAST, CAND, LASTIN, LSTOUT, DROP, REMPOS, BL, RPOS, LAST1,
X + LAST2, FIRST1, FIRST2, POS1, POS2, BLOCK1, BLOCK2, CASE1,
X + CASE2, POSI, POSJ, BESTB1, BESTB2, BESTP1, BESTP2, RANK,
X + MXRANK, INC
X DOUBLE PRECISION ONE, ZERO, MINUS1, TEMP, EPS, DMAX, ABOVE1,
X + SUM, SMALL, HUNDRD
X DOUBLE PRECISION DELTA
X REAL RAND
X LOGICAL CHANGE
X DATA ONE/1.D0/, ZERO/0.D0/, ABOVE1/1.0001D0/, EPS/1.D-15/,
X + MINUS1/-1.D0/, SMALL/1.D-04/, HUNDRD/100.D0/
C
C Argument checks
C
X IFAULT = 0
X IF (DIM1 .LT. NCAND) IFAULT = 1
X IF (K .GT. N) IFAULT = IFAULT + 2
X IF (NRBAR .LT. K* (K-1)/2) IFAULT = IFAULT + 4
X IF (K .NE. KIN+NBLOCK) IFAULT = IFAULT + 8
X IF (NBLOCK .GT. 1) THEN
X L = 0
X DO 10 BLOCK = 1, NBLOCK
X 10 L = L + BLKSIZ(BLOCK)
X IF (N .NE. L) IFAULT = IFAULT + 16
X ELSE
X IF (N .NE. BLKSIZ(1)) IFAULT = IFAULT + 16
X END IF
C
C NB = max(1, NBLOCK) so that we can force it to go through DO-loops
C once. NIN = no. of design points forced into the design.
C
X NB = MAX(1, NBLOCK)
X NIN = 0
X DO 20 I = 1, NB
X IF (IN(I) .LT. 0) GO TO 30
X IF (IN(I) .GT. 0) THEN
X IF (IN(I) .GT. BLKSIZ(I)) GO TO 30
X NIN = NIN + IN(I)
X END IF
X 20 CONTINUE
X IF (NIN .LE. N) GO TO 40
X 30 IFAULT = IFAULT + 32
X 40 CONTINUE
X IF (IFAULT .NE. 0) RETURN
X CALL CLEAR(K, NRBAR, D, RBAR, IFAULT)
C
C Set up an array of tolerances.
C
X DO 50 I = 1, K
X 50 TOL(I) = ZERO
X BLOCK = 1
X DO 70 CASE = 1, NCAND
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, CASE)
X DO 60 I = 1, K
X 60 TOL(I) = TOL(I) + ABS(XX(I))
X 70 CONTINUE
X TEMP = FLOAT(N) * EPS / NCAND
X DO 80 I = 1, K
X IF (I .LE. NBLOCK) THEN
X TOL(I) = EPS
X ELSE
X TOL(I) = TOL(I) * TEMP
X END IF
X 80 CONTINUE
C
C Form initial Cholesky factorization.
C
X POS = 1
X DO 120 BLOCK = 1, NB
X IF (RSTART) THEN
X LAST1 = (IN(BLOCK) + BLKSIZ(BLOCK))/2
X INC = SQRT(FLOAT(NCAND) + SMALL)
X END IF
X DO 110 I = 1, BLKSIZ(BLOCK)
X IF (RSTART .AND. I .GT. IN(BLOCK)) THEN
X POINT = 1 + NCAND * RAND()
C
C If I <= LAST1, use a random point, otherwise find the candidate
C which maximizes the rank, and then maximizes the subspace
C determinant for that rank.
C
X IF (I .GT. LAST1) THEN
X MXRANK = 0
X LNDET = -HUNDRD
X DO 100 CAND = 1, NCAND, INC
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, POINT)
X CALL MODTR2(K, NRBAR, XX, D, RBAR, TOL, RANK, SUM)
X IF (RANK .LT. MXRANK) GO TO 90
X IF (RANK .EQ. MXRANK .AND. SUM .LT. LNDET) GO TO 90
X BEST = POINT
X MXRANK = RANK
X LNDET = SUM * ABOVE1
X 90 POINT = POINT + INC
X IF (POINT .GT. NCAND) POINT = POINT - NCAND
X 100 CONTINUE
X POINT = BEST
X END IF
X PICKED(POS) = POINT
X ELSE
C
C Case in which a full design has been input, or points are to be
C forced into the design.
C
X POINT = PICKED(POS)
X END IF
C
C Augment the Cholesky factorization.
C
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, POINT)
X CALL MODTRI(K, NRBAR, ONE, XX, D, RBAR, TOL)
X POS = POS + 1
X 110 CONTINUE
X 120 CONTINUE
C
C Adjust factorization in case of singular matrix.
C
X CALL SINGM(K, NRBAR, D, RBAR, TOL, WK, IFAULT)
C
C If rank of input design < K, try replacing points.
C
X IF (IFAULT .EQ. 0) GO TO 280
C
C Find first row of Cholesky factorization with a zero multiplier.
C
X 180 DO 190 POS = 1, K
X IF (D(POS) .LT. TOL(POS)) GO TO 200
X 190 CONTINUE
X GO TO 280
C
C Find linear relationship between variable in position POS and the
C previous variables.
C
X 200 L = POS - 1
X DO 210 I = 1, POS-1
X WK(I) = RBAR(L)
X L = L + K - I - 1
X 210 CONTINUE
X CALL REGCF(K, NRBAR, D, RBAR, WK, TOL, WK, POS-1, IFAULT)
C
C Find a candidate point which does not satisfy this linear
C relationship. Use a random start.
C
X BL = 1
X CASE = 1 + NCAND * RAND()
X DO 230 CAND = 1, NCAND
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BL, XX, CASE)
X SUM = XX(POS)
X DO 220 I = 1, POS-1
X 220 SUM = SUM - WK(I) * XX(I)
X IF (ABS(SUM) .GT. HUNDRD * TOL(POS)) GO TO 240
X CASE = CASE + 1
X IF (CASE .GT. NCAND) CASE = 1
X 230 CONTINUE
C
C Failed to find any candidate point which would make the design
C of higher rank.
C
X IFAULT = -1
X RETURN
C
C Before adding the point, find one which it can replace without
C lowering the rank.
C
X 240 BL = 0
X TEMP = ONE - SMALL
X POS = IN(1) + 1
X DO 270 BLOCK = 1, NB
X DO 260 J = IN(BLOCK)+1, BLKSIZ(BLOCK)
X L = PICKED(POS)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, L)
X CALL BKSUB2(RBAR, NRBAR, K, XX, WK)
X SUM = ZERO
X DO 250 I = 1, K
X 250 IF (D(I) .GT. TOL(I)) SUM = SUM + WK(I)**2 / D(I)
X IF (SUM .LT. TEMP) THEN
X TEMP = SUM
X REMPOS = POS
X BL = BLOCK
X END IF
X POS = POS + 1
X 260 CONTINUE
X IF (BLOCK .LT. NBLOCK) POS = POS + IN(BLOCK+1)
X 270 CONTINUE
C
C If BL = 0 it means that any point removed from the existing design
C would reduce the rank.
C
X IF (BL .EQ. 0) THEN
X IFAULT = -1
X RETURN
X END IF
C
C Add candidate CASE in block BL, then delete the design point
C already in that position.
C
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BL, XX, CASE)
X CALL MODTRI(K, NRBAR, ONE, XX, D, RBAR, TOL)
X L = PICKED(REMPOS)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BL, XX, L)
X CALL MODTRI(K, NRBAR, MINUS1, XX, D, RBAR, TOL)
X PICKED(REMPOS) = CASE
X
X GO TO 180
C----------------------------------------------------------------------
C
C Design is now of full rank.
C
C Calculate z'z for all candidate points.
C z is the solution of R'z = x, so that z'z = x'.inv(X'X).x
C WK holds sqrt(D) times vector z on return from BKSUB2.
C
X 280 DO 310 BLOCK = 1, NB
X DO 300 CASE = 1, NCAND
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, CASE)
X CALL BKSUB2(RBAR, NRBAR, K, XX, WK)
X TEMP = ZERO
X DO 290 I = 1, K
X 290 TEMP = TEMP + WK(I)**2 / D(I)
X ZPZ(CASE, BLOCK) = TEMP
X 300 CONTINUE
X 310 CONTINUE
X
C----------------------------------------------------------------------
C
C Start of Fedorov exchange algorithm
C
X LASTIN = 0
X LSTOUT = 0
X 320 CHANGE = .FALSE.
X LAST = 0
X DO 420 BLOCK = 1, NB
X FIRST = LAST + 1 + IN(BLOCK)
X LAST = LAST + BLKSIZ(BLOCK)
X DMAX = SMALL
X BEST = 0
C
C Start at a random position within the block.
C I = no. of point being considered for deletion.
C
X POS = FIRST + (BLKSIZ(BLOCK) - IN(BLOCK)) * RAND()
X DO 350 CASE = IN(BLOCK)+1, BLKSIZ(BLOCK)
X POS = POS + 1
X IF (POS .GT. LAST) POS = FIRST
X I = PICKED(POS)
X IF (I .EQ. LASTIN) GO TO 350
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, I)
X CALL BKSUB2(RBAR, NRBAR, K, XX, WK)
X CALL BKSUB1(RBAR, NRBAR, K, WK, WK, TOL, D)
C
C Cycle through the candidates for exchange, using a random start.
C J = no. of point being considered for addition.
C
X J = 1 + NCAND * RAND()
X DO 340 CAND = 1, NCAND
X J = J + 1
X IF (J .GT. NCAND) J = 1
X IF (J .EQ. I .OR. J .EQ. LSTOUT) GO TO 340
C... The Cauchy-Schwarz test.
X TEMP = ZPZ(J,BLOCK) - ZPZ(I,BLOCK)
X IF (TEMP .LT. DMAX) GO TO 340
C
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, J)
X SUM = ZERO
X DO 330 L = 1, K
X 330 SUM = SUM + XX(L) * WK(L)
X TEMP = TEMP + SUM**2 - ZPZ(I,BLOCK)*ZPZ(J,BLOCK)
X IF (TEMP .GT. DMAX) THEN
X DMAX = TEMP * ABOVE1
X BEST = J
X REMPOS = POS
X DROP = I
X END IF
X 340 CONTINUE
X 350 CONTINUE
C
C Exchange points BEST & DROP in position REMPOS, if the determinant
C is increased.
C
X IF (BEST .NE. 0) THEN
X CHANGE = .TRUE.
X IF (NB .EQ. 1) THEN
X LASTIN = BEST
X LSTOUT = DROP
X END IF
C
C Add the new point, BEST, first to avoid ill-conditioning.
C Update z'z.
C
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, BEST)
X CALL BKSUB2(RBAR, NRBAR, K, XX, WK)
X CALL BKSUB1(RBAR, NRBAR, K, WK, WK, TOL, D)
X CALL MODTRI(K, NRBAR, ONE, XX, D, RBAR, TOL)
X TEMP = ONE + ZPZ(BEST,BLOCK)
X DO 360 BL = 1, NB
X DO 380 CASE = 1, NCAND
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BL, XX, CASE)
X SUM = ZERO
X DO 370 L = 1, K
X 370 SUM = SUM + XX(L)*WK(L)
X ZPZ(CASE,BL) = ZPZ(CASE,BL) - SUM**2 / TEMP
X 380 CONTINUE
X 360 CONTINUE
C
C Remove the point DROP, and update z'z.
C
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, DROP)
X CALL BKSUB2(RBAR, NRBAR, K, XX, WK)
X CALL BKSUB1(RBAR, NRBAR, K, WK, WK, TOL, D)
X CALL MODTRI(K, NRBAR, MINUS1, XX, D, RBAR, TOL)
X TEMP = ONE - ZPZ(DROP,BLOCK)
X DO 390 BL = 1, NB
X DO 410 CASE = 1, NCAND
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BL, XX, CASE)
X SUM = ZERO
X DO 400 L = 1, K
X 400 SUM = SUM + XX(L)*WK(L)
X ZPZ(CASE,BL) = ZPZ(CASE,BL) + SUM**2 / TEMP
X 410 CONTINUE
X 390 CONTINUE
C
X PICKED(REMPOS) = BEST
X END IF
X 420 CONTINUE
C
C Repeat until there is no further improvement
C
X IF (CHANGE) GO TO 320
C---------------------------------------------------------------------
C
C If there is more than one block, try swapping treatments between
C blocks. This is the Cook & Nachtsheim (1989) algorithm.
C
X IF (NBLOCK .LE. 1) GO TO 500
C
C
C... RPOS is the position of the first element in RBAR after the rows
C for the block constants.
X RPOS = NBLOCK * K - NBLOCK * (NBLOCK + 1)/2 + 1
C
X 430 LAST1 = 0
C... POS1 & POS2 will hold the positions of the start of the means of
C the X-variables in the two blocks being considered in RBAR.
X POS1 = NBLOCK
X DMAX = SMALL
X CHANGE = .FALSE.
X DO 490 BLOCK1 = 1, NBLOCK-1
X FIRST1 = LAST1 + 1 + IN(BLOCK1)
X LAST1 = LAST1 + BLKSIZ(BLOCK1)
X LAST2 = LAST1
X POS2 = POS1 + K - 1 - BLOCK1
X DO 480 BLOCK2 = BLOCK1+1, NBLOCK
X FIRST2 = LAST2 + 1 + IN(BLOCK2)
X LAST2 = LAST2 + BLKSIZ(BLOCK2)
X DO 470 CASE1 = IN(BLOCK1)+1, BLKSIZ(BLOCK1)
X POSI = FIRST1 - 1 + CASE1
X I = PICKED(POSI)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, I)
X DO 440 L = 1, KIN
X 440 ZPZ(L,1) = XX(L+NBLOCK)
X DO 460 CASE2 = IN(BLOCK2)+1, BLKSIZ(BLOCK2)
X POSJ = FIRST2 - 1 + CASE2
X J = PICKED(POSJ)
X IF (I .EQ. J) GO TO 460
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, J)
X DO 450 L = 1, KIN
X 450 ZPZ(L,2) = XX(L+NBLOCK)
C... Pass the orthogonal factorization to DELTA with the top NBLOCK
C rows removed, i.e. without that part relating to the blocks.
X TEMP = DELTA(KIN, ZPZ(1,1), ZPZ(1,2), RBAR(POS1),
X + RBAR(POS2), BLKSIZ(BLOCK1), BLKSIZ(BLOCK2), NRBAR,
X + D(NBLOCK+1), RBAR(RPOS), ZPZ(1,3), WK, XX, ZPZ(1,2))
X IF (TEMP .GT. DMAX) THEN
X DMAX = TEMP * ABOVE1
X BESTB1 = BLOCK1
X BESTB2 = BLOCK2
X BESTP1 = POSI
X BESTP2 = POSJ
X CHANGE = .TRUE.
X END IF
X 460 CONTINUE
X 470 CONTINUE
X POS2 = POS2 + K - 1 - BLOCK2
X 480 CONTINUE
X POS1 = POS1 + K - 1 - BLOCK1
X 490 CONTINUE
C
C If CHANGE = .TRUE. then make the swap, otherwise the search ends.
C
X IF (CHANGE) THEN
X I = PICKED(BESTP1)
X J = PICKED(BESTP2)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BESTB2, XX, I)
X CALL MODTRI(K, NRBAR, ONE, XX, D, RBAR, TOL)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BESTB1, XX, J)
X CALL MODTRI(K, NRBAR, ONE, XX, D, RBAR, TOL)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BESTB1, XX, I)
X CALL MODTRI(K, NRBAR, MINUS1, XX, D, RBAR, TOL)
X CALL GETX(X, DIM1, KIN, NBLOCK, K, BESTB2, XX, J)
X CALL MODTRI(K, NRBAR, MINUS1, XX, D, RBAR, TOL)
X PICKED(BESTP1) = J
X PICKED(BESTP2) = I
X GO TO 430
X END IF
C
C
C Calculate log of determinant.
C
X 500 LNDET = ZERO
X DO 510 I = 1, K
X 510 LNDET = LNDET + LOG(D(I))
X RETURN
X END
C
C
X SUBROUTINE MODTRI(NP, NRBAR, WEIGHT, XROW, D, RBAR, TOL)
C
C ALGORITHM AS75.1 with minor modifications
C
C Modify a triangular (Cholesky) decomposition.
C Calling this routine updates d and rbar by adding another design
C point with weight = WEIGHT, which may be negative.
C
C *** WARNING Array XROW is overwritten.
C
X INTEGER NP, NRBAR
X DOUBLE PRECISION WEIGHT, XROW(NP), D(NP), RBAR(*), TOL(NP)
C
C Local variables
C
X INTEGER I, K, NEXTR
X DOUBLE PRECISION ZERO, W, XI, DI, WXI, DPI, CBAR, SBAR, XK
C
X DATA ZERO/0.D0/
C
X W = WEIGHT
X NEXTR = 1
X DO 30 I = 1, NP
C
C Skip unnecessary transformations. Test on exact zeroes must be
C used or stability can be destroyed.
C
X IF (W .EQ. ZERO) RETURN
C
X XI = XROW(I)
X IF (ABS(XI) .LT. TOL(I)) THEN
X NEXTR = NEXTR + NP - I
X GO TO 30
X END IF
C
X DI = D(I)
X WXI = W*XI
X DPI = DI + WXI*XI
C
C Test for new singularity.
C
X IF (DPI .LT. TOL(I)) THEN
X DPI = ZERO
X CBAR = ZERO
X SBAR = ZERO
X W = ZERO
X ELSE
X CBAR = DI/DPI
X SBAR = WXI/DPI
X W = CBAR*W
X END IF
C
X D(I) = DPI
X DO 20 K = I + 1, NP
X XK = XROW(K)
X XROW(K) = XK - XI*RBAR(NEXTR)
X RBAR(NEXTR) = CBAR*RBAR(NEXTR) + SBAR*XK
X NEXTR = NEXTR + 1
X 20 CONTINUE
X 30 CONTINUE
C
X RETURN
X END
C
C
X SUBROUTINE MODTR2(NP, NRBAR, XROW, D, RBAR, TOL, RANK, LNDET)
C
C ALGORITHM AS75.1 with modifications.
C Calculate the effect of an update of a QR factorization upon the
C rank and determinant, without changing D or RBAR.
C
C *** WARNING Array XROW is overwritten.
C
X INTEGER NP, NRBAR, RANK
X DOUBLE PRECISION XROW(NP), D(NP), RBAR(*), TOL(NP), LNDET
C
C Local variables
C
X INTEGER I, K, NEXTR, J
X DOUBLE PRECISION ZERO, ONE, W, XI, DI, WXI, DPI, CBAR, XK
C
X DATA ZERO/0.D0/, ONE/1.D0/
C
X W = ONE
X RANK = 0
X LNDET = ZERO
X NEXTR = 1
X DO 30 I = 1, NP
C
C Skip unnecessary transformations. Test on exact zeroes must be
C used or stability can be destroyed.
C
X IF (W .EQ. ZERO) THEN
X DO 10 J = I, NP
X IF (D(J) .GT. TOL(J)) THEN
X RANK = RANK + 1
X LNDET = LNDET + LOG(D(J))
X END IF
X 10 CONTINUE
X RETURN
X END IF
C
X XI = XROW(I)
X IF (ABS(XI) .LT. TOL(I)) THEN
X IF (D(I) .GT. TOL(I)) THEN
X RANK = RANK + 1
X LNDET = LNDET + LOG(D(I))
X END IF
X NEXTR = NEXTR + NP - I
X GO TO 30
X END IF
C
X DI = D(I)
X WXI = W*XI
X DPI = DI + WXI*XI
C
C Test for new singularity.
C
X IF (DPI .LT. TOL(I)) THEN
X DPI = ZERO
X CBAR = ZERO
X W = ZERO
X ELSE
X CBAR = DI/DPI
X W = CBAR*W
X LNDET = LNDET + LOG(DPI)
X RANK = RANK + 1
X END IF
C
X DO 20 K = I + 1, NP
X XK = XROW(K)
X XROW(K) = XK - XI*RBAR(NEXTR)
X NEXTR = NEXTR + 1
X 20 CONTINUE
X 30 CONTINUE
C
X RETURN
X END
C
X SUBROUTINE GETX(X, DIM1, KIN, NBLOCK, K, BLOCK, XX, CASE)
C
C Copy one case from X to XX.
C
X INTEGER DIM1, KIN, NBLOCK, K, BLOCK, CASE
X REAL X(DIM1, KIN)
X DOUBLE PRECISION XX(K)
C
C Local variables.
C
X INTEGER I, J
X DOUBLE PRECISION ZERO, ONE
X DATA ZERO/0.D0/, ONE/1.D0/
C
X DO 10 I = 1, NBLOCK
X IF (I .NE. BLOCK) THEN
X XX(I) = ZERO
X ELSE
X XX(I) = ONE
X END IF
X 10 CONTINUE
X J = NBLOCK + 1
X DO 20 I = 1, KIN
X XX(J) = X(CASE, I)
X J = J + 1
X 20 CONTINUE
C
X RETURN
X END
C
X SUBROUTINE BKSUB1(RBAR, NRBAR, K, RHS, SOLN, TOL, D)
C
C Solves D R y = z for y (SOLN), where z = RHS.
C RBAR is an upper-triangular matrix with implicit 1's on it's
C diagonal, stored by rows.
C
X INTEGER NRBAR, K
X DOUBLE PRECISION RBAR(NRBAR), RHS(K), SOLN(K), TOL(K), D(K)
C
C Local variables
C
X INTEGER ROW, COL, POS
X DOUBLE PRECISION ZERO, TEMP
X DATA ZERO /0.D0/
C
X POS = K * (K - 1) / 2
X DO 20 ROW = K, 1, -1
X IF (D(ROW) .GT. TOL(ROW)) THEN
X TEMP = RHS(ROW) / D(ROW)
X DO 10 COL = K, ROW+1, -1
X TEMP = TEMP - RBAR(POS) * SOLN(COL)
X POS = POS - 1
X 10 CONTINUE
X SOLN(ROW) = TEMP
X ELSE
X POS = POS - K + ROW
X SOLN(ROW) = ZERO
X END IF
X 20 CONTINUE
C
X RETURN
X END
C
C
X SUBROUTINE BKSUB2(RBAR, NRBAR, K, RHS, SOLN)
C
C Solves R'(sqrt(D).z) = x where (sqrt(D).z) = SOLN and x = RHS.
C RBAR is an upper-triangular matrix with implicit 1's on it's
C diagonal, stored by rows.
C
X INTEGER NRBAR, K
X DOUBLE PRECISION RBAR(NRBAR), RHS(K), SOLN(K)
C
C Local variables
C
X INTEGER ROW, COL, POS
X DOUBLE PRECISION TEMP
C
X SOLN(1) = RHS(1)
X DO 20 ROW = 2, K
X TEMP = RHS(ROW)
X POS = ROW - 1
X DO 10 COL = 1, ROW-1
X TEMP = TEMP - RBAR(POS) * SOLN(COL)
X POS = POS + K - COL - 1
X 10 CONTINUE
X SOLN(ROW) = TEMP
X 20 CONTINUE
C
X RETURN
X END
C
X DOUBLE PRECISION FUNCTION DELTA(K, XJ, XL, XBARI, XBARK, NI, NK,
X + NRBAR, D, RBAR, Z, A, B, DIFF)
C
C Calculate the delta function for the swap of case J in block I
C with case L in block K. Uses the method of Cook & Nachtsheim.
C
X INTEGER K, NI, NK, NRBAR
X DOUBLE PRECISION XJ(K), XL(K), XBARI(K), XBARK(K), D(K),
X + RBAR(NRBAR), Z(K,3), A(K), B(K), DIFF(K)
X DOUBLE PRECISION DOTP
C
C Local variables
C
X INTEGER I
X DOUBLE PRECISION CONST, ONE, TWO, TEMP, E11, E12, E21, E22
X DATA ONE/1.D0/, TWO/2.D0/
C
C Calculate vectors DIFF, A & B.
C
X CONST = TWO - ONE/NI - ONE/NK
X DO 10 I = 1, K
X TEMP = XJ(I) - XL(I)
X DIFF(I) = -TEMP
X A(I) = TEMP - XBARI(I) + XBARK(I)
X B(I) = A(I) - CONST * TEMP
X 10 CONTINUE
C
C Calculate the z-vectors by back-substitution.
C Z1 for A, Z2 for B and Z3 for DIFF.
C N.B. The solutions returned from BKSUB2 have the I-th element
C multiplied by sqrt(D(I)).
C
X CALL BKSUB2(RBAR, NRBAR, K, A, Z(1,1))
X CALL BKSUB2(RBAR, NRBAR, K, B, Z(1,2))
X CALL BKSUB2(RBAR, NRBAR, K, DIFF, Z(1,3))
C
C Calculate the elements E11, E12, E21 & E22 as dot-products of the
C appropriate z-vectors.
C
X E11 = DOTP(K, Z(1,3), Z(1,1), D)
X E12 = DOTP(K, Z(1,3), Z(1,3), D)
X E21 = DOTP(K, Z(1,2), Z(1,1), D)
X E22 = DOTP(K, Z(1,2), Z(1,3), D)
C
C Return the determinant of the matrix: E11+1 E12
C (minus 1) E21 E22+1
C
X DELTA = (E11+ONE) * (E22+ONE) - E12 * E21 - ONE
X RETURN
X END
C
C
X DOUBLE PRECISION FUNCTION DOTP(K, X, Y, D)
C
C [Formerly named DOTPRD, but renamed here to avoid loader table
C conflicts. W. Venables.]
C
C Dot-product scaled by vector D.
C
X INTEGER K
X DOUBLE PRECISION X(K), Y(K), D(K)
C
C Local variables
C
X INTEGER I
X DOUBLE PRECISION ZERO
X DATA ZERO/0.D0/
C
X DOTP = ZERO
X DO 10 I = 1, K
X 10 DOTP = DOTP + X(I) * Y(I) / D(I)
X RETURN
X END
C
C
X SUBROUTINE SINGM(NP, NRBAR, D, RBAR, TOL, WORK, IFAULT)
C
C Modified from:
C ALGORITHM AS274.5 APPL. STATIST. (1992) VOL.41, NO.2
C
C Checks for singularities, and adjusts orthogonal
C reductions produced by AS75.1.
C
C Auxiliary routine called: MODTRI
C
X INTEGER NP, NRBAR, IFAULT
X DOUBLE PRECISION D(NP), RBAR(NRBAR), TOL(NP), WORK(NP)
C
C Local variables
C
X DOUBLE PRECISION ZERO, TEMP
X INTEGER COL, POS, ROW, NP2, POS2, J
C
X DATA ZERO/0.D0/
C
C Check input parameters
C
X IFAULT = 0
X IF (NP .LE. 0) IFAULT = 1
X IF (NRBAR .LT. NP*(NP-1)/2) IFAULT = IFAULT + 2
X IF (IFAULT .NE. 0) RETURN
C
X DO 10 COL = 1, NP
X 10 WORK(COL) = SQRT(D(COL))
C
X DO 40 COL = 1, NP
C
C Set elements within RBAR to zero if they are less than TOL(COL) in
C absolute value after being scaled by the square root of their row
C multiplier.
C
X TEMP = TOL(COL)
X POS = COL - 1
X DO 20 ROW = 1, COL-1
X IF (ABS(RBAR(POS)) * WORK(ROW) .LT. TEMP) RBAR(POS) = ZERO
X POS = POS + NP - ROW - 1
X 20 CONTINUE
C
C If diagonal element is near zero, set it to zero, and use MODTRI
C to augment the projections in the lower rows of the factorization.
C
X IF (WORK(COL) .LE. TEMP) THEN
X IFAULT = IFAULT - 1
X IF (COL .LT. NP) THEN
X NP2 = NP - COL
X POS2 = POS + NP - COL + 1
X IF (NP2 .GT. 1) THEN
X CALL MODTRI(NP2, NP2*(NP2-1)/2, D(COL), RBAR(POS+1),
X + D(COL+1), RBAR(POS2), TOL)
X ELSE
X CALL MODTRI(1, 0, D(COL), RBAR(POS+1), D(COL+1), RBAR(1),
X + TOL)
X END IF
X DO 30 J = POS+1, POS2-1
X 30 RBAR(J) = ZERO
X END IF
X D(COL) = ZERO
X END IF
X 40 CONTINUE
X RETURN
X END
C
X SUBROUTINE XXTR(NP, NRBAR, D, RBAR, NREQ, TRACE, RINV)
X
C
C Calculate the trace of the inverse of X'X (= R'R).
C
X INTEGER NP, NRBAR, NREQ
X DOUBLE PRECISION D(NP), RBAR(NRBAR), TRACE, RINV(*)
C
C Local variables
C
X INTEGER POS, ROW, COL
X DOUBLE PRECISION ONE, ZERO
X DATA ONE/1.D0/, ZERO/0.D0/
C
C Get the inverse of R
C
X CALL INV(NP, NRBAR, RBAR, NREQ, RINV)
C
C Trace = the sum of the diagonal elements of RINV * (1/D) * (RINV)'
C
X TRACE = ZERO
X POS = 1
X DO 20 ROW = 1, NREQ
X TRACE = TRACE + ONE / D(ROW)
X DO 10 COL = ROW+1, NREQ
X TRACE = TRACE + RINV(POS)**2 / D(COL)
X POS = POS + 1
X 10 CONTINUE
X 20 CONTINUE
C
X RETURN
X END
C
C
X SUBROUTINE INV(NP, NRBAR, RBAR, NREQ, RINV)
C
C ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2
C
C Invert first NREQ rows and columns of Cholesky factorization
C produced by AS 75.1.
C
X INTEGER NP, NRBAR, NREQ
X REAL RBAR(*), RINV(*)
C
C Local variables.
C
X INTEGER POS, ROW, COL, START, K, POS1, POS2
X REAL SUM, ZERO
C
X DATA ZERO/0.0/
C
C Invert RBAR ignoring row multipliers, from the bottom up.
C
X POS = NREQ * (NREQ-1)/2
X DO 30 ROW = NREQ-1, 1, -1
X START = (ROW-1) * (NP+NP-ROW)/2 + 1
X DO 20 COL = NREQ, ROW+1, -1
X POS1 = START
X POS2 = POS
X SUM = ZERO
X DO 10 K = ROW+1, COL-1
X POS2 = POS2 + NREQ - K
X SUM = SUM - RBAR(POS1) * RINV(POS2)
X POS1 = POS1 + 1
X 10 CONTINUE
X RINV(POS) = SUM - RBAR(POS1)
X POS = POS - 1
X 20 CONTINUE
X 30 CONTINUE
C
X RETURN
X END
X
C-------------------------------------------------------------------------
X
X SUBROUTINE CLEAR(NP, NRBAR, D, RBAR, IER)
C
C Modified version of:
C ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2
C
C Sets arrays to zero prior to calling INCLUD
C
X INTEGER NP, NRBAR, IER
X DOUBLE PRECISION D(NP), RBAR(*)
C
C Local variables
C
X INTEGER I
X DOUBLE PRECISION ZERO
C
X DATA ZERO/0.D0/
C
C Some checks.
C
X IER = 0
X IF (NP .LT. 1) IER = 1
X IF (NRBAR .LT. NP*(NP-1)/2) IER = IER + 2
X IF (IER .NE. 0) RETURN
C
X DO 10 I = 1, NP
X D(I) = ZERO
X 10 CONTINUE
X DO 20 I = 1, NRBAR
X 20 RBAR(I) = ZERO
X RETURN
X END
C
C
X SUBROUTINE REGCF(NP, NRBAR, D, RBAR, THETAB, TOL, BETA,
X + NREQ, IER)
C
C ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2
C
C Modified version of AS75.4 to calculate regression coefficients
C for the first NREQ variables, given an orthogonal reduction from
C AS75.1.
C
X INTEGER NP, NRBAR, NREQ, IER
X DOUBLE PRECISION D(NP), RBAR(*), THETAB(NP), TOL(NP),
X + BETA(NP)
C
C Local variables
C
X INTEGER I, J, NEXTR
X DOUBLE PRECISION ZERO
C
X DATA ZERO/0.D0/
C
C Some checks.
C
X IER = 0
X IF (NP .LT. 1) IER = 1
X IF (NRBAR .LT. NP*(NP-1)/2) IER = IER + 2
X IF (NREQ .LT. 1 .OR. NREQ .GT. NP) IER = IER + 4
X IF (IER .NE. 0) RETURN
C
X DO 20 I = NREQ, 1, -1
X IF (SQRT(D(I)) .LT. TOL(I)) THEN
X BETA(I) = ZERO
X D(I) = ZERO
X GO TO 20
X END IF
X BETA(I) = THETAB(I)
X NEXTR = (I-1) * (NP+NP-I)/2 + 1
X DO 10 J = I+1, NREQ
X BETA(I) = BETA(I) - RBAR(NEXTR) * BETA(J)
X NEXTR = NEXTR + 1
X 10 CONTINUE
X 20 CONTINUE
C
X RETURN
X END
C
C
C
C
X REAL FUNCTION RAND()
C
C The Wichmann & Hill random number generator
C Algorithm AS183, Appl. Statist., 31, 188-190, 1982.
C The cycle length is 6.95E+12.
C This random number generator is very slow compared with most
C others, but it is dependable, and the results are reproducible.
C
X INTEGER*4 IX, IY, IZ
X COMMON /RANDC/ IX, IY, IZ
C
X IX = MOD(171*IX, 30269)
X IY = MOD(172*IY, 30307)
X IZ = MOD(170*IZ, 30323)
X RAND = MOD(FLOAT(IX)/30269. + FLOAT(IY)/30307. +
X 1 FLOAT(IZ)/30323. , 1.0)
X RETURN
X END
X
X SUBROUTINE SETRAN(INIT)
C
C Simple initializer for use in S-PLUS routines - WNV April, 1995
C
X INTEGER*4 INIT(3), IX, IY, IZ
X COMMON /RANDC/ IX, IY, IZ
X IX = INIT(1)
X IY = INIT(2)
X IZ = INIT(3)
X RETURN
X END
X
C subroutine tesran(n, x)
C integer n
C real rand
C double precision x(*)
C do i = 1,n
C x(i) = rand()
C end do
C return
C end
X
X
X
SHAR_EOF
chmod 0644 Dopt.f ||
echo 'restore of Dopt.f failed'
Wc_c="`wc -c < 'Dopt.f'`"
test 28206 -eq "$Wc_c" ||
echo 'Dopt.f: original size 28206, current size' "$Wc_c"
fi
# ============= Dopt.q ==============
if test -f 'Dopt.q' -a X"$1" != X"-c"; then
echo 'x - skipping Dopt.q (File already exists)'
else
echo 'x - extracting Dopt.q (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'Dopt.q' &&
Dopt <- function(treatments,
X blocks,
X data.treatments = sys.parent(),
X data.blocks = sys.parent(),
X initial.points = rep(0, npoints),
X fixed.points = rep(F, npoints),
X RNG = round(runif(3) * 30269),
X eps = .Machine$single.eps,
X runs)
{
X if(missing(blocks)) {
X if(missing(runs))
X stop("Either blocks or runs must be specified")
X else blocks <- runs
X }
X if(inherits(blocks, "formula")) {
X vars <- attr(terms(blocks, data = data.blocks), "variables")
X block.data.frame <- as.data.frame(lapply(as.list(vars),
X function(x, dat)
X as.factor(eval(as.expression(x), local = dat)),
X data.blocks))
X names(block.data.frame) <- as.character(vars)
X vars <- lapply(unclass(block.data.frame),
X function(s) format(as.character(s)))
X bf <- factor(do.call("paste", c(vars, list(sep = ""))))
X } else {
X if(length(blocks) == 1) {
X bf <- rep(1, as.numeric(blocks))
X block.data.frame <- data.frame("(Intercept)" = bf)
X }
X else {
X bf <- as.factor(blocks)
X block.data.frame <- data.frame(Blocks = bf)
X }
X }
X npoints <- length(bf)
X if(missing(treatments))
X stop("treatments must be specified")
X if(!inherits(treatments, "formula")) {
X if(length(treatments) == 1)
X Treatments <- as.factor(1:as.numeric(treatments))
X else {
X Treatments <- eval(treatments, data.treatments)
X }
X treatments <- ~ Treatments
X data.treatments <- sys.frame()
X }
X X <- model.matrix(treatments, data.treatments)
X intercept <- any(as.logical(match(dimnames(X)[[2]],
X "(Intercept)", nomatch = 0)))
X if(intercept) {
X X <- scale(X, scale = F)
X bf <- factor(as.numeric(bf))
X nblocks <- length(levels(bf))
X if(nblocks == npoints)
X stop("All blocks have only one plot.")
X blocksizes <- table(bf)
X if(nblocks == 1)
X names(block.data.frame) <- "(Intercept)"
X }
X else {
X nblocks <- 0
X blocksizes <- npoints
X block.data.frame <- data.frame(.Intercept. = bf)
X o <- 1:npoints
X }
X s <- svd(X)
X X <- s$u[, s$d > eps]
X ncand <- dim(X)[1]
X kin <- dim(X)[2]
X k <- kin + nblocks
X nrbar <- (k * (k - 1))/2
X X <- X * sqrt(ncand/npoints)
X storage.mode(X) <- "single"
X rstart <- missing(initial.points)
X if(!rstart) {
X initial.points <- round(as.numeric(initial.points))
X fixed.points <- as.logical(fixed.points)
X length(initial.points) <- npoints
X length(fixed.points) <- npoints
X initial.points[is.na(initial.points)] <- 0
X fixed.points[is.na(fixed.points)] <- F
X if(any(initial.points < 0 | initial.points > ncand))
X stop("initial points invalidly specified")
X }
X in.nbl <- 0
X m <- sum(where <- (initial.points == 0))
X initial.points[where] <- sample(1:ncand, m, rep = T)
X if(nblocks > 0) {
X o <- order(bf, !fixed.points)
X block.data.frame <- block.data.frame[o, , drop = F]
X initial.points <- initial.points[o]
X fixed.points <- fixed.points[o]
X in.nbl <- tapply(fixed.points, bf, sum)
X }
X .Fortran("setran",
X as.integer(RNG))
X fo <- .Fortran("fedrov",
X X = X,
X dim1 = as.integer(ncand),
X ncand = as.integer(ncand),
X kin = as.integer(kin),
X n = as.integer(npoints),
X nblock = as.integer(nblocks),
X in.nbl = as.integer(in.nbl),
X blocksiz = as.integer(blocksizes),
X k = as.integer(k),
X rstart = as.logical(rstart),
X nrbar = as.integer(nrbar),
X d = double(k),
X rbar = double(nrbar),
X picked = as.integer(initial.points),
X lndet = double(1),
X xx = double(k),
X tol = double(k),
X zpz = double(ncand * max(1, nblocks)),
X xk = double(k),
X ifault = integer(1))
X if(fo$ifault < 0)
X stop("No full rank starting design found.")
X if(fo$ifault != 0)
X stop("Esoteric fault in fedorov. Complain!")
X X <- X[fo$picked, ]
X if(intercept)
X if(nblocks > 1)
X X <- qr.resid(qr(model.matrix( ~ bf,
X data = sys.frame())), X)
X else scale(X, scale = F)
X cre <- sort(svd(X, 0, 0)$d^2)
X vars <- attr(terms(treatments, data = data.treatments),
X "variables")
X simple <- all(sapply(as.character(vars), function(x, dat)
X exists(x, where = dat) || exists(x), data.treatments))
X if(simple) {
X treatment.data.frame <- as.data.frame(lapply(
X as.list(vars), function(x, dat)
X eval(as.expression(x), local = dat),
X data.treatments))[fo$picked, , drop = F]
X names(treatment.data.frame) <- as.character(vars)
X }
X else {
X treatment.data.frame <- as.data.frame(
X model.matrix(treatments,
X data.treatments)[fo$picked, , drop = F])
X }
X dat <- cbind(block.data.frame, treatment.data.frame,
X Point.No = fo$picked)
X if(any(i <- match("(Intercept)", names(dat), nomatch = 0)))
X dat <- dat[, - i]
X if(any(i <- match(".Intercept.", names(dat), nomatch = 0)))
X dat <- dat[, - i]
X if(any(fixed.points))
X dat$Fixed <- fixed.points
X attr(dat, "log(det)") <- fo$lndet
X attr(dat, "RNG") <- RNG
X attr(dat, "CRE") <- cre[cre > eps]
X dat[order(o), ]
}
SHAR_EOF
chmod 0644 Dopt.q ||
echo 'restore of Dopt.q failed'
Wc_c="`wc -c < 'Dopt.q'`"
test 5660 -eq "$Wc_c" ||
echo 'Dopt.q: original size 5660, current size' "$Wc_c"
fi
# ============= First.lib.q ==============
if test -f 'First.lib.q' -a X"$1" != X"-c"; then
echo 'x - skipping First.lib.q (File already exists)'
else
echo 'x - extracting First.lib.q (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'First.lib.q' &&
".First.lib"<-
function(library, section)
{
X obj <- paste(library, "/", section, "/", "Dopt.o", sep = "")
X invisible(dyn.load(obj))
}
SHAR_EOF
chmod 0644 First.lib.q ||
echo 'restore of First.lib.q failed'
Wc_c="`wc -c < 'First.lib.q'`"
test 134 -eq "$Wc_c" ||
echo 'First.lib.q: original size 134, current size' "$Wc_c"
fi
# ============= First.lib2.q ==============
if test -f 'First.lib2.q' -a X"$1" != X"-c"; then
echo 'x - skipping First.lib2.q (File already exists)'
else
echo 'x - extracting First.lib2.q (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'First.lib2.q' &&
".First.lib"<-
function(library, section)
{
X obj <- paste(library, "/", section, "/", "Dopt.o", sep = "")
X invisible(dyn.load2(obj))
}
SHAR_EOF
chmod 0644 First.lib2.q ||
echo 'restore of First.lib2.q failed'
Wc_c="`wc -c < 'First.lib2.q'`"
test 135 -eq "$Wc_c" ||
echo 'First.lib2.q: original size 135, current size' "$Wc_c"
fi
# ============= INDEX ==============
if test -f 'INDEX' -a X"$1" != X"-c"; then
echo 'x - skipping INDEX (File already exists)'
else
echo 'x - extracting INDEX (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'INDEX' &&
Dopt --> Function to construct approximately D-optimal
X experimental designs from candidate points. Only simple
X blocking is allowed but arbitrary treatment structures
X are accommodated.
SHAR_EOF
chmod 0644 INDEX ||
echo 'restore of INDEX failed'
Wc_c="`wc -c < 'INDEX'`"
test 189 -eq "$Wc_c" ||
echo 'INDEX: original size 189, current size' "$Wc_c"
fi
# ============= INS ==============
if test -f 'INS' -a X"$1" != X"-c"; then
echo 'x - skipping INS (File already exists)'
else
echo 'x - extracting INS (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'INS' &&
#! /bin/sh
###
### Adjust the following as necessary
X
S=Splus # command to run S-PLUS
X
First=First.lib.q # if using dyn.load
#First=First.lib2.q # if using dyn.load2
X
### No need to adjust below here on most UNIX systems - I hope!
X
PATH=/bin:$PATH; export PATH
X
if test ! -d 'Dopt'
then
X echo making directory Dopt...
X mkdir Dopt
fi
X
cd Dopt
X
if test ! -d '.Data'
then
X echo creating directory "'Dopt/.Data'"...
X mkdir '.Data'
fi
X
cd '.Data'
X
if test ! -d '.Help'
then
X echo creating directory "'Dopt/.Data/.Help'"...
X mkdir '.Help'
fi
cd ..
X
echo copying files to Dopt...
cp ../Dopt.[fdq] ../$First ../README ../NOTES ../INS .
X
echo compiling Fortran...
X
$S COMPILE Dopt.f
X
echo installing functions...
X
cat Dopt.q $First | $S 2>&1 > install.out || exit 1
cp Dopt.d .Data/.Help/Dopt
Splus help.findsum .Data
X
echo you may now test the library.
X
Xexit 0
SHAR_EOF
chmod 0755 INS ||
echo 'restore of INS failed'
Wc_c="`wc -c < 'INS'`"
test 859 -eq "$Wc_c" ||
echo 'INS: original size 859, current size' "$Wc_c"
fi
# ============= INSTALLATION ==============
if test -f 'INSTALLATION' -a X"$1" != X"-c"; then
echo 'x - skipping INSTALLATION (File already exists)'
else
echo 'x - extracting INSTALLATION (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'INSTALLATION' &&
REQUIREMENTS:
X
To install this function you will need
X
1. The dynamic loading facility, either dyn.load() or dyn.load2()
X
2. Access to a Fortran compiler that generates dynamically
X loadable binaries for your version of S-PLUS
X
3. Something at least functinally equivalent to the UNIX Bourne
X shell, sh, to decouple the shell archive file (but if you are
X reading this, I suppose you have got that far!)
X
-----------------------------------------------------------------
X
INSTALLATION and TESTING:
X
The recommended (UNIX) proceedure is as follows:
X
A. Make a new, empty directory, move the shell archive file to it
X and expand it there. For example
X
X $ mkdir dopt
X $ mv Dopt.shar dopt
X $ cd dopt
X $ sh Dopt.shar
X
B. If your S-PLUS command is `Splus' and your system uses
X `dyn.load()', fine, otherwise you will need to adjust the two
X lines at the top of the shell script INS where these are
X defined.
X
X The test shell script TEST also needs to have the S-PLUS
X command corrected.
X
C. To install the Dopt library as a subdirectory of the current
X directory:
X
X $ ./INS
X
D. If everything goes fine you should test the library with
X
X $ ./TEST
X
X This will generate an output file `test.out' that should
X only differ with the comparison file `master.out' in the
X banner at the top.
X
X There can be major differences in the files, though, and the
X function still be correct, simply because of differences in
X floating point arithmetic on your machine and on mine.
X
E. If everything seems OK, move the Dopt subdirectory of your
X temporary working directory (and all its contents, of course)
X to $SHOME/library, or to some convenient place where users may
X have access to it.
X
-----------------------------------------------------------------
X
PROBLEMS and POTENTIAL Problems:
X
1. I have only tested this on Sun SPARC IPCs using SunOS 4.1.1
X and 4.1.3, and with S-PLUS Rel. 3.1 and 3.2.
X
2. It is easy to underestimate the size of the computational
X problem for many designs. The last part of the test procedure
X is deliberately a fairly large problem, but much larger ones
X occur regularly in practice. This can generate huge storage
X and virtual memory requirements and may cause some irritation
X if you work on a shared memory system.
X
3. The routine purports to use its own internal random number
X generator and should therefore have its results reproducible
X if the seed is reset beforehand. This appears not to be so;
X for some odd reason to make the results reproducible you need
X both to set the RNG local to Dopt and the .Random.seed vector
X for the S-PLUS RNG. I do not know if this is a bug in my
X front end or some problem with the Fortran, but so far it has
X caused no problems in getting good designs.
X
_________________________________________________________________
Bill Venables, Department of Statistics, Tel: +61 8 303 3026
The University of Adelaide, Fax: +61 8 232 5670
South AUSTRALIA. 5005. Email: Bill.Venables@adelaide.edu.au
SHAR_EOF
chmod 0644 INSTALLATION ||
echo 'restore of INSTALLATION failed'
Wc_c="`wc -c < 'INSTALLATION'`"
test 3054 -eq "$Wc_c" ||
echo 'INSTALLATION: original size 3054, current size' "$Wc_c"
fi
# ============= NOTES ==============
if test -f 'NOTES' -a X"$1" != X"-c"; then
echo 'x - skipping NOTES (File already exists)'
else
echo 'x - extracting NOTES (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'NOTES' &&
Extract from notes supplied by the authors of the original fortran program.
X
X
[NOTE: The fortran function DOPT has been renamed FEDROV in the version for
X use in dynamic loading in S-PLUS to avoid name conflicts.
X W. Venables.]
X
X D O P T
X _______
X
X by Alan Miller & Nam-Ky Nguyen
X
DOPT is an algorithm for generating experimental designs which are close to
being D-optimal designs, using Fedorov's exchange algorithm. The algorithm
allows for blocking of designs in blocks which may vary in size. It is
particularly suitable for producing fractional factorial designs when some
of the factors have more than two levels. Such cases usually require a
large number of experimental runs and it is rarely possible to complete
them all in one batch, and hence the design often needs to be blocked.
Where there is more than one block, an algorithm due to Cook & Nachtsheim
is used to swap design points between blocks to attempt to increase the
determinant.
X
X
Use of DOPT
___________
X
To use DOPT in any other program, the user must supply a 2-D array, X, of
candidate points from which DOPT will select n points, possibly with
duplication. The candidate points could be read in from a file, or
generated in the program. See the programs DRIVER, IBQUAD and DYKSTRA for
examples.
X
The optimum design depends upon the model to be fitted. DOPT attempts to
maximize the determinant of X'X for the selected design points. Here X
consists of only those rows corresponding to the selected points. If say a
quadratic surface:
X 2 2
X z = a + a x + a y + a x + a xy + a y
X 0 1 2 3 4 5
X
is to be fitted, then the input (or generated) matrix X of candidate points
might be:
X 2 2
X (x) (y) (x ) (xy) (y )
X
X -1 -1 +1 +1 +1
X -1 0 +1 0 0
X -1 +1 +1 -1 +1
X 0 -1 0 0 +1
X 0 0 0 0 0
X 0 +1 0 0 +1
X +1 -1 +1 -1 +1
X +1 0 +1 0 0
X +1 +1 +1 +1 +1
X
If you want to fit a cubic term in x then an extra column for the cube of x
can be added at the right-hand side (or anywhere else!), but at least four
different values of x would be needed for a cubic.
X
In simple designs, such as for fitting quadratic curves, or for estimating
effects of factors which have only two levels, values such as -1, 0 and +1
can be used for X and then translated into say low, medium and high actual
values for the real life variables.
X
The arguments for DOPT are as follows:
X
XX(DIM1, NCAND) Real, input Array of candidate points from which
X the design points will be chosen.
X
DIM1 Integer, input 1st dimension of X in the calling program.
X
NCAND Integer, input Number of candidate points.
X
KIN Integer, input The number of columns of the design matrix.
X
N Integer, input The number of points required in the design.
X
NBLOCK Integer, input The number of blocks in the design. If no
X constant is to be fitted in the model, set
X NBLOCK = 0.
X
IN(NBL) Integer, input The number of points, if any, to be forced
X into each block. N.B. If NBLOCK = 0, the
** NBL = max(NBLOCK, 1) ** value of IN(1) must be set.
X
BLKSIZ(NBL) Integer, input The number of design points in each block.
X
K Integer, input KIN + the number of blocks.
X
RSTART Logical, input = .TRUE. for DOPT to generate a random
X starting design,
X = .FALSE. if a starting design is being
X input in array PICKED (see below). This
X option would be used if you want to see if
X any improvement can be made to an existing
X design.
X
NRBAR Integer, input The dimension of RBAR. Must be at least
X K(K-1)/2.
X
D(K) D.prec., output Array of row multipliers for the Cholesky
X factorization of the final X'X matrix.
X
RBAR(NRBAR) D.prec., output Array containing the upper triangle of the
X Cholesky factor stored by rows, excluding
X the diagonal elements which are all 1's.
X
PICKED(N) Integer, I/O 1. INPUT. PICKED(I) contains the position
X of the I-th design point in the input
X array of candidate points, X. Used
X when one or more values of IN(I) is > 0
X or when RSTART = .FALSE.
X 2. OUTPUT. PICKED(I) contains the positions
X of the selected design points in the input
X array of candidate points, X.
X
LNDET D.prec., output Logarithm of the determinant of X'X.
X
XXX(K) D.prec., workspace
X
TOL(K) D.prec., output Array of tolerances used for testing for
X singularities.
X
ZPZ(NCAND,NBL) D.prec., workspace
X
WK(K) D.prec., workspace
X
IFAULT Integer, output Error indicator
X = 0 if no errors detected
X = -1 if no full-rank starting design is
X found. Check your design specification
X if this error is reported.
X = 1 if DIM1 < NCAND
X = 2 if K < N
X = 4 if NRBAR < K(K-1)/2
X = 8 if K not equal to KIN + NBLOCK
X = 16 if sum of block sizes in array BLKSIZ
X is not equal to N
X = 32 if any IN(I) < 0 or IN(I) > BLKSIZ(I)
X
N.B. If more than one error is detected with a positive error indicator, then
the values above are summed. If a non-zero value for IFAULT is returned, it
could be caused by calling DOPT with the wrong number or types of arguments.
X
X
X
A note on memory requirements
_____________________________
X
The array X can be very large. e.g. for a (2^8) x (3^3) factorial, there
are 256 x 27 = 6912 candidate points, and the second dimension of this
array needs to be at least 58 to fit all 2-factor interactions and a
quadratic surface. If 4 bytes are required for each REAL value stored, the
minimum memory requirement for this array is 6912 x 58 x 4 = 1603584 bytes,
or just over 1.5Mbytes.
X
This space requirement can be reduced in three ways:
1. X can be declared as an INTEGER*2 or INTEGER*1 array if the compiler
X allows it. If so, the declaration for X must be changed in routines
X DOPT and GETX. This is of course only appropriate when X takes integer
X values.
X
2. Subroutine GETX is used to supply a requested candidate point to DOPT in
X vector XX. In the case above, the 2nd dimension of X could be reduced
X to 11 if routine GETX is re-written to calculate the interaction and
X quadratic terms. Notice that GETX puts 0 and 1 values into the first
X NBL positions of XX to indicate which block a design point is in.
X
3. The use of array X can be eliminated, and its elements can be generated
X in GETX as required.
X
Option 2 will slow down execution, and option 3 will make it very much slower.
X
------------------------------------------------------------------------------
X
For a reference on the methods used, prior to the publication in Applied
Statistics, see:
X
Nguyen, N-K. and Miller, A.J. (1992). A review of some exchange algorithms
for constructing discrete D-optimal designs. Comput. Statist. & Data Anal.,
vol.14, 489-498.
X
X
Please report errors or difficulties to:
X Alan Miller
X CSIRO Division of Mathematics & Statistics
X Private Mail Bag 10
X Rosebank MDC
X Clayton, Victoria 3169
X Australia
X
Telephone: +61 3 542-2266 Fax: +61 3 542-2474
e-mail : alan @ mel.dms.CSIRO.AU
X
SHAR_EOF
chmod 0644 NOTES ||
echo 'restore of NOTES failed'
Wc_c="`wc -c < 'NOTES'`"
test 8603 -eq "$Wc_c" ||
echo 'NOTES: original size 8603, current size' "$Wc_c"
fi
# ============= README ==============
if test -f 'README' -a X"$1" != X"-c"; then
echo 'x - skipping README (File already exists)'
else
echo 'x - extracting README (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'README' &&
Dopt: Generate D-optimal designs using the Fedorov exchange algorithm
SHAR_EOF
chmod 0644 README ||
echo 'restore of README failed'
Wc_c="`wc -c < 'README'`"
test 70 -eq "$Wc_c" ||
echo 'README: original size 70, current size' "$Wc_c"
fi
# ============= STATEMENT ==============
if test -f 'STATEMENT' -a X"$1" != X"-c"; then
echo 'x - skipping STATEMENT (File already exists)'
else
echo 'x - extracting STATEMENT (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'STATEMENT' &&
Redistribution:
X
Permission is given to re-distribute this software freely for any
non-commercial purpose.
X
The S front end was was written by Bill Venables, who retains the
copyright.
X
The Fortran program was written by Alan Miller and Nam-Ky Nguyen,
who have both given their acquiescence and support to the free
distribution of this material, however the copyright for the
fortran program itself is held by the Royal Statistical Society.
Any use of this software in commercial products must be done with
the permission of the Royal Statistical Society. A contact
address is given at the top of the fortran file. The Fortran
file is also available on statlib, with the permission of the
Royal Statistical Society.
X
Support:
X
I would be very pleased to hear of any improvements, suggestions
or bugs in this software. Anyone is free to contact me at the
email address given below, but I cannot undertake to provide a
support service.
X
Disclaimer:
X
No warranty is given on the results of this software. No claims
are made for its exactness or suitability for any purpose.
X
Bill Venables (Bill.Venables@Adelaide.edu.au)
21 April, 1995
SHAR_EOF
chmod 0644 STATEMENT ||
echo 'restore of STATEMENT failed'
Wc_c="`wc -c < 'STATEMENT'`"
test 1138 -eq "$Wc_c" ||
echo 'STATEMENT: original size 1138, current size' "$Wc_c"
fi
# ============= TEST ==============
if test -f 'TEST' -a X"$1" != X"-c"; then
echo 'x - skipping TEST (File already exists)'
else
echo 'x - extracting TEST (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'TEST' &&
#! /bin/sh
###
### Adjust the following as necessary
X
S=Splus # command to run S-PLUS
X
### No need to adjust below here on most UNIX systems - I hope!
X
PATH=/bin:$PATH; export PATH
X
if test ! -d '.Data'
then
X echo creating directory "'Dopt/.Data'"...
X mkdir '.Data'
fi
X
$S < test.q 2>&1 > test.out
X
echo Now compare '"test.out"' with '"master.out"'
diff test.out master.out
X
Xexit 0
SHAR_EOF
chmod 0755 TEST ||
echo 'restore of TEST failed'
Wc_c="`wc -c < 'TEST'`"
test 385 -eq "$Wc_c" ||
echo 'TEST: original size 385, current size' "$Wc_c"
fi
# ============= master.out ==============
if test -f 'master.out' -a X"$1" != X"-c"; then
echo 'x - skipping master.out (File already exists)'
else
echo 'x - extracting master.out (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'master.out' &&
S-PLUS : Copyright (c) 1988, 1993 Statistical Sciences, Inc.
S : Copyright AT&T.
Version 3.2 Release 1 for Sun SPARC, SunOS 4.x : 1993
Working data will be in .Data
X
X
Dopt library successfully loaded.
X
X
Cubic regression 0 <= x <= 100 design points:
X [1] 0 0 0 0 0 27 28 28 28 28 72 72 72 73 73 100 100
[18] 100 100 100
X
X
Quartic regression through origin, 0 <= x <= 100 design points
X [1] 17 17 17 17 17 50 50 50 50 50 82 83 83 83 83 100 100
[18] 100 100 100
X
X
Block design: v = b = 7, k = r = 4, concurrency matrix:
X 1 2 3 4 5 6 7
1 4 2 2 2 2 2 2
2 2 4 2 2 2 2 2
3 2 2 4 2 2 2 2
4 2 2 2 4 2 2 2
5 2 2 2 2 4 2 2
6 2 2 2 2 2 4 2
7 2 2 2 2 2 2 4
X
Canonical efficiencies:
[1] 0.875 0.875 0.875 0.875 0.875 0.875
X
X
2^14 main effect plan in 16 runs, with restrictions:
X A B C D E G H I J K L M N O Point.No Fixed
X 1 - - - - - - - - - - - - - - 1 T
X 2 - - - + + + - + + + - - - + 17 F
X 3 + + + - + + + - - - + - - - 122 F
X 4 + - + + - - - - + - - + + + 91 F
X 5 + - - - + + - - + + + - - + 73 F
X 6 - + + - + - - + - - + - + + 55 F
X 7 + + - + - + - - - + - + - + 109 F
X 8 - - - + + - + + + + - - + - 16 F
X 9 - + - - - + + + - + + + - - 38 F
10 - + + + + + + + - - - - - - 66 F
11 - + - - - - - + - + + + + + 35 F
12 + - - - - + + - + + + + - - 70 F
13 - - + - - - + + + - + + + - 20 F
14 - - + + - + - + + - - + - + 29 F
15 + + - + + - + - - + - - + - 112 F
16 + + + + + + + + + + + + + + 2 T
X
Crude factor replications:
X A B C D E G H I J K L M N O
+ 7 8 7 8 8 9 8 9 8 9 8 8 7 8
- 9 8 9 8 8 7 8 7 8 7 8 8 9 8
SHAR_EOF
chmod 0644 master.out ||
echo 'restore of master.out failed'
Wc_c="`wc -c < 'master.out'`"
test 1705 -eq "$Wc_c" ||
echo 'master.out: original size 1705, current size' "$Wc_c"
fi
#
touch -am 1231235999 $$.touch >/dev/null 2>&1
if test ! -f 1231235999 && test -f $$.touch; then
shar_touch=touch
else
shar_touch=:
echo 'WARNING: not restoring timestamps'
fi
rm -f 1231235999 $$.touch
#
# ============= test.q ==============
if test -f 'test.q' && test X"$1" != X"-c"; then
echo 'x - skipping test.q (File already exists)'
else
echo 'x - extracting test.q (text)'
sed 's/^X//' << 'SHAR_EOF' > 'test.q' &&
{if(ers <- exists(".Random.seed"))
X RS <- .Random.seed
X.Random.seed <- c(21, 14, 49, 32, 43, 1, 32, 22, 36, 23, 28, 3)
options(digits=7, width=75)
invisible(NULL)}
X
library(Dopt, lib.loc = unix("pwd"))
cat("\n\nDopt library successfully loaded.\n")
X
RNG <- c(10332, 7168, 20256)
X
x <- 0:100
dsn <- Dopt(~x + x^2 + x^3, 20, RNG = RNG)
cat("\n\nCubic regression 0 <= x <= 100 design points:\n")
sort(dsn$x)
X
x <- x/100
dsn <- Dopt(~x + x^2 + x^3 + x^4 - 1, 20, RNG = RNG + 1)
cat("\n\nQuartic regression through origin, 0 <= x <= 100 design points\n")
sort(dsn$Point.No) - 1
X
Bl <- factor(rep(1:7, 4))
Tr <- factor(1:7)
dsn <- Dopt(~Tr, ~Bl, RNG = RNG + 2)
CM <- crossprod(table(dsn$Bl, dsn$Tr))
cat("\n\nBlock design: v = b = 7, k = r = 4, concurrency matrix:\n")
CM
e <- sort(eigen(diag(7) - CM/(4 * 4), only.values = T)$va)[-1]
cat("\nCanonical efficiencies:\n")
e
X
dtr <- 0:1
for(i in 1:6) dtr <- rbind(cbind(0, dtr), cbind(1, dtr))
dtr <- rbind(0, rbind(1, cbind(dtr, 1-dtr)))
dtr <- as.data.frame(matrix(c("-","+")[dtr+1], nrow = nrow(dtr)))
names(dtr) <- LETTERS[-6][1:ncol(dtr)]
init <- c(1, rep(0, 14), 2)
dsn <- Dopt(~., 16, dtr, initial = init, fixed =
X as.logical(init), RNG = RNG + 3)
cat("\n\n2^14 main effect plan in 16 runs, with restrictions:\n")
dsn
cat("\nCrude factor replications:\n")
as.data.frame(lapply(dsn[,1:14], table))
X
{if(ers) .Random.seed <- RS
rm(dtr, dsn, x, Bl, Tr, e, CM, RNG, init, ers, RS)
invisible(NULL)}
SHAR_EOF
$shar_touch -am 0422102395 'test.q' &&
chmod 0644 'test.q' ||
echo 'restore of test.q failed'
shar_count="`wc -c < 'test.q'`"
test 1447 -eq "$shar_count" ||
echo "test.q: original size 1447, current size $shar_count"
fi
exit 0