From dab15318276dccee228a9f75d75b40f5bb441286 Mon Sep 17 00:00:00 2001 From: RickardSjogren Date: Thu, 22 Sep 2016 07:58:32 +0200 Subject: [PATCH 1/4] Added fractional factorial design specified by resolution. --- pyDOE/doe_factorial.py | 107 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 1 deletion(-) diff --git a/pyDOE/doe_factorial.py b/pyDOE/doe_factorial.py index b49bb3a..6aad2ad 100644 --- a/pyDOE/doe_factorial.py +++ b/pyDOE/doe_factorial.py @@ -14,9 +14,14 @@ """ import re +import string +from itertools import dropwhile, combinations, islice + import numpy as np +from scipy.special import binom + -__all__ = ['np', 'fullfact', 'ff2n', 'fracfact'] +__all__ = ['np', 'fullfact', 'ff2n', 'fracfact', 'fracfact_by_res'] def fullfact(levels): """ @@ -229,6 +234,100 @@ def fracfact(gen): # Return the fractional factorial design return H + + +def fracfact_by_res(n, res): + """ + Create a 2-level fractional factorial design with `n` factors + and resolution `res`. + + Parameters + ---------- + n : int + The number of factors in the design. + res : int + Desired design resolution + + Returns + ------- + H : 2d-array + A m-by-`n` matrix, the fractional factorial design. m is the + minimal amount of rows possible for creating a fractional + factorial design matrix at resolution `res` + + Raises + ------ + ValueError + If the current design is not possible to construct. + + Notes + ----- + The resolution of a design is defined as the length of the shortest + word in the defining relation. The resolution describes the level of + confounding between factors and interaction effects, where higher + resolution indicates lower degree of confounding. + + For example, consider the 2^4-1-design defined by + + gen = "a b c ab" + + The factor "d" is defined by "ab" with defining relation I="abd", where + I is the unit vector. In this simple example the shortest word is "abd" + meaning that this is a resolution III-design. + + In practice resolution III-, IV- and V-designs are most commonly applied. + + * III: Main effects may be confounded with two-factor interactions. + * IV: Main effects are unconfounded by two-factor interactions, but + two-factor interactions may be confounded with each other. + * V: Main effects unconfounded with up to four-factor interactions, + two-factor interactions unconfounded with up to three-factor + interactions. Three-factor interactions may be confounded with + each other. + + Examples + -------- + :: + >>> fracfact_by_res(6, 3) + array([[-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.]]) + + >>> fracfact_by_res(5, 5) + Traceback (most recent call last): + ... + ValueError: design not possible + """ + # Determine minimum required number of base-factors. + min_fac = next(dropwhile(lambda n_: _n_fac_at_res(n_, res) < n, + range(res - 1, n)), None) + + if min_fac is None: + raise ValueError('design not possible') + elif min_fac > len(string.ascii_lowercase): + # This check needs to be done to make sure that the number + # of available are enough since `fracfact` parses design generator + # characters. In practice, this is highly theoretical and it is + # much more likely to run into memory-issues. + raise ValueError('design requires too many base-factors.') + + # Get base factors. + factors = list(string.ascii_lowercase[:min_fac]) + + # Fill out with factor combinations until `n` factors. + factor_combs = (''.join(c) for r in range(res - 1, len(factors)) + for c in combinations(factors, r)) + extra_factors = list(islice(factor_combs, n - len(factors))) + + # Concatenate `gen` string for `fracfact`. + gen = ' '.join(factors + extra_factors) + return fracfact(gen) + def _grep(haystack, needle): try: @@ -241,3 +340,9 @@ def _grep(haystack, needle): if needle in item: locs += [idx] return locs + +def _n_fac_at_res(n, res): + """ Calculate number of possible factors for fractional factorial + design with `n` base factors at resolution `res`. + """ + return sum(binom(n, r) for r in range(res - 1, n)) + n From f3ea4b89fa91897fe44ba8c5fc1e5946805b02b1 Mon Sep 17 00:00:00 2001 From: RickardSjogren Date: Thu, 22 Sep 2016 08:18:20 +0200 Subject: [PATCH 2/4] Added doc-section on fractional factorial design specified by resolution. --- doc/factorial.rst | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/doc/factorial.rst b/doc/factorial.rst index bdc700a..4f4655c 100644 --- a/doc/factorial.rst +++ b/doc/factorial.rst @@ -176,6 +176,27 @@ to the keyword ``columns``:: Care should be taken to decide the appropriate alias structure for your design and the effects that folding has on it. +2-Level Fractional-Factorial specified by resolution (``fracfact_by_res``) +-------------------------------------------------------------------------- + +This function constructs a minimal design at given resolution. It does so +by constructing a generator string with a minimal number of base factors +and passes it to ``fracfact``. This approach favors convenience over +fine-grained control over which factors that are confounded. + +To construct a 6-factor, resolution III-design, ``fractfact_by_res`` +is used like this:: + + >>> fracfact_by_res(6, 3) + array([[-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.]]) + .. index:: Plackett-Burman .. _plackett_burman: From 038efe235edf1e885d70ff1c3e3b28bcb3f9dff8 Mon Sep 17 00:00:00 2001 From: RickardSjogren Date: Fri, 26 Jan 2018 08:36:33 +0100 Subject: [PATCH 3/4] Renamed repo after fork. Prepare for release of fork. --- .gitignore | 1 + README => README.rst | 4 +- doc/conf.py | 6 +- {pyDOE => pyDOE2}/__init__.py | 64 +- {pyDOE => pyDOE2}/build_regression_matrix.py | 194 +++--- {pyDOE => pyDOE2}/doe_box_behnken.py | 182 ++--- {pyDOE => pyDOE2}/doe_composite.py | 320 ++++----- {pyDOE => pyDOE2}/doe_factorial.py | 696 +++++++++---------- {pyDOE => pyDOE2}/doe_fold.py | 122 ++-- {pyDOE => pyDOE2}/doe_lhs.py | 488 ++++++------- {pyDOE => pyDOE2}/doe_plackett_burman.py | 192 ++--- {pyDOE => pyDOE2}/doe_repeat_center.py | 86 +-- {pyDOE => pyDOE2}/doe_star.py | 158 ++--- {pyDOE => pyDOE2}/doe_union.py | 94 +-- {pyDOE => pyDOE2}/var_regression_matrix.py | 102 +-- setup.py | 14 +- 16 files changed, 1362 insertions(+), 1361 deletions(-) rename README => README.rst (96%) rename {pyDOE => pyDOE2}/__init__.py (74%) rename {pyDOE => pyDOE2}/build_regression_matrix.py (96%) rename {pyDOE => pyDOE2}/doe_box_behnken.py (93%) rename {pyDOE => pyDOE2}/doe_composite.py (95%) rename {pyDOE => pyDOE2}/doe_factorial.py (96%) rename {pyDOE => pyDOE2}/doe_fold.py (96%) rename {pyDOE => pyDOE2}/doe_lhs.py (96%) rename {pyDOE => pyDOE2}/doe_plackett_burman.py (96%) rename {pyDOE => pyDOE2}/doe_repeat_center.py (96%) rename {pyDOE => pyDOE2}/doe_star.py (96%) rename {pyDOE => pyDOE2}/doe_union.py (96%) rename {pyDOE => pyDOE2}/var_regression_matrix.py (96%) diff --git a/.gitignore b/.gitignore index 5fa8f85..e6eaf30 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ dist build MANIFEST *.egg-info +.idea diff --git a/README b/README.rst similarity index 96% rename from README rename to README.rst index 95d0b59..beca9a9 100644 --- a/README +++ b/README.rst @@ -1,8 +1,8 @@ ===================================================== -``pyDOE``: The experimental design package for python +``pyDOE2``: The experimental design package for python ===================================================== -The ``pyDOE`` package is designed to help the +The ``pyDOE`` package is designed to help the **scientist, engineer, statistician,** etc., to construct appropriate **experimental designs**. diff --git a/doc/conf.py b/doc/conf.py index 233ae20..d47af67 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -17,7 +17,7 @@ sys.path.insert(0, os.path.abspath('..')) -import pyDOE +import pyDOE2 # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the @@ -43,7 +43,7 @@ master_doc = 'index_TOC' # General information about the project. -project = u'pyDOE' +project = u'pyDOE2' if date.today().year!=2013: copyright = u'2013–%d, Abraham Lee' % date.today().year else: @@ -56,7 +56,7 @@ # The short X.Y version. version = '1' # The full version, including alpha/beta/rc tags. -release = pyDOE.__version__ +release = pyDOE2.__version__ # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/pyDOE/__init__.py b/pyDOE2/__init__.py similarity index 74% rename from pyDOE/__init__.py rename to pyDOE2/__init__.py index 5bbd357..cfd4aec 100644 --- a/pyDOE/__init__.py +++ b/pyDOE2/__init__.py @@ -1,32 +1,32 @@ -""" -================================================================================ -pyDOE: Design of Experiments for Python -================================================================================ - -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -# from __future__ import absolute_import - -__author__ = 'Abraham Lee' -__version__ = '0.3.8' - -from pyDOE.doe_box_behnken import * -from pyDOE.doe_composite import * -from pyDOE.doe_factorial import * -from pyDOE.doe_lhs import * -from pyDOE.doe_fold import * -from pyDOE.doe_plackett_burman import * -from pyDOE.var_regression_matrix import * - +""" +================================================================================ +pyDOE: Design of Experiments for Python +================================================================================ + +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +# from __future__ import absolute_import + +__author__ = 'Abraham Lee' +__version__ = '0.3.8' + +from pyDOE2.doe_box_behnken import * +from pyDOE2.doe_composite import * +from pyDOE2.doe_factorial import * +from pyDOE2.doe_lhs import * +from pyDOE2.doe_fold import * +from pyDOE2.doe_plackett_burman import * +from pyDOE2.var_regression_matrix import * + diff --git a/pyDOE/build_regression_matrix.py b/pyDOE2/build_regression_matrix.py similarity index 96% rename from pyDOE/build_regression_matrix.py rename to pyDOE2/build_regression_matrix.py index cb2fe4d..2bcb4e4 100644 --- a/pyDOE/build_regression_matrix.py +++ b/pyDOE2/build_regression_matrix.py @@ -1,97 +1,97 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np - -def grep(haystack, needle): - start = 0 - while True: - start = haystack.find(needle, start) - if start==-1: - return - yield start - start += len(needle) - -def build_regression_matrix(H, model, build=None): - """ - Build a regression matrix using a DOE matrix and a list of monomials. - - Parameters - ---------- - H : 2d-array - model : str - build : bool-array - - Returns - ------- - R : 2d-array - - """ - ListOfTokens = model.split(' ') - if H.shape[1]==1: - size_index = len(str(H.shape[0])) - else: - size_index = len(str(H.shape[1])) - - if build is None: - build = [True]*len(ListOfTokens) - - # Test if the vector has the wrong direction (lines instead of columns) - if H.shape[0]==1: - H = H.T - - # Collect the list of monomials - Monom_Index = [] - for i in range(len(ListOfTokens)): - if build[i]: - Monom_Index += [grep(ListOfTokens, 'x' + str(0)*(size_index - \ - len(str(i))) + str(i))] - - Monom_Index = -np.sort(-Monom_Index) - Monom_Index = np.unique(Monom_Index) - - if H.shape[1]==1: - nb_var = H.shape[0] # vector "mode": the number of vars is equal to the number of lines of H - VectorMode = True - - for i in range(nb_var): - for j in range(ListOfTokens.shape[0]): - ListOfTokens[j] = ListOfTokens[j].replace( - 'x' + str(0)*(size_index - len(str(i))) + str(i), - 'H(' + str(i) + ')') - else: - nb_var = H.shape[0] # matrix "mode": the number of vars is equal to the number of columns of H - VectorMode = False - - for i in range(nb_var): - for j in range(ListOfTokens.shape[0]): - ListOfTokens[j] = ListOfTokens[j].replace( - 'x' + str(0)*(size_index - len(str(i))) + str(i), - 'H[i,' + str(i) + ')') - - # Now build the regression matrix - if VectorMode: - R = np.zeros((len(ListOfTokens), 1)) - for j in range(len(ListOfTokens)): - R[j, 0] = eval(ListOfTokens[j]) - else: - R = np.zeros((H.shape[0], len(ListOfTokens))) - for i in range(H.shape[0]): - for j in range(len(ListOfTokens)): - R[i, j] = eval(ListOfTokens[j]) - - return R - - +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np + +def grep(haystack, needle): + start = 0 + while True: + start = haystack.find(needle, start) + if start==-1: + return + yield start + start += len(needle) + +def build_regression_matrix(H, model, build=None): + """ + Build a regression matrix using a DOE matrix and a list of monomials. + + Parameters + ---------- + H : 2d-array + model : str + build : bool-array + + Returns + ------- + R : 2d-array + + """ + ListOfTokens = model.split(' ') + if H.shape[1]==1: + size_index = len(str(H.shape[0])) + else: + size_index = len(str(H.shape[1])) + + if build is None: + build = [True]*len(ListOfTokens) + + # Test if the vector has the wrong direction (lines instead of columns) + if H.shape[0]==1: + H = H.T + + # Collect the list of monomials + Monom_Index = [] + for i in range(len(ListOfTokens)): + if build[i]: + Monom_Index += [grep(ListOfTokens, 'x' + str(0)*(size_index - \ + len(str(i))) + str(i))] + + Monom_Index = -np.sort(-Monom_Index) + Monom_Index = np.unique(Monom_Index) + + if H.shape[1]==1: + nb_var = H.shape[0] # vector "mode": the number of vars is equal to the number of lines of H + VectorMode = True + + for i in range(nb_var): + for j in range(ListOfTokens.shape[0]): + ListOfTokens[j] = ListOfTokens[j].replace( + 'x' + str(0)*(size_index - len(str(i))) + str(i), + 'H(' + str(i) + ')') + else: + nb_var = H.shape[0] # matrix "mode": the number of vars is equal to the number of columns of H + VectorMode = False + + for i in range(nb_var): + for j in range(ListOfTokens.shape[0]): + ListOfTokens[j] = ListOfTokens[j].replace( + 'x' + str(0)*(size_index - len(str(i))) + str(i), + 'H[i,' + str(i) + ')') + + # Now build the regression matrix + if VectorMode: + R = np.zeros((len(ListOfTokens), 1)) + for j in range(len(ListOfTokens)): + R[j, 0] = eval(ListOfTokens[j]) + else: + R = np.zeros((H.shape[0], len(ListOfTokens))) + for i in range(H.shape[0]): + for j in range(len(ListOfTokens)): + R[i, j] = eval(ListOfTokens[j]) + + return R + + diff --git a/pyDOE/doe_box_behnken.py b/pyDOE2/doe_box_behnken.py similarity index 93% rename from pyDOE/doe_box_behnken.py rename to pyDOE2/doe_box_behnken.py index b3816e6..86646d8 100644 --- a/pyDOE/doe_box_behnken.py +++ b/pyDOE2/doe_box_behnken.py @@ -1,91 +1,91 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np -from pyDOE.doe_factorial import ff2n -from pyDOE.doe_repeat_center import repeat_center - -__all__ = ['bbdesign'] - -def bbdesign(n, center=None): - """ - Create a Box-Behnken design - - Parameters - ---------- - n : int - The number of factors in the design - - Optional - -------- - center : int - The number of center points to include (default = 1). - - Returns - ------- - mat : 2d-array - The design matrix - - Example - ------- - :: - - >>> bbdesign(3) - array([[-1., -1., 0.], - [ 1., -1., 0.], - [-1., 1., 0.], - [ 1., 1., 0.], - [-1., 0., -1.], - [ 1., 0., -1.], - [-1., 0., 1.], - [ 1., 0., 1.], - [ 0., -1., -1.], - [ 0., 1., -1.], - [ 0., -1., 1.], - [ 0., 1., 1.], - [ 0., 0., 0.], - [ 0., 0., 0.], - [ 0., 0., 0.]]) - - """ - assert n>=3, 'Number of variables must be at least 3' - - # First, compute a factorial DOE with 2 parameters - H_fact = ff2n(2) - # Now we populate the real DOE with this DOE - - # We made a factorial design on each pair of dimensions - # - So, we created a factorial design with two factors - # - Make two loops - Index = 0 - nb_lines = (n*(n-1)/2)*H_fact.shape[0] - H = repeat_center(n, nb_lines) - - for i in range(n - 1): - for j in range(i + 1, n): - Index = Index + 1 - H[max([0, (Index - 1)*H_fact.shape[0]]):Index*H_fact.shape[0], i] = H_fact[:, 0] - H[max([0, (Index - 1)*H_fact.shape[0]]):Index*H_fact.shape[0], j] = H_fact[:, 1] - - if center is None: - if n<=16: - points= [0, 0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] - center = points[n] - else: - center = n - - H = np.c_[H.T, repeat_center(n, center).T].T - - return H +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np +from pyDOE2.doe_factorial import ff2n +from pyDOE2.doe_repeat_center import repeat_center + +__all__ = ['bbdesign'] + +def bbdesign(n, center=None): + """ + Create a Box-Behnken design + + Parameters + ---------- + n : int + The number of factors in the design + + Optional + -------- + center : int + The number of center points to include (default = 1). + + Returns + ------- + mat : 2d-array + The design matrix + + Example + ------- + :: + + >>> bbdesign(3) + array([[-1., -1., 0.], + [ 1., -1., 0.], + [-1., 1., 0.], + [ 1., 1., 0.], + [-1., 0., -1.], + [ 1., 0., -1.], + [-1., 0., 1.], + [ 1., 0., 1.], + [ 0., -1., -1.], + [ 0., 1., -1.], + [ 0., -1., 1.], + [ 0., 1., 1.], + [ 0., 0., 0.], + [ 0., 0., 0.], + [ 0., 0., 0.]]) + + """ + assert n>=3, 'Number of variables must be at least 3' + + # First, compute a factorial DOE with 2 parameters + H_fact = ff2n(2) + # Now we populate the real DOE with this DOE + + # We made a factorial design on each pair of dimensions + # - So, we created a factorial design with two factors + # - Make two loops + Index = 0 + nb_lines = (n*(n-1)/2)*H_fact.shape[0] + H = repeat_center(n, nb_lines) + + for i in range(n - 1): + for j in range(i + 1, n): + Index = Index + 1 + H[max([0, (Index - 1)*H_fact.shape[0]]):Index*H_fact.shape[0], i] = H_fact[:, 0] + H[max([0, (Index - 1)*H_fact.shape[0]]):Index*H_fact.shape[0], j] = H_fact[:, 1] + + if center is None: + if n<=16: + points= [0, 0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] + center = points[n] + else: + center = n + + H = np.c_[H.T, repeat_center(n, center).T].T + + return H diff --git a/pyDOE/doe_composite.py b/pyDOE2/doe_composite.py similarity index 95% rename from pyDOE/doe_composite.py rename to pyDOE2/doe_composite.py index 8d475f2..038b6da 100644 --- a/pyDOE/doe_composite.py +++ b/pyDOE2/doe_composite.py @@ -1,160 +1,160 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np -from pyDOE.doe_factorial import ff2n -from pyDOE.doe_star import star -from pyDOE.doe_union import union -from pyDOE.doe_repeat_center import repeat_center - -__all__ = ['ccdesign'] - -def ccdesign(n, center=(4, 4), alpha='orthogonal', face='circumscribed'): - """ - Central composite design - - Parameters - ---------- - n : int - The number of factors in the design. - - Optional - -------- - center : int array - A 1-by-2 array of integers, the number of center points in each block - of the design. (Default: (4, 4)). - alpha : str - A string describing the effect of alpha has on the variance. ``alpha`` - can take on the following values: - - 1. 'orthogonal' or 'o' (Default) - - 2. 'rotatable' or 'r' - - face : str - The relation between the start points and the corner (factorial) points. - There are three options for this input: - - 1. 'circumscribed' or 'ccc': This is the original form of the central - composite design. The star points are at some distance ``alpha`` - from the center, based on the properties desired for the design. - The start points establish new extremes for the low and high - settings for all factors. These designs have circular, spherical, - or hyperspherical symmetry and require 5 levels for each factor. - Augmenting an existing factorial or resolution V fractional - factorial design with star points can produce this design. - - 2. 'inscribed' or 'cci': For those situations in which the limits - specified for factor settings are truly limits, the CCI design - uses the factors settings as the star points and creates a factorial - or fractional factorial design within those limits (in other words, - a CCI design is a scaled down CCC design with each factor level of - the CCC design divided by ``alpha`` to generate the CCI design). - This design also requires 5 levels of each factor. - - 3. 'faced' or 'ccf': In this design, the star points are at the center - of each face of the factorial space, so ``alpha`` = 1. This - variety requires 3 levels of each factor. Augmenting an existing - factorial or resolution V design with appropriate star points can - also produce this design. - - Notes - ----- - - Fractional factorial designs are not (yet) available here. - - 'ccc' and 'cci' can be rotatable design, but 'ccf' cannot. - - If ``face`` is specified, while ``alpha`` is not, then the default value - of ``alpha`` is 'orthogonal'. - - Returns - ------- - mat : 2d-array - The design matrix with coded levels -1 and 1 - - Example - ------- - :: - - >>> ccdesign(3) - array([[-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. , 0. , 0. ], - [ 0. , 0. , 0. ], - [-1.82574186, 0. , 0. ], - [ 1.82574186, 0. , 0. ], - [ 0. , -1.82574186, 0. ], - [ 0. , 1.82574186, 0. ], - [ 0. , 0. , -1.82574186], - [ 0. , 0. , 1.82574186], - [ 0. , 0. , 0. ], - [ 0. , 0. , 0. ], - [ 0. , 0. , 0. ], - [ 0. , 0. , 0. ]]) - - - """ - # Check inputs - assert isinstance(n, int) and n>1, '"n" must be an integer greater than 1.' - assert alpha.lower() in ('orthogonal', 'o', 'rotatable', - 'r'), 'Invalid value for "alpha": {:}'.format(alpha) - assert face.lower() in ('circumscribed', 'ccc', 'inscribed', 'cci', - 'faced', 'ccf'), 'Invalid value for "face": {:}'.format(face) - - try: - nc = len(center) - except: - raise TypeError('Invalid value for "center": {:}. Expected a 1-by-2 array.'.format(center)) - else: - if nc!=2: - raise ValueError('Invalid number of values for "center" (expected 2, but got {:})'.format(nc)) - - # Orthogonal Design - if alpha.lower() in ('orthogonal', 'o'): - H2, a = star(n, alpha='orthogonal', center=center) - - # Rotatable Design - if alpha.lower() in ('rotatable', 'r'): - H2, a = star(n, alpha='rotatable') - - # Inscribed CCD - if face.lower() in ('inscribed', 'cci'): - H1 = ff2n(n) - H1 = H1/a # Scale down the factorial points - H2, a = star(n) - - # Faced CCD - if face.lower() in ('faced', 'ccf'): - H2, a = star(n) # Value of alpha is always 1 in Faced CCD - H1 = ff2n(n) - - # Circumscribed CCD - if face.lower() in ('circumscribed', 'ccc'): - H1 = ff2n(n) - - C1 = repeat_center(n, center[0]) - C2 = repeat_center(n, center[1]) - - H1 = union(H1, C1) - H2 = union(H2, C2) - H = union(H1, H2) - - return H +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np +from pyDOE2.doe_factorial import ff2n +from pyDOE2.doe_star import star +from pyDOE2.doe_union import union +from pyDOE2.doe_repeat_center import repeat_center + +__all__ = ['ccdesign'] + +def ccdesign(n, center=(4, 4), alpha='orthogonal', face='circumscribed'): + """ + Central composite design + + Parameters + ---------- + n : int + The number of factors in the design. + + Optional + -------- + center : int array + A 1-by-2 array of integers, the number of center points in each block + of the design. (Default: (4, 4)). + alpha : str + A string describing the effect of alpha has on the variance. ``alpha`` + can take on the following values: + + 1. 'orthogonal' or 'o' (Default) + + 2. 'rotatable' or 'r' + + face : str + The relation between the start points and the corner (factorial) points. + There are three options for this input: + + 1. 'circumscribed' or 'ccc': This is the original form of the central + composite design. The star points are at some distance ``alpha`` + from the center, based on the properties desired for the design. + The start points establish new extremes for the low and high + settings for all factors. These designs have circular, spherical, + or hyperspherical symmetry and require 5 levels for each factor. + Augmenting an existing factorial or resolution V fractional + factorial design with star points can produce this design. + + 2. 'inscribed' or 'cci': For those situations in which the limits + specified for factor settings are truly limits, the CCI design + uses the factors settings as the star points and creates a factorial + or fractional factorial design within those limits (in other words, + a CCI design is a scaled down CCC design with each factor level of + the CCC design divided by ``alpha`` to generate the CCI design). + This design also requires 5 levels of each factor. + + 3. 'faced' or 'ccf': In this design, the star points are at the center + of each face of the factorial space, so ``alpha`` = 1. This + variety requires 3 levels of each factor. Augmenting an existing + factorial or resolution V design with appropriate star points can + also produce this design. + + Notes + ----- + - Fractional factorial designs are not (yet) available here. + - 'ccc' and 'cci' can be rotatable design, but 'ccf' cannot. + - If ``face`` is specified, while ``alpha`` is not, then the default value + of ``alpha`` is 'orthogonal'. + + Returns + ------- + mat : 2d-array + The design matrix with coded levels -1 and 1 + + Example + ------- + :: + + >>> ccdesign(3) + array([[-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. , 0. , 0. ], + [ 0. , 0. , 0. ], + [-1.82574186, 0. , 0. ], + [ 1.82574186, 0. , 0. ], + [ 0. , -1.82574186, 0. ], + [ 0. , 1.82574186, 0. ], + [ 0. , 0. , -1.82574186], + [ 0. , 0. , 1.82574186], + [ 0. , 0. , 0. ], + [ 0. , 0. , 0. ], + [ 0. , 0. , 0. ], + [ 0. , 0. , 0. ]]) + + + """ + # Check inputs + assert isinstance(n, int) and n>1, '"n" must be an integer greater than 1.' + assert alpha.lower() in ('orthogonal', 'o', 'rotatable', + 'r'), 'Invalid value for "alpha": {:}'.format(alpha) + assert face.lower() in ('circumscribed', 'ccc', 'inscribed', 'cci', + 'faced', 'ccf'), 'Invalid value for "face": {:}'.format(face) + + try: + nc = len(center) + except: + raise TypeError('Invalid value for "center": {:}. Expected a 1-by-2 array.'.format(center)) + else: + if nc!=2: + raise ValueError('Invalid number of values for "center" (expected 2, but got {:})'.format(nc)) + + # Orthogonal Design + if alpha.lower() in ('orthogonal', 'o'): + H2, a = star(n, alpha='orthogonal', center=center) + + # Rotatable Design + if alpha.lower() in ('rotatable', 'r'): + H2, a = star(n, alpha='rotatable') + + # Inscribed CCD + if face.lower() in ('inscribed', 'cci'): + H1 = ff2n(n) + H1 = H1/a # Scale down the factorial points + H2, a = star(n) + + # Faced CCD + if face.lower() in ('faced', 'ccf'): + H2, a = star(n) # Value of alpha is always 1 in Faced CCD + H1 = ff2n(n) + + # Circumscribed CCD + if face.lower() in ('circumscribed', 'ccc'): + H1 = ff2n(n) + + C1 = repeat_center(n, center[0]) + C2 = repeat_center(n, center[1]) + + H1 = union(H1, C1) + H2 = union(H2, C2) + H = union(H1, H2) + + return H diff --git a/pyDOE/doe_factorial.py b/pyDOE2/doe_factorial.py similarity index 96% rename from pyDOE/doe_factorial.py rename to pyDOE2/doe_factorial.py index 6aad2ad..153e7db 100644 --- a/pyDOE/doe_factorial.py +++ b/pyDOE2/doe_factorial.py @@ -1,348 +1,348 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import re -import string -from itertools import dropwhile, combinations, islice - -import numpy as np -from scipy.special import binom - - -__all__ = ['np', 'fullfact', 'ff2n', 'fracfact', 'fracfact_by_res'] - -def fullfact(levels): - """ - Create a general full-factorial design - - Parameters - ---------- - levels : array-like - An array of integers that indicate the number of levels of each input - design factor. - - Returns - ------- - mat : 2d-array - The design matrix with coded levels 0 to k-1 for a k-level factor - - Example - ------- - :: - - >>> fullfact([2, 4, 3]) - array([[ 0., 0., 0.], - [ 1., 0., 0.], - [ 0., 1., 0.], - [ 1., 1., 0.], - [ 0., 2., 0.], - [ 1., 2., 0.], - [ 0., 3., 0.], - [ 1., 3., 0.], - [ 0., 0., 1.], - [ 1., 0., 1.], - [ 0., 1., 1.], - [ 1., 1., 1.], - [ 0., 2., 1.], - [ 1., 2., 1.], - [ 0., 3., 1.], - [ 1., 3., 1.], - [ 0., 0., 2.], - [ 1., 0., 2.], - [ 0., 1., 2.], - [ 1., 1., 2.], - [ 0., 2., 2.], - [ 1., 2., 2.], - [ 0., 3., 2.], - [ 1., 3., 2.]]) - - """ - n = len(levels) # number of factors - nb_lines = np.prod(levels) # number of trial conditions - H = np.zeros((nb_lines, n)) - - level_repeat = 1 - range_repeat = np.prod(levels) - for i in range(n): - range_repeat /= levels[i] - lvl = [] - for j in range(levels[i]): - lvl += [j]*level_repeat - rng = lvl*range_repeat - level_repeat *= levels[i] - H[:, i] = rng - - return H - -################################################################################ - -def ff2n(n): - """ - Create a 2-Level full-factorial design - - Parameters - ---------- - n : int - The number of factors in the design. - - Returns - ------- - mat : 2d-array - The design matrix with coded levels -1 and 1 - - Example - ------- - :: - - >>> ff2n(3) - array([[-1., -1., -1.], - [ 1., -1., -1.], - [-1., 1., -1.], - [ 1., 1., -1.], - [-1., -1., 1.], - [ 1., -1., 1.], - [-1., 1., 1.], - [ 1., 1., 1.]]) - - """ - return 2*fullfact([2]*n) - 1 - -################################################################################ - -def fracfact(gen): - """ - Create a 2-level fractional-factorial design with a generator string. - - Parameters - ---------- - gen : str - A string, consisting of lowercase, uppercase letters or operators "-" - and "+", indicating the factors of the experiment - - Returns - ------- - H : 2d-array - A m-by-n matrix, the fractional factorial design. m is 2^k, where k - is the number of letters in ``gen``, and n is the total number of - entries in ``gen``. - - Notes - ----- - In ``gen`` we define the main factors of the experiment and the factors - whose levels are the products of the main factors. For example, if - - gen = "a b ab" - - then "a" and "b" are the main factors, while the 3rd factor is the product - of the first two. If we input uppercase letters in ``gen``, we get the same - result. We can also use the operators "+" and "-" in ``gen``. - - For example, if - - gen = "a b -ab" - - then the 3rd factor is the opposite of the product of "a" and "b". - - The output matrix includes the two level full factorial design, built by - the main factors of ``gen``, and the products of the main factors. The - columns of ``H`` follow the sequence of ``gen``. - - For example, if - - gen = "a b ab c" - - then columns H[:, 0], H[:, 1], and H[:, 3] include the two level full - factorial design and H[:, 2] includes the products of the main factors. - - Examples - -------- - :: - - >>> fracfact("a b ab") - array([[-1., -1., 1.], - [ 1., -1., -1.], - [-1., 1., -1.], - [ 1., 1., 1.]]) - - >>> fracfact("A B AB") - array([[-1., -1., 1.], - [ 1., -1., -1.], - [-1., 1., -1.], - [ 1., 1., 1.]]) - - >>> fracfact("a b -ab c +abc") - array([[-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.]]) - - """ - # Recognize letters and combinations - A = [item for item in re.split('\-?\s?\+?', gen) if item] # remove empty strings - C = [len(item) for item in A] - - # Indices of single letters (main factors) - I = [i for i, item in enumerate(C) if item==1] - - # Indices of letter combinations (we need them to fill out H2 properly). - J = [i for i, item in enumerate(C) if item!=1] - - # Check if there are "-" or "+" operators in gen - U = [item for item in gen.split(' ') if item] # remove empty strings - - # If R1 is either None or not, the result is not changed, since it is a - # multiplication of 1. - R1 = _grep(U, '+') - R2 = _grep(U, '-') - - # Fill in design with two level factorial design - H1 = ff2n(len(I)) - H = np.zeros((H1.shape[0], len(C))) - H[:, I] = H1 - - # Recognize combinations and fill in the rest of matrix H2 with the proper - # products - for k in J: - # For lowercase letters - xx = np.array([ord(c) for c in A[k]]) - 97 - - # For uppercase letters - if np.any(xx<0): - xx = np.array([ord(c) for c in A[k]]) - 65 - - H[:, k] = np.prod(H1[:, xx], axis=1) - - # Update design if gen includes "-" operator - if R2: - H[:, R2] *= -1 - - # Return the fractional factorial design - return H - - -def fracfact_by_res(n, res): - """ - Create a 2-level fractional factorial design with `n` factors - and resolution `res`. - - Parameters - ---------- - n : int - The number of factors in the design. - res : int - Desired design resolution - - Returns - ------- - H : 2d-array - A m-by-`n` matrix, the fractional factorial design. m is the - minimal amount of rows possible for creating a fractional - factorial design matrix at resolution `res` - - Raises - ------ - ValueError - If the current design is not possible to construct. - - Notes - ----- - The resolution of a design is defined as the length of the shortest - word in the defining relation. The resolution describes the level of - confounding between factors and interaction effects, where higher - resolution indicates lower degree of confounding. - - For example, consider the 2^4-1-design defined by - - gen = "a b c ab" - - The factor "d" is defined by "ab" with defining relation I="abd", where - I is the unit vector. In this simple example the shortest word is "abd" - meaning that this is a resolution III-design. - - In practice resolution III-, IV- and V-designs are most commonly applied. - - * III: Main effects may be confounded with two-factor interactions. - * IV: Main effects are unconfounded by two-factor interactions, but - two-factor interactions may be confounded with each other. - * V: Main effects unconfounded with up to four-factor interactions, - two-factor interactions unconfounded with up to three-factor - interactions. Three-factor interactions may be confounded with - each other. - - Examples - -------- - :: - >>> fracfact_by_res(6, 3) - array([[-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.]]) - - >>> fracfact_by_res(5, 5) - Traceback (most recent call last): - ... - ValueError: design not possible - """ - # Determine minimum required number of base-factors. - min_fac = next(dropwhile(lambda n_: _n_fac_at_res(n_, res) < n, - range(res - 1, n)), None) - - if min_fac is None: - raise ValueError('design not possible') - elif min_fac > len(string.ascii_lowercase): - # This check needs to be done to make sure that the number - # of available are enough since `fracfact` parses design generator - # characters. In practice, this is highly theoretical and it is - # much more likely to run into memory-issues. - raise ValueError('design requires too many base-factors.') - - # Get base factors. - factors = list(string.ascii_lowercase[:min_fac]) - - # Fill out with factor combinations until `n` factors. - factor_combs = (''.join(c) for r in range(res - 1, len(factors)) - for c in combinations(factors, r)) - extra_factors = list(islice(factor_combs, n - len(factors))) - - # Concatenate `gen` string for `fracfact`. - gen = ' '.join(factors + extra_factors) - return fracfact(gen) - - -def _grep(haystack, needle): - try: - haystack[0] - except (TypeError, AttributeError): - return [0] if needle in haystack else [] - else: - locs = [] - for idx, item in enumerate(haystack): - if needle in item: - locs += [idx] - return locs - -def _n_fac_at_res(n, res): - """ Calculate number of possible factors for fractional factorial - design with `n` base factors at resolution `res`. - """ - return sum(binom(n, r) for r in range(res - 1, n)) + n +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import re +import string +from itertools import dropwhile, combinations, islice + +import numpy as np +from scipy.special import binom + + +__all__ = ['np', 'fullfact', 'ff2n', 'fracfact', 'fracfact_by_res'] + +def fullfact(levels): + """ + Create a general full-factorial design + + Parameters + ---------- + levels : array-like + An array of integers that indicate the number of levels of each input + design factor. + + Returns + ------- + mat : 2d-array + The design matrix with coded levels 0 to k-1 for a k-level factor + + Example + ------- + :: + + >>> fullfact([2, 4, 3]) + array([[ 0., 0., 0.], + [ 1., 0., 0.], + [ 0., 1., 0.], + [ 1., 1., 0.], + [ 0., 2., 0.], + [ 1., 2., 0.], + [ 0., 3., 0.], + [ 1., 3., 0.], + [ 0., 0., 1.], + [ 1., 0., 1.], + [ 0., 1., 1.], + [ 1., 1., 1.], + [ 0., 2., 1.], + [ 1., 2., 1.], + [ 0., 3., 1.], + [ 1., 3., 1.], + [ 0., 0., 2.], + [ 1., 0., 2.], + [ 0., 1., 2.], + [ 1., 1., 2.], + [ 0., 2., 2.], + [ 1., 2., 2.], + [ 0., 3., 2.], + [ 1., 3., 2.]]) + + """ + n = len(levels) # number of factors + nb_lines = np.prod(levels) # number of trial conditions + H = np.zeros((nb_lines, n)) + + level_repeat = 1 + range_repeat = np.prod(levels) + for i in range(n): + range_repeat /= levels[i] + lvl = [] + for j in range(levels[i]): + lvl += [j]*level_repeat + rng = lvl*range_repeat + level_repeat *= levels[i] + H[:, i] = rng + + return H + +################################################################################ + +def ff2n(n): + """ + Create a 2-Level full-factorial design + + Parameters + ---------- + n : int + The number of factors in the design. + + Returns + ------- + mat : 2d-array + The design matrix with coded levels -1 and 1 + + Example + ------- + :: + + >>> ff2n(3) + array([[-1., -1., -1.], + [ 1., -1., -1.], + [-1., 1., -1.], + [ 1., 1., -1.], + [-1., -1., 1.], + [ 1., -1., 1.], + [-1., 1., 1.], + [ 1., 1., 1.]]) + + """ + return 2*fullfact([2]*n) - 1 + +################################################################################ + +def fracfact(gen): + """ + Create a 2-level fractional-factorial design with a generator string. + + Parameters + ---------- + gen : str + A string, consisting of lowercase, uppercase letters or operators "-" + and "+", indicating the factors of the experiment + + Returns + ------- + H : 2d-array + A m-by-n matrix, the fractional factorial design. m is 2^k, where k + is the number of letters in ``gen``, and n is the total number of + entries in ``gen``. + + Notes + ----- + In ``gen`` we define the main factors of the experiment and the factors + whose levels are the products of the main factors. For example, if + + gen = "a b ab" + + then "a" and "b" are the main factors, while the 3rd factor is the product + of the first two. If we input uppercase letters in ``gen``, we get the same + result. We can also use the operators "+" and "-" in ``gen``. + + For example, if + + gen = "a b -ab" + + then the 3rd factor is the opposite of the product of "a" and "b". + + The output matrix includes the two level full factorial design, built by + the main factors of ``gen``, and the products of the main factors. The + columns of ``H`` follow the sequence of ``gen``. + + For example, if + + gen = "a b ab c" + + then columns H[:, 0], H[:, 1], and H[:, 3] include the two level full + factorial design and H[:, 2] includes the products of the main factors. + + Examples + -------- + :: + + >>> fracfact("a b ab") + array([[-1., -1., 1.], + [ 1., -1., -1.], + [-1., 1., -1.], + [ 1., 1., 1.]]) + + >>> fracfact("A B AB") + array([[-1., -1., 1.], + [ 1., -1., -1.], + [-1., 1., -1.], + [ 1., 1., 1.]]) + + >>> fracfact("a b -ab c +abc") + array([[-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.]]) + + """ + # Recognize letters and combinations + A = [item for item in re.split('\-?\s?\+?', gen) if item] # remove empty strings + C = [len(item) for item in A] + + # Indices of single letters (main factors) + I = [i for i, item in enumerate(C) if item==1] + + # Indices of letter combinations (we need them to fill out H2 properly). + J = [i for i, item in enumerate(C) if item!=1] + + # Check if there are "-" or "+" operators in gen + U = [item for item in gen.split(' ') if item] # remove empty strings + + # If R1 is either None or not, the result is not changed, since it is a + # multiplication of 1. + R1 = _grep(U, '+') + R2 = _grep(U, '-') + + # Fill in design with two level factorial design + H1 = ff2n(len(I)) + H = np.zeros((H1.shape[0], len(C))) + H[:, I] = H1 + + # Recognize combinations and fill in the rest of matrix H2 with the proper + # products + for k in J: + # For lowercase letters + xx = np.array([ord(c) for c in A[k]]) - 97 + + # For uppercase letters + if np.any(xx<0): + xx = np.array([ord(c) for c in A[k]]) - 65 + + H[:, k] = np.prod(H1[:, xx], axis=1) + + # Update design if gen includes "-" operator + if R2: + H[:, R2] *= -1 + + # Return the fractional factorial design + return H + + +def fracfact_by_res(n, res): + """ + Create a 2-level fractional factorial design with `n` factors + and resolution `res`. + + Parameters + ---------- + n : int + The number of factors in the design. + res : int + Desired design resolution + + Returns + ------- + H : 2d-array + A m-by-`n` matrix, the fractional factorial design. m is the + minimal amount of rows possible for creating a fractional + factorial design matrix at resolution `res` + + Raises + ------ + ValueError + If the current design is not possible to construct. + + Notes + ----- + The resolution of a design is defined as the length of the shortest + word in the defining relation. The resolution describes the level of + confounding between factors and interaction effects, where higher + resolution indicates lower degree of confounding. + + For example, consider the 2^4-1-design defined by + + gen = "a b c ab" + + The factor "d" is defined by "ab" with defining relation I="abd", where + I is the unit vector. In this simple example the shortest word is "abd" + meaning that this is a resolution III-design. + + In practice resolution III-, IV- and V-designs are most commonly applied. + + * III: Main effects may be confounded with two-factor interactions. + * IV: Main effects are unconfounded by two-factor interactions, but + two-factor interactions may be confounded with each other. + * V: Main effects unconfounded with up to four-factor interactions, + two-factor interactions unconfounded with up to three-factor + interactions. Three-factor interactions may be confounded with + each other. + + Examples + -------- + :: + >>> fracfact_by_res(6, 3) + array([[-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.]]) + + >>> fracfact_by_res(5, 5) + Traceback (most recent call last): + ... + ValueError: design not possible + """ + # Determine minimum required number of base-factors. + min_fac = next(dropwhile(lambda n_: _n_fac_at_res(n_, res) < n, + range(res - 1, n)), None) + + if min_fac is None: + raise ValueError('design not possible') + elif min_fac > len(string.ascii_lowercase): + # This check needs to be done to make sure that the number + # of available are enough since `fracfact` parses design generator + # characters. In practice, this is highly theoretical and it is + # much more likely to run into memory-issues. + raise ValueError('design requires too many base-factors.') + + # Get base factors. + factors = list(string.ascii_lowercase[:min_fac]) + + # Fill out with factor combinations until `n` factors. + factor_combs = (''.join(c) for r in range(res - 1, len(factors)) + for c in combinations(factors, r)) + extra_factors = list(islice(factor_combs, n - len(factors))) + + # Concatenate `gen` string for `fracfact`. + gen = ' '.join(factors + extra_factors) + return fracfact(gen) + + +def _grep(haystack, needle): + try: + haystack[0] + except (TypeError, AttributeError): + return [0] if needle in haystack else [] + else: + locs = [] + for idx, item in enumerate(haystack): + if needle in item: + locs += [idx] + return locs + +def _n_fac_at_res(n, res): + """ Calculate number of possible factors for fractional factorial + design with `n` base factors at resolution `res`. + """ + return sum(binom(n, r) for r in range(res - 1, n)) + n diff --git a/pyDOE/doe_fold.py b/pyDOE2/doe_fold.py similarity index 96% rename from pyDOE/doe_fold.py rename to pyDOE2/doe_fold.py index 27054e0..b8d64ac 100644 --- a/pyDOE/doe_fold.py +++ b/pyDOE2/doe_fold.py @@ -1,61 +1,61 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np - -__all__ = ['fold'] - -def fold(H, columns=None): - """ - Fold a design to reduce confounding effects. - - Parameters - ---------- - H : 2d-array - The design matrix to be folded. - columns : array - Indices of of columns to fold (Default: None). If ``columns=None`` is - used, then all columns will be folded. - - Returns - ------- - Hf : 2d-array - The folded design matrix. - - Examples - -------- - :: - - """ - H = np.array(H) - assert len(H.shape)==2, 'Input design matrix must be 2d.' - - if columns is None: - columns = range(H.shape[1]) - - Hf = H.copy() - - for col in columns: - vals = np.unique(H[:, col]) - assert len(vals)==2, 'Input design matrix must be 2-level factors only.' - - for i in range(H.shape[0]): - Hf[i, col] = vals[0] if H[i, col]==vals[1] else vals[1] - - Hf = np.vstack((H, Hf)) - - return Hf - - +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np + +__all__ = ['fold'] + +def fold(H, columns=None): + """ + Fold a design to reduce confounding effects. + + Parameters + ---------- + H : 2d-array + The design matrix to be folded. + columns : array + Indices of of columns to fold (Default: None). If ``columns=None`` is + used, then all columns will be folded. + + Returns + ------- + Hf : 2d-array + The folded design matrix. + + Examples + -------- + :: + + """ + H = np.array(H) + assert len(H.shape)==2, 'Input design matrix must be 2d.' + + if columns is None: + columns = range(H.shape[1]) + + Hf = H.copy() + + for col in columns: + vals = np.unique(H[:, col]) + assert len(vals)==2, 'Input design matrix must be 2-level factors only.' + + for i in range(H.shape[0]): + Hf[i, col] = vals[0] if H[i, col]==vals[1] else vals[1] + + Hf = np.vstack((H, Hf)) + + return Hf + + diff --git a/pyDOE/doe_lhs.py b/pyDOE2/doe_lhs.py similarity index 96% rename from pyDOE/doe_lhs.py rename to pyDOE2/doe_lhs.py index 3e3dccd..a293673 100644 --- a/pyDOE/doe_lhs.py +++ b/pyDOE2/doe_lhs.py @@ -1,244 +1,244 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np -from math import factorial - -__all__ = ['lhs'] - -def lhs(n, samples=None, criterion=None, iterations=None): - """ - Generate a latin-hypercube design - - Parameters - ---------- - n : int - The number of factors to generate samples for - - Optional - -------- - samples : int - The number of samples to generate for each factor (Default: n) - criterion : str - Allowable values are "center" or "c", "maximin" or "m", - "centermaximin" or "cm", and "correlation" or "corr". If no value - given, the design is simply randomized. - iterations : int - The number of iterations in the maximin and correlations algorithms - (Default: 5). - - Returns - ------- - H : 2d-array - An n-by-samples design matrix that has been normalized so factor values - are uniformly spaced between zero and one. - - Example - ------- - A 3-factor design (defaults to 3 samples):: - - >>> lhs(3) - array([[ 0.40069325, 0.08118402, 0.69763298], - [ 0.19524568, 0.41383587, 0.29947106], - [ 0.85341601, 0.75460699, 0.360024 ]]) - - A 4-factor design with 6 samples:: - - >>> lhs(4, samples=6) - array([[ 0.27226812, 0.02811327, 0.62792445, 0.91988196], - [ 0.76945538, 0.43501682, 0.01107457, 0.09583358], - [ 0.45702981, 0.76073773, 0.90245401, 0.18773015], - [ 0.99342115, 0.85814198, 0.16996665, 0.65069309], - [ 0.63092013, 0.22148567, 0.33616859, 0.36332478], - [ 0.05276917, 0.5819198 , 0.67194243, 0.78703262]]) - - A 2-factor design with 5 centered samples:: - - >>> lhs(2, samples=5, criterion='center') - array([[ 0.3, 0.5], - [ 0.7, 0.9], - [ 0.1, 0.3], - [ 0.9, 0.1], - [ 0.5, 0.7]]) - - A 3-factor design with 4 samples where the minimum distance between - all samples has been maximized:: - - >>> lhs(3, samples=4, criterion='maximin') - array([[ 0.02642564, 0.55576963, 0.50261649], - [ 0.51606589, 0.88933259, 0.34040838], - [ 0.98431735, 0.0380364 , 0.01621717], - [ 0.40414671, 0.33339132, 0.84845707]]) - - A 4-factor design with 5 samples where the samples are as uncorrelated - as possible (within 10 iterations):: - - >>> lhs(4, samples=5, criterion='correlate', iterations=10) - - """ - H = None - - if samples is None: - samples = n - - if criterion is not None: - assert criterion.lower() in ('center', 'c', 'maximin', 'm', - 'centermaximin', 'cm', 'correlation', - 'corr'), 'Invalid value for "criterion": {}'.format(criterion) - else: - H = _lhsclassic(n, samples) - - if criterion is None: - criterion = 'center' - - if iterations is None: - iterations = 5 - - if H is None: - if criterion.lower() in ('center', 'c'): - H = _lhscentered(n, samples) - elif criterion.lower() in ('maximin', 'm'): - H = _lhsmaximin(n, samples, iterations, 'maximin') - elif criterion.lower() in ('centermaximin', 'cm'): - H = _lhsmaximin(n, samples, iterations, 'centermaximin') - elif criterion.lower() in ('correlate', 'corr'): - H = _lhscorrelate(n, samples, iterations) - - return H - -################################################################################ - -def _lhsclassic(n, samples): - # Generate the intervals - cut = np.linspace(0, 1, samples + 1) - - # Fill points uniformly in each interval - u = np.random.rand(samples, n) - a = cut[:samples] - b = cut[1:samples + 1] - rdpoints = np.zeros_like(u) - for j in range(n): - rdpoints[:, j] = u[:, j]*(b-a) + a - - # Make the random pairings - H = np.zeros_like(rdpoints) - for j in range(n): - order = np.random.permutation(range(samples)) - H[:, j] = rdpoints[order, j] - - return H - -################################################################################ - -def _lhscentered(n, samples): - # Generate the intervals - cut = np.linspace(0, 1, samples + 1) - - # Fill points uniformly in each interval - u = np.random.rand(samples, n) - a = cut[:samples] - b = cut[1:samples + 1] - _center = (a + b)/2 - - # Make the random pairings - H = np.zeros_like(u) - for j in range(n): - H[:, j] = np.random.permutation(_center) - - return H - -################################################################################ - -def _lhsmaximin(n, samples, iterations, lhstype): - maxdist = 0 - - # Maximize the minimum distance between points - for i in range(iterations): - if lhstype=='maximin': - Hcandidate = _lhsclassic(n, samples) - else: - Hcandidate = _lhscentered(n, samples) - - d = _pdist(Hcandidate) - if maxdist>> x = np.array([[0.1629447, 0.8616334], - ... [0.5811584, 0.3826752], - ... [0.2270954, 0.4442068], - ... [0.7670017, 0.7264718], - ... [0.8253975, 0.1937736]]) - >>> _pdist(x) - array([ 0.6358488, 0.4223272, 0.6189940, 0.9406808, 0.3593699, - 0.3908118, 0.3087661, 0.6092392, 0.6486001, 0.5358894]) - - """ - - x = np.atleast_2d(x) - assert len(x.shape)==2, 'Input array must be 2d-dimensional' - - m, n = x.shape - if m<2: - return [] - - d = [] - for i in range(m - 1): - for j in range(i + 1, m): - d.append((sum((x[j, :] - x[i, :])**2))**0.5) - - return np.array(d) - +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np +from math import factorial + +__all__ = ['lhs'] + +def lhs(n, samples=None, criterion=None, iterations=None): + """ + Generate a latin-hypercube design + + Parameters + ---------- + n : int + The number of factors to generate samples for + + Optional + -------- + samples : int + The number of samples to generate for each factor (Default: n) + criterion : str + Allowable values are "center" or "c", "maximin" or "m", + "centermaximin" or "cm", and "correlation" or "corr". If no value + given, the design is simply randomized. + iterations : int + The number of iterations in the maximin and correlations algorithms + (Default: 5). + + Returns + ------- + H : 2d-array + An n-by-samples design matrix that has been normalized so factor values + are uniformly spaced between zero and one. + + Example + ------- + A 3-factor design (defaults to 3 samples):: + + >>> lhs(3) + array([[ 0.40069325, 0.08118402, 0.69763298], + [ 0.19524568, 0.41383587, 0.29947106], + [ 0.85341601, 0.75460699, 0.360024 ]]) + + A 4-factor design with 6 samples:: + + >>> lhs(4, samples=6) + array([[ 0.27226812, 0.02811327, 0.62792445, 0.91988196], + [ 0.76945538, 0.43501682, 0.01107457, 0.09583358], + [ 0.45702981, 0.76073773, 0.90245401, 0.18773015], + [ 0.99342115, 0.85814198, 0.16996665, 0.65069309], + [ 0.63092013, 0.22148567, 0.33616859, 0.36332478], + [ 0.05276917, 0.5819198 , 0.67194243, 0.78703262]]) + + A 2-factor design with 5 centered samples:: + + >>> lhs(2, samples=5, criterion='center') + array([[ 0.3, 0.5], + [ 0.7, 0.9], + [ 0.1, 0.3], + [ 0.9, 0.1], + [ 0.5, 0.7]]) + + A 3-factor design with 4 samples where the minimum distance between + all samples has been maximized:: + + >>> lhs(3, samples=4, criterion='maximin') + array([[ 0.02642564, 0.55576963, 0.50261649], + [ 0.51606589, 0.88933259, 0.34040838], + [ 0.98431735, 0.0380364 , 0.01621717], + [ 0.40414671, 0.33339132, 0.84845707]]) + + A 4-factor design with 5 samples where the samples are as uncorrelated + as possible (within 10 iterations):: + + >>> lhs(4, samples=5, criterion='correlate', iterations=10) + + """ + H = None + + if samples is None: + samples = n + + if criterion is not None: + assert criterion.lower() in ('center', 'c', 'maximin', 'm', + 'centermaximin', 'cm', 'correlation', + 'corr'), 'Invalid value for "criterion": {}'.format(criterion) + else: + H = _lhsclassic(n, samples) + + if criterion is None: + criterion = 'center' + + if iterations is None: + iterations = 5 + + if H is None: + if criterion.lower() in ('center', 'c'): + H = _lhscentered(n, samples) + elif criterion.lower() in ('maximin', 'm'): + H = _lhsmaximin(n, samples, iterations, 'maximin') + elif criterion.lower() in ('centermaximin', 'cm'): + H = _lhsmaximin(n, samples, iterations, 'centermaximin') + elif criterion.lower() in ('correlate', 'corr'): + H = _lhscorrelate(n, samples, iterations) + + return H + +################################################################################ + +def _lhsclassic(n, samples): + # Generate the intervals + cut = np.linspace(0, 1, samples + 1) + + # Fill points uniformly in each interval + u = np.random.rand(samples, n) + a = cut[:samples] + b = cut[1:samples + 1] + rdpoints = np.zeros_like(u) + for j in range(n): + rdpoints[:, j] = u[:, j]*(b-a) + a + + # Make the random pairings + H = np.zeros_like(rdpoints) + for j in range(n): + order = np.random.permutation(range(samples)) + H[:, j] = rdpoints[order, j] + + return H + +################################################################################ + +def _lhscentered(n, samples): + # Generate the intervals + cut = np.linspace(0, 1, samples + 1) + + # Fill points uniformly in each interval + u = np.random.rand(samples, n) + a = cut[:samples] + b = cut[1:samples + 1] + _center = (a + b)/2 + + # Make the random pairings + H = np.zeros_like(u) + for j in range(n): + H[:, j] = np.random.permutation(_center) + + return H + +################################################################################ + +def _lhsmaximin(n, samples, iterations, lhstype): + maxdist = 0 + + # Maximize the minimum distance between points + for i in range(iterations): + if lhstype=='maximin': + Hcandidate = _lhsclassic(n, samples) + else: + Hcandidate = _lhscentered(n, samples) + + d = _pdist(Hcandidate) + if maxdist>> x = np.array([[0.1629447, 0.8616334], + ... [0.5811584, 0.3826752], + ... [0.2270954, 0.4442068], + ... [0.7670017, 0.7264718], + ... [0.8253975, 0.1937736]]) + >>> _pdist(x) + array([ 0.6358488, 0.4223272, 0.6189940, 0.9406808, 0.3593699, + 0.3908118, 0.3087661, 0.6092392, 0.6486001, 0.5358894]) + + """ + + x = np.atleast_2d(x) + assert len(x.shape)==2, 'Input array must be 2d-dimensional' + + m, n = x.shape + if m<2: + return [] + + d = [] + for i in range(m - 1): + for j in range(i + 1, m): + d.append((sum((x[j, :] - x[i, :])**2))**0.5) + + return np.array(d) + diff --git a/pyDOE/doe_plackett_burman.py b/pyDOE2/doe_plackett_burman.py similarity index 96% rename from pyDOE/doe_plackett_burman.py rename to pyDOE2/doe_plackett_burman.py index fb24f2d..2c3e43b 100644 --- a/pyDOE/doe_plackett_burman.py +++ b/pyDOE2/doe_plackett_burman.py @@ -1,96 +1,96 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import math -import numpy as np -from scipy.linalg import toeplitz, hankel - -__all__ = ['pbdesign'] - -def pbdesign(n): - """ - Generate a Plackett-Burman design - - Parameter - --------- - n : int - The number of factors to create a matrix for. - - Returns - ------- - H : 2d-array - An orthogonal design matrix with n columns, one for each factor, and - the number of rows being the next multiple of 4 higher than n (e.g., - for 1-3 factors there are 4 rows, for 4-7 factors there are 8 rows, - etc.) - - Example - ------- - - A 3-factor design:: - - >>> pbdesign(3) - array([[-1., -1., 1.], - [ 1., -1., -1.], - [-1., 1., -1.], - [ 1., 1., 1.]]) - - A 5-factor design:: - - >>> pbdesign(5) - array([[-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.]]) - - """ - assert n>0, 'Number of factors must be a positive integer' - keep = int(n) - n = 4*(int(n/4) + 1) # calculate the correct number of rows (multiple of 4) - f, e = np.frexp([n, n/12., n/20.]) - k = [idx for idx, val in enumerate(np.logical_and(f==0.5, e>0)) if val] - - assert isinstance(n, int) and k!=[], 'Invalid inputs. n must be a multiple of 4.' - - k = k[0] - e = e[k] - 1 - - if k==0: # N = 1*2**e - H = np.ones((1, 1)) - elif k==1: # N = 12*2**e - H = np.vstack((np.ones((1, 12)), np.hstack((np.ones((11, 1)), - toeplitz([-1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1], - [-1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1]))))) - elif k==2: # N = 20*2**e - H = np.vstack((np.ones((1, 20)), np.hstack((np.ones((19, 1)), - hankel( - [-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]) - )))) - - # Kronecker product construction - for i in range(e): - H = np.vstack((np.hstack((H, H)), np.hstack((H, -H)))) - - # Reduce the size of the matrix as needed - H = H[:, 1:(keep + 1)] - - return np.flipud(H) - - +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import math +import numpy as np +from scipy.linalg import toeplitz, hankel + +__all__ = ['pbdesign'] + +def pbdesign(n): + """ + Generate a Plackett-Burman design + + Parameter + --------- + n : int + The number of factors to create a matrix for. + + Returns + ------- + H : 2d-array + An orthogonal design matrix with n columns, one for each factor, and + the number of rows being the next multiple of 4 higher than n (e.g., + for 1-3 factors there are 4 rows, for 4-7 factors there are 8 rows, + etc.) + + Example + ------- + + A 3-factor design:: + + >>> pbdesign(3) + array([[-1., -1., 1.], + [ 1., -1., -1.], + [-1., 1., -1.], + [ 1., 1., 1.]]) + + A 5-factor design:: + + >>> pbdesign(5) + array([[-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.]]) + + """ + assert n>0, 'Number of factors must be a positive integer' + keep = int(n) + n = 4*(int(n/4) + 1) # calculate the correct number of rows (multiple of 4) + f, e = np.frexp([n, n/12., n/20.]) + k = [idx for idx, val in enumerate(np.logical_and(f==0.5, e>0)) if val] + + assert isinstance(n, int) and k!=[], 'Invalid inputs. n must be a multiple of 4.' + + k = k[0] + e = e[k] - 1 + + if k==0: # N = 1*2**e + H = np.ones((1, 1)) + elif k==1: # N = 12*2**e + H = np.vstack((np.ones((1, 12)), np.hstack((np.ones((11, 1)), + toeplitz([-1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1], + [-1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1]))))) + elif k==2: # N = 20*2**e + H = np.vstack((np.ones((1, 20)), np.hstack((np.ones((19, 1)), + hankel( + [-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]) + )))) + + # Kronecker product construction + for i in range(e): + H = np.vstack((np.hstack((H, H)), np.hstack((H, -H)))) + + # Reduce the size of the matrix as needed + H = H[:, 1:(keep + 1)] + + return np.flipud(H) + + diff --git a/pyDOE/doe_repeat_center.py b/pyDOE2/doe_repeat_center.py similarity index 96% rename from pyDOE/doe_repeat_center.py rename to pyDOE2/doe_repeat_center.py index d6d7bc4..c69b314 100644 --- a/pyDOE/doe_repeat_center.py +++ b/pyDOE2/doe_repeat_center.py @@ -1,43 +1,43 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np - -def repeat_center(n, repeat): - """ - Create the center-point portion of a design matrix - - Parameters - ---------- - n : int - The number of factors in the original design - repeat : int - The number of center points to repeat - - Returns - ------- - mat : 2d-array - The center-point portion of a design matrix (elements all zero). - - Example - ------- - :: - - >>> repeat_center(3, 2) - array([[ 0., 0., 0.], - [ 0., 0., 0.]]) - - """ - return np.zeros((repeat, n)) +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np + +def repeat_center(n, repeat): + """ + Create the center-point portion of a design matrix + + Parameters + ---------- + n : int + The number of factors in the original design + repeat : int + The number of center points to repeat + + Returns + ------- + mat : 2d-array + The center-point portion of a design matrix (elements all zero). + + Example + ------- + :: + + >>> repeat_center(3, 2) + array([[ 0., 0., 0.], + [ 0., 0., 0.]]) + + """ + return np.zeros((repeat, n)) diff --git a/pyDOE/doe_star.py b/pyDOE2/doe_star.py similarity index 96% rename from pyDOE/doe_star.py rename to pyDOE2/doe_star.py index ba72920..788a1ec 100644 --- a/pyDOE/doe_star.py +++ b/pyDOE2/doe_star.py @@ -1,79 +1,79 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np - -def star(n, alpha='faced', center=(1, 1)): - """ - Create the star points of various design matrices - - Parameters - ---------- - n : int - The number of variables in the design - - Optional - -------- - alpha : str - Available values are 'faced' (default), 'orthogonal', or 'rotatable' - center : array - A 1-by-2 array of integers indicating the number of center points - assigned in each block of the response surface design. Default is - (1, 1). - - Returns - ------- - H : 2d-array - The star-point portion of the design matrix (i.e. at +/- alpha) - a : scalar - The alpha value to scale the star points with. - - Example - ------- - :: - - >>> star(3) - array([[-1., 0., 0.], - [ 1., 0., 0.], - [ 0., -1., 0.], - [ 0., 1., 0.], - [ 0., 0., -1.], - [ 0., 0., 1.]]) - - """ - # Star points at the center of each face of the factorial - if alpha=='faced': - a = 1 - elif alpha=='orthogonal': - nc = 2**n # factorial points - nco = center[0] # center points to factorial - na = 2*n # axial points - nao = center[1] # center points to axial design - # value of alpha in orthogonal design - a = (n*(1 + nao/float(na))/(1 + nco/float(nc)))**0.5 - elif alpha=='rotatable': - nc = 2**n # number of factorial points - a = nc**(0.25) # value of alpha in rotatable design - else: - raise ValueError('Invalid value for "alpha": {:}'.format(alpha)) - - # Create the actual matrix now. - H = np.zeros((2*n, n)) - for i in range(n): - H[2*i:2*i+2, i] = [-1, 1] - - H *= a - - return H, a +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np + +def star(n, alpha='faced', center=(1, 1)): + """ + Create the star points of various design matrices + + Parameters + ---------- + n : int + The number of variables in the design + + Optional + -------- + alpha : str + Available values are 'faced' (default), 'orthogonal', or 'rotatable' + center : array + A 1-by-2 array of integers indicating the number of center points + assigned in each block of the response surface design. Default is + (1, 1). + + Returns + ------- + H : 2d-array + The star-point portion of the design matrix (i.e. at +/- alpha) + a : scalar + The alpha value to scale the star points with. + + Example + ------- + :: + + >>> star(3) + array([[-1., 0., 0.], + [ 1., 0., 0.], + [ 0., -1., 0.], + [ 0., 1., 0.], + [ 0., 0., -1.], + [ 0., 0., 1.]]) + + """ + # Star points at the center of each face of the factorial + if alpha=='faced': + a = 1 + elif alpha=='orthogonal': + nc = 2**n # factorial points + nco = center[0] # center points to factorial + na = 2*n # axial points + nao = center[1] # center points to axial design + # value of alpha in orthogonal design + a = (n*(1 + nao/float(na))/(1 + nco/float(nc)))**0.5 + elif alpha=='rotatable': + nc = 2**n # number of factorial points + a = nc**(0.25) # value of alpha in rotatable design + else: + raise ValueError('Invalid value for "alpha": {:}'.format(alpha)) + + # Create the actual matrix now. + H = np.zeros((2*n, n)) + for i in range(n): + H[2*i:2*i+2, i] = [-1, 1] + + H *= a + + return H, a diff --git a/pyDOE/doe_union.py b/pyDOE2/doe_union.py similarity index 96% rename from pyDOE/doe_union.py rename to pyDOE2/doe_union.py index 71bf60a..4922b5b 100644 --- a/pyDOE/doe_union.py +++ b/pyDOE2/doe_union.py @@ -1,47 +1,47 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np - -def union(H1, H2): - """ - Join two matrices by stacking them on top of each other. - - Parameters - ---------- - H1 : 2d-array - The matrix that goes on top of the new matrix - H2 : 2d-array - The matrix that goes on bottom of the new matrix - - Returns - ------- - mat : 2d-array - The new matrix that contains the rows of ``H1`` on top of the rows of - ``H2``. - - Example - ------- - :: - - >>> union(np.eye(2), -np.eye(2)) - array([[ 1., 0.], - [ 0., 1.], - [-1., 0.], - [ 0., -1.]]) - - """ - H = np.r_[H1, H2] - return H +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np + +def union(H1, H2): + """ + Join two matrices by stacking them on top of each other. + + Parameters + ---------- + H1 : 2d-array + The matrix that goes on top of the new matrix + H2 : 2d-array + The matrix that goes on bottom of the new matrix + + Returns + ------- + mat : 2d-array + The new matrix that contains the rows of ``H1`` on top of the rows of + ``H2``. + + Example + ------- + :: + + >>> union(np.eye(2), -np.eye(2)) + array([[ 1., 0.], + [ 0., 1.], + [-1., 0.], + [ 0., -1.]]) + + """ + H = np.r_[H1, H2] + return H diff --git a/pyDOE/var_regression_matrix.py b/pyDOE2/var_regression_matrix.py similarity index 96% rename from pyDOE/var_regression_matrix.py rename to pyDOE2/var_regression_matrix.py index 7ae5549..93d17d0 100644 --- a/pyDOE/var_regression_matrix.py +++ b/pyDOE2/var_regression_matrix.py @@ -1,51 +1,51 @@ -""" -This code was originally published by the following individuals for use with -Scilab: - Copyright (C) 2012 - 2013 - Michael Baudin - Copyright (C) 2012 - Maria Christopoulou - Copyright (C) 2010 - 2011 - INRIA - Michael Baudin - Copyright (C) 2009 - Yann Collette - Copyright (C) 2009 - CEA - Jean-Marc Martinez - - website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. It has been converted to Python by -Abraham Lee. -""" - -import numpy as np - -def var_regression_matrix(H, x, model, sigma=1): - """ - Compute the variance of the 'regression error'. - - Parameters - ---------- - H : 2d-array - The regression matrix - x : 2d-array - The coordinates to calculate the regression error variance at. - model : str - A string of tokens that define the regression model (e.g. - '1 x1 x2 x1*x2') - sigma : scalar - An estimate of the variance (default: 1). - - Returns - ------- - var : scalar - The variance of the regression error, evaluated at ``x``. - - """ - x = np.atleast_2d(x) - H = np.atleast_2d(H) - - if x.shape[0]==1: - x = x.T - - if np.rank(H)<(np.dot(H.T, H)).shape[0]: - raise ValueError("model and DOE don't suit together") - - x_mod = build_regression_matrix(x, model) - var = sigma**2*np.dot(np.dot(x_mod.T, np.linalg.inv(np.dot(H.T, H))), x_mod) - return var +""" +This code was originally published by the following individuals for use with +Scilab: + Copyright (C) 2012 - 2013 - Michael Baudin + Copyright (C) 2012 - Maria Christopoulou + Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + Copyright (C) 2009 - Yann Collette + Copyright (C) 2009 - CEA - Jean-Marc Martinez + + website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. It has been converted to Python by +Abraham Lee. +""" + +import numpy as np + +def var_regression_matrix(H, x, model, sigma=1): + """ + Compute the variance of the 'regression error'. + + Parameters + ---------- + H : 2d-array + The regression matrix + x : 2d-array + The coordinates to calculate the regression error variance at. + model : str + A string of tokens that define the regression model (e.g. + '1 x1 x2 x1*x2') + sigma : scalar + An estimate of the variance (default: 1). + + Returns + ------- + var : scalar + The variance of the regression error, evaluated at ``x``. + + """ + x = np.atleast_2d(x) + H = np.atleast_2d(H) + + if x.shape[0]==1: + x = x.T + + if np.rank(H)<(np.dot(H.T, H)).shape[0]: + raise ValueError("model and DOE don't suit together") + + x_mod = build_regression_matrix(x, model) + var = sigma**2*np.dot(np.dot(x_mod.T, np.linalg.inv(np.dot(H.T, H))), x_mod) + return var diff --git a/setup.py b/setup.py index abf56af..274e429 100644 --- a/setup.py +++ b/setup.py @@ -5,15 +5,15 @@ def read(fname): return open(os.path.join(os.path.dirname(__file__), fname)).read() setup( - name='pyDOE', + name='pyDOE2', version="0.3.8", author='Abraham Lee', author_email='tisimst@gmail.com', description='Design of experiments for Python', - url='https://github.com/tisimst/pyDOE', + url='https://github.com/clicumu/pyDOE2', license='BSD License (3-Clause)', - long_description=read('README'), - packages=['pyDOE'], + long_description=read('README.rst.rst'), + packages=['pyDOE2'], install_requires=['numpy', 'scipy'], keywords=[ 'DOE', @@ -30,12 +30,12 @@ def read(fname): 'License :: OSI Approved :: BSD License', 'Operating System :: OS Independent', 'Programming Language :: Python', - 'Programming Language :: Python :: 2.6', 'Programming Language :: Python :: 2.7', - 'Programming Language :: Python :: 3.0', - 'Programming Language :: Python :: 3.1', 'Programming Language :: Python :: 3.2', 'Programming Language :: Python :: 3.3', + 'Programming Language :: Python :: 3.4', + 'Programming Language :: Python :: 3.5', + 'Programming Language :: Python :: 3.6', 'Topic :: Education', 'Topic :: Scientific/Engineering', 'Topic :: Scientific/Engineering :: Mathematics', From e32f7bc5e2ca2cdab83f350834f3561ec017e0b7 Mon Sep 17 00:00:00 2001 From: RickardSjogren Date: Fri, 26 Jan 2018 09:01:10 +0100 Subject: [PATCH 4/4] Updated README --- README.md | 88 ++++++++++++++++++++++++++++++++++++++++ README.rst | 115 ----------------------------------------------------- setup.py | 2 +- 3 files changed, 89 insertions(+), 116 deletions(-) create mode 100644 README.md delete mode 100644 README.rst diff --git a/README.md b/README.md new file mode 100644 index 0000000..e7d839d --- /dev/null +++ b/README.md @@ -0,0 +1,88 @@ +pyDOE2: An experimental design package for python +===================================================== + +`pyDOE2` is a fork of the [`pyDOE`](https://github.com/tisimst/pyDOE) package +that is designed to help the scientist, engineer, statistician, etc., to +construct appropriate experimental designs. + +This fork came to life to solve bugs and issues that remained unsolved in the +original package. + +Capabilities +------------ + +The package currently includes functions for creating designs for any +number of factors: + +- Factorial Designs + - General Full-Factorial (``fullfact``) + - 2-level Full-Factorial (``ff2n``) + - 2-level Fractional Factorial (``fracfact``) + - Plackett-Burman (``pbdesign``) +- Response-Surface Designs + - Box-Behnken (``bbdesign``) + - Central-Composite (``ccdesign``) +- Randomized Designs + - Latin-Hypercube (``lhs``) + +See the original [package homepage](http://pythonhosted.org/pyDOE) for details +on usage and other notes. + +Requirements +------------ + +- NumPy +- SciPy + +Installation and download +------------------------- + +Through pip: + +``` +pip install pyDOE2 +``` + + +Credits +------- + +`pyDOE` was originally published by the following individuals for use with +Scilab: + +- Copyright (C) 2012 - 2013 - Michael Baudin +- Copyright (C) 2012 - Maria Christopoulou +- Copyright (C) 2010 - 2011 - INRIA - Michael Baudin +- Copyright (C) 2009 - Yann Collette +- Copyright (C) 2009 - CEA - Jean-Marc Martinez + +- Website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros + +Much thanks goes to these individuals. + +And thanks goes out to the following for finding and offering solutions for +bugs: + +- Ashmeet Singh + +And thanks to the following individuals for forking and working with `pyDOE2`: + +- Copyright (C) 2018 - Rickard Sjögren and Daniel Svensson + + +License +------- + +This package is provided under two licenses: + +1. The *BSD License* (3-clause) +2. Any other that the author approves (just ask!) + +References +---------- + +- [Factorial designs](http://en.wikipedia.org/wiki/Factorial_experiment) +- [Plackett-Burman designs](http://en.wikipedia.org/wiki/Plackett-Burman_design) +- [Box-Behnken designs](http://en.wikipedia.org/wiki/Box-Behnken_design) +- [Central composite designs](http://en.wikipedia.org/wiki/Central_composite_design) +- [Latin-Hypercube designs](http://en.wikipedia.org/wiki/Latin_hypercube_sampling) diff --git a/README.rst b/README.rst deleted file mode 100644 index beca9a9..0000000 --- a/README.rst +++ /dev/null @@ -1,115 +0,0 @@ -===================================================== -``pyDOE2``: The experimental design package for python -===================================================== - -The ``pyDOE`` package is designed to help the -**scientist, engineer, statistician,** etc., to construct appropriate -**experimental designs**. - -Capabilities -============ - -The package currently includes functions for creating designs for any -number of factors: - -- *Factorial Designs* - - #. **General Full-Factorial** (``fullfact``) - - #. **2-level Full-Factorial** (``ff2n``) - - #. **2-level Fractional Factorial** (``fracfact``) - - #. **Plackett-Burman** (``pbdesign``) - -- *Response-Surface Designs* - - #. **Box-Behnken** (``bbdesign``) - - #. **Central-Composite** (``ccdesign``) - -- *Randomized Designs* - - #. **Latin-Hypercube** (``lhs``) - -*See the* `package homepage`_ *for details on usage and other notes* - -What's New -========== - -In this release, an incorrect indexing variable in the internal LHS function -``_pdist`` has been corrected so point-distances are now calculated accurately. - -Requirements -============ - -- NumPy -- SciPy - -Installation and download -========================= - -See the `package homepage`_ for helpful hints relating to downloading -and installing pyDOE. - -Source Code -=========== - -The latest, bleeding-edge but working `code -`_ -and `documentation source -`_ are -available `on GitHub `_. - -Contact -======= - -Any feedback, questions, bug reports, or success stores should -be sent to the `author`_. I'd love to hear from you! - -Credits -======= - -This code was originally published by the following individuals for use with -Scilab: - -- Copyright (C) 2012 - 2013 - Michael Baudin -- Copyright (C) 2012 - Maria Christopoulou -- Copyright (C) 2010 - 2011 - INRIA - Michael Baudin -- Copyright (C) 2009 - Yann Collette -- Copyright (C) 2009 - CEA - Jean-Marc Martinez - -- Website: forge.scilab.org/index.php/p/scidoe/sourcetree/master/macros - -Much thanks goes to these individuals. - -And thanks goes out to the following for finding and offering solutions for -bugs: - -- Ashmeet Singh - -License -======= - -This package is provided under two licenses: - -1. The *BSD License* (3-clause) -2. Any other that the author approves (just ask!) - -References -========== - -- `Factorial designs`_ -- `Plackett-Burman designs`_ -- `Box-Behnken designs`_ -- `Central composite designs`_ -- `Latin-Hypercube designs`_ - -.. _author: mailto:tisimst@gmail.com -.. _Factorial designs: http://en.wikipedia.org/wiki/Factorial_experiment -.. _Box-Behnken designs: http://en.wikipedia.org/wiki/Box-Behnken_design -.. _Central composite designs: http://en.wikipedia.org/wiki/Central_composite_design -.. _Plackett-Burman designs: http://en.wikipedia.org/wiki/Plackett-Burman_design -.. _Latin-Hypercube designs: http://en.wikipedia.org/wiki/Latin_hypercube_sampling -.. _package homepage: http://pythonhosted.org/pyDOE -.. _lhs documentation: http://pythonhosted.org/pyDOE/randomized.html#latin-hypercube diff --git a/setup.py b/setup.py index 274e429..148499d 100644 --- a/setup.py +++ b/setup.py @@ -12,7 +12,7 @@ def read(fname): description='Design of experiments for Python', url='https://github.com/clicumu/pyDOE2', license='BSD License (3-Clause)', - long_description=read('README.rst.rst'), + long_description=read('README.md'), packages=['pyDOE2'], install_requires=['numpy', 'scipy'], keywords=[