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

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:: 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

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