8
8
* it under the terms of the GNU Lesser General Public License as published by
9
9
* the Free Software Foundation, either version 3 of the License, or
10
10
* (at your option) any later version.
11
- *
11
+ *
12
12
* This program is distributed in the hope that it will be useful,
13
13
* but WITHOUT ANY WARRANTY; without even the implied warranty of
14
14
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15
15
* GNU Lesser General Public License for more details.
16
- *
16
+ *
17
17
* You should have received a copy of the GNU Lesser General Public License
18
18
* along with this program. If not, see <http://www.gnu.org/licenses/>.
19
- *
19
+ *
20
20
* Reference:
21
21
* <ul><li><a href="http://lib.stat.cmu.edu/designs/">Statlib Designs</a></li>
22
22
* <li><a href="http://lib.stat.cmu.edu/designs/oa.c">Owen's Orthogonal Array Algorithms</a></li></ul>
36
36
#include < sstream>
37
37
#include < iostream>
38
38
#include < numeric>
39
+ #include < array>
39
40
40
41
#ifdef RCOMPILE
41
42
#include < Rcpp.h>
@@ -99,153 +100,153 @@ namespace oacpp {
99
100
100
101
/* *
101
102
* @page oa_main_page Orthogonal Array Library
102
- *
103
+ *
103
104
* From the original documentation by Owen:
104
- *
105
+ *
105
106
*<blockquote>
106
107
* From: owen@stat.stanford.edu
107
- *
108
- * These programs construct and manipulate orthogonal
108
+ *
109
+ * These programs construct and manipulate orthogonal
109
110
* arrays. They were prepared by
110
111
* - Art Owen
111
112
* - Department of Statistics
112
113
* - Sequoia Hall
113
114
* - Stanford CA 94305
114
- *
115
+ *
115
116
* They may be freely used and shared. This code comes
116
117
* with no warranty of any kind. Use it at your own
117
118
* risk.
118
- *
119
+ *
119
120
* I thank the Semiconductor Research Corporation and
120
121
* the National Science Foundation for supporting this
121
122
* work.
122
- *
123
+ *
123
124
* I thank Randall Tobias of SAS Inc. for many helpful
124
125
* electronic discussions that lead to improvements in
125
126
* these programs.
126
127
* </blockquote>
127
- *
128
+ *
128
129
* @tableofcontents
129
- *
130
+ *
130
131
* @section orthogonal_arrays_sec Orthogonal Arrays
131
132
* <blockquote>
132
- * An orthogonal array <code>A</code> is a matrix of <code>n</code> rows, <code>k</code>
133
+ * An orthogonal array <code>A</code> is a matrix of <code>n</code> rows, <code>k</code>
133
134
* columns with every element being one of <code>q</code> symbols
134
135
* <code>0,...,q-1</code>. The array has strength <code>t</code> if, in every <code>n</code> by <code>t</code>
135
136
* submatrix, the <code>q^t</code> possible distinct rows, all appear
136
137
* the same number of times. This number is the index
137
- * of the array, commonly denoted <code>lambda</code>. Clearly,
138
- * <code>lambda*q^t = n</code>. Geometrically, if one were to "plot" the
138
+ * of the array, commonly denoted <code>lambda</code>. Clearly,
139
+ * <code>lambda*q^t = n</code>. Geometrically, if one were to "plot" the
139
140
* submatrix with one plotting axis for each of the <code>t</code> columns
140
141
* and one point in <code>t</code> dimensional space for each row, the
141
142
* result would be a grid of <code>q^t</code> distinct points. There would
142
143
* be <code>lambda</code> "overstrikes" at each point of the grid.
143
- *
144
+ *
144
145
* The notation for such an array is <code>OA( n, k, q, t )</code>.
145
- *
146
+ *
146
147
* If <code>n <= q^(t+1)</code>, then the <code>n</code> rows "should" plot as
147
148
* <code>n</code> distinct points in every <code>n</code> by <code>t+1</code> dimensional subarray.
148
149
* When this fails to hold, the array has the "coincidence
149
150
* defect".
150
- *
151
+ *
151
152
* Owen (1992,199?) describes some uses for randomized
152
153
* orthogonal arrays, in numerical integration, computer
153
154
* experiments and visualization of functions. Those
154
155
* references contain further references to the literature,
155
- * that provide further explanations. A strength 1 randomized
156
+ * that provide further explanations. A strength 1 randomized
156
157
* orthogonal array is a Latin hypercube
157
158
* sample, essentially so or exactly so, depending on
158
159
* the definition used for Latin hypercube sampling.
159
160
* The arrays constructed here have strength 2 or more, it
160
161
* being much easier to construct arrays of strength 1.
161
- *
162
+ *
162
163
* The randomization is achieved by independent
163
164
* uniform permutation of the symbols in each column.
164
- *
165
+ *
165
166
* To investigate a function <code>f</code> of <code>d</code> variables, one
166
167
* has to have an array with <code>k >= d</code>. One may also
167
168
* have a maximum value of <code>n</code> in mind and a minimum value
168
169
* for the number <code>q</code> of distinct levels to investigate.
169
- *
170
+ *
170
171
* It is entirely possible that no array of strength <code>t > 1</code>
171
172
* is compatible with these conditions. The programs
172
173
* below provide some choices to pick from, hopefully
173
174
* without too much of a compromise.
174
- *
175
+ *
175
176
* The constructions used are based on published
176
177
* algorithms that exploit properties of Galois fields.
177
178
* Because of this the number of levels <code>q</code> must be
178
179
* a prime power. That is <code>q = p^r</code> where <code>p</code> is prime
179
180
* and <code>r >= 1</code> is an integer.
180
- *
181
+ *
181
182
* The Galois field arithmetic for the prime powers is
182
183
* based on tables published by Knuth and Alanen (1964)
183
184
* below. The resulting fields have been tested by the
184
185
* methods described in Appendix 2 of that paper and
185
186
* they passed. This is more a test of the accuracy of
186
187
* my transcription than of the original tables.
187
188
* </blockquote>
188
- *
189
+ *
189
190
* @section avail_prime_sec Available Prime Powers
190
- *
191
+ *
191
192
* <blockquote>
192
193
* The designs given here require a prime power for
193
- * the number of levels. They presently work for the
194
+ * the number of levels. They presently work for the
194
195
* following prime powers:
195
- *
196
+ *
196
197
* All Primes
197
198
* All prime powers <code>q = p^r</code> where <code>p < 50</code> and <code>q < 10^9</code>
198
- *
199
+ *
199
200
* Here are some of the smaller prime powers:
200
- *
201
+ *
201
202
* - Powers of 2: 4 8 16 32 64 128 256 512
202
203
* - Powers of 3: 9 27 81 243 729
203
204
* - Powers of 5: 25 125 625
204
- * - Powers of 7: 49 343
205
+ * - Powers of 7: 49 343
205
206
* - Square of 11: 121
206
207
* - Square of 13: 169
207
- *
208
+ *
208
209
* Here are some useful primes:
209
- *
210
+ *
210
211
* - 2,3,5,7,11,13,17,19,23,29,31,37,101,251,401
211
- *
212
+ *
212
213
* The first row are small primes, the second row are
213
214
* primes that are 1 more than a "round number". The small
214
215
* primes lead to small arrays. An array with 101 levels
215
216
* is useful for exploring a function at levels 0.00 0.01
216
217
* through 1.00. Keep in mind that a strength 2 array on
217
218
* 101 levels requires <code>101^2 = 10201</code> experimental runs,
218
219
* so it is only useful where large experiments are possible.
219
- *
220
+ *
220
221
* Note that some of these will require more
221
222
* memory than your computer has. For example,
222
223
* with a large prime like 10663, the program knows
223
224
* the Galois field, but can't allocate enough
224
225
* memory:
225
- *
226
+ *
226
227
* <code>bose 10663</code>
227
228
* - Unable to allocate 1927'th row in an integer matrix.
228
229
* - Unable to allocate space for Galois field on 10663 elements.
229
230
* - Construction failed for <code>GF(10663)</code>.
230
231
* - Could not construct Galois field needed for Bose design.
231
- *
232
+ *
232
233
* The smallest prime power not covered is <code>53^2 = 2809</code>.
233
234
* The smallest strength 2 array with 2809 symbols has
234
235
* <code>2809^2 = 7890481</code> rows. Therefore the missing prime powers
235
236
* are only needed in certain enormous arrays, not in the
236
237
* small ones of most practical use. In any event there
237
238
* are some large primes and prime powers in the program
238
239
* if an enormous array is needed.
239
- *
240
+ *
240
241
* To add <code>GF(p^r)</code> for some new prime power p^r,
241
242
* consult Alanen and Knuth for instructions on how
242
243
* to search for an appropriate indexing polynomial,
243
244
* and for how to translate that polynomial into a
244
- * replacement rule for <code>x^r</code>.
245
+ * replacement rule for <code>x^r</code>.
245
246
* </blockquote>
246
- *
247
+ *
247
248
* @section methods Methods
248
- *
249
+ *
249
250
* <ul>
250
251
* <li>@ref oacpp::COrthogonalArray::bose</li>
251
252
* <li>@ref oacpp::COrthogonalArray::bush</li>
@@ -264,9 +265,9 @@ namespace oacpp {
264
265
* <li>@ref oacpp::COrthogonalArray::oaagree</li>
265
266
* <li>@ref oacpp::COrthogonalArray::oadimen</li>
266
267
* </ul>
267
- *
268
+ *
268
269
* @section tips Tips On Use
269
- *
270
+ *
270
271
* <blockquote>
271
272
* It is faster to generate only the columns you need.
272
273
* For example
@@ -275,11 +276,11 @@ namespace oacpp {
275
276
* <code>bose 101</code>
276
277
* generates 102 columns. If you only want 4 columns the
277
278
* former saves a lot of time.
278
- *
279
+ *
279
280
* Passing the <code>q n k</code> on the command line is more difficult
280
281
* than letting the computer figure them out, but it
281
282
* allows more error checking.
282
- *
283
+ *
283
284
* In practical use, I would try first to use a Bose
284
285
* design. Then I would consider either an Addelman-
285
286
* Kempthorne or Bose-Bush design to see whether it
@@ -289,9 +290,9 @@ namespace oacpp {
289
290
* large number of runs is possible a Bush design may
290
291
* work well, since it can have high strength.
291
292
* </blockquote>
292
- *
293
+ *
293
294
* @section references References
294
- *
295
+ *
295
296
* <blockquote>
296
297
* Here are the references for the constructions used:
297
298
* <ul>
@@ -303,7 +304,7 @@ namespace oacpp {
303
304
* </ul>
304
305
* This book provides a large list of orthogonal array constructions:
305
306
* <ul><li>Aloke Dey (1985) "Orthogonal Fractional Factorial Designs" Halstead Press</li></ul>
306
- *
307
+ *
307
308
* These papers discuss randomized orthogonal arrays, the second
308
309
* is being revised in parallel with development of the software
309
310
* described here:
@@ -320,9 +321,9 @@ namespace oacpp {
320
321
* <li>M. Stein (1987) Technometrics 29, 143-151</li>
321
322
* </ul>
322
323
* </blockquote>
323
- *
324
+ *
324
325
* @section implement Implementation Details
325
- *
326
+ *
326
327
* <blockquote>
327
328
* Galois fields are implemented through arrays that
328
329
* store their addition and multiplication tables. Some
@@ -333,15 +334,15 @@ namespace oacpp {
333
334
* was not considered to be a burden. Subtraction and
334
335
* division are implemented through vectors of additive
335
336
* and multiplicative inverses, derived from the tables.
336
- * The tables for <code>GF(p^r)</code> are constructed using a
337
+ * The tables for <code>GF(p^r)</code> are constructed using a
337
338
* representation of the field elements as polynomials in <code>x</code>
338
339
* with coefficients as integers modulo <code>p</code> and a special
339
340
* rule (derived from minimal polynomials) for handling
340
341
* products involving <code>x^r</code>. These rules are taken from
341
342
* published references. The rules have not all
342
- * been checked for accuracy, because some of the fields are
343
+ * been checked for accuracy, because some of the fields are
343
344
* very large (e.g. 16807 elements).
344
- *
345
+ *
345
346
* The functions that manipulate orthogonal arrays
346
347
* keep the arrays in integer matrices. This might be
347
348
* a problem for applications that require enormous
@@ -355,11 +356,11 @@ namespace oacpp {
355
356
* row by row, so it is not too hard to change the program
356
357
* to place the elements on an output stream as they
357
358
* are computed and do away with the storage.
358
- *
359
+ *
359
360
* The functions that test the strength of the
360
361
* arrays may be very far from optimally fast.
361
362
* </blockquote>
362
- *
363
+ *
363
364
* @section compile_oa Compiling oalib
364
365
* When compiling <code>oalib</code> these preprocessor directives are used:
365
366
* - NDEBUG defined for a release build
0 commit comments