# 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 – Number of rows

• factor_levels_for_groups – Factor levels for the groups

• verbose – Verbosity level

• 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

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:

array – 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:
• array – 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

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

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

Calculate second order interaction model for 2-level array

Parameters:

array – Array to calculate second order interaction model from

Returns:

Array interaction effects

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

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

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

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.

convert array_link to Eigen 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)

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

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

calculate the rank of the second order interaction matrix of an array directly

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