Welcome to oapackage’s documentation!

The Orthogonal Array package contains functionality to generate and analyse orthogonal arrays, optimal designs and conference designs.

Introduction

Orthogonal arrays, optimal designs and conference designs are important tools for the design of experiments [EN95] [HSS99] [WH09]. The Orthogonal Array package contains functionality to generate and analyse these types of designs. To generate the arrays and designs, the package uses the exhaustive enumeration algorithm of [SEN10] and the optimization algorithm of [ES17]. To analyze the arrays and designs, the package includes a wide variety of relevant statistical and combinatorial criteria. A large collection of orthogonal arrays, optimal designs and conference designs generated with the package are available in the Orthogonal Array package website [Een18].

Example usage

The Orthogonal Array package can be used to generate and manipulate arrays and designs. Additionally, it can calculate some of their statistical properties. The following example shows how to generate an orthogonal array with 8 runs and 2 factors, and calculate three relevant statistical properties:

Calculate D-efficiency

>>> import oapackage
>>> array=oapackage.exampleArray(0) # define an orthogonal array
>>> array.showarray()
array:
  0   0
  0   0
  0   1
  0   1
  1   0
  1   0
  1   1
  1   1
>>> D = array.Defficiency() # calculate the D-efficiency for estimating the interaction effects model
>>> array_rank = array.rank() # calculate the rank of the design
>>> print('D-efficiency %f, rank %d' % (D, array_rank) )
D-efficiency 1.000000, rank 2
>>> gwlp = array.GWLP() # calculate the generalized word length pattern
>>> print('Generalized wordlength pattern: %s' % (gwlp,) )
Generalized wordlength pattern: (1.0, 0.0, 0.0)

The statistical properties of the arrays and designs are introduced in Properties of designs.

Interfaces

The Orthogonal Array package has interfaces in C++ and Python for generating, manipulating and analyzing all the types of arrays and designs. In this documentation, you will find references to both the Python and the C++ interface. The package also includes several command line tools.

For the generation of optimal designs [ES17], the Orthogonal Array package has also a Matlab interface; see the documentation README.Matlab.md.

License

The code is available under a BSD style license; see the file LICENSE for details. If you use this code or any of the results, please cite this program as follows:

OApackage: A Python package for generation and analysis of orthogonal arrays, optimal designs and conference designs, P.T. Eendebak, A.R. Vazquez, Journal of Open Source Software, 2019
Complete Enumeration of Pure-Level and Mixed-Level Orthogonal Arrays, E.D. Schoen, P.T. Eendebak, M.V.M. Nguyen, Volume 18, Issue 2, pages 123-140, 2010.

Acknowledgements

The code and ideas for this package have been contributed by Eric Schoen, Ruben Snepvangers, Vincent Brouerius van Nidek, Alan Roberto Vazquez and Pieter Thijs Eendebak.

Installation

The packge is continously tested on Linux and Windows. The Python interface is available from the Python Package Index. The package can be installed from the command line using pip:

$ pip install OApackage

The source code for the package is available on https://github.com/eendebakpt/oapackage. The command line tools use a cmake build system. From the command line, type the following:

$ mkdir -p build
$ cd build
$ cmake ..
$ make
$ make install

This creates the command line utilities and a C++ library.

To compile the Python interface use

$ python setup.py build
$ python setup.py install --user

The Python interface requires Numpy [TheScipycommunity12], Matplotlib [Hun07] and Swig. The code has been tested with Python 3.6, 3.7 and 3.8.

The R interface to the optimal design functionality of the package is available from CRAN. For the Matlab and Octave interface to the optimal design functionality see the file README.Matlab.md.

Example notebooks

This section contains several examples for generating, manipulating and analyzing arrays and designs using the Orthogonal Arrays package. The examples are shown using Jupyter notebooks, which can be found on github https://github.com/eendebakpt/oapackage/tree/master/docs/examples.

Example script for Python interface to Orthogonal Array package

Pieter Eendebak pieter.eendebak@gmail.com

[1]:

import numpy as np
import oapackage

print("oapackage version: %s" % oapackage.version())

oapackage version: 2.6.0

Load an example array.

[2]:

array = oapackage.exampleArray(0)
array.showarray()

Calculate properties of the array such as the D-efficiency for the main effects model, the generalized word length pattern and the rank.

[3]:

print("D-efficiency %f, rank %d" % (array.Defficiency(), array.rank()))
print("Generalized wordlength pattern: %s" % str(array.GWLP()))

D-efficiency 1.000000, rank 2
Generalized wordlength pattern: (1.0, 0.0, 0.0)

Calculate the generalized word length pattern for another example array.

[11]:

array = oapackage.exampleArray(11)
print("Generalized wordlength pattern: %s" % oapackage.oahelper.gwlp2str(array.GWLP()))

Generalized wordlength pattern: 1.00,0.02273,0.03926,0.5434,0.8244,2.217,1.043,0.126,0.002066

Indexing

The array_link object can be indexed as a normal array.

[12]:

array[0:5, 2:3]

[12]:

array_link: 5, 1

[14]:

array[0:5, 2:4].showarray()

[15]:

print(array[0, 2])

Numpy

We can convert between Numpy arrays and array_link objects. Note that an array_link is always integer valued.

[19]:

X = (4 * np.random.rand(20, 10)).astype(int)
array = oapackage.array_link(X)
array.showarraycompact()

Convert from array_link back to Numpy array.

[20]:

X2 = np.array(array)

Example arrays

The package contains several example arrays and designs. They can be retrieved using the method exampleArray.

[4]:

import oapackage
al = oapackage.exampleArray(-1, 1)

exampleArray 0: array in OA(8,2, 2^2)
exampleArray 1: array 3 in OA(16, 2, 2^5)
exampleArray 2: array 6 in OA(16, 2, 2^6)
exampleArray 3: array ? in OA(32, 3, 2^7)
exampleArray 4: array 4 in OA(16, 2, 2^7)
exampleArray 5: array 0 in OA(24, 2, 4 3 2^a)
exampleArray 6: array in OA(4, 2, 2^a)
exampleArray 7: array 0 in OA(4, 2, 2^a)?
exampleArray 8: array in OA(40, 3, 2^7)
exampleArray 9: array in A(40, 2^7), D-optimal
exampleArray 10: array in OA(9, 3^3)
exampleArray 11: D-optimal array in OA(44, 2^8)
exampleArray 12: even-odd array OA(64, 2^13)
exampleArray 13: array in OA(25, 2^5)
exampleArray 14: design in D(28, 2^5), D-efficiency is low
exampleArray 15: design in D(56, 2^10), D-efficiency is low
exampleArray 16: array in OA(32, 2, 2^5)
exampleArray 17: unique array in OA(64, 4, 2^7)
exampleArray 18: conference matrix of size 16, 7
exampleArray 19: conference matrix of size 4, 3
exampleArray 20: first LMC-0 double conference matrix in DC(24,3)
exampleArray 21: second LMC-0 double conference matrix in DC(16,4)
exampleArray 22: LMC-0 double conference matrix in DC(32,4)
exampleArray 23: LMC-0 double conference matrix in DC(32,6)
exampleArray 24: design in OA(64, 3, 2^16) (even-odd)
exampleArray 25: design in OA(64, 3, 2^16) (even-odd)
exampleArray 26: design in OA(64, 3, 2^16) (even-odd)
exampleArray 27: design in OA(64, 3, 2^16) (even-odd)
exampleArray 28: conference design in C(4, 3) in LMC0 form
exampleArray 29: conference design in C(4, 3)
exampleArray 30: conference design in C(8,4)
exampleArray 31: conference design in C(8,4)
exampleArray 32: first double conference design in DC(18,4)
exampleArray 33: second double conference design in DC(18,4)
exampleArray 34: third double conference design in DC(18,4)
exampleArray 35: first double conference design in DC(20,4)
exampleArray 36: second double conference design in DC(20,4)
exampleArray 37: third double conference design in DC(20,4)
exampleArray 38: LMC0 conference design in C(30,3)
exampleArray 39: first LMC0 conference design in C(8,6)
exampleArray 40: first conference design in C(14, 5)
exampleArray 41: second conference design in C(14, 5)
exampleArray 42: third conference design in C(14, 5)
exampleArray 43: 2x2 array with zeros and a singe value -1
exampleArray 44: D-optimal strength 3 ortogonal array in OA(40,3, 2^7)
exampleArray 45: first conference design in C(20,8)
exampleArray 46: second conference design in C(20,8)
exampleArray 47: third conference design in C(20,8)
exampleArray 48: last conference design in C(20,8)
exampleArray 49: array 4347 C(20,8)
exampleArray 50: array 4506 C(20,8)
exampleArray 51: first array in C(12,4)
exampleArray 52: second array in C(12,4)
exampleArray 53: third array in C(12,4)
exampleArray 54: root array in OA(6,2,3^1 2^1)
exampleArray: no example array with index 55 exists

Select an example array and show it.

[3]:

al=oapackage.exampleArray(2, 1)
al.showarray()

exampleArray 2: array 6 in OA(16, 2, 2^6)
array:
  0   0   0   0   0   0
  0   0   0   0   0   0
  0   0   0   1   1   1
  0   0   0   1   1   1
  0   1   1   0   0   1
  0   1   1   0   1   0
  0   1   1   1   0   1
  0   1   1   1   1   0
  1   0   1   0   0   1
  1   0   1   0   1   1
  1   0   1   1   0   0
  1   0   1   1   1   0
  1   1   0   0   1   0
  1   1   0   0   1   1
  1   1   0   1   0   0
  1   1   0   1   0   1

Example to write and read files with arrays from disk

The package can write arrays and designs to a custom file format. There is a text based format (specified by oapackage.ATEXT) and a binary format (specified by oapackage.ABINARY). See the online documentation for details on the file formats.

[1]:

import oapackage
import tempfile

Create a list of two example arrays.

[2]:

lst = [oapackage.exampleArray(32), oapackage.exampleArray(33)]
print(lst)

[array_link: 18, 4, array_link: 18, 4]

Write the two arrays to a file on disk.

[3]:

filename = tempfile.mktemp(".oa")
_ = oapackage.writearrayfile(filename, lst, oapackage.ATEXT)

Display information about the file written.

[4]:

oapackage.oainfo(filename)

file C:\Users\EENDEB~1\AppData\Local\Temp\tmpxc4dnsl6.oa: 18 rows, 4 columns, 2 arrays, mode text, nbits 0

The arrays can be read from disk again using the command readarrayfile.

[6]:

lst2 = oapackage.readarrayfile(filename)
print(lst2)

(array_link: 18, 4, array_link: 18, 4)

The first array in the list that was read from disk is:

[7]:

lst2[0].showarray()

array:
  0   0   0   0
  0   0   0   0
  1   1   1   1
  1   1   1   1
  1   1  -1  -1
  1   1  -1  -1
  1  -1   1  -1
  1  -1   1  -1
  1  -1  -1   1
  1  -1  -1   1
 -1   1   1  -1
 -1   1   1  -1
 -1   1  -1   1
 -1   1  -1   1
 -1  -1   1   1
 -1  -1   1   1
 -1  -1  -1  -1
 -1  -1  -1  -1

Enumerate orthogonal arrays

The Orthogonal Array package can completely enumerate all orthogonal arrays of a specified class. In this notebook, we enumerate specific classes of three-level orthogonal arrays and mixel-level orthogonal arrays.

First, we specify the class of three-level orthogonal arrays to enumerate. For example, we consider three-level orthogonal arrays of strength 2 with 27 runs and 8 factors.

[1]:

import oapackage
run_size = 27
strength=2
number_of_factors=8
factor_levels = 3
arrayclass=oapackage.arraydata_t(factor_levels, run_size, strength, number_of_factors)
print(arrayclass)

arrayclass: N 27, k 8, strength 2, s {3,3,3,3,3,3,3,3}, order 0

Second, we create the root array as the starting point of our enumeration.

[2]:

ll2=[arrayclass.create_root()]
ll2[0].showarraycompact()

Third, extend the root array. It is also possible to extend a list of arrays.

[3]:

list3columns = oapackage.extend_arraylist(ll2, arrayclass)
print('extended to %d arrays with 3 columns' % len(list3columns))
list4columns = oapackage.extend_arraylist(list3columns, arrayclass)
print('extended to %d arrays with 4 columns' % len(list4columns))

extended to 9 arrays with 3 columns
extended to 711 arrays with 4 columns

It is also possible to extend selected arrays from a list.

[4]:

ll = oapackage.extend_arraylist(list4columns[0:8], arrayclass)
print('extended first 2 arrays to %d arrays with 5 columns' % len(ll))

extended first 2 arrays to 189 arrays with 5 columns

By adding one column at a time we can enumerate all three-level orthogonal arrays of strength 2 with 27 runs and 8 factors. The total computation time for this would be a couple of hours.

Mixed-level orthogonal arrays

The package can also enumerate mixed-level orthogonal arrays. For instance, consider enumerating all 16-run strength-2 orthogonal arrays with one four-level factor and nine two-level factors.

[9]:

run_size = 16
strength=2
number_of_factors = 10
factor_levels=[4,2,2,2,2,2,2,2,2,2]
arrayclass=oapackage.arraydata_t(factor_levels, run_size, strength, number_of_factors)
print(arrayclass)

arrayclass: N 16, k 10, strength 2, s {4,2,2,2,2,2,2,2,2,2}, order 0

Create the root array as the starting point of our enumeration.

[10]:

al=arrayclass.create_root()
al.showarraycompact()

For these arrays, we can extend a single array or lists of arrays.

[14]:

array_list=[arrayclass.create_root()]
array_list_3columns=oapackage.extend_arraylist(array_list, arrayclass)
array_list_4columns=oapackage.extend_arraylist(array_list_3columns, arrayclass)
print('extended to %d arrays with 3 columns' % len(array_list_3columns))
print('extended to %d arrays with 4 columns' % len(array_list_4columns))

extended to 3 arrays with 3 columns
extended to 10 arrays with 4 columns

Finally, enumerate all 16-run strength-2 orthogonal arrays with one four-level factor and nine two-level factors.

[15]:

arrays = array_list_4columns
for extension_column in range(5, number_of_factors+1):
    extensions=oapackage.extend_arraylist(arrays, arrayclass)
    print('extended to %d arrays with %d columns' % (len(extensions), extension_column))
    arrays=extensions

extended to 28 arrays with 5 columns
extended to 65 arrays with 6 columns
extended to 110 arrays with 7 columns
extended to 123 arrays with 8 columns
extended to 110 arrays with 9 columns
extended to 72 arrays with 10 columns

Notes

The numbers of isomorphism classes for various types of classes can be found at the webpage series of orthogonal arrays.
For larger number of arrays the command line tools are more convenient and more memory efficient.

Generate orthogonal arrays with high D-efficiency

This notebook contains example code from the article Two-level designs to estimate all main effects and two-factor interactions by Eendebak, P. T. and Schoen, E. D. This example shows how to generate orthogonal arrays with a high \(D\)-efficiency in a reasonable amount of time (< 1 minute). For more results and details, see the paper.

Generate a D-optimal orthogonal array of strength 2 with 32 runs and 7 factors.

[12]:

import numpy as np
import matplotlib.pyplot as plt
import oapackage
%matplotlib inline

[2]:

run_size = 32
number_of_factors=7
factor_levels=2
strength=2
nkeep=24 # Number of designs to keep at each stage

arrayclass=oapackage.arraydata_t(factor_levels, run_size, strength, number_of_factors)
print('In this example we generate orthogonal arrays in the class: %s' % arrayclass)

In this example we generate orthogonal arrays in the class: arrayclass: N 32, k 7, strength 2, s {2,2,2,2,2,2,2}, order 0

First, generate orthogonal arrays with the function extend_arraylist. Next, keep the arrays with the best \(D\)-efficiency.

[7]:

arraylist=[arrayclass.create_root()]

#%% Extend arrays and filter based on D-efficiency
options=oapackage.OAextend()
options.setAlgorithmAuto(arrayclass)

for extension_column in range(strength+1, number_of_factors+1):
    print('extend %d arrays with %d columns with a single column' % (len(arraylist), arraylist[0].n_columns) )
    arraylist_extensions = oapackage.extend_arraylist(arraylist, arrayclass, options)

    # Select the best arrays based on the D-efficiency
    dd = np.array([a.Defficiency() for a in arraylist_extensions])
    ind = np.argsort(dd)[::-1][0:nkeep]
    selection = [ arraylist_extensions[ii] for ii in ind]
    dd=dd[ind]
    print('  generated %d arrays, selected %d arrays with D-efficiency %.4f to %.4f' % (len(arraylist_extensions), len(ind), dd.min(), dd.max() ) )

    arraylist = selection

extend 1 arrays with 2 columns with a single column
  generated 5 arrays, selected 5 arrays with D-efficiency 0.0000 to 1.0000
extend 5 arrays with 3 columns with a single column
  generated 19 arrays, selected 19 arrays with D-efficiency 0.0000 to 1.0000
extend 19 arrays with 4 columns with a single column
  generated 491 arrays, selected 24 arrays with D-efficiency 0.9183 to 1.0000
extend 24 arrays with 5 columns with a single column
  generated 2475 arrays, selected 24 arrays with D-efficiency 0.9196 to 1.0000
extend 24 arrays with 6 columns with a single column
  generated 94 arrays, selected 24 arrays with D-efficiency 0.7844 to 0.8360

Show the best array from the list of D-optimal orthogonal arrays.

[20]:

selected_array = selection[0]
print('Generated a design in OA(%d, %d, 2^%d) with D-efficiency %.4f' % (selected_array.n_rows, arrayclass.strength, selected_array.n_columns, dd[0] ) )
print('The array is (in transposed form):\n')
selected_array.transposed().showarraycompact()

Generated a design in OA(32, 2, 2^7) with D-efficiency 0.8360
The array is (in transposed form):

00000000000000001111111111111111
00000000111111110000000011111111
00000111000111110001111100000111
00011001011001110110011100011001
00101010101010111010101100101010
01001011010110011101000101101001
01110010001101010111001010100011

We calculate the \(D\)-, \(D_s\)- and \(D_1\)-efficiencies.

[21]:

efficiencies = np.array([array.Defficiencies() for array in arraylist])
print(efficiencies)

[[0.8360354  0.72763273 1.        ]
 [0.83481532 0.78995185 1.        ]
 [0.82829522 0.70849985 1.        ]
 [0.82829522 0.65037401 1.        ]
 [0.82545286 0.75411444 1.        ]
 [0.8185915  0.65157148 1.        ]
 [0.81676643 0.68965719 1.        ]
 [0.81559279 0.72218391 1.        ]
 [0.81557799 0.66966068 1.        ]
 [0.81302642 0.74299714 1.        ]
 [0.80973102 0.73354611 1.        ]
 [0.80286424 0.67455544 1.        ]
 [0.79907785 0.72697978 1.        ]
 [0.79762174 0.58515515 1.        ]
 [0.79514152 0.65199837 1.        ]
 [0.79348419 0.74167976 1.        ]
 [0.7928208  0.70038649 1.        ]
 [0.79225237 0.71961467 1.        ]
 [0.79171981 0.62994172 1.        ]
 [0.79160277 0.66978711 1.        ]
 [0.7883109  0.62733236 1.        ]
 [0.78683679 0.60799604 1.        ]
 [0.78477264 0.69310129 1.        ]
 [0.7843891  0.7533355  1.        ]]

Visualize the \(D\)-efficiencies using a scatter plot.

[19]:

plt.plot(efficiencies[:,0], efficiencies[:,1], '.b')
plt.title('Scatterplot of effiencies of designs with %d columns' % arraylist[0].n_columns)
plt.xlabel('D-efficiency')
_=plt.ylabel('Ds-efficiency')

_images/examples_example_optimal_orthogonal_arrays_11_0.png

Generate D-efficient designs

This notebook contains example code from the article Two-level designs to estimate all main effects and two-factor interactions by Eendebak, P. T. and Schoen, E. D. This example shows how to generate D-efficient designs with a user-specified optimization function.

[1]:

import numpy as np
import oapackage
import oapackage.Doptim

%matplotlib inline

Define the class of designs to generate.

[2]:

run_size = 40
number_of_factors = 7
factor_levels = 2
strength = 0
arrayclass = oapackage.arraydata_t(factor_levels, run_size, strength, number_of_factors)
print("We generate D-efficient designs with %d rows and %d columns\n" % (run_size, number_of_factors))

We generate D-efficient designs with 40 rows and 7 columns

Generate a single D-efficient design using \(\alpha=(1,2,0)\) as the parameters for the optimization function. For details on this parameter and its corresponding optimization function, see Two-Level Designs to Estimate All Main Effects and Two-Factor Interactions.

[3]:

alpha = [1, 2, 0]
scores, design_efficiencies, designs, ngenerated = oapackage.Doptim.Doptimize(
    arrayclass, nrestarts=30, optimfunc=alpha, selectpareto=True
)

Doptim: optimization class 40.2-2-2-2-2-2-2
Doptimize: iteration 0/30
Doptimize: iteration 29/30
Doptim: done (11 arrays, 0.4 [s])

[4]:

print("\nGenerated %d designs, the efficiencies for these designs are:" % len(designs))
for ii, d in enumerate(designs):
    dd = d.Defficiencies()
    print("array %d: D-efficiency %.4f, Ds-efficiency %.4f" % (ii, dd[0], dd[1]))

D = [d.Defficiency() for d in designs]
best = np.argmax(D)
print("\nThe design with the highest D-efficiency (%.4f) is:\n" % D[best])

designs[best].transposed().showarraycompact()


Generated 11 designs, the efficiencies for these designs are:
array 0: D-efficiency 0.8815, Ds-efficiency 0.9807
array 1: D-efficiency 0.9076, Ds-efficiency 0.9628
array 2: D-efficiency 0.8670, Ds-efficiency 0.9827
array 3: D-efficiency 0.8945, Ds-efficiency 0.9669
array 4: D-efficiency 0.9027, Ds-efficiency 0.9540
array 5: D-efficiency 0.8851, Ds-efficiency 0.9549
array 6: D-efficiency 0.8737, Ds-efficiency 0.9581
array 7: D-efficiency 0.9036, Ds-efficiency 0.9400
array 8: D-efficiency 0.8614, Ds-efficiency 0.9595
array 9: D-efficiency 0.8897, Ds-efficiency 0.9418
array 10: D-efficiency 0.9046, Ds-efficiency 0.9203

The design with the highest D-efficiency (0.9076) is:

0010011001011001010101011101101100011101
0011010111100101111100010111001100100010
1011011111001010100010111001000010011010
0101100001100111110011010001001100111110
0011100010001101001011101101110100010001
0011110000011011101001110011101010101000
1000110011101000101110001001101111001101

Optimizing with a different optimization target leads to different D-efficient designs. Below we compare the sets of designs generated with optimization target [1,0,0] and [1,2,0].

[5]:

scores0, design_efficiencies0, designs0, _ = oapackage.Doptim.Doptimize(
    arrayclass, nrestarts=30, optimfunc=[1, 0, 0], selectpareto=True
)

Doptim: optimization class 40.2-2-2-2-2-2-2
Doptimize: iteration 0/30
Doptimize: iteration 29/30
Doptim: done (13 arrays, 0.4 [s])

[7]:

def combineEfficiencyData(lst):
    data = np.zeros((0, 4))

    for jj, dds in enumerate(lst):
        dds_index = np.hstack((dds, jj * np.ones((len(dds), 1))))
        data = np.vstack((data, dds_index))
    return data


design_efficiencies_combined = combineEfficiencyData([design_efficiencies0, design_efficiencies])
plot_handles = oapackage.generateDscatter(
    design_efficiencies_combined, ndata=3, lbls=["Optimize $D$", "Optimize $D+2D_s$"], verbose=0
)

Pareto: 23 optimal values, 24 objects

_images/examples_example_generation_optimal_designs_10_1.png

<Figure size 432x288 with 0 Axes>

[ ]:

Example script to use Nauty from Python

Nauty can be used to reduce any graph into a normal form. In this notebook, we show how to use the Nauty functionality from Python.

[1]:

import numpy as np
import oapackage

Define a function to invert a permutation.

[2]:

def inverse_permutation(perm):
    inverse = [0] * len(perm)
    for i, p in enumerate(perm):
        inverse[p] = i
    return inverse

We define a graph with 5 nodes. The graph is defined by the incidence matrix of size \(5\times 5\) and a coloring with two colors.

[3]:

graph = np.zeros((5, 5), dtype=int)
graph[0, 1] = graph[0, 2] = graph[0, 3] = graph[1, 3] = 1
graph = np.maximum(graph, graph.T)  # make array symmetric
colors = [0, 0, 0, 1, 1]

Reduce the graph to normal form using Nauty.

[4]:

help(oapackage.reduceGraphNauty)

Help on function reduceGraphNauty in module oalib:

reduceGraphNauty(G, colors=None, verbose=1)
    Return vertex transformation reducing array to normal form

    The reduction is calculated using `Nauty <http://users.cecs.anu.edu.au/~bdm/nauty/>`_

    Args:
        G (numpy array or array_link) :   the graph in incidence matrix form
        colors (list or None): an optional vertex coloring
    Returns:
        v: relabelling of the vertices

[4]:

def reduce(graph, colors):
    tr = oapackage.reduceGraphNauty(graph, colors=colors, verbose=0)
    tri = inverse_permutation(tr)

    graph_reduced = oapackage.transformGraphMatrix(graph, tri)
    colors_reduced = [colors[idx] for idx in tr]
    return graph_reduced, colors_reduced, tr


graph_reduced, colors_reduced, tr = reduce(graph, colors)

print("input graph: ")
print(graph)

print("normal form reduction: %s" % (tr,))
print("reduced graph: ")
print(graph_reduced)
print("colors reduced: %s" % (colors_reduced,))

input graph:
[[0 1 1 1 0]
 [1 0 0 1 0]
 [1 0 0 0 0]
 [1 1 0 0 0]
 [0 0 0 0 0]]
normal form reduction: (1, 2, 0, 4, 3)
reduced graph:
[[0 0 1 0 1]
 [0 0 1 0 0]
 [1 1 0 0 1]
 [0 0 0 0 0]
 [1 0 1 0 0]]
colors reduced: [0, 0, 0, 1, 1]

Apply a random permutation to the graph and reduce the graph again.

[5]:

perm = np.random.permutation(5)
iperm = inverse_permutation(perm)
print("random permutation: %s" % (perm,))
graph2 = graph[perm, :][:, perm]
colors2 = [colors[idx] for idx in perm]

random permutation: [3 2 0 4 1]

Show the transformed matrix and color vector.

[8]:

print(graph2)
print(colors2)

[[0 0 0 0 0]
 [0 0 0 1 1]
 [0 0 0 0 1]
 [0 1 0 0 1]
 [0 1 1 1 0]]
[1, 1, 0, 0, 0]

[6]:

graph2_reduced, colors2_reduced, tr2 = reduce(graph2, colors2)

print("input graph: ")
print(graph2)

print("tr2: %s" % (tr2,))
print("reduced graph: ")
print(graph2_reduced)

print("colors2_reduced: %s" % (colors2_reduced,))

input graph:
[[0 0 1 0 1]
 [0 0 1 0 0]
 [1 1 0 0 1]
 [0 0 0 0 0]
 [1 0 1 0 0]]
tr2: (4, 1, 2, 3, 0)
reduced graph:
[[0 0 1 0 1]
 [0 0 1 0 0]
 [1 1 0 0 1]
 [0 0 0 0 0]
 [1 0 1 0 0]]
colors2_reduced: [0, 0, 0, 1, 1]

Check that the two reduced graphs are equal.

[9]:

if np.all(graph_reduced == graph2_reduced):
    print("reduced arrays are equal!")
if np.all(colors_reduced == colors2_reduced):
    print("reduced colors are equal!")

reduced arrays are equal!
reduced colors are equal!

Analyse isomorphisms of a set of orthogonal arrays with N=56

In this example, we show how to identify isomorphic arrays from a list of two-level strength-2 orthogonal arrays with 56 runs and 28 factors.

[1]:

import numpy as np
import oapackage
import oapackage.graphtools
from oapackage.graphtools import selectIsomorphismClasses

Read the arrays from a file and determine their class.

[2]:

sols=oapackage.readarrayfile('OAN56K28.oa')
arrayclass=oapackage.arraylink2arraydata(sols[0])
print(arrayclass)
print('loaded %d arrays' % len(sols))

arrayclass: N 56, k 28, strength 2, s {2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, order 0
loaded 9 arrays

Determine the unique isomorphism classes using Nauty.

[3]:

b,mm=selectIsomorphismClasses(sols, verbose=0)

print('from %d arrays selected %d unique isomorphism classes' % (len(mm),np.unique(b).size))
print('indices: %s' % str(b) )

from 9 arrays selected 6 unique isomorphism classes
indices: [0 4 5 3 1 5 3 2 4]

Verify that the first two arrays are indeed non-isomorphic:

[4]:

jj=np.abs(sols[0].Jcharacteristics(4))
n0, _=np.histogram(jj, [0,8,16,24])
print(n0)
jj=np.abs(sols[1].Jcharacteristics(4))
n1, _=np.histogram(jj, [0,8,16,24])
print(n1)

[    0 18018  2457]
[    0 18032  2436]

Since the first two arrays have different \(J\)-characteristics (see the section Statistical properties of orthogonal arrays for details), they are non-isomorphic. For the isomorphic arrays, it is possible to obtain the array transformation to make the arrays identical using the function reduceConferenceTransformation.

Generation and analysis of conference designs

In this notebook, we show how to generate conference designs and calculate properties of these designs. For details on conference designs and their properties, see Properties of conference designs and A Classification Criterion for Definitive Screening Designs, Schoen et al., 2018 and [Enumeration and Classification of Definitive Screening Designs] (in preparation).

Load required libraries and define the class of conference designs to enumerate.

[3]:

import oapackage
conference_class=oapackage.conference_t(12, 6, 0)
print(conference_class)

conference class: number of rows 12, number of columns 6

Define the root array and extend the lists of conference designs.

[4]:

conference_designs=[[conference_class.create_root_three_columns()]]

for ii, ncols in enumerate(range(4, 8)):
    arrays = oapackage.extend_conference (conference_designs[ii], conference_class, verbose=0)
    conference_designs.append(arrays)
    print('extension resulted in %d designs with %d columns' % (len(arrays), ncols))

extension resulted in 9 designs with 4 columns
extension resulted in 42 designs with 5 columns
extension resulted in 123 designs with 6 columns
extension resulted in 184 designs with 7 columns

Calculate properties of conference designs

Here, we show how to calculate relevant properties of conference designs. Select a 12-run 7-factor conference design generated previously.

[5]:

design = conference_designs[4][0]
design.showarray()

array:
  0   1   1   1   1   1   1
  1   0  -1  -1  -1  -1  -1
  1   1   0  -1  -1   1   1
  1   1   1   0   1  -1  -1
  1   1   1  -1   0  -1   1
  1   1  -1   1   1   0  -1
  1   1  -1   1  -1   1   0
  1  -1   1   1  -1   1  -1
  1  -1   1   1  -1  -1   1
  1  -1   1  -1   1   1  -1
  1  -1  -1   1   1  -1   1
  1  -1  -1  -1   1   1   1

A sensible criterion to evaluate conference designs is the so-called \(F_4\)-vector (Schoen et al., 2019). We can compute the \(F_4\)-vector of a conference design as follows.

[6]:

design.FvaluesConference(4)

[6]:

(0, 25, 10)

Schoen et al., 2019 showed that conference designs are good building blocks for definitive screening designs (Xiao et al. 2012). The Orthogonal Array package can construct a definitive screening design from a conference design.

[9]:

dsd = oapackage.conference2DSD(design)
dsd.showarray()

array:
  0   1   1   1   1   1   1
  1   0  -1  -1  -1  -1  -1
  1   1   0  -1  -1   1   1
  1   1   1   0   1  -1  -1
  1   1   1  -1   0  -1   1
  1   1  -1   1   1   0  -1
  1   1  -1   1  -1   1   0
  1  -1   1   1  -1   1  -1
  1  -1   1   1  -1  -1   1
  1  -1   1  -1   1   1  -1
  1  -1  -1   1   1  -1   1
  1  -1  -1  -1   1   1   1
  0  -1  -1  -1  -1  -1  -1
 -1   0   1   1   1   1   1
 -1  -1   0   1   1  -1  -1
 -1  -1  -1   0  -1   1   1
 -1  -1  -1   1   0   1  -1
 -1  -1   1  -1  -1   0   1
 -1  -1   1  -1   1  -1   0
 -1   1  -1  -1   1  -1   1
 -1   1  -1  -1   1   1  -1
 -1   1  -1   1  -1  -1   1
 -1   1   1  -1  -1   1  -1
 -1   1   1   1  -1  -1  -1
  0   0   0   0   0   0   0

For the resulting definitive screening design, the Orthogonal Array package can compute some statistical properties based on projections into a smaller number of factors.

[7]:

PEC4, PIC4, PPC4 = oapackage.conference.conferenceProjectionStatistics(design, ncolumns=4, verbose=1)

conferenceProjectionStatistics: projection to 4 columns: PEC 0.286 PIC 3.111 PPC 0.458

Example code for delete-one-factor projections

Any orthogonal array can be reduced to delete-one-factor projection form using reduceDOPform. The method is described in the article A canonical form for non-regular arrays based on generalized wordlength pattern values of delete-one-factor projections by Eendebak, P. T. In this notebook, we reduce an example array to its delete-one-factor projection form.

Load required libraries and select an example orthogonal array with 16 runs and 7 factors.

[4]:

import oapackage

al = oapackage.exampleArray(4)
al = oapackage.reduceDOPform(al)
al.showarray()

array:
  0   0   0   0   0   0   0
  0   0   0   0   0   0   1
  0   0   0   1   1   1   0
  0   0   0   1   1   1   1
  0   1   1   0   0   1   0
  0   1   1   0   0   1   1
  0   1   1   1   1   0   0
  0   1   1   1   1   0   1
  1   0   1   0   1   0   0
  1   0   1   0   1   1   1
  1   0   1   1   0   0   0
  1   0   1   1   0   1   1
  1   1   0   0   1   0   1
  1   1   0   0   1   1   0
  1   1   0   1   0   0   1
  1   1   0   1   0   1   0

A key property of the delete-of-factor projection form is that the generalized word length patterns (GWLPs) of the projections are ordered. For details on the GWLP, see Statistical properties of orthogonal arrays.

[5]:

print("GWLP %s" % str(al.GWLP()))
for ii in range(0, al.n_columns):
    bl = al.deleteColumn(ii)
    print("Delete column %d: GWLP %s" % (ii, str(bl.GWLP())))

GWLP (1.0, 0.0, 0.0, 3.5, 2.5, 0.5, 0.5, 0.0)
Delete column 0: GWLP (1.0, 0.0, 0.0, 1.5, 1.0, 0.5, 0.0)
Delete column 1: GWLP (1.0, 0.0, 0.0, 1.75, 0.75, 0.25, 0.25)
Delete column 2: GWLP (1.0, 0.0, 0.0, 1.75, 0.75, 0.25, 0.25)
Delete column 3: GWLP (1.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0)
Delete column 4: GWLP (1.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0)
Delete column 5: GWLP (1.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0)
Delete column 6: GWLP (1.0, 0.0, 0.0, 3.0, 2.0, 0.0, 0.0)

The symmetry group of the projection GWLPs can be calculated. This symmetry group determines how fast an array can be reduced to normal form.

[6]:

dof_values = oapackage.projectionDOFvalues(al)
sg = oapackage.symmetry_group(dof_values, False)
sg.show(1)

symmetry group: 7 elements, 4 subgroups: 1 2 3 1

It is also possible to reduce mixed-level arrays to their delete-of-factor projection forms. We now show an example involving an orthogonal array with 24 runs, one four-level factor, one three-level factor and three two-level factors.

[8]:

al = oapackage.exampleArray(5)
al.showarray()
dopgwlp = oapackage.projectionGWLPs(al)
print("delete-one-factor GWLP values %s" % ([x.raw_values() for x in dopgwlp],))

array:
  0   0   0   0   0
  0   0   0   0   0
  0   1   0   0   0
  0   1   1   1   1
  0   2   1   1   1
  0   2   1   1   1
  1   0   0   1   1
  1   0   0   1   1
  1   1   0   0   1
  1   1   1   1   0
  1   2   1   0   0
  1   2   1   0   0
  2   0   1   0   1
  2   0   1   0   1
  2   1   0   1   0
  2   1   1   0   1
  2   2   0   1   0
  2   2   0   1   0
  3   0   1   1   0
  3   0   1   1   0
  3   1   0   1   1
  3   1   1   0   0
  3   2   0   0   1
  3   2   0   0   1
delete-one-factor GWLP values [(1.0, 0.0, 0.0, 0.0, 0.6666666666666666), (1.0, 0.0, 0.0, 2.111111111111111, 0.0), (1.0, 0.0, 0.0, 1.8888888888888888, 0.4444444444444444), (1.0, 0.0, 0.0, 2.3333333333333335, 0.0), (1.0, 0.0, 0.0, 1.8888888888888888, 0.4444444444444444)]

The delete-of-factor values for mixel-level arrays consists of the factor level of the deleted column and the GWLP of the projection.

[9]:

oapackage.projectionGWLPs
dopgwp = oapackage.projectionGWLPs(al)
arrayclass = oapackage.arraylink2arraydata(al)
dofvalues = oapackage.projectionDOFvalues(al)
factor_levels = arrayclass.factor_levels()

for column, dof_value in enumerate(dofvalues):
    print("column %d: factor level %s" % (column, factor_levels[column]))
    print("   delete-of-factor value: %s" % (list(dof_value.raw_values()),))

column 0: factor level 4
   delete-of-factor value: [-4.0, 1.0, 0.0, 0.0, 0.0, 0.6666666666666666]
column 1: factor level 3
   delete-of-factor value: [-3.0, 1.0, 0.0, 0.0, 2.111111111111111, 0.0]
column 2: factor level 2
   delete-of-factor value: [-2.0, 1.0, 0.0, 0.0, 1.8888888888888888, 0.4444444444444444]
column 3: factor level 2
   delete-of-factor value: [-2.0, 1.0, 0.0, 0.0, 2.3333333333333335, 0.0]
column 4: factor level 2
   delete-of-factor value: [-2.0, 1.0, 0.0, 0.0, 1.8888888888888888, 0.4444444444444444]

Show the reduced array.

[10]:

reduced_array = al.reduceDOP()
al.showarray()

array:
  0   0   0   0   0
  0   0   0   0   0
  0   1   0   0   0
  0   1   1   1   1
  0   2   1   1   1
  0   2   1   1   1
  1   0   0   1   1
  1   0   0   1   1
  1   1   0   0   1
  1   1   1   1   0
  1   2   1   0   0
  1   2   1   0   0
  2   0   1   0   1
  2   0   1   0   1
  2   1   0   1   0
  2   1   1   0   1
  2   2   0   1   0
  2   2   0   1   0
  3   0   1   1   0
  3   0   1   1   0
  3   1   0   1   1
  3   1   1   0   0
  3   2   0   0   1
  3   2   0   0   1

Minimal number of runs for an orthognal array

In this example we calculate the minimum number of runs required for an orthogonal array. For more details, see Schoen et al.

[6]:

import itertools
import numpy as np


def minimum_number_of_runs(factor_levels, strength):
    """Calculate the minimum number of runs for an orthogonal array

    The minimum number of runs is based on the strength conditions. IWhether a design actually exists
    Args:
        factor_levels: Factor levels of the design
        strength: Strength of the array
    Returns:
        Minimum number of runs

    """
    runs = [np.prod(tt) for tt in itertools.combinations(factor_levels, strength)]
    N = np.lcm.reduce(runs)
    return N

We run the method on several examples.

[7]:

strength = 3
factor_levels = [2, 3, 3, 4, 5]
N = minimum_number_of_runs(factor_levels, strength)

print(f"for a design of strength {strength} and factor levels {factor_levels} we require (a multiple of) {N} runs")

for a design of strength 3 and factor levels [2, 3, 3, 4, 5] we require (a multiple of) 360 runs

[12]:

strength = 2
factor_levels = [3, 2, 2, 2, 2]
N = minimum_number_of_runs(factor_levels, strength)

print(f"for a design of strength {strength} and factor levels {factor_levels} we require (a multiple of) {N} runs")

for a design of strength 2 and factor levels [3, 2, 2, 2, 2] we require (a multiple of) 12 runs

[ ]:

Example of GWLP calculation for mixed-level designs

This notebook contains example code to compute the type-specific generalized word-length pattern (GWLP) for mixed-level designs, more specifically for regular four-and-two-level designs.

In four-and-two-level designs, the four-level factors are constructed us the grouping scheme of Wu & Zhang (1989) where the levels of a four-level factor \(A\) are based on the levels of three two-level factors \(a_1\), \(a_2\) and \(a_3\), called the pseudo-factors, where \(I=a_1a_2a_3\):

\[\begin{split}\begin{array}{ccccc} a_1 & a_2 & a_3 & & A \\ 1 & 1 & 0 & \rightarrow & 0 \\ 1 & 0 & 1 & \rightarrow & 1 \\ 0 & 1 & 1 & \rightarrow & 2 \\ 0 & 0 & 0 & \rightarrow & 3 \\ \end{array}\end{split}\]

The two-level factors can either be main factors or be entirely aliased with a combination of the main factors. Such aliased factors are called added factors. If one of the main factors, used in an added factor, is also used as pseudo-factor in a four-level factor, then the added factor has type \(I\). If it is used as pseudo-factor in two distinct four-level factor, it has type \(II\), etc \(\ldots\).

[1]:

import numpy as np

import oapackage

Create a \(4^{1}2^{7}\) four-and-two-level regular design in 32 runs with 1 four-level factor and 7 two-level factors.

[2]:

array = oapackage.exampleArray(56, 1)
array = array.selectFirstColumns(7)
arrayclass = oapackage.arraylink2arraydata(array)
array.showarraycompact()

exampleArray 56: design in OA(32, 42^{18})
0000000
0100110
0101010
0001100
0110011
0010101
0011001
0111111
1010100
1110010
1111110
1011000
1100111
1000001
1001101
1101011
2101101
2001011
2000111
2100001
2011110
2111000
2110100
2010010
3111001
3011111
3010011
3110101
3001010
3101100
3100000
3000110

The GWLP for mixed-level designs is computed using the adapted forms of equations (7) and (8) in Xu and Wu (2001).

Distance distribution

First, the distance distribution is first computed for all combinations of \((i,j)\) with \(i \in \{0,1\}\) and \(j=0,\ldots,7\) and \(i\neq j\). In this example with a \(4^12^{7}\) design, that is a \(2 \times 8\) matrix \(B\), where \(B_{i,j}\) is the number of rows that have \((1-i)\) different elements if their four-level parts and \((7-j)\) different elements in their two-level parts, divided by the number of runs.

[3]:

Dm = oapackage.distance_distribution_mixed(array, 0)
D = np.array(Dm).astype(int)
print(f"distance distribution for mixed-level design\n{D}")

distance distribution for mixed-level design
[[1 0 1 4 1 0 1]
 [0 2 6 8 6 2 0]]

This is verified since the sum of the matrix is 32 which is equal to \(\binom{32}{2} 32^{-1}\).

McWilliams Transforms

Now, the second step in the computation of the GWLP, is to use the McWilliams Transforms to obtain the GWLP from the distance distribution matrix. The McWilliams Transform will create another \(2 \times 8\) matrix \(B^{\prime}\), obtained using the following formula

\[B_{j_{1}, j_{2}}^{\prime}=N^{-1} \sum_{i_{1}=0}^{m} \sum_{i_{2}=0}^{n} B_{i_{1}, i_{2}} P_{j_{1}}\left(i_{1} ; 1, 4\right) P_{j_{2}}\left(i_{2} ; 7, 2\right)\]

where \(P_j(x,n,s)\) is the Krawtchouck polynomials for a total of \(n\) factors with \(s\) levels.

[4]:

N = array.n_rows
factor_levels_for_groups = arrayclass.factor_levels_column_groups()
Bprime = oapackage.macwilliams_transform_mixed(Dm, N, factor_levels_for_groups, verbose=0)
print(f"MacWilliams transform:")
print(np.array(Bprime).astype(int))

MacWilliams transform:
[[1. 0. 0. 0. 1. 0. 0.]
 [0. 0. 2. 0. 4. 0. 0.]]

This matrix is equivalent to the type specific generalized word-length pattern, where the row index indicates the type of words (0 to 1) and the column index indicates the length of the words (0 to 7).

The generalized word-length pattern (\(A\)) can be obtained by summing the rows of the \(B^{\prime}\) matrix anti-diagonally. That is:

\[A_j(D) = \sum_{i^{\prime}+j^{\prime}=j} B^{\prime}_{i^{\prime},j^{\prime}}\]

[4]:

gwlp_mixed = oapackage.GWLPmixed(array, 0)
print(f"GWLP: {gwlp_mixed}")

GWLP: (1.0, 0.0, 0.0, 2.0, 1.0, 4.0, 0.0, 0.0)

Calculate a Pareto optimal set

Pareto optimality (or multi-objective optimization) allows one to search for optimal solutions for an optimization problem with multiple objectives. The Pareto class in the Orthogonal Array package allows one to calculate the Pareto optimal elements (called the Pareto frontier).

[1]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import oapackage

First, define a dataset of 50 random vectors. The vectors have length 2, so there are 2 objectives to be optimized.

[2]:

datapoints=np.random.rand(2, 50)

for ii in range(0, datapoints.shape[1]):
    w=datapoints[:,ii]
    fac=.6+.4*np.linalg.norm(w)
    datapoints[:,ii]=(1/fac)*w

h=plt.plot(datapoints[0,:], datapoints[1,:], '.b', markersize=16, label='Non Pareto-optimal')
_=plt.title('The input data', fontsize=15)
plt.xlabel('Objective 1', fontsize=16)
plt.ylabel('Objective 2', fontsize=16)

Create a structure (ParetoDoubleLong) to keep track of the data.

[3]:

pareto=oapackage.ParetoDoubleLong()

for ii in range(0, datapoints.shape[1]):
    w=oapackage.doubleVector( (datapoints[0,ii], datapoints[1,ii]))
    pareto.addvalue(w, ii)

pareto.show(verbose=1)

Pareto: 12 optimal values, 12 objects

Plot the results.

[4]:

lst=pareto.allindices() # the indices of the Pareto optimal designs

optimal_datapoints=datapoints[:,lst]

h=plt.plot(datapoints[0,:], datapoints[1,:], '.b', markersize=16, label='Non Pareto-optimal')
hp=plt.plot(optimal_datapoints[0,:], optimal_datapoints[1,:], '.r', markersize=16, label='Pareto optimal')
plt.xlabel('Objective 1', fontsize=16)
plt.ylabel('Objective 2', fontsize=16)
plt.xticks([])
plt.yticks([])
_=plt.legend(loc=3, numpoints=1)

Isomorphism reduction for conference designs

In this example, we show how to test if two conference designs are isomorphic. We consider conference designs with 10 runs and 3 factors, and calculate a reduction to their normal form. Using the reduction, we determine if the two designs are isomorphic.

[1]:

import oapackage
import numpy as np

A = np.array(
    [0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1, 1, -1, 1, 1, -1, -1, 1, -1, -1]
).reshape(10, 3)
B = np.array(
    [0, 1, 1, 1, 0, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 0, 1, -1, 1, 1, -1, 1, 1, -1, -1]
).reshape(10, 3)
array1 = oapackage.makearraylink(A)
array2 = oapackage.makearraylink(B)

array1.showarray()
array2.showarray()

We calculate the normal forms of the conference designs using the function reduceConference or reduceConferenceTransformation. The result of the former is the reduced design, while the result of the latter is an object describing the tranformation to normal form. The normal form is calculated using Nauty.

[2]:

help(oapackage.reduceConference)

Help on function reduceConference in module oalib:

reduceConference(arg1, verbose=0)
    reduceConference(array_link arg1, int verbose=0) -> array_link
    reduceConference(array_link arg1) -> array_link



    Reduce conference matrix to normal form using Nauty

    See also: reduceConferenceTransformation

[3]:

T1 = oapackage.reduceConferenceTransformation(array1, verbose=1)
T2 = oapackage.reduceConferenceTransformation(array2, verbose=1)
T1.show()

reduceConferenceTransformation: reduce design with 10 rows, 3 columns
reduceConferenceTransformation: reduce design with 10 rows, 3 columns
row permutation: {2,0,1,9,7,8,5,6,3,4}
  row flips: {1,1,-1,1,1,1,1,1,1,1}
column permutation: {0,1,2}
  col flips: {1,-1,-1}

We can check whether the designs are isomorphic by comparing the normal forms.

[4]:

design_equal = T1.apply(array1) == T2.apply(array2)
print("designs isomorphic? %s" % design_equal)

designs isomorphic? 1

The designs are isomorphic. So, it is possible to calculate a reduction of the second design into the first design.

[5]:

TT = T1.inverse() * T2
TT.show()

row permutation: {0,1,8,9,6,7,2,4,5,3}
  row flips: {-1,1,1,1,1,1,1,1,1,1}
column permutation: {0,1,2}
  col flips: {1,-1,-1}

[6]:

r1 = T1.apply(array1)
r1.showarray()
r2 = T2.apply(array2)
r2.showarray()

Calculate some properties of a conference design; see the section Properties of conference designs for details.

[7]:

print("array 1: is_conference() %d" % array1.is_conference())
print("array 1: J2-characteristics %s" % (oapackage.Jcharacteristics_conference(array1, number_of_columns=2),))

array 1: is_conference() 1
array 1: J2-characteristics (0, 0, 0)

[8]:

help(oapackage.Jcharacteristics_conference)

Help on function Jcharacteristics_conference in module oalib:

Jcharacteristics_conference(array, number_of_columns, verbose=0)
    Jcharacteristics_conference(array_link array, int number_of_columns, int verbose=0) -> intVector
    Jcharacteristics_conference(array_link array, int number_of_columns) -> intVector

Data representation

All designs handled by the OApackage are integer valued. The designs (whether these are orthogonal arrays, optimal designs or conferences designs) are stored in an array_link() object. The definitions of orthogonal arrays, optimal designs and conference designs are included in the section Definitions of arrays and designs.

Data structures

The package contains several data structures. Here, we describe the main structures and their use.

array_link()
The structure containing an orthogonal array is called the array_link() structure. Lists of arrays are stored in the arraylist_t() object, which is implemented as a std::deque container.

arrayfile_t()
This object allows for reading and writing of arrays to disk.

arraydata_t()
The structure describing a certain class of orthogonal arrays or optimal designs.

conference_t()
The structure describing a certain class of conference designs.

array_transformation_t()
The structure describing a transformation of an orthogonal array, which includes row and column permutations, as well as permutations to the symbols in one or more columns.

conference_transformation_t()
The structure describing a transformation of conference design or double conference design, which includes row and column permutations, as well as sign switches to the elements in one or more rows and columns.

Representing arrays

The structure containing an orthogonal array is called the array_link structure. It consists of a specified number of rows and columns, the data (integer valued) and an index. In the Python interface, the array_link() object can be indexed just as normal arrays.

It is also possible to convert to a Numpy array. The array_link object implements the Python array interface, so most operations from packages such as Numpy work on the array_link() object.

Array representation and indexing in Python

>>> import oapackage; import numpy as np
>>> al=oapackage.exampleArray(0)
>>> al.showarray() 
array:
  0   0
  0   0
  0   1
  0   1
  1   0
  1   0
  1   1
  1   1

>>> al[2,1]
1
>>> X=np.array(al)
>>> X
array([[0, 0],
       [0, 0],
       [0, 1],
       [0, 1],
       [1, 0],
       [1, 0],
       [1, 1],
       [1, 1]], dtype=int16)

The C++ class is array_link.

Classes of arrays

The arraydata_t object represents data about a class of orthogonal arrays, e.g. the class \({\operatorname{OA}(N; t; s^k)}\); see Definitions of arrays and designs. The conference_t object represents data about a class of conference designs.

Array transformations

Transformations of (orthogonal) arrays consist of row, column and level permutations. A transformation is represented by the array_transformation_t object.

For a given transformation, the column permutations are applied first, then the level permutations and finally the row permutations. The level and column permutations are not commutative.

The conference design transformations also allow for row sign switches and are described by the class conference_transformation_t.

Reading and writing arrays

Reading and writing arrays to disk can be done with the arrayfile_t class.

Write an array or a list of arrays to disk

>>> import oapackage
>>> list_of_arrays = [oapackage.exampleArray(24), oapackage.exampleArray(25)]
>>> _ = oapackage.writearrayfile('test.oa', list_of_arrays)
>>> oapackage.oainfo('test.oa')
file test.oa: 64 rows, 16 columns, 2 arrays, mode text, nbits 0
>>> al=oapackage.exampleArray()
>>> af=oapackage.arrayfile_t('test.oa', al.n_rows, al.n_columns)
>>> af.append_array(al)
>>> print(af)
file test.oa: 8 rows, 2 columns, 1 arrays, mode text, nbits 8
>>> af.closefile()

The arrays can be written in text or binary format. For more details on the file format, see the section File formats.

The Python interface is oalib.arrayfile_t() and the C++ interface is arrayfile_t.

File formats

The Orthogonal Array package stores arrays in a custom file format. There is a text format which is easily readable by humans and a binary format which is faster to process and memory efficient.

Plain text array files

Arrays are stored in plain text files with extension .oa. The first line contains the number of columns, the number of rows and the number of arrays (or -1 if the number of arrays is not specified). Then, for each array, a single line with the index of the array, followed by \({N}\) lines containing the array.

A typical example of a text file is the following:

This file contains exactly 1 array with 8 rows and 5 columns.

Binary array files

Every binary file starts with a header, which has the following format:

[INT32] 65 (magic identifier)
[INT32] b: Format: number of bits per number. Currently supported are 1 and 8
[INT32] N: number of rows
[INT32] k: kumber of columns
[INT32] Number of arrays (can be -1 if unknown)
[INT32] Binary format number: 1001: normal, 1002: binary diff, 1003: binary diff zero
[INT32] Reserved integer
[INT32] Reserved integer

The format of the remainder of the binary file depends on the binary format specified. For the normal binary format, the format is as follows. For each array, the number is specified in the header:

[INT32] Index
[Nxk elements] The elements contain b bits

If the number of bits per number is 1 (e.g. a 2-level array), then the data is padded with zeros to a multiple of 64 bits. The data of the array is stored in column-major order. The binary file format allows for random access reading and writing. The binary diff and binary diff zero formats are special formats.

A binary array file can be compressed using gzip. Most tools in the Orthogonal Array package can read these compressed files transparently. Writing to compressed array files is not supported at the moment.

Data files

The analysis tool (oaanalyse) writes data to disk in binary format. The format consists of a binary header:

[FLOAT64] Magic number 30397995;
[FLOAT64] Magic number 12224883;
[FLOAT64] nc: Number of rows
[FLOAT64] nr: Number of columns

After the header there follow nc*nr [FLOAT64] values.

MD5 sums

To check data integrity on disk, the packages includes functions to generate MD5 sums of designs.

Calculate md5 sum of a design

>>> import oapackage; al=oapackage.exampleArray(0)
>>> al.md5()
'6454c492239a8e01e3c01a864583abf2'

The C++ functions are array_link::md5() and md5().

Command line interface

Several command line tools are included in the Orthogonal Array package. For each tool, help can be obtained from the command line by using the switch -h. The tools include the following:

oainfo

This program reads Orthogonal Array package data files and reports the contents of the files. For example:

$ oainfo result-8.2-2-2-2.oa
Orthogonal Array package 1.8.7
oainfo: reading 1 file(s)
file result-8.2-2-2.oa: 8 rows, 3 columns, 2 arrays, mode text, nbits 0
$

oacat

Shows the contents of a file with orthogonal arrays for a data file.

oacheck

Checks or reduces an array to canonical form.

oaextendsingle

Extends a set of arrays in LMC form with one or more columns.

oacat

Shows the contents of an array file or data file.

Usage: oacat [OPTIONS] [FILES]

oajoin

Reads one or more files from disk and join all the array files into a single list.

oasplit

Takes a single array file as input and splits the arrays into a specified number of output files.

oapareto

Calculates the set of Pareto optimal arrays in a file with arrays.

oaanalyse

Calculates various statistical properties of arrays in a file. The properties are described in section Properties of designs.

alternate text — Orthogonal array in \(\mathrm{OA}(18, 2 3^a, 2)\).

Generation of designs

The Orthogonal Array package can be used to generate several classes of arrays and designs. A collection of arrays and designs generated by the package is available on the website http://www.pietereendebak.nl/oapackage/index.html.

Generation of orthogonal arrays

A list of arrays in LMC form (i.e., lexicographically minimum in columns) can be extended to a list of arrays in LMC form with one additional column. Details about the algorithm are described in [SEN10].

The main function for array extension is the function extend_arraylist(). The arguments for this function are the list of arrays to extend, a specification of the class of arrays in arraydata_t and the options OAextend for the algorithm.

An example of a session that extends an array is:

>>> import oapackage
>>> nrows=8; ncols=3;
>>> arrayclass=oapackage.arraydata_t(2, nrows, 2, ncols)
>>> root_array=arrayclass.create_root()
>>> root_array.showarraycompact()
00
00
01
01
10
10
11
11
>>> array_list=oapackage.extend_array(root_array, arrayclass)
>>> print('found %d extensions of the root array' % len(array_list))
found 2 extensions of the root array

A more detailed example is included in Enumerate orthogonal arrays.

Conference designs

A conference design is an \(N\times k\) matrix with entries 0, -1, +1 such that i) in each column the symbol 0 occurs exactly one time and ii) all columns are orthogonal to each other. For details on conference designs, see the section Properties of conference designs and [SEG19]. An example of a session to generate conference designs is the following:

Generate conference designs with 8 rows

>>> import oapackage
>>> conference_class=oapackage.conference_t(8, 6, 0)
>>> array = conference_class.create_root_three_columns()
>>> array.showarray()
array:
  0   1   1
  1   0  -1
  1   1   0
  1   1   1
  1   1  -1
  1  -1   1
  1  -1   1
  1  -1  -1
>>> l4=oapackage.extend_conference ([array], conference_class, verbose=0)
>>> l5=oapackage.extend_conference ( l4, conference_class, verbose=0)
>>> l6=oapackage.extend_conference ( l5, conference_class, verbose=0)
>>> print('number of non-isomorphic conference designs with 6 columns: %d' % len(l6) )
number of non-isomorphic conference designs with 6 columns: 11

An example notebook with more functionality is included in Generation and analysis of conference designs. The full interface for conference designs is available in the Interface for conference designs.

The main functions to extend conference and double conference designs are extend_conference() and extend_double_conference(), respectively. The low-level functions for generating candidate extension columns of conference and double conference designs are generateConferenceExtensions() and generateDoubleConferenceExtensions(), respectively.

The conference designs are generated in LMC0 form.

Calculation of D-efficient designs

D-efficient designs (sometimes called D-optimal designs) can be calculated with the function oapackage.Doptim.Doptimize(). This function uses a coordinate-exchange algorithm to generate designs with good properties for the \(D\)-efficiency. With the coordinate-exchange algorithm, the following target function \(T\) is optimized:

\[T = \alpha_1 D_{\text{eff}} + \alpha_2 D_{s, \text{eff}} + \alpha_3 D_{1, \text{eff}}\]

Here, \(\alpha\) is a weight vector specified by the user. Details on the \(D_{\text{eff}}\), \(D_{s, \text{eff}}\) and \(D_{1, \text{eff}}\) can be found in the section Optimality criteria for D-efficient designs.

A Python script to generate D-efficient designs with 40 runs and 7 factors is shown below.

Example of Doptimize usage

>>> N=40; s=2; k=7;
>>> arrayclass=oapackage.arraydata_t(s, N, 0, k)
>>> print('We generate optimal designs with: %s' % arrayclass)
We generate optimal designs with: arrayclass: N 40, k 7, strength 0, s {2,2,2,2,2,2,2}, order 0
>>> alpha=[1,2,0]
>>> scores, dds, designs, ngenerated = oapackage.Doptimize(arrayclass, nrestarts=40, optimfunc=alpha, selectpareto=True, verbose=0)
Doptimize: iteration 0/40
Doptimize: iteration 39/40
>>> print('Generated %d designs, the best D-efficiency is %.4f' % (len(designs), dds[:,0].max() ))
Generated 10 designs, the best D-efficiency is 0.9198

The parameters of the Doptimize() function are documented in the code.

To calculate the \(D\)-, \(D_s\)- and \(D_1\)-efficiencies of the designs, we can use the method Defficiencies(). For details on these efficiencies, see the section Optimality criteria for D-efficient designs and [ES17].

In [ES17], it is shown that one can optimize a linear combination of the \(D\)-efficiency and \(D_s\)-efficiency to generate a rich set of good compromise designs. From the generated designs, the optimal ones according to Pareto optimality can be selected.

_images/motivating-40-d-2-2-2-2-2-2-2-scatterplot-ndata2.png — Scatterplot for the \(D\)-efficiency and \(D_s\)-efficiency for generated designs in \({\operatorname{OA}(40; 2; 2^7)}\). The Pareto optimal designs are colored, while the non-Pareto optimal designs are grey. For reference the strength-3 orthogonal array with highest D-efficiency is also included in the plot.

Even-odd arrays

The even-odd arrays are a special class of orthognal arrays with at least one of the odd \(J_k\)-characteristics unequal to zero. More information on this class of designs will appear later.

Normal form of arrays

The Orthogonal Array package contains functions to reduce arrays and designs to canonical form with respect to some ordering. The default ordering for orthogonal arrays is the lexicographic ordering in columns [SEN10]. The default ordering for conference designs is the LMC0 ordering [SEG19]. Alternative orderings include the delete-one-factor projection ordering introduced in [Een13] or the even-odd ordering. For a given ordering of a set of arrays, the minimal element of all arrays in an isomorphism class defines a unique representative of that isomorphism class.

Specialized packages such as Nauty [McK81], [MP13] can also reduce arrays to their canonical form using state-of-the-art, graph-isomorphism methods. However, these methods do not take into account the special structure of the arrays and so, they cannot be tailored to create normal forms of a specific form.

Reduction to LMC normal form

The Orthogonal Array package implements theory and methods from the article Complete enumeration of pure-level and mixed-level orthogonal arrays, Schoen et al. to reduce orthogonal arrays to their LMC normal form. The C++ function to perform the reduction is reduceLMCform(). An example on how to use this function is shown below.

>>> import oapackage
>>> oapackage.set_srand(1)
>>> array = oapackage.exampleArray(1, 0).selectFirstColumns(3)
>>> array = array.randomperm()
>>> print('input array:'); array.transposed().showarraycompact()
input array:
1100010111010001
0101100110100011
0000001111101101
>>> reduced_array = oapackage.reduceLMCform(array)
>>> print('reduced array:'); reduced_array.transposed().showarraycompact()
reduced array:
0000000011111111
0000111100001111
0001011101110001

It is also possible to check whether an array is in normal form with the LMCcheck() method:

>>> import oapackage
>>> array = oapackage.exampleArray(1)
>>> lmc_type = oapackage.LMCcheck(array)
>>> if lmc_type == oapackage.LMC_MORE:
...      print('array is in minimal form')
... elif lmc_type == oapackage.LMC_LESS:
...      print('array is not in minimal form')
array is in minimal form

Reduction to delete-one-factor projection form

The article A canonical form for non-regular arrays based on generalized word length pattern values of delete-one-factor projections [Een13] describes a canonical form of an orthogonal array based on delete-one-factor projections. The C++ interface to delete-one-factor projection form is reduceDOPform(). The reduction method works well for large arrays with a large variation in the projection values.

An example on how to use this reduction is shown in Example code for delete-one-factor projections, which can be found in the example notebooks section.

Reduction using graph isomorphisms

The function reduceOAnauty() reduces an orthogonal array to Nauty canonical form. To reduce general graphs to Nauty canonical form, the Orthogonal Array package includes the function reduceGraphNauty().

Reduce a design to normal form using Nauty

>>> oapackage.set_srand(1)
>>> al = oapackage.exampleArray(0).randomperm()
>>> al.showarray()
array:
  0   0
  0   1
  1   1
  0   1
  1   0
  0   0
  1   0
  1   1
>>> transformation=oapackage.reduceOAnauty(al, 0)
>>> transformation.show()
array transformation: N 8
column permutation: {0,1}
level perms:
{0,1}
{0,1}
row permutation: {0,5,1,3,4,6,2,7}
>>> alr=transformation.apply(al)
>>> alr.showarray()
array:
  0   0
  0   0
  0   1
  0   1
  1   0
  1   0
  1   1
  1   1

Normal forms for conference designs

For conference designs, a convenient normal form is the LMC0 ordering (sometimes also called L0 ordering) [SEG19].

LMC0 ordering

The LMC0 ordering for conference designs is defined in three steps:

Definition LMC0 i: Order of elements
The LMC0 order of the factor levels -1, 0 and +1 is 0 < +1 < -1.

Definition LMC0 ii: Order of columns
A column a is smaller than a column b according to LMC0 ordering, notated as a < b, if either of the following conditions hold:

If we replace the values -1 by +1 in both columns, then the first element where the columns differ is smaller in a than in b according to Definition 1.

The zeros in column a are in the same position as the zeros in column b, and the first element where the columns differ is smaller in a than in b according to Definition i.

Definition LMC0 iii: Order of designs
Conference design A is smaller than conference design B according to LMC0 ordering, notated as A < B, if the first column where the designs differ is smaller in A than in B.

The definition implies that the ordering of designs is column-by-column and that the position of zeros in the columns is dominant over the values +1, -1. To check whether a design is in LMC0 form, we can use LMC0check():

Conference design in normal form

>>> array = oapackage.exampleArray(53,1)
exampleArray 53: third array in C(12,4)
>>> array.showarray()
array:
  0   1   1   1
  1   0  -1   1
  1   1   0  -1
  1   1   1  -1
  1   1   1  -1
  1   1  -1   1
  1   1  -1   1
  1  -1   1   0
  1  -1   1   1
  1  -1   1   1
  1  -1  -1  -1
  1  -1  -1  -1
>>> oapackage.LMC0check(array) == oapackage.LMC_LESS
True

Properties of designs

This section shows the structural and statistical properties of the orthogonal arrays, conference designs and D-efficient designs generated by the Orthogonal Array package. The properties of the arrays and designs are calculated using the array_link object or functions from the package.

Definitions of arrays and designs

Before introducing the structural and statistical properties, we define orthogonal arrays, conference designs and D-efficient designs:

An orthogonal array (OA) of strength \({t}\), \({N}\) runs and \({n}\) factors at \({s}\) levels is an \({N}\times {n}\) array of symbols \(0, \ldots,({s}-1)\), such that for every subset of \({t}\) columns, every \({t}\)-tuple occurs equally often [Rao47]. The set of all strength-\({t}\) OAs with \({N}\) runs and \({n}\) factors at \({s}\) levels is denoted by \({\operatorname{OA}({N}; {t}; {s}^{n})}\). If \({s=2}\), the OA is called a two-level OA and the set of all strength-\({t}\) two-level OAs with \({N}\) runs and \({n}\) factors is denoted as \({\operatorname{OA}({N}; {t}; {2}^{n})}\).

For \({N}\) even, a conference design [SEG19] \(C\) is an \({N}\times {n}\) array which satisfies \({C}^{T}C = (n-1) I_{n}\), with \({C}_{ii} = 0\) and \({C}_{ij} \in \{-1,1\}\), for \({i} \neq {j}\) and \({i}, {j} = 1, \ldots, n\). A \({N}\times {N}\) conference design \(E\) such that \(E{E}^{T} = (n-1) I_{n}\) is called a conference matrix; see [EN95], [CD06] and [XLB12].

A D-optimal design [DAT07] (\(X\)) is an \({N}\times {n}\) array which maximizes the \(D\)-efficiency, defined as \({(\operatorname{det}({X}^{T}_{M}{X}^{\phantom{T}}_{M})^{1/p})/N}\), for a given \({N}\times {p}\) model matrix \({X}_{M}\) (for details see Model matrices). The Orthogonal Array package uses a coordinate-exchange algorithm to generate designs that optimize the \(D\)-efficiency. Since there is no guarantee that the resulting designs have the largest possible \(D\)-efficiency, we refer to them as D-efficient designs in this documentation. An orthogonal array is called D-optimal orthogonal array if it provides the largest \(D\)-efficiency among all comparable orthogonal arrays.

Structural properties of an array

The OApackage can calculate the rank of an array, defined as the maximum number of linearly independent column or row vectors in the array. The rank of an array is useful for several other functions in the package. For two-level arrays, the OApackage can also check if the arrays are foldover arrays. A two-level array is called a foldover array if half of its runs are mirror images of the other half, in the sense that the factor levels are changed from 0 to 1 and from 1 to 0.

For example, to calculate the rank of a two-level orthogonal array and determine whether the array is a foldover array, one can use the methods array_link::rank() and array_link::foldover():

Calculate rank of array and test for foldover

>>> array = oapackage.exampleArray(1) # Select an example two-level orthogonal array
>>> array.showarray() # Show the two-level orthogonal array
array:
  0   0   0   0   0
  0   0   0   0   0
  0   0   0   1   1
  0   0   1   0   1
  0   1   0   1   0
  0   1   1   0   0
  0   1   1   1   1
  0   1   1   1   1
  1   0   0   1   1
  1   0   1   0   1
  1   0   1   1   0
  1   0   1   1   0
  1   1   0   0   1
  1   1   0   0   1
  1   1   0   1   0
  1   1   1   0   0
 >>> print(array.rank()) # Calculate the rank of the array
 5
 >>> print(array.foldover()) # Determine if the array is foldover
 False

Other structural properties such as whether an array involves two levels or is symetric can be found in the documentation of array_link, which shows the full set of methods available.

Model matrices

For orthogonal arrays and conference designs, the OApackage can calculate different model matrices. The model matrices available depend on the type of array and design.

Model matrices for two-level orthogonal arrays

For two-level orthogonal arrays, the levels of the array are first coded according to the map \({0 \rightarrow -1}\) and \({1 \rightarrow +1}\). The coded matrix is referred to as the design matrix. The main effect contrast vectors are given by the columns in the design matrix. The contrast vectors associated to the two-factor interactions are calculated by taking products between two different columns in the design matrix. The model matrix consists of the intercept column (i.e. a columns of ones) and the contrast vectors associated to the main effects and, optionally, the two-factor interactions.

Model matrices for conference designs

The model matrix for a conference design consists of the intercept column (i.e. a columns of ones) and the contrast vectors associated to the main effects and, optionally, the second-order effects (two-factor interactions and quadratic effects). The main effect contrast vectors are given by the columns in the conference design. The contrast vectors associated to the second-order effects are calculated by taking products between two columns in the conference design.

Model matrices for mixed-level orthogonal arrays

For mixed-level orthogonal arrays, the main effect contrast vectors are defined by the Helmert contrasts. The contrast vectors associated to the two-factor interactions are calculated by taking products between two different columns in the matrix containing the Helmert contrasts of the array; see Model matrices for mixed-level orthogonal arrays for details. The model matrix consists of the intercept column (i.e. a columns of ones) and the contrast vectors associated to the main effects and, optionally, the two-factor interactions.

An example on how to generate an interaction model matrix for a two-level orthogonal array is shown below.

Calculate interaction effects model matrix

>>> array=oapackage.exampleArray(0,1)
exampleArray 0: array in OA(8,2, 2^2)
>>> array.showarray()
array:
  0   0
  0   0
  0   1
  0   1
  1   0
  1   0
  1   1
  1   1
>>> M=oapackage.array2modelmatrix(array, 'i')
>>> print(M)
[[ 1. -1. -1.  1.]
 [ 1. -1. -1.  1.]
 [ 1. -1.  1. -1.]
 [ 1. -1.  1. -1.]
 [ 1.  1. -1. -1.]
 [ 1.  1. -1. -1.]
 [ 1.  1.  1.  1.]
 [ 1.  1.  1.  1.]]

Statistical properties of orthogonal arrays

Orthogonal arrays are commonly evaluated in terms of their generalized wordlength pattern [XW01] (GWLP). Two-level OAs are also commonly evaluated in terms of their \({J}_{k}\)-characteristics and \(F\)-vectors [DT99]. The OApackage can calculate all these statistical criteria: array_link::GWLP(), array_link::Fvalues(), array_link::Jcharacteristics().

The following example shows how to calculate the GWLP, \({F}_{k}\)-values and \({J}_{k}\)-characteristics from an array_link object:

Calculate GWLP and F-values

>>> al=oapackage.exampleArray(1,1) # Select an example array
exampleArray 1: array 3 in OA(16, 2, 2^5)
>>> gwlp = al.GWLP() # Calculate its generalized word length pattern
>>> print('GWLP: %s'% str(gwlp) )
GWLP: (1.0, 0.0, 0.0, 1.0, 1.0, 0.0)
>>> print('F3-value: %s' % str(al.Fvalues(3))) # Calculate the F_3 values
F3-value: (4, 6)
>>> print('F4-value: %s' % str(al.Fvalues(4))) # Calculate the F_3 values
F4-value: (1, 4)
>>> print('J3-characteristics: %s' % str(al.Jcharacteristics(3))) # Calculate the J_3-characteristics
J3-characteristics: (-8, -8, 0, 0, 0, -8, 0, -8, 0, 0)

We now briefly mention some technical details of the \({J}_{k}\)-characteristics, the \({F}_{k}\)-values and the GWLP.

\({J}_{k}\)-characteristics

To calculate \({J}_{k}\)-characteristics of a two-level OA, the OApackage codes the levels of the array as \(-1\) and \(+1\). To this end, the package uses the mapping \({0 \rightarrow -1}\) and \({1 \rightarrow +1}\). Let \(D\) be an \({N}\times {n}\) with coded levels \(-1\) and \(+1\). For \({S} = \{l_1, \ldots, l_k\}\), a subset of \(k\) different factors of \(D = (d_{il})\), define

\[j_k (S; D) = \sum_{i = 1}^{N} d_{i l_1} \cdots d_{i l_k}.\]

The \({|{j}_{k} (S; D)|}\) values are called the \({J}_{k}\)-characteristics, which necessarily equal \(N - 4q\) [DT02], where \({q} \leq N/4\) is a non-negative integer.

\({F}_{k}\)-values

The \({F}_{k}\)-vector collects the frequencies of all the \({J}_{k}\)-characteristics. More specifically, the vector \({F}_{k} = (f_{k1}, \ldots, f_{kv})\), where \(v = N/4\) and \(f_{ku}\) denotes the frequency of the \({J}_{k}\)-characteristics which are equal to \(4(v + 1 - u)\). When calculating an \({F}_{k}\)-vector, the OApackage shows only the vector \((f_{k1}, \ldots, f_{kv})\), whose elements are referred to as the \({F}_{k}\)-values.

Generalized word length pattern

Consider an OA, \({D}\), of strength \({t}\) with \({N}\) runs and \({n}\) factors at \({s}\) levels. Let \({X_0}\) be a column of ones, \({X_1}\) the matrix involving the contrast vectors associated with the main effects, and \({X_j}\) the matrix involving the contrast vectors associated with the \({j}\)-factor interactions, \({j \geq 2}\). We assume that the column vectors in \({X_1}\) are normalized so that they have the same length \({\sqrt{N}}\). For \({j = 0, \ldots, n}\), let

\[A_j (D) = N^{-2} 1_{N}^{T} X_{j} X_{j}^{T} 1_{N}^{\phantom{T}},\]

where \(1_{N}\) denotes the \(N \times 1\) column of ones. The value of \({A}_{j}(D)\) is invariant to the choice of the orthonormal contrasts used; see [XW01] for details. The vector \({(A_0(D), \ldots, A_n (D) )}\) is called the generalized word length pattern (GWLP). To increase the speed of the computations for the GWLP, the OApackage uses the distance distribution and the MacWilliams identities as in [XW01] and [Xu09].

Optimality criteria for D-efficient designs

In [ES17], D-efficient designs for the model including the intercept, all main effects and all two-factor interactions are generated. The OApackage provides functionality to compute the optimality criteria used to generate the D-efficient designs in [ES17]. Moreover, the package can calculate the well-known \(A\)- and \(E\)-optimality criteria from the literature on Optimal Experimental Design [DAT07]. The functions to perform the calulcations are array_link::Defficiency(), array_link::DsEfficiency(), array_link::Aefficiency(), array_link::Eefficiency().

The following example shows how to calculate the \(D\)-, \({D}_{s}\)-, \(A\)- and \(E\)-efficiency for a design that permits the estimation of the interaction model.

Calculate optimality criteria for D-efficient designs

# Select an array that can estimate the interaction model
>>> al = oapackage.exampleArray(11, 1)
exampleArray 11: D-optimal array in OA(44, 2^8)
>>> print('D-efficiency: %.4f' % al.Defficiency())
D-efficiency: 0.8879
>>> print('Ds-efficiency (Eendebak and Schoen, 2017): %.4f' % al.DsEfficiency())
Ds-efficiency (Eendebak and Schoen, 2017): 0.8059
>>> print('A-efficiency for the interaction model: %.4f' % al.Aefficiency())
A-efficiency for the interaction model: 0.7906
>>> print('E-efficiency for the interaction model: %.4f' % al.Eefficiency())
E-efficiency for the interaction model: 0.3602

Calculation of \(D\)-, \(A\)- and \(E\)-efficiency

Let \({X}\) be again the \({N}\times {p}\) interaction model matrix (see section Model matrices) consisting of a column of ones and the contrast vectors associated to the main and two-factor interactions of \({n}\) factors, where \({p = 1 + n + (n)(n-1)/2}\). The \(D\)-, \(A\)- and \(E\)-efficiency are calculated using the eigenvalues of the singular-value decomposition (SVD) of \({X}\). To calculate the rank of a matrix, the lower-upper (LU) decomposition, as implemented in the Eigen package [GJ+10], is used.

Let \(\lambda_1, \ldots, \lambda_p\) be the eigenvalues of the SVD of \({X}\). The OApackage calculates the \(D\)-, \(A\)- and \(E\)-efficiency of a design \(D\) as follows:

\[\begin{split}{D_{\text{eff}}(D)} = (\prod_j \lambda_j)^{1/p} / N \label{formula:Defficiency}, \\ {A_{\text{eff}}(D)} = N (\sum_j \lambda_j^{-1})/m \label{formula:VIF}, \\ {E_{\text{eff}}(D)} = \min_j \lambda_j. \label{formula:E-efficiency}\end{split}\]

\(D_s\)-efficiency and \(D_1\)-efficiency

In [ES17], the \(D_s\)-efficiency is used to assess the joint precision of the main effects in the interaction model. Let the interaction model matrix \({X}\) be split into \({X_{1}}\), containing the contrast vectors associated with the main effects only, and \({X_{02}}\), containing the intercept column and the contrast vectors associated to the two-factor interactions. The \(D_{s}\)-criterion of a design \(D\) is defined as

\[{D_{s,\text{crit}}(D)} = \operatorname{det}(X^{T}X) / \operatorname{det}(X_{02}^{T} X_{02}^{\phantom{T}}), \label{formula:Dsefficiency}\]

where \({X_{02}}\) is necessarily of full rank. Similar to the calculations of the \(D\)-efficiency, the OApackage calculates the \(D_{s}\)-criterion using the eigenvalues of the SVD of the matrices \({X}\) and \({X_{01}}\). Finally, the package calculates the \(D_{s}\)-efficiency of \(D\) as \(D_{s,\text{eff}}(A) = D_{s,\text{crit}}(A)^{1/m}\), where \(m\) is the number of factors.

In a similar way the \(D_{1}\)-efficiency of a design \({A}\) with \(n\) factors and model matrix of intercept and main effects \({X_{01}}\), is defined as

\[D_{1,\text{eff}}(A) = ( \operatorname{det}((X_{01})^{T}(X_{01}) )^{1/(n+1)} \label{formula:D1efficiency}\]

Projection capacities

Other relevant statistical criteria to evaluate a two-level design with \(N\) runs and \(k\) factors include the so-called projection estimation capacity (PEC) and projection information capacity (PIC) [LST07]. These criteria focus on the projections of the two-level design onto a smaller number of factors. More specifically, the PEC and PIC summarize the performance of all the \(N\)-run subdesigns with \(l \leq k\) factors in terms of the capacity to estimate the interaction model and the \(D\)-efficiency for this model, respectively.

The PEC and PIC are based on the so-called PEC and PIC sequences, which are formally defined as follows. Let \(PEC_{l}\) denote the proportion of \(N\)-run \(l\)-factor subdesigns that permit the estimation of the interaction model in \(l\) factors, that is, the model including the intercept, all \(l\) main effects and all \(l(l-1)/2\) two-factor interactions. The PEC sequence is the vector \((PEC_{1}, PEC_{2}, \ldots, PEC_{k})\). Now, let \(PIC_{l}\) denote the average \(D\)-efficiency for the interaction model in \(l\) factors accross all \(N\)-run \(l\)-factor subdesigns. The PIC sequence is the vector \((PIC_{1}, PIC_{2}, \ldots, PIC_{k})\). The OApackage can calculate the PEC and PIC sequences of two-level designs with PECsequence() and PICsequence(), respectively.

The following example shows how to compute the PEC and PIC sequences of a two-level orthogonal array using the OApackage.

Calculate the PEC and PIC sequences

>>> al=oapackage.exampleArray(1,1)
exampleArray 1: array 3 in OA(16, 2, 2^5)
>>> PEC = al.PECsequence()
>>> print('PEC sequence: %s'% ','.join(['%.2f' % x for x in PEC]) )
PEC sequence: 1.00,1.00,1.00,0.80,0.00
>>> PIC = al.PICsequence()
>>> print('PIC sequence: %s'% ','.join(['%.2f' % x for x in PIC]) )
PIC sequence: 1.00,1.00,0.95,0.66,0.00

Properties of conference designs

In [SEG19], it is shown that the \(F_4\) vector is useful for classifying definitive screening designs [XLB12] that are generated by folding over a conference design. To calculate the \(F_4\) vector, we first need to compute the \(J_4\)-characteristics of the conference design. The calculations for the \(J_k\)-characteristics of conference designs are similar as for orthogonal arrays; see Statistical properties of orthogonal arrays. Consider a definitive screening design constructed from an \(N\)-run conference design with at least four factors. The \(F_4\) vector of this design collects the frequencies of the \(J_4\)-characteristics of \({2N − 8\lambda}\) for \({\lambda = 1,\ldots,N/4}\) when \(N\) is a multiple of 4, or \({\lambda = 1, \ldots, (N − 2)/4}\) when \(N\) is an odd multiple of 2.

Note: the \(J_4\)-values of a definitive screening design generated by folding over a conference design are twice the value of the \(J_4\)-values of the conference design. The \(F_4\) vector of a conference design and the corresponding definitive screening design are equal.

Calculate the F4 vector for a conference design

>>> import oapackage
>>> array=oapackage.exampleArray(47,1)
exampleArray 47: third conference design in C(20,8)
>>> F4=array.FvaluesConference(4)
>>> print(F4)
(0, 2, 4, 51, 13)
>>> definitive_screening_design = oapackage.conference2DSD(array)

The individual \(J_k\)-characteristics can be calculated with the method Jcharacteristics_conference(). For conference designs, we can calculate the projection statistics using conferenceProjectionStatistics().

Calculate projection statistics for conference designs

>>> array = oapackage.exampleArray(46, 1)
exampleArray 46: second conference design in C(20,8)
>>> pec, pic, ppc = oapackage.conference.conferenceProjectionStatistics(array)
>>> print('Projection estimation capacity for 4 columns: %.3f'  % pec)
Projection estimation capacity for 4 columns: 0.986
>>> J3 = oapackage.Jcharacteristics_conference(array, number_of_columns = 3)

C++ library

The full documentation for the C++ library can be build from the source code using doxygen. For convenience we provide links the main functions on this page.

Interface for D-optimal designs

Contains functions to generate optimal designs.

For more information see “Two-Level Designs to Estimate All Main Effects and Two-Factor Interactions”, P.T. Eendebak and E.D. Schoen, 2017

Enums

enum coordinate_exchange_method_t

Different methods for the optimization. The default method DOPTIM_SWAP is a coordinate-exchange algorithms.

Values:

enumerator DOPTIM_UPDATE: replace a random element with a random value

enumerator DOPTIM_SWAP: swap two elements at random

enumerator DOPTIM_FLIP: randomly flip an element between 0 and 1

enumerator DOPTIM_AUTOMATIC: automatically select one of the methods

enumerator DOPTIM_NONE: perform no optimization

Functions

double scoreD(const std::vector<double> efficiencies, const std::vector<double> weights)

Calculate score from a set of efficiencies

The score is the weighted sum of the efficiencies.

Parameters

efficiencies – Vector with calculated efficiencies
weights – Weights for the efficiencies

Returns

Weighted sum of the efficiencies

DoptimReturn Doptimize(const arraydata_t &arrayclass, int nrestarts, const std::vector<double> alpha, int verbose, coordinate_exchange_method_t method = DOPTIM_AUTOMATIC, int niter = 300000, double maxtime = 100000, int nabort = 5000)

Generates optimal designs for the specified class of designs

The method uses a coordinate-exchange algorithm to optimze a target function defined by the optimziation paramaters. The optimization is performed multiple times to prevent finding a design in a local minmum of the target function.

The method is described in more detail in “Two-Level Designs to Estimate All Main Effects and Two-Factor Interactions”, Eendebak et al., 2015, Technometrics, https://doi.org/10.1080/00401706.2016.1142903.

Parameters

arrayclass – Class of designs to optimize
nrestarts – Number of restarts to perform
alpha – Optimization parameters. The target function is alpha_1 D + alpha_2 D_s + alpha D_1
verbose – Verbosity level
method – Method for optimization algorithm
niter – Maximum number of iterations for each restart
maxtime – Maximum calculation time. If this time is exceeded, the function is aborted
nabort – Maximum number of iterations when no improvement is found

Returns

A structure with the generated optimal designs

DoptimReturn DoptimizeMixed(const arraylist_t &sols, const arraydata_t &arrayclass, const std::vector<double> alpha, int verbose = 1, int nabort = -1)

Function to generate optimal designs with mixed optimization approach

This function is beta code. See Doptimize for detauls of the parameters.

array_link optimDeff(const array_link &array, const arraydata_t &arrayclass, std::vector<double> alpha, int verbose = 1, coordinate_exchange_method_t optimmethod = DOPTIM_AUTOMATIC, int niter = 100000, int nabort = 0)

Optimize a design according to the optimization function specified.

Arguments:

Parameters

array – Array to be optimized
arrayclass – Structure describing the design class
alpha – 3x1 array with optimization parameters
verbose – Verbosity level
optimmethod – Optimization method to use
niter – Number of iterations
nabort – Number of iterations after which to abort when no improvements are found

Returns

Optimized designs

struct DoptimReturn

#include <Deff.h>

Structure containing results of the Doptimize function

Public Members

std::vector<std::vector<double>> dds: calculated efficiencies for the generated designs

arraylist_t designs: designs generated

int nrestarts: number of restarts performed

int _nimproved

Interface for array properties

Contains functions to calculate properties of arrays.

Author: Pieter Eendebak pieter.eendebak@gmail.com Copyright: See LICENSE.txt file that comes with this distribution

Defines

GWLPvalue

Typedefs

typedef mvalue_t<double> DOFvalue: delete-one-factor projection value

Enums

enum model_matrix_t

Values:

enumerator MODEL_CONSTANT: only the intercept

enumerator MODEL_MAIN: intercept and main effects

enumerator MODEL_INTERACTION: intercept, main effects and second order interactions

enumerator MODEL_SECONDORDER: intercept, main effects and second order effects(interactions and quadratic effects)

enumerator MODEL_INVALID: invalid model

enum paretomethod_t

Values:

enumerator PARETOFUNCTION_DEFAULT

enumerator PARETOFUNCTION_J5

Functions

void DAEefficiencyWithSVD(const Eigen::MatrixXd &secondorder_interaction_matrix, double &Deff, double &vif, double &Eeff, int &rank, int verbose): Calculate D-efficiency and VIF-efficiency and E-efficiency values using SVD.

int array2rank_Deff_Beff(const array_link &al, std::vector<double> *ret = 0, int verbose = 0)

Calculate the rank of the second order interaction matrix of an orthogonal array

The model is the intercept, main effects and interaction effects The rank, D-efficiency, VIF-efficiency and E-efficiency are appended to the second argument

The return vector is filled with the rank, Defficiency, VIF efficiency and Eefficiency

double Defficiency(const array_link &orthogonal_array, int verbose = 0): Calculate D-efficiency for a 2-level array using symmetric eigenvalue decomposition.

std::vector<double> Defficiencies(const array_link &array, const arraydata_t &arrayclass, int verbose = 0, int addDs0 = 0)

Calculate efficiencies for an array

Parameters

array – Array to use in calculation
arrayclass – Specification of the array class
verbose – Verbosity level
addDs0 – If True, then add the Ds0-efficiency to the output

Returns

Vector with the calculate D-efficiency, the main effect robustness (or Ds-optimality) and D1-efficiency for an orthogonal array

double VIFefficiency(const array_link &orthogonal_array, int verbose = 0): Calculate VIF-efficiency of matrix.

double Aefficiency(const array_link &orthogonal_array, int verbose = 0): Calculate A-efficiency of matrix.

double Eefficiency(const array_link &orthogonal_array, int verbose = 0): Calculate E-efficiency of matrix (1 over the VIF-efficiency)

std::vector<double> Aefficiencies(const array_link &orthogonal_array, int verbose = 0): calculate various A-efficiencies

std::vector<double> projDeff(const array_link &array, int number_of_factors, int verbose = 0)

Calculate D-efficiencies for all projection designs

Parameters

array – Design to calculate D-efficiencies for
number_of_factors – Number of factors into which to project
verbose – Verbosity level

Returns

Vector with calculated D-efficiencies

std::vector<double> PECsequence(const array_link &array, int verbose = 0)

Calculate the projection estimation capacity sequence for a design

The PECk of a design is the fraction of estimable second-order models in k factors. The vector (PEC1, PEC2, …, ) is called the projection estimation capacity sequence. See “Ranking Non-regular Designs”, J.L. Loeppky, 2004.

Parameters

array – Input array
verbose – Verbosity level

Returns

Vector with the caculated PEC sequence

std::vector<double> PICsequence(const array_link &array, int verbose = 0)

Calculate the projection information capacity sequence for a design.

The PICk of a design is the average D-efficiency of estimable second-order models in k factors. The vector (PIC1, PIC2, …, ) is called the PIC sequence.

Parameters

array – Input array
verbose – Verbosity level

Returns

Vector with the caculated PIC sequence

ndarray<double> macwilliams_transform_mixed(const ndarray<double> &B, int N, const std::vector<int> &factor_levels_for_groups, int verbose = 0)

Calculate MacWilliams transform.

Calculate MacWilliams transform for mixed level data.

Parameters

B – Input array
N –
verbose – Verbosity level
factor_levels_for_groups – Factor levels for the groups
B – Input array
N – Number of runs
factor_levels_for_groups – Factor levels
verbose – Verbosity level

Returns

MacWilliams transform of the input array

Returns

Transform if the input array

std::vector<double> distance_distribution(const array_link &array)

Return the distance distribution of a design

The distance distribution is described in “Generalized minimum aberration for asymmetrical fractional factorial designs”, Wu and Xu, 2001

Parameters: al – Array for which to calculate the distribution
Returns: Distance distribution

std::vector<int> distance_distribution_shape(const arraydata_t arrayclass)

Return shape of distance distribution for mixed level design

Parameters: arrayclass – Specification of the array class
Returns: Shape of the distance distribution

ndarray<double> distance_distribution_mixed(const array_link &array, int verbose = 0)

Return the distance distribution of a mixed-level design

The distance distribution is described in “Generalized minimum aberration for asymmetrical fractional factorial designs”, Wu and Xu, 2001. For mixed-level designs more details can be found in “A canonical form for non-regular arrays based on generalized wordlength pattern values of delete-one-factor projections”, Eendebak, 2014.

Parameters

al – Array for which to calculate the distribution
verbose – Verbosity level

Returns

Distance distribution

void distance_distribution_mixed_inplace(const array_link &al, ndarray<double> &B, int verbose = 0)

template<class Type> std::vector<double> macwilliams_transform(std::vector<Type> B, int N, int s): Calculate MacWilliams transform.

std::vector<int> Jcharacteristics(const array_link &array, int number_of_columns = 4, int verbose = 0)

Calculate Jk-characteristics of a matrix

The calculated Jk-values are signed.

Parameters

array – Array to calculate Jk-characteristics for
number_of_columns – Number of columns
verbose – Verbosity level

Returns

Vector with calculated Jk-characteristics

std::vector<double> GWLP(const array_link &array, int verbose = 0, int truncate = 1)

Calculate GWLP (generalized wordlength pattern)

The method used for calculation is from Xu and Wu (2001), “Generalized minimum aberration for asymmetrical

fractional factorial desings”. For non-symmetric arrays see “Algorithmic Construction of Efficient Fractional Factorial Designs With Large Run

Sizes”, Xu, Technometrics, 2009.

A more detailed description of the generalized wordlength pattern can also be found in the documentation at https://oapackage.readthedocs.io/.

Parameters

array – Array to calculate the GWLP value for
verbose – Verbosity level
truncate – If True then round values near zero to solve double precision errors

Returns

Vector with calculated generalized wordlength pattern

std::vector<double> GWLPmixed(const array_link &array, int verbose = 0, int truncate = 1)

Calculate GWLP (generalized wordlength pattern) for mixed-level arrays.

The method used for calculation is from “Algorithmic Construction of Efficient Fractional Factorial Designs With Large Run

Sizes”, Xu, Technometrics, 2009.

Parameters

array – Array to calculate the GWLP value for
verbose – Verbosity level
truncate – If True then round values near zero to solve double precision errors

Returns

Vector with calculated generalized wordlength pattern

std::vector<GWLPvalue> projectionGWLPs(const array_link &al): calculate delete-one-factor GWLP (generalized wordlength pattern) projections

std::vector<GWLPvalue> sortGWLP(std::vector<GWLPvalue>): sort a list of GWLP values and return the sorted list

double CL2discrepancy(const array_link &array)

Calculate centered L2-discrepancy of a design

The method is from “A connection between uniformity and aberration in regular fractions of two-level factorials”, Fang and Mukerjee, 2000

array_link array2secondorder(const array_link &array)

Calculate second order interaction model for 2-level array

Parameters: array – Array to calculate second order interaction model from
Returns: Array interaction effects

array_link array2xf(const array_link &array)

calculate second order interaction model for 2-level array

Parameters: array – Array to calculate second order interaction model from
Returns: Array with intercept, main effects and interaction effects

array_link conference_design2modelmatrix(const array_link &conference_design, const char *mode, int verbose = 0)

Calculate model matrix for a conference design

Parameters

conference_design – Conference design
mode – Can be ‘m’ for main effects, ‘i’ for interaction effects or ‘q’ for quadratic effects
verbose – Verbosity level

Returns

Calculated model matrix

Eigen::MatrixXd array2modelmatrix(const array_link &array, const char *mode, int verbose = 0)

Convert orthogonal array or conference design to model matrix

The model matrix consists of the intercept, main effects and (optionally) the interaction effects and quadratic effects. The order in the interaction effects is (c1, c2)=(0,0), (1,0), (2,0), (2,1), … with c2<c1 for columns c1, c2. The size of the model matrix calculated by this function is given by array2modelmatrix_sizes.

For conference designs the method conference_design2modelmatrix is used. For orthogonal array the calculated is performed with array2eigenModelMatrixMixed.

Parameters

array – Orthogonal array or conference design
mode – Type of model matrix to calculate. Can be ‘m’ for main effects, ‘i’ for interaction effects or ‘q’ for quadratic effects
verbose – Verbosity level

Returns

Calculated model matrix

std::vector<int> array2modelmatrix_sizes(const array_link &array)

Return the sizes of the model matrices calculated

Parameters: array – Orthogonal array or conference designs
Returns: List with the sizes of the model matrix for: only intercept; intercept, main; intercept, main, and iteraction terms, intercept, main and full second order

Eigen::MatrixXd array2xfeigen(const array_link &array)

calculate second order interaction model for 2-level array

Parameters: array – Array to calculate second order interaction model from
Returns: Array with intercept, main effects and interaction effects

int arrayrankFullPivQR(const array_link &al, double threshold = -1): return rank of an array based on Eigen::FullPivHouseholderQR

int arrayrankColPivQR(const array_link &al, double threshold = -1): return rank of an array based on Eigen::ColPivHouseholderQR

int arrayrankFullPivLU(const array_link &al, double threshold = -1): return rank of an array based on Eigen::FullPivLU

int arrayrankSVD(const array_link &al, double threshold = -1): return rank of an array based on Eigen::JacobiSVD

int arrayrank(const array_link &array): calculate the rank of an array

int arrayrankInfo(const Eigen::MatrixXd&, int verbose = 1): Return rank of an array. Information about the different methods for rank calculation is printed to stdout.

int arrayrankInfo(const array_link &array, int verbose = 1): Return rank of an array. Information about the different methods for rank calculation is printed to stdout.

Eigen::MatrixXd arraylink2eigen(const array_link &array): convert array_link to Eigen matrix

double conditionNumber(const array_link &matrix): Return the condition number of a matrix.

void calculateParetoEvenOdd(const std::vector<std::string> infiles, const char *outfile, int verbose = 1, arrayfilemode_t afmode = ABINARY, int nrows = -1, int ncols = -1, paretomethod_t paretomethod = PARETOFUNCTION_DEFAULT)

Calculate the Pareto optimal arrays from a list of array files

Pareto optimality is calculated according to (rank; A3,A4; F4)

Pareto<mvalue_t<long>, long> parsePareto(const arraylist_t &arraylist, int verbose, paretomethod_t paretomethod = PARETOFUNCTION_DEFAULT)

mvalue_t<long> A3A4(const array_link &al)

calculate A3 and A4 value for array

Parameters: al – Array for which to calculate A3 and A4
Returns: Object with A3 and A4

inline mvalue_t<long> F4(const array_link &al, int verbose = 1): calculate F4 value for 2-level array

template<class IndexType> Pareto<mvalue_t<long>, IndexType>::pValue calculateArrayParetoRankFA(const array_link &array, int verbose)

Calculate properties of an array and create a Pareto element

The values calculated are:

1) Rank (higher is better) 2) A3, A4 (lower is better) 3) F4 (lower is better, sum of elements is constant)

Valid for 2-level arrays of strength at least 3

template<class IndexType> void addJmax(const array_link &al, typename Pareto<mvalue_t<long>, IndexType>::pValue &p, int verbose = 1): add Jmax criterium to Pareto set

template<class IndexType> Pareto<mvalue_t<long>, IndexType>::pValue calculateArrayParetoJ5(const array_link &al, int verbose): Calculate Pareto element with J5 criterium.

template<class IndexType> inline void parseArrayPareto(const array_link &array, IndexType i, Pareto<mvalue_t<long>, IndexType> &pset, int verbose)

Add array to list of Pareto optimal arrays

The values to be optimized are:

1) Rank (higher is better) 2) A3, A4 (lower is better) 3) F4 (lower is better, sum of elements is constant)

inline double Cvalue2Dvalue(double Cvalue, int number_of_columns): convert C value to D-efficiency value

inline double Dvalue2Cvalue(double Defficiency, int number_of_columns): convert D-efficiency value to C value

template<class Type> class ndarray

#include <arrayproperties.h>

Class representing an n-dimensional array

The data is stored in a flat array. The dimensions are stored in a vector dims.

Public Functions

inline ndarray(const std::vector<int> dims)

Class represensing an n-dimensional array.

Parameters: dims – Dimension of the array

inline ndarray(const ndarray<Type> &rhs): Copy constructor Copies the internal data

inline ~ndarray()

inline void initialize(const Type value): Initialize array with specified value.

inline int sizeof_type() const: Return size of ndarray template type.

inline bool type_is_floating_point() const: Return True is the data type is of floating point type.

inline void info() const

inline std::string idxstring(int linear_idx) const: Convert linear index to string representing the index.

inline long totalsize() const: size of the array (product of all dimensions)

inline void show() const: print the array to stdout

inline void linear2idx(int ndx, int *nidx = 0) const: convert a linear index to normal indices

inline void linear2idx(int ndx, std::vector<int> &nidx) const: convert a linear index to normal indices

inline int getlinearidx(int *idx) const: From an n-dimensional index return the linear index in the data.

inline void *data_pointer() const: Return pointer to data.

inline void setconstant(Type val): set all values of the array to specified value

inline void set(int *idx, Type val): set value at position

inline void setlinear(int idx, Type val): set value using linear index

inline Type getlinear(int idx) const: get value using linear index

inline Type get(int *idx) const: get value using n-dimensional index

Public Members

Type *data

std::vector<int> dims

int k: dimensions of the array

int n

std::vector<int> cumdims

std::vector<int> cumprod

Private Functions

inline void initialize_internal_structures(const std::vector<int> dimsx)

Initialize internal structures

Data pointer is created, but not set with data

Parameters: dimsx – Dimensions of the array

class rankStructure

#include <arrayproperties.h>

Structure to efficiently calculate the rank of the second order interaction matrix of many arrays

The efficiency is obtained if the arrays share a common subarray. The theory is described in “Efficient rank calculation for matrices with a common submatrix”, Eendebak, 2016

Public Types

typedef Eigen::FullPivHouseholderQR<Eigen::MatrixXd> EigenDecomp

Public Functions

inline rankStructure(const array_link &al, int nsub = 3, int verbose = 0): constructor

inline rankStructure(int nsub = 3, int id = -1): constructor

void info() const: print information about the rank structure

void updateStructure(const array_link &al): update the structure cache with a new array

int rankdirect(const Eigen::MatrixXd &array) const: calculate the rank of an array directly, uses special threshold

int rankxfdirect(const array_link &array) const: calculate the rank of the second order interaction matrix of an array directly

int rankxf(const array_link &array): calculate the rank of the second order interaction matrix of an array using the cache system

Public Members

array_link alsub

int r

int verbose: verbosity level

int ks: number of columns of subarray in cache

int nsub: number of columns to subtract from array when updating cache

int id: used for debugging

Private Members

EigenDecomp decomp: decomposition of subarray

Eigen::MatrixXd Qi

long ncalc: internal structure

long nupdate

Interface for array tools

Contains the array_link class and related classes.

This file contains method and classes to work with (orthogonal) arrays.

Author: Pieter Eendebak pieter.eendebak@gmail.com Copyright: See LICENSE.txt file that comes with this distribution

Defines

MPI_ARRAY_T

Typedefs

typedef Eigen::MatrixXd MatrixFloat

typedef Eigen::ArrayXd ArrayFloat

typedef Eigen::VectorXd VectorFloat

typedef double eigenFloat

typedef short int array_t: data type for elements of orthogonal arrays

typedef const short int carray_t: constant version of array_t

typedef short int rowindex_t

typedef int colindex_t: type used for row indexing

typedef const int const_colindex_t: type used for column indexing

typedef array_t *array_p

pointer to array

constant version of type used for column indexing

typedef carray_t *carray_p: pointer to constant array

typedef rowindex_t *rowperm_t

typedef colindex_t *colperm_t: type of row permutation

typedef array_t *levelperm_t: type of column permutation

typedef int vindex_t: type of level permutation

typedef signed char conf_t: data type for elements of conference designs

typedef std::vector<conf_t> conference_column: data type for column of a conference design

typedef std::vector<conference_column> conference_column_list: list of columns of conference designs

typedef std::deque<array_link> arraylist_t: container with arrays

Enums

enum ordering_t

Values:

enumerator ORDER_LEX: lexicograph minimal by columns ordering

enumerator ORDER_J5: J5 based ordering.

Functions

void throw_runtime_exception(const std::string exception_message)

void eigenInfo(const MatrixFloat m, const char *str = "eigen", int verbose = 1)

Print information about an Eigen matrix

Parameters

m – Matrix about which to print information
str – String to prepend in output
verbose – Verbosity level

void print_eigen_matrix(const MatrixFloat matrix): Print Eigen matrix to stdout

void eigen2numpyHelper(double *pymat1, int n, const MatrixFloat &m)

int sizeof_array_t(): return size in bytes of array_t type

int sizeof_double(): return size in bytes of double type

inline std::vector<int> possible_F_values(int N, int strength): possible values for J-values of 2-level design

bool file_exists(const std::string filename): return true if the specified file exists

bool file_exists(const char *filename): return true if the specified file exists

bool oa_file_exists(const char *filename): return true if the specified oa file exists

bool oa_file_exists(const std::string filename): return true if the specified oa file exists

arraydata_t *readConfigFile(const char *file): Read array configuration from file.

inline void copy_array(const array_t *src, array_t *const dst, const int nrows, const int ncols): Make a copy of an array.

inline int destroy_array(array_t *array)

Delete an array.

Parameters: array –
Returns

static inline array_t *create_array(const int nrows, const int ncols)

Create an array.

Parameters

nrows – Number of rows
ncols – Number of columns

Returns

inline array_t *create_array(const arraydata_t *ad): Create an array from an arraydata_t structure.

inline array_t *clone_array(const array_t *const array, const rowindex_t nrows, const colindex_t ncols): Clone an array.

int compareLMC(const array_link &lhs, const array_link &rhs): Return -1 if the first array is smaller in LMC ordering than the second array, 0 if equal and 1 otherwise

array_link exampleArray(int idx = 0, int verbose = 0)

Return example array

Parameters

idx – Index of example array to return
verbose – If True, then print information about the array to stdout

std::vector<int> Jcharacteristics_conference(const array_link &array, int number_of_columns, int verbose = 0)

Calculate Jk-characteristics for a conference design

Parameters

array – Conference design
number_of_columns – Specifies the number of columns to use
verbose – Verbosity level

Returns

A vector of calculated inner products between all combinations of k columns.

array_link hstack(const array_link &array1, const array_link &array2)

concatenate 2 arrays in vertical direction

concatenate 2 arrays in horizontal direction

array_link hstack(const array_link &array, const conference_column &column): concatenate array and conference_column

array_link hstacklastcol(const array_link &A, const array_link &B): concatenate the last column of array B to array A

conference_column vstack(const conference_column &column_top, const conference_column &column_bottom): concatenate two columns

void perform_column_permutation(const array_link source, array_link &target, const std::vector<int> perm): perform column permutation for an array

void perform_row_permutation(const array_link source, array_link &target, const std::vector<int> perm): perform row permutation for an array

arraydata_t arraylink2arraydata(const array_link &array, int extracols = 0, int strength = 2)

create arraydata_t structure from array

Parameters

array – Array to use as input specifiction for array class
extracols – Number of extra columns to add to the number of columns of the array
strength – Strength to set in the array class. If -1, then use the strength of the array

arraylist_t addConstant(const arraylist_t &lst, int value): add a constant value to all arrays in a list

std::vector<int> getJcounts(arraylist_t *arraylist, int N, int k, int verbose = 1): Return number of arrays with j_{2n+1}=0 for number_of_arrays<m

void create_root(array_t *array, const arraydata_t *arrayclass): set first columns of an array to root form

void create_root(const arraydata_t *arrayclass, arraylist_t &solutions): Creates the root of an orthogonal array. The root is appended to the list of arrays.

int array_diff(carray_p A, carray_p B, const rowindex_t r, const colindex_t c, rowindex_t &rpos, colindex_t &cpos): Compare 2 arrays and return position of first difference.

inline void fastJupdate(const array_t *array, rowindex_t N, const int J, const colindex_t *column_indices, array_t *tmp): helper function to calculate J-values

int jvalue(const array_link &array, const int J, const int *column_indices): Calculate J-value for a 2-level array

int jvaluefast(const array_t *array, rowindex_t N, const int J, const colindex_t *column_indices)

Calculate J-value for a column combination of a 2-level array

We assume the array has values 0 and 1. No boundary checks are performed.

std::vector<jstruct_t> analyseArrays(const arraylist_t &arraylist, const int verbose, const int jj = 4): Analyse a list of arrays.

void showArrayList(const arraylist_t &lst): print a list of arrays to stdout

long nArrays(const char *fname): return number of arrays in an array file

void arrayfileinfo(const char *filename, int &number_of_arrays, int &number_of_rows, int &number_of_columns)

return information about file with arrays

Parameters

filename – Filename of array file
number_of_arrays – Variable is set with number of arrays
number_of_rows – Variable is set with number of rows
number_of_columns – Variable is set with number of columns

arraylist_t readarrayfile(const char *fname, int verbose = 1, int *setcols = 0)

Read all arrays in a file

Parameters

fname – Filename to read from
verbose – Verbosity level
setcols – Pointer to return number of columns from array file

Returns

List of arrays

int readarrayfile(const char *filename, arraylist_t *arraylist, int verbose = 1, int *setcols = 0, int *setrows = 0, int *setbits = 0)

Read all arrays in a file and append then to an array list

Parameters

filename – Filename to read from
arraylist – Pointer to list of arrays
verbose – Verbosity level
setcols – Reference that is set with the number of columns from the file
setrows – Reference that is set with the number of rows from the file
setbits – Reference that is set with the number of bits from the file

Returns

int writearrayfile(const char *filename, const arraylist_t &arraylist, arrayfile::arrayfilemode_t mode = arrayfile::ATEXT, int nrows = NRAUTO, int ncols = NRAUTO)

Write a list of arrays to file on disk

Parameters

filename – Filename to use
arraylist – List of arrays to write
mode – Mode for the file with designs
nrows – If the list of arrays is empty, use this number of rows for the design file
ncols – If the list of arrays is empty, use this number of rows for the design file

Returns

Value zero if succesfull

int writearrayfile(const char *filename, const array_link &array, arrayfile::arrayfilemode_t mode = arrayfile::ATEXT): Write a single array to file.

int append_arrayfile(const char *filename, const array_link array): Append a single array to an array file. creates a new file if no file exists.

void selectArrays(const std::string filename, std::vector<int> &idx, arraylist_t &fl, int verbose = 0): Make a selection of arrays from binary array file, append to list.

array_link selectArrays(std::string filename, int index): Select a single array from a file.

arraylist_t selectArrays(const arraylist_t &input_list, std::vector<int> &idx): Make a selection of arrays.

arraylist_t selectArrays(const arraylist_t &input_list, std::vector<long> &idx): Make a selection of arrays.

void selectArrays(const arraylist_t &input_list, std::vector<int> &idx, arraylist_t &output_list): Make a selection of arrays, append to list.

void selectArrays(const arraylist_t &input_list, std::vector<long> &idx, arraylist_t &output_list): Make a selection of arrays, append to list.

template<class Container, class IntType> void keepElements(Container &al, std::vector<IntType> &idx): From a container keep all elements with specified indices.

template<class Container, class IntType> void removeElements(Container &al, std::vector<IntType> &idx): From a container remove all elements with specified indices.

template<class MType> void selectArraysMask(const arraylist_t &al, std::vector<MType> &mask, arraylist_t &rl): Make a selection of arrays from a list, append to list.

template<class IndexType> void appendArrays(const arraylist_t &al, const typename std::vector<IndexType> &idx, arraylist_t &lst): Append selection of arrays to existing list.

void appendArrays(const arraylist_t &arrays_to_append, arraylist_t &dst): Append set of arrays to existing list.

template<class atype> void write_array_format(const atype *array, const int nrows, const int ncols, int width = 3): Write a formatted array

template<class atype> void write_array_format(FILE *fid, const atype *array, const int nrows, const int ncols): Write an array to a file pointer.

template<class atype> void write_array_latex(std::ostream &ss, const atype *array, const int nrows, const int ncols): write an array in latex style

void convert_array_file(std::string input_filename, std::string output_filename, arrayfile::arrayfilemode_t output_format, int verbose = 0): Convert a file with arrays to a different format

bool readbinheader(FILE *fid, int &nr, int &nc)

Read header for binary data file. Return true if valid header file

The header consists of 4 integers: 2 magic numbers, then the number of rows and columns

void writebinheader(FILE *fid, int number_rows, int number_columns): Write header for binary data file.

template<class Type> void vector2doublebinfile(const std::string fname, std::vector<Type> vals, int writeheader = 1): Write a vector of numeric elements to binary file as double values.

void vectorvector2binfile(const std::string fname, const std::vector<std::vector<double>> vals, int writeheader, int na): Write a vector of vector elements to binary file.

MatrixFloat array2eigenX1(const array_link &array, int intercept = 1)

Convert 2-level array to main effects in Eigen format

Parameters

array – Array to convert
intercept – If True, then include the intercept

Returns

The main effects model

MatrixFloat array2eigenX2(const array_link &array)

Convert 2-level array to second order interaction matrix in Eigen format

The intercept and main effects are not included.

Parameters: array – Array to convert
Returns: The second order interaction model

MatrixFloat array2eigenModelMatrix(const array_link &array)

Convert 2-level array to second order interaction model matrix (intercept, main effects, interaction effects)

Parameters: array – Design of which to calculate the model matrix
Returns: Eigen matrix with the model matrix

std::pair<MatrixFloat, MatrixFloat> array2eigenModelMatrixMixed(const array_link &array, int verbose = 1)

Create first and second order model matrix for mixed-level orthogonal array

For 2-level arrays a direct calculation is used. For mixel-level arrays Helmert contrasts are used.

Parameters

array – Input array
verbose – Verbosity level

Returns

Pair with main effects and two-factor interaction model

std::vector<int> numberModelParams(const array_link &array, int order = -1)

Calculate number of parameters in the model matrix

A list of integers is returned, with the number of columns in:

The intercept (always 1)
The main effects
The interaction effects (second order interaction terms without quadratics)
The quadratic effects

Parameters

array – Orthogonal array or conference design
order – Not used any more

Returns

List of sizes

int arrayInFile(const array_link &array, const char *array_file, int verbose = 1)

return index of specified array in a file. returns -1 if array is not found

Parameters

array – Array to find
array_file – Location if file with arrays
verbose – Verbosity level

Returns

Position of array in list

int arrayInList(const array_link &array, const arraylist_t &arrays, int verbose = 1)

return index of specified array in a list. returns -1 if array is not found

Parameters

array – Array to find
arrays – List of arrays
verbose – Verbosity level

Returns

Position of array in list

Variables

const int NRAUTO = 0

struct arraydata_t

#include <arraytools.h>

Specifies a class of arrays.