import numpy as np
from scipy.linalg import toeplitz, hankel, cholesky
from scipy import spatial
from scipy import stats
from numpy import ma
import re
import itertools

"""
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.

https://github.com/tirthajyoti/Design-of-experiment-Python
https://pythonhosted.org/pyDOE/

GSD Copyright (C) 2018 - Rickard Sjoegren
https://github.com/clicumu/pyDOE2/blob/master/pyDOE2
"""


# Function for constructing a DataFrame from a numpy array generated by PyDOE function and individual lists
def construct_df(x, factor_lists):
    out = []
    for col in x:
        row = []
        for index in range(len(col)):
            row.append(factor_lists[index][int(col[index])])
        out.append(row)

    return out


# Function for constructing a DataFrame from a matrix with floating point numbers between -1 and +1
def construct_df_from_matrix(x, factor_array):
    assert 0

    """
    row_num = x.shape[0]  # Number of rows in the matrix x
    col_num = x.shape[1]  # Number of columns in the matrix x

    empty = np.zeros((row_num, col_num))

    def simple_substitution(idx, factor_list):
        if idx == -1:
            return factor_list[0]
        elif idx == 0:
            return factor_list[1]
        elif idx == 1:
            return factor_list[2]
        else:
            alpha = np.abs(factor_list[2] - factor_list[0]) / 2
            if idx < 0:
                beta = np.abs(idx) - 1
                return factor_list[0] - (beta * alpha)
            else:
                beta = idx - 1
                return factor_list[2] + (beta * alpha)

    for i in range(row_num):
        for j in range(col_num):
            empty[i, j] = simple_substitution(x[i, j], factor_array[j])

    return pd.DataFrame(data=empty)
    """

# Function for constructing a DataFrame from a matrix with floating point numbers between 0 and 1
def construct_df_from_random_matrix(x, factor_lists):
    out = []
    for w in x:
        row = []
        for index in range(len(w)):
            row.append(factor_lists[index][0] + (w[index] * np.fabs(factor_lists[index][1] - factor_lists[index][0])))

        out.append(row)

    return out


# Function for building Box-Behnken designs from a dictionary of process variables
def build_box_behnken(factor_level_ranges):
    """
    Builds a Box-Behnken design dataframe from a dictionary of factor/level ranges.
    Note 3 levels of factors are necessary. If not given, the function will automatically create 3 levels by linear mid-section method.
    Example of the dictionary:
    {'Pressure':[50, 70],'Temperature':[290, 350],'Flow rate':[0.9, 1.1]}

	In statistics, Box–Behnken designs are experimental designs for response surface methodology, devised by George E. P. Box and Donald Behnken in 1960, to achieve the following goals:
		* Each factor, or independent variable, is placed at one of three equally spaced values, usually coded as −1, 0, +1. (At least three levels are needed for the following goal.)
		* The design should be sufficient to fit a quadratic model, that is, one containing squared terms, products of two factors, linear terms and an intercept.
		* The ratio of the number of experimental points to the number of coefficients in the quadratic model should be reasonable (in fact, their designs kept it in the range of 1.5 to 2.6).*estimation variance should more or less depend only on the distance from the centre (this is achieved exactly for the designs with 4 and 7 factors), and should not vary too much inside the smallest (hyper)cube containing the experimental points.
	"""

    # creating a mid-point by averaging them
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) == 2:
            factor_level_ranges[key].append((factor_level_ranges[key][0] + factor_level_ranges[key][1]) / 2)
            factor_level_ranges[key].sort()


    factor_count = len(factor_level_ranges)
    factor_lists = []

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = bbdesign(factor_count, center=1)
    x = x + 1  # Adjusting the index up by 1

    df = construct_df(x, factor_lists)

    return df


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 = int((0.5 * n * (n - 1)) * 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


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)


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


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


# Function for building full factorial DataFrame from a dictionary of process variables
# ======================================================================================

def build_full_fact(factor_level_ranges):
    """
    Builds a full factorial design dataframe from a dictionary of factor/level ranges
    Example of the process variable dictionary:
    {'Pressure':[50,60,70],'Temperature':[290, 320, 350],'Flow rate':[0.9,1.0]}
    """

    factor_lvl_count = []
    factor_lists = []

    for key in factor_level_ranges:
        factor_lvl_count.append(len(factor_level_ranges[key]))
        factor_lists.append(factor_level_ranges[key])

    x = fullfact(factor_lvl_count)
    df = construct_df(x, factor_lists)

    return df


# ================================================================================================================================================================
# Function for building 2-level fractional factorial DataFrame from a dictionary of process variables and a generator string representing the confounded variable
# ================================================================================================================================================================

def build_frac_fact(factor_level_ranges, gen_string):
    """
    Builds a full factorial design dataframe from a dictionary of factor/level ranges.
    Only min and max values of the range are required.
	Example of the dictionary:
    {'Pressure':[50,70],'Temperature':[290, 350],'Flow rate':[0.9,1.0]}

	This function requires a little more knowledge of how the confounding will be allowed.
	This means that some factor effects get muddled with other interaction effects, so it’s harder to distinguish between them).

	Let’s assume that we just can’t afford (for whatever reason) the number of runs in a full-factorial design. We can systematically decide on a fraction of the full-factorial by allowing some of the factor main effects to be confounded with other factor interaction effects.
	This is done by defining an alias structure that defines, symbolically, these interactions. These alias structures are written like “C = AB” or “I = ABC”, or “AB = CD”, etc.
	These define how one column is related to the others.

	EXAMPLE
    ------------
    For example, the alias “C = AB” or “I = ABC” indicate that there are three factors (A, B, and C) and that the main effect of factor C is confounded with the interaction effect of the product AB, and by extension, A is confounded with BC and B is confounded with AC.
	A full- factorial design with these three factors results in a design matrix with 8 runs, but we will assume that we can only afford 4 of those runs.
	To create this fractional design, we need a matrix with three columns, one for A, B, and C, only now where the levels in the C column is created by the product of the A and B columns.
    """

    factor_count = len(factor_level_ranges)
    factor_lists = []

    # Assigning the end point to the high level.
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) != 2:
            factor_level_ranges[key][1] = factor_level_ranges[key][-1]
            factor_level_ranges[key] = factor_level_ranges[key][:2]

    if factor_count != len(gen_string.split(' ')):
        print("Length of the generator string for the fractional factorial build does not match the length of the process variables dictionary")
        return None

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = fracfact(gen_string)

    def index_change(x):
        if x == -1:
            return 0
        else:
            return x

    vfunc = np.vectorize(index_change)
    x = vfunc(x)

    df = construct_df(x, factor_lists)

    return df


# =====================================================================================
# Function for building Plackett-Burman designs from a dictionary of process variables
# =====================================================================================

def build_plackett_burman(factor_level_ranges):
    """
    Builds a Plackett-Burman dataframe from a dictionary of factor/level ranges.
    Only min and max values of the range are required.
    Example of the dictionary:
    {'Pressure':[50,70],'Temperature':[290, 350],'Flow rate':[0.9,1.0]}

	Plackett–Burman designs are experimental designs presented in 1946 by Robin L. Plackett and J. P. Burman while working in the British Ministry of Supply.(Their goal was to find experimental designs for investigating the dependence of some measured quantity on a number of independent variables (factors), each taking L levels, in such a way as to minimize the variance of the estimates of these dependencies using a limited number of experiments.

    Interactions between the factors were considered negligible. The solution to this problem is to find an experimental design where each combination of levels for any pair of factors appears the same number of times, throughout all the experimental runs (refer to table).
	A complete factorial design would satisfy this criterion, but the idea was to find smaller designs.

	These designs are unique in that the number of trial conditions (rows) expands by multiples of four (e.g. 4, 8, 12, etc.).
	The max number of columns allowed before a design increases the number of rows is always one less than the next higher multiple of four.
    """

    # Assigning the end point to the high level.
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) != 2:
            factor_level_ranges[key][1] = factor_level_ranges[key][-1]
            factor_level_ranges[key] = factor_level_ranges[key][:2]

    factor_count = len(factor_level_ranges)
    factor_lists = []

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = pbdesign(factor_count)

    def index_change(x):
        if x == -1:
            return 0
        else:
            return x

    vfunc = np.vectorize(index_change)
    x = vfunc(x)

    df = construct_df(x, factor_lists)
    

    return df


# ===================================================================================
# Function for building Box-Behnken designs from a dictionary of process variables
# ===================================================================================

def build_box_behnken(factor_level_ranges):
    """
    Builds a Box-Behnken design dataframe from a dictionary of factor/level ranges.
    Note 3 levels of factors are necessary. If not given, the function will automatically create 3 levels by linear mid-section method.
    Example of the dictionary:
    {'Pressure':[50,60,70],'Temperature':[290, 320, 350],'Flow rate':[0.9,1.0,1.1]}

	In statistics, Box–Behnken designs are experimental designs for response surface methodology, devised by George E. P. Box and Donald Behnken in 1960, to achieve the following goals:
		* Each factor, or independent variable, is placed at one of three equally spaced values, usually coded as −1, 0, +1. (At least three levels are needed for the following goal.)
		* The design should be sufficient to fit a quadratic model, that is, one containing squared terms, products of two factors, linear terms and an intercept.
		* The ratio of the number of experimental points to the number of coefficients in the quadratic model should be reasonable (in fact, their designs kept it in the range of 1.5 to 2.6).*estimation variance should more or less depend only on the distance from the centre (this is achieved exactly for the designs with 4 and 7 factors), and should not vary too much inside the smallest (hyper)cube containing the experimental points.
	"""

    # Creating a mid-point by averaging them
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) == 2:
            factor_level_ranges[key].append((factor_level_ranges[key][0] + factor_level_ranges[key][1]) / 2)
            factor_level_ranges[key].sort()

    factor_count = len(factor_level_ranges)
    factor_lists = []

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = bbdesign(factor_count, center=1)
    x = x + 1  # Adjusting the index up by 1

    df = construct_df(x, factor_lists)

    return df


# Function for building central-composite (Box-Wilson) designs from a dictionary of process variables
def build_central_composite(factor_level_ranges, center=(2, 2), alpha='o', face='ccc'):
    """
    Builds a central-composite design dataframe from a dictionary of factor/level ranges.
    Only min and max values of the range are required.
    Example of the dictionary:
    {'Pressure':[50,70],'Temperature':[290, 350],'Flow rate':[0.9,1.0]}

	In statistics, a central composite design is an experimental design, useful in response surface methodology, for building a second order (quadratic) model for the response variable without needing to use a complete three-level factorial experiment.
	The design consists of three distinct sets of experimental runs:
		* A factorial (perhaps fractional) design in the factors studied, each having two levels;
		* A set of center points, experimental runs whose values of each factor are the medians of the values used in the factorial portion. This point is often replicated in order to improve the precision of the experiment;
		* A set of axial points, experimental runs identical to the centre points except for one factor, which will take on values both below and above the median of the two factorial levels, and typically both outside their range. All factors are varied in this way.
    """
    # Assigning the end point to the high level.
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) != 2:
            factor_level_ranges[key][1] = factor_level_ranges[key][-1]
            factor_level_ranges[key] = factor_level_ranges[key][:2]

    # Creates the mid-points by averaging the low and high levels
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) == 2:
            factor_level_ranges[key].append((factor_level_ranges[key][0] + factor_level_ranges[key][1]) / 2)
            factor_level_ranges[key].sort()

    factor_count = len(factor_level_ranges)
    factor_lists = []

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = ccdesign(factor_count, center=center, alpha=alpha, face=face)
    factor_lists = np.array(factor_lists)

    # rescale x from (-1, 1) to (0, 1)

    df = construct_df_from_matrix(x, factor_lists)
    
    return df


def build_lhs(factor_level_ranges, num_samples=None):
    """
    Builds a Latin Hypercube design dataframe from a dictionary of factor/level ranges.
    Only min and max values of the range are required.
    Example of the dictionary:
    {'Pressure':[50,70],'Temperature':[290, 350],'Flow rate':[0.9,1.0]}
    num_samples: Number of samples to be generated

	Latin hypercube sampling (LHS) is a form of stratified sampling that can be applied to multiple variables. The method commonly used to reduce the number or runs necessary for a Monte Carlo simulation to achieve a reasonably accurate random distribution. LHS can be incorporated into an existing Monte Carlo model fairly easily, and work with variables following any analytical probability distribution.
    """

    # Assigning the end point to the high level.
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) != 2:
            factor_level_ranges[key][1] = factor_level_ranges[key][-1]
            factor_level_ranges[key] = factor_level_ranges[key][:2]

    factor_count = len(factor_level_ranges)
    factor_lists = []

    if num_samples == None:
        num_samples = factor_count

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = lhs(n=factor_count, samples=num_samples)
    factor_lists = np.array(factor_lists)

    df = construct_df_from_random_matrix(x, factor_lists)

    return df


def lhs(n, samples=None, criterion=None, iterations=None, random_state=None, correlation_matrix=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).
    randomstate : np.random.RandomState, int
         Random state (or seed-number) which controls the seed and random draws
    correlation_matrix : ndarray
         Enforce correlation between factors (only used in lhsmu)
    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, random_state=42)
        array([[ 0.12484671,  0.95539205,  0.24399798],
               [ 0.53288616,  0.38533955,  0.86703834],
               [ 0.68602787,  0.31690477,  0.38533151]])
    A 4-factor design with 6 samples::
        >>> lhs(4, samples=6, random_state=42)
        array([[ 0.06242335,  0.19266575,  0.88202411,  0.89439364],
               [ 0.19266977,  0.53538985,  0.53030416,  0.49498498],
               [ 0.71737371,  0.75412607,  0.17634727,  0.71520486],
               [ 0.63874044,  0.85658231,  0.33676408,  0.31102936],
               [ 0.43351917,  0.45134543,  0.12199899,  0.53056742],
               [ 0.93530882,  0.15845238,  0.7386575 ,  0.09977641]])
    A 2-factor design with 5 centered samples::
        >>> lhs(2, samples=5, criterion='center', random_state=42)
        array([[ 0.1,  0.9],
               [ 0.5,  0.5],
               [ 0.7,  0.1],
               [ 0.3,  0.7],
               [ 0.9,  0.3]])
    A 3-factor design with 4 samples where the minimum distance between
    all samples has been maximized::
        >>> lhs(3, samples=4, criterion='maximin', random_state=42)
        array([[ 0.69754389,  0.2997106 ,  0.96250964],
               [ 0.10585037,  0.09872038,  0.73157522],
               [ 0.25351996,  0.65148999,  0.07337204],
               [ 0.91276926,  0.97873992,  0.42783549]])
    A 4-factor design with 5 samples where the samples are as uncorrelated
    as possible (within 10 iterations)::
        >>> lhs(4, samples=5, criterion='correlation', iterations=10, random_state=42)
        array([[ 0.72088348,  0.05121366,  0.97609357,  0.92487081],
               [ 0.49507404,  0.51265511,  0.00808672,  0.37915272],
               [ 0.22217816,  0.2878673 ,  0.24034384,  0.42786629],
               [ 0.91977309,  0.93895699,  0.64061224,  0.14213258],
               [ 0.04719698,  0.70796822,  0.53910322,  0.78857071]])
    """
    H = None

    if random_state is None:
        random_state = np.random.RandomState()
    elif not isinstance(random_state, np.random.RandomState):
        random_state = np.random.RandomState(random_state)

    if samples is None:
        samples = n

    if criterion is not None:
        if not criterion.lower() in ('center', 'c', 'maximin', 'm',
                                     'centermaximin', 'cm', 'correlation',
                                     'corr', 'lhsmu'):
            raise ValueError('Invalid value for "criterion": {}'.format(criterion))

    else:
        H = _lhsclassic(n, samples, random_state)

    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, random_state)
        elif criterion.lower() in ('maximin', 'm'):
            H = _lhsmaximin(n, samples, iterations, 'maximin', random_state)
        elif criterion.lower() in ('centermaximin', 'cm'):
            H = _lhsmaximin(n, samples, iterations, 'centermaximin', random_state)
        elif criterion.lower() in ('correlation', 'corr'):
            H = _lhscorrelate(n, samples, iterations, random_state)
        elif criterion.lower() in ('lhsmu'):
            # as specified by the paper. M is set to 5
            H = _lhsmu(n, samples, correlation_matrix, random_state, M=5)

    return H


################################################################################

def _lhsclassic(n, samples, randomstate):
    # Generate the intervals
    cut = np.linspace(0, 1, samples + 1)

    # Fill points uniformly in each interval
    u = randomstate.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 = randomstate.permutation(range(samples))
        H[:, j] = rdpoints[order, j]

    return H


################################################################################

def _lhscentered(n, samples, randomstate):
    # Generate the intervals
    cut = np.linspace(0, 1, samples + 1)

    # Fill points uniformly in each interval
    u = randomstate.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] = randomstate.permutation(_center)

    return H


################################################################################

def _lhsmaximin(n, samples, iterations, lhstype, randomstate):
    maxdist = 0

    # Maximize the minimum distance between points
    for i in range(iterations):
        if lhstype == 'maximin':
            Hcandidate = _lhsclassic(n, samples, randomstate)
        else:
            Hcandidate = _lhscentered(n, samples, randomstate)

        d = spatial.distance.pdist(Hcandidate, 'euclidean')
        if maxdist < np.min(d):
            maxdist = np.min(d)
            H = Hcandidate.copy()

    return H


################################################################################

def _lhscorrelate(n, samples, iterations, randomstate):
    mincorr = np.inf

    # Minimize the components correlation coefficients
    for i in range(iterations):
        # Generate a random LHS
        Hcandidate = _lhsclassic(n, samples, randomstate)
        R = np.corrcoef(Hcandidate.T)
        if np.max(np.abs(R[R != 1])) < mincorr:
            mincorr = np.max(np.abs(R - np.eye(R.shape[0])))
            H = Hcandidate.copy()

    return H


################################################################################

def _lhsmu(N, samples=None, corr=None, random_state=None, M=5):
    if random_state is None:
        random_state = np.random.RandomState()
    elif not isinstance(random_state, np.random.RandomState):
        random_state = np.random.RandomState(random_state)

    if samples is None:
        samples = N

    I = M * samples

    rdpoints = random_state.uniform(size=(I, N))

    dist = spatial.distance.cdist(rdpoints, rdpoints, metric='euclidean')
    D_ij = ma.masked_array(dist, mask=np.identity(I))

    index_rm = np.zeros(I - samples, dtype=int)
    i = 0
    while i < I - samples:
        order = ma.sort(D_ij, axis=1)

        avg_dist = ma.mean(order[:, 0:2], axis=1)
        min_l = ma.argmin(avg_dist)

        D_ij[min_l, :] = ma.masked
        D_ij[:, min_l] = ma.masked

        index_rm[i] = min_l
        i += 1

    rdpoints = np.delete(rdpoints, index_rm, axis=0)

    if (corr is not None):
        # check if covariance matrix is valid
        assert type(corr) == np.ndarray
        assert corr.ndim == 2
        assert corr.shape[0] == corr.shape[1]
        assert corr.shape[0] == N

        norm_u = stats.norm().ppf(rdpoints)
        L = cholesky(corr, lower=True)

        norm_u = np.matmul(norm_u, L)

        H = stats.norm().cdf(norm_u)
    else:
        H = np.zeros_like(rdpoints, dtype=float)
        rank = np.argsort(rdpoints, axis=0)

        for l in range(samples):
            low = float(l) / samples
            high = float(l + 1) / samples

            l_pos = rank == l
            H[l_pos] = random_state.uniform(low, high, size=N)
    return H


################################################################################

def _pdist(x):
    """
    Calculate the pair-wise point distances of a matrix

    Parameters
    ----------
    x : 2d-array
        An m-by-n array of scalars, where there are m points in n dimensions.

    Returns
    -------
    d : array
        A 1-by-b array of scalars, where b = m*(m - 1)/2. This array contains
        all the pair-wise point distances, arranged in the order (1, 0),
        (2, 0), ..., (m-1, 0), (2, 1), ..., (m-1, 1), ..., (m-1, m-2).

    Examples
    --------
    ::

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


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


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


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


def build_gsd(levels, reduction, n=1):
    """
    Create a Generalized Subset Design (GSD).
    Parameters
    ----------
    levels : array-like
        Number of factor levels per factor in design.
    reduction : int
        Reduction factor (bigger than 1). Larger `reduction` means fewer
        experiments in the design and more possible complementary designs.
    n : int
        Number of complementary GSD-designs (default 1). The complementary
        designs are balanced analogous to fold-over in two-level fractional
        factorial designs.
    Returns
    -------
    H : 2d-array | list of 2d-arrays
        `n` m-by-k matrices where k is the number of factors (equal
        to the length of `factor_levels`. The number of rows, m, will
        be approximately equal to the grand product of the factor levels
        divided by `reduction`.
    Raises
    ------
    ValueError
        If input is valid or if design construction fails. Design can fail
        if `reduction` is too large compared to values of `levels`.
    Notes
    -----
    The Generalized Subset Design (GSD) [1]_ or generalized factorial design is
    a generalization of traditional fractional factorial designs to problems
    where factors can have more than two levels.
    In many application problems factors can have categorical or quantitative
    factors on more than two levels. Previous reduced designs have not been
    able to deal with such types of problems. Full multi-level factorial
    designs can handle such problems but are however not economical regarding
    the number of experiments.
    Note for commercial users, the application of GSD to testing of product
    characteristics in a processing facility is patented [2]_
    Examples
    --------
    An example with three factors using three, four and
    six levels respectively reduced with a factor 4 ::
        >>> gsd([3, 4, 6], 4)
        array([[0, 0, 0],
               [0, 0, 4],
               [0, 1, 1],
               [0, 1, 5],
               [0, 2, 2],
               [0, 3, 3],
               [1, 0, 1],
               [1, 0, 5],
               [1, 1, 2],
               [1, 2, 3],
               [1, 3, 0],
               [1, 3, 4],
               [2, 0, 2],
               [2, 1, 3],
               [2, 2, 0],
               [2, 2, 4],
               [2, 3, 1],
               [2, 3, 5]])
    Two complementary designs with two factors using three and
    four levels reduced with a factor 2 ::
        >>> gsd([3, 4], 2, n=2)[0]
        array([[0, 0],
               [0, 2],
               [2, 0],
               [2, 2],
               [1, 1],
               [1, 3]])
        >>> gsd([3, 4], 2, n=2)[1]
        array([[0, 1],
               [0, 3],
               [2, 1],
               [2, 3],
               [1, 0],
               [1, 2]])
    If design fails ValueError is raised ::
        >>> gsd([2, 3], 5)
        Traceback (most recent call last):
         ...
        ValueError: reduction too large compared to factor levels
    References
    ----------
    .. [1] Surowiec, Izabella, Ludvig Vikstrom, Gustaf Hector, Erik Johansson,
       Conny Vikstrom, and Johan Trygg. "Generalized Subset Designs in
       Analytical Chemistry." Analytical Chemistry 89, no. 12 (June 20, 2017):
       6491-97. https://doi.org/10.1021/acs.analchem.7b00506.
    .. [2] Vikstrom, Ludvig, Conny Vikstrom, Erik Johansson, and Gustaf Hector.
       Computer-implemented systems and methods for generating
       generalized fractional designs. US9746850 B2, filed May 9,
       2014, and issued August 29, 2017. http://www.google.se/patents/US9746850.
    """
    try:
        assert all(isinstance(v, int) for v in levels), \
            'levels has to be sequence of integers'
        assert isinstance(reduction, int) and reduction > 1, \
            'reduction has to be integer larger than 1'
        assert isinstance(n, int) and n > 0, \
            'n has to be positive integer'
    except AssertionError as e:
        raise ValueError(e)

    partitions = _make_partitions(levels, reduction)
    latin_square = _make_latin_square(reduction)
    ortogonal_arrays = _make_orthogonal_arrays(latin_square, len(levels))

    try:
        designs = [_map_partitions_to_design(partitions, oa) - 1 for oa in
                   ortogonal_arrays]
    except ValueError:
        raise ValueError('reduction too large compared to factor levels')

    if n == 1:
        return designs[0]
    else:
        return designs[:n]


def _make_orthogonal_arrays(latin_square, n_cols):
    """
    Augment latin-square to the specified number of columns to produce
    an orthogonal array.
    """
    p = len(latin_square)

    first_row = latin_square[0]
    A_matrices = [np.array([[v]]) for v in first_row]

    while A_matrices[0].shape[1] < n_cols:
        new_A_matrices = list()

        for i, A_matrix in enumerate(A_matrices):
            sub_a = list()
            for constant, other_A in zip(first_row,
                                         np.array(A_matrices)[latin_square[i]]):
                constant_vec = np.repeat(constant, len(other_A))[:, np.newaxis]
                combined = np.hstack([constant_vec, other_A])
                sub_a.append(combined)

            new_A_matrices.append(np.vstack(sub_a))

        A_matrices = new_A_matrices

        if A_matrices[0].shape[1] == n_cols:
            break

    return A_matrices


def _map_partitions_to_design(partitions, ortogonal_array):
    """
    Map partitioned factor to final design using orthogonal-array produced
    by augmenting latin square.
    """
    assert len(
        partitions) == ortogonal_array.max() + 1 and ortogonal_array.min() == 0, \
        'Orthogonal array indexing does not match partition structure'

    mappings = list()
    for row in ortogonal_array:
        if any(not partitions[p][factor] for factor, p in enumerate(row)):
            continue

        partition_sets = [partitions[p][factor] for factor, p in enumerate(row)]
        mapping = list(itertools.product(*partition_sets))
        mappings.append(mapping)

    return np.vstack(mappings)


def _make_partitions(factor_levels, num_partitions):
    """
    Balanced partitioning of factors.
    """
    partitions = list()
    for partition_i in range(1, num_partitions + 1):
        partition = list()

        for num_levels in factor_levels:
            part = list()
            for level_i in range(1, num_levels):
                index = partition_i + (level_i - 1) * num_partitions
                if index <= num_levels:
                    part.append(index)

            partition.append(part)

        partitions.append(partition)

    return partitions


def _make_latin_square(n):
    numbers = np.arange(n)
    latin_square = np.vstack([np.roll(numbers, -i) for i in range(n)])
    return latin_square


def build_halton(factor_level_ranges, num_samples=None):
    """
    Builds a quasirandom dataframe from a dictionary of factor/level ranges using prime numbers as seed.
    Only min and max values of the range are required.
    Example of the dictionary which is needed as the input:
    {'Pressure':[50,70],'Temperature':[290, 350],'Flow rate':[0.9,1.0]}
    num_samples: Number of samples to be generated
    Quasirandom sequence using the default initialization with first n prime numbers equal to the number of factors/variables.
    """
    for key in factor_level_ranges:
        if len(factor_level_ranges[key]) != 2:
            factor_level_ranges[key][1] = factor_level_ranges[key][-1]
            factor_level_ranges[key] = factor_level_ranges[key][:2]
            print(
                f"{key} had more than two levels. Assigning the end point to the high level."
            )

    factor_count = len(factor_level_ranges)
    factor_lists = []

    if num_samples == None:
        num_samples = factor_count

    for key in factor_level_ranges:
        factor_lists.append(factor_level_ranges[key])

    x = halton(
        num_points=num_samples, dimension=factor_count
    )  # create Halton matrix design
    factor_lists = np.array(factor_lists)

    df = construct_df_from_random_matrix(x, factor_lists)
    # df.columns = factor_level_ranges.keys()
    return df


def halton(num_points, dimension):
    """Halton sequence.
    :param int dim: dimension
    :param int n_sample: number of samples.
    :return: sequence of Halton.
    :rtype: array_like (n_samples, n_features)
    """
    big_number = 10
    while 'Not enought primes':
        base = _primes_from_2_to(big_number)[:dimension]
        if len(base) == dimension:
            break
        big_number += 1000

    # Generate a sample using a Van der Corput sequence per dimension.
    sample = [_van_der_corput(num_points + 1, dimension) for dimension in base]
    sample = np.stack(sample, axis=-1)[1:]

    return sample


def _primes_from_2_to(n):
    """Prime number from 2 to n.
    From `StackOverflow <https://stackoverflow.com/questions/2068372>`_.
    :param int n: sup bound with ``n >= 6``.
    :return: primes in 2 <= p < n.
    :rtype: list
    """
    sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool)
    for i in range(1, int(n ** 0.5) // 3 + 1):
        if sieve[i]:
            k = 3 * i + 1 | 1
            sieve[k * k // 3::2 * k] = False
            sieve[k * (k - 2 * (i & 1) + 4) // 3::2 * k] = False
    return np.r_[2, 3, ((3 * np.nonzero(sieve)[0][1:] + 1) | 1)]


def _van_der_corput(n_sample, base=2):
    """Van der Corput sequence.
    :param int n_sample: number of element of the sequence.
    :param int base: base of the sequence.
    :return: sequence of Van der Corput.
    :rtype: list (n_samples,)
    """
    sequence = []
    for i in range(n_sample):
        n_th_number, denom = 0., 1.
        while i > 0:
            i, remainder = divmod(i, base)
            denom *= base
            n_th_number += remainder / denom
        sequence.append(n_th_number)

    return sequence
