Latin Hypercube Samples (lhs)
1.0
R, C++, and Rcpp code to generate Latin hypercube samples
|
#include <COrthogonalArray.h>
Public Member Functions | |
COrthogonalArray () | |
~COrthogonalArray () | |
void | addelkemp (int q, int k, int *n) |
void | addelkemp3 (int q, int k, int *n) |
void | addelkempn (int akn, int q, int k, int *n) |
void | bose (int q, int k, int *n) |
void | bosebush (int q, int k, int *n) |
void | bosebushl (int lambda, int q, int k, int *n) |
void | bush (int q, int k, int *n) |
void | busht (int str, int q, int k, int *n) |
int | oaagree (bool verbose) |
int | oatriple (bool verbose) |
void | oadimen () |
void | oarand (int is, int js, int ks, int ls) |
int | oastr (bool verbose) |
bool | oastr1 (bool verbose) |
bool | oastr2 (bool verbose) |
bool | oastr3 (bool verbose) |
bool | oastr4 (bool verbose) |
bool | oastrt (int t, bool verbose) |
int | getnrows () |
int | getncols () |
int | getq () |
const bclib::matrix< int > & | getoa () |
const std::string | getMessage () |
int | getReturnCode () |
Orthogonal Array Class
A collection of functions used as an API for Art Owen's oa library.
oacpp::COrthogonalArray::COrthogonalArray | ( | ) |
Default Constructor
|
inline |
Default Destructor
void oacpp::COrthogonalArray::addelkemp | ( | int | q, |
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Addelman Kempthorne algorithm
From the original documentation:
The addelkemp program produces
OA( 2q^2, k, q, 2 ), k <= 2q+1
, for odd prime powersq
. Even prime powers may be produced using bosebush above. This construction is based on:S. Addelman and O. Kempthorne (1961) Annals of Mathematical Statistics, Vol 32 pp 1167-1176.
using
n=2
in their notation.
2q
columns can be constructed without a coincidence defect. Settingk=2q+1
leads to an array with the coincidence defect. Some triples of columns contain duplicate rows. (The lack of a coincidence defect has been verified forq = 2,3,4,5,7,9,11,13,17,19,23,25
andk = 2q
.)This construction should work for all prime powers
q
, but it failed to do so for even powers greater than 4. This may have been a programming error, or it may have stemmed from misunderstanding of the description of the algorithm. The program rejects requests withq = 2^r
forr > 2
. The Bose Bush construction handles these cases.The description of the construction for odd prime powers calls for some arithmetic involving the number 4. In Galois fields with
3^r
elements, there is no 4. Replacing 4 by 1 for these fields works whenq = 3,9,27
(brute force verification).
q | the number of symbols (0,...,q-1) | |
k | the number of columns in the array. k <= q+1 | |
[out] | n | the number of rows in the array, n = 2q^2 |
std::runtime_error |
void oacpp::COrthogonalArray::addelkemp3 | ( | int | q, |
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Addelman Kempthorne algorithm
From the original documentation:
The addelkemp3 program produces
OA( 2*q^3, k, q, 2 ), k <= 2q^2+2q+1
, for prime powersq
.q
may be an odd prime power, orq
may be 2 or 4.This construction is based on:
S. Addelman and O. Kempthorne (1961) Annals of Mathematical Statistics, Vol 32 pp 1167-1176.
using
n=3
in their notation.Coincidences are much harder to understand with these designs. For example
addelkemp3 3 9
does lead to a number of triple coincidences, that is pairs of rows in which 3 columns agree, but no quadruple coincidences.addelkemp3 9 28
produces an extra column that figures in some quadruple coincidences.As for addelkemp above, 4 is replaced by 1 in fields that do not have an element 4. Also powers of 2 larger than 4 are not allowed, as described above for addelkemp.
The article is quite vague on this. Page 1173 states "When n>2 the same procedure will yield the desired plans if Lemma 5a is used in place of Lemma 5." Page 1175 provides the example n=3,q=3 which is OA( 54,25,3,2 ). Based on this example it is possible to make an educated guess as to how the construction generalizes to n=3. The resulting OA's are seen, by brute force to be of strength 2 for q=2,3,4,5,7,11. These OAs are:
- OA( 16, 13, 2, 2 )
- OA( 54, 25, 3, 2 )
- OA( 128, 41, 4, 2 )
- OA( 250, 61, 5, 2 )
- OA( 686, 113, 7, 2 )
- OA( 1458, 181, 9, 2 )
- OA( 2662, 265, 11, 2 )
The one with q=7 required 212709392 comparisons to determine that it really is of strength 2. This took roughly 11.5 minutes on a DEC 5000/240 workstation (real and elapsed in this case). The array with q=11 took 1.12671e+10 comparisons to verify its strength. This took roughly 10 1/2 hours.
For even q, only q= 2 or 4 are available. The prescription given in Addelman and Kempthorne (1961) does not appear to work. Commented out code below attempts to implement that prescription. It seemed to be impossible to find a constant b[1],c[1] pair.
q | the number of symbols (0,...,q-1) | |
k | the number of columns in the array. k <= 2q^2+2q+1 | |
[out] | n | the number of rows in the array, n = 2q^3 |
std::runtime_error |
void oacpp::COrthogonalArray::addelkempn | ( | int | akn, |
int | q, | ||
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Addelman Kempthorne algorithm
This method is not included by default in Art Owens's project, but it is in the code. It is not compiled in the makefile. Adding it as a target to the makefile creates a successful build.
From the original documentation:
The article is quite vague on this. Page 1173 states "When n>2 the same procedure will yield the desired plans if Lemma 5a is used in place of Lemma 5." Page 1175 provides the example n=3,q=3 which is OA( 54,25,3,2 ). Based on this example it is possible to make an educated guess as to how the construction generalizes.
akn | the exponent on q for the number of rows n = 2q^akn | |
q | the number of symbols (0,...,q-1) | |
k | the number of columns in the array. k <= 2(q^akn-1)/(q-1) - 1 | |
[out] | n | the number of rows in the array, n = 2q^akn |
std::runtime_error |
void oacpp::COrthogonalArray::bose | ( | int | q, |
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Bose algorithm
From the original documentation:
The bose program produces
OA( q^2, k, q, 2 ), k <= q+1
for prime powersq
. This is based on:R.C. Bose (1938) Sankhya Vol 3 pp 323-338
q | the number of symbols (0,...,q-1) | |
k | the number of columns in the array. k <= q+1 | |
[out] | n | the number of rows in the array, n = q^2 |
std::runtime_error |
void oacpp::COrthogonalArray::bosebush | ( | int | q, |
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Bose-Bush algorithm
From the original documentation:
The bosebush program produces
OA( 2q^2, k, q, 2 ), k <= 2q+1
, for powers of 2,q = 2^r
. This construction is based on:R.C. Bose and K.A. Bush (1952) Annals of Mathematical Statistics, Vol 23 pp 508-524.
2q
columns can be constructed without a coincidence defect. Settingk = 2q+1
leads to an array with the coincidence defect. Some triples of columns contain duplicate rows. (The lack of a coincidence defect has been verified forq = 2,4,8,16,32
andk = 2q
.)
q | the number of symbols (0,...,q-1) | |
k | the number of columns in the array. k <= q+1 | |
[out] | n | the number of rows in the array, n = q^2 |
std::runtime_error |
void oacpp::COrthogonalArray::bosebushl | ( | int | lambda, |
int | q, | ||
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Bose-Bush algorithm
From the original documentation:
The bosebushl program produces
OA( lambda*q^2, k, q, 2 ),
, for prime powers
k <= lambda*q+1q
andlambda > 1
. Bothq
andlambda
must be powers of the same prime. This construction is based on:R.C. Bose and K.A. Bush (1952) Annals of Mathematical Statistics, Vol 23 pp 508-524.
Coincidences are harder to understand with these designs. For example
bosebushl 3 9
does lead to a number of triple coincidences, that is pairs of rows in which 3 columns agree, but no quadruple coincidences.bosebush 3 9 28
produces an extra column that figures in some quadruple coincidences.The arrays produced by this program are not always the largest possible. The article by Bose and Bush cited above describes ways of adjoining some extra columns.
When
k <= lambda*q
, the program produces an array that is "completely resolvable". What this means is that the rows of the array may be split into lambda*q consecutive nonoverlapping sets of rows each of which isOA( q,k,q,1 )
.
lambda | |
q | |
k | |
n |
std::runtime_error |
void oacpp::COrthogonalArray::bush | ( | int | q, |
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Bush algorithm
From the original documentation:
The bush program produces
OA( q^3, k, q, 3 ), k <= q+1
for prime powersq
. This strength 3 construction is based on:K.A. Bush (1952) Annals of Mathematical Statistics, Vol 23 pp 426-434
This construction is the most commonly used special case of busht given below.
q | the number of symbols (0,...,q-1) | |
k | the number of columns in the array. k <= q+1 | |
[out] | n | the number of rows in the array, n = q^3 |
std::runtime_error |
void oacpp::COrthogonalArray::busht | ( | int | str, |
int | q, | ||
int | k, | ||
int * | n | ||
) |
Construct an orthogonal array using the Bush algorithm
From the original documentation:
The bush program produces
OA( q^t, k, q, t ), k <= q+1, t>=3
, for prime powersq
. This strengtht
construction is based on:K.A. Bush (1952) Annals of Mathematical Statistics, Vol 23 pp 426-434
str | |
q | |
k | |
n |
std::runtime_error |
|
inline |
Get available warning message
|
inline |
column accessor
|
inline |
row accessor
|
inline |
orthogonal array accessor
|
inline |
symbol accessor
|
inline |
Get the method return code
int oacpp::COrthogonalArray::oaagree | ( | bool | verbose | ) |
Count the number of columns for which each pair of rows agree
From the original documentation:
This program counts the number of columns in which each pair of distinct rows agree.
Input is described above under OA input conventions.
Examples:
COrthogonalArray coa; int n; coa.addelkemp3(3, 25, &n); coa.oagree(true);
This example finds that in OA( 54, 25, 3, 2 ) produced by addelkemp3 there exist pairs of rows agreeing in 9 columns. The first rows to attain this are rows 0 and 9, the 1st and 10th rows.
COrthogonalArray coa; int n; coa.addelkemp3(3, 24, &n); coa.agree(true);
The second example finds that in OA( 54, 24, 3, 2 ) produced by addelkemp3 there exist pairs of rows agreeing in 8 columns. No pairs of rows agree in 9 columns.
verbose | Should messages be printed about the findings? |
|
inline |
Print the dimension of the orthogonal array
void oacpp::COrthogonalArray::oarand | ( | int | is, |
int | js, | ||
int | ks, | ||
int | ls | ||
) |
Randomize an orthogonal array
From the original documentation
This program permutes the symbols in each column. The permutations are uniformly distributed (all
q!
permutations have the same probability) and all columns are permuted independently.Input is described above under OA input conventions, with exceptions noted below to allow passing a random seed. If oarand is called twice with the same input array, the same permuted output will result both times, unless different seeds are given.
The random number generator is a version of the Marsaglia-Zaman random number generator, transliterated into C from FORTRAN. The seed must be four integers between 1 and 168 inclusive, with not all values equal to 1.
is | seed |
js | seed |
ks | seed |
ls | seed |
int oacpp::COrthogonalArray::oastr | ( | bool | verbose | ) |
Find the strength of an orthogonal array
This program reads an orthogonal array strength by brute force computation. In addition to the strength t
described above under the heading orthogonal arrays, strength 0 is taken to mean that the array indeed has all its elements in the range 0..q-1.
An array of strength t > 0
is also of strength s
for all 0 <= s < t
. The program starts testing t = 0
and increases t
until it finds t
for which the array is not strength t
.
Finding the strength of an array by brute force is lightning fast for small arrays but very slow for larger arrays. When the job is large enough, intermediate results are printed so the user can decide whether or not to kill the job, based on how much progress is being made.
The function that calculates strength has an argument verbose. In oastr the array strength function is called with verbose=2. This prints to standard output a description of progress as the strength check proceeds. If one wants to use this function in other settings, calling it with verbose=1 shuts off standard output but leaves the warnings to standard error, and verbose=0 shuts off all output.
verbose | should diagnostic message be printed? |
bool oacpp::COrthogonalArray::oastr1 | ( | bool | verbose | ) |
Similar to oastr, but only checking strength 1
From the original documentation:
Check whether the array in standard input is really of strength 1. Use brute force. For OA( nrow, ncol, q, ? ) it takes work roughly proportional to ncol * nrow * q to decide if ? >= 1. The user is warned if this is likely to be too much work.
The program calls exit(0) if the input array has strength
- It calls exit(1) if the array is not of strength 1, or if the input is invalid, or if it is impossible to allocate enough memory to find out.
Note that an array of strength larger than 1 is a fortiori of strength 1 and will pass this test.
verbose | should diagnostic message be printed? |
bool oacpp::COrthogonalArray::oastr2 | ( | bool | verbose | ) |
Similar to oastr, but only checking strength 2
From the original documentation
Check whether the array in standard input is really of strength 2. Use brute force. For OA( nrow, ncol, q, ? ) it takes work roughly proportional to ncol^2 * nrow * q^2/2 to decide if ? >= 2. The user is warned if this is likely to be too much work.
The program calls exit(0) if the input array has strength
- It calls exit(1) if the array is not of strength 2, or if the input is invalid, or if it is impossible to allocate enough memory to find out.
The program exits at the first sign that the array is not of strength 2. This can save lots of work if the problem shows up early, but it doesn't give a complete list of the array's shortcomings. Such a list could be very large.
Note that an array of strength larger than 2 is a fortiori of strength 2 and will pass this test.
verbose | should diagnostic message be printed? |
bool oacpp::COrthogonalArray::oastr3 | ( | bool | verbose | ) |
Similar to oastr, but only checking strength 3
verbose | should diagnostic message be printed? |
bool oacpp::COrthogonalArray::oastr4 | ( | bool | verbose | ) |
Similar to oastr, but only checking strength 4
verbose | should diagnostic message be printed? |
bool oacpp::COrthogonalArray::oastrt | ( | int | t, |
bool | verbose | ||
) |
Similar to oastr, but only checking for strength t
t | the strength to check for |
verbose | should diagnostic messages be printed? |
int oacpp::COrthogonalArray::oatriple | ( | bool | verbose | ) |
Count the number of columns for which each three rows agree
From the original documentation:
This program reports on triple coincidences. For all triples of distinct columns, it counts the number of distinct pairs of rows in which the triple of columns agree.
Input is described above under OA input conventions.
Examples:
COrthogonalArray coa; int n; coa.bosebush(8, 16, &n); coa.agree(true);
There are 0 distinct triples of columns that agree in at least two distinct rows.
COrthogonalArray coa; int n; coa.bosebush(8, 17, &n); coa.agree(true);
Warning: The Bose-Bush construction with ncol = 2q+1 has a defect. While it is still an OA(2q^2,2q+1,q,2), there exist some pairs of rows that agree in three columns.
There are 8 distinct triples of columns that agree in at least two distinct rows.
The warning above is generated by bosebush. The rest shows that there are triple coincidences. Notice that they all involve the 17th column (which is column 16, since the first one is column 0). The other 16 columns can be organized into 8 pairs with each pair forming a triple with column 16 and no other triples agreeing in any row.
verbose | Should messages be printed about the findings? |