def lhd(<EOL>dist=None,<EOL>size=None,<EOL>dims=<NUM_LIT:1>,<EOL>form="<STR_LIT>",<EOL>iterations=<NUM_LIT:100>,<EOL>showcorrelations=False,<EOL>):
assert dims > <NUM_LIT:0>, '<STR_LIT>'<EOL>if not size or not dist:<EOL><INDENT>return None<EOL><DEDENT>def _lhs(x, samples=<NUM_LIT:20>):<EOL><INDENT>"""<STR_LIT>"""<EOL>segmentSize = <NUM_LIT:1.0> / samples<EOL>numVars = x.shape[<NUM_LIT:1>]<EOL>out = np.zeros((samples, numVars))<EOL>pointValue = np.zeros(samples)<EOL>for n in range(numVars):<EOL><INDENT>for i in range(samples):<EOL><INDENT>segmentMin = i * segmentSize<EOL>point = segmentMin + (np.random.random() * segmentSize)<EOL>pointValue[i] = (point * (x[<NUM_LIT:1>, n] - x[<NUM_LIT:0>, n])) + x[<NUM_LIT:0>, n]<EOL><DEDENT>out[:, n] = pointValue<EOL><DEDENT>return _mix(out)<EOL><DEDENT>def _mix(data, dim="<STR_LIT>"):<EOL><INDENT>"""<STR_LIT>"""<EOL>data = np.atleast_2d(data)<EOL>n = data.shape[<NUM_LIT:0>]<EOL>if dim == "<STR_LIT>":<EOL><INDENT>data = data.T<EOL><DEDENT>data_rank = list(range(n))<EOL>for i in range(data.shape[<NUM_LIT:1>]):<EOL><INDENT>new_data_rank = np.random.permutation(data_rank)<EOL>vals, order = np.unique(<EOL>np.hstack((data_rank, new_data_rank)), return_inverse=True<EOL>)<EOL>old_order = order[:n]<EOL>new_order = order[-n:]<EOL>tmp = data[np.argsort(old_order), i][new_order]<EOL>data[:, i] = tmp[:]<EOL><DEDENT>if dim == "<STR_LIT>":<EOL><INDENT>data = data.T<EOL><DEDENT>return data<EOL><DEDENT>if form is "<STR_LIT>":<EOL><INDENT>if hasattr(dist, "<STR_LIT>"):  <EOL><INDENT>nvars = len(dist)<EOL>x = np.vstack((np.zeros(nvars), np.ones(nvars)))<EOL>unif_data = _lhs(x, samples=size)<EOL>dist_data = np.empty_like(unif_data)<EOL>for i, d in enumerate(dist):<EOL><INDENT>dist_data[:, i] = d.ppf(unif_data[:, i])<EOL><DEDENT><DEDENT>else:  <EOL><INDENT>nvars = dims<EOL>x = np.vstack((np.zeros(nvars), np.ones(nvars)))<EOL>unif_data = _lhs(x, samples=size)<EOL>dist_data = np.empty_like(unif_data)<EOL>for i in range(nvars):<EOL><INDENT>dist_data[:, i] = dist.ppf(unif_data[:, i])<EOL><DEDENT><DEDENT><DEDENT>elif form is "<STR_LIT>":<EOL><INDENT>def euclid_distance(arr):<EOL><INDENT>n = arr.shape[<NUM_LIT:0>]<EOL>ans = <NUM_LIT:0.0><EOL>for i in range(n - <NUM_LIT:1>):<EOL><INDENT>for j in range(i + <NUM_LIT:1>, n):<EOL><INDENT>d = np.sqrt(<EOL>np.sum(<EOL>[(arr[i, k] - arr[j, k]) ** <NUM_LIT:2> for k in range(arr.shape[<NUM_LIT:1>])]<EOL>)<EOL>)<EOL>ans += <NUM_LIT:1.0> / d ** <NUM_LIT:2><EOL><DEDENT><DEDENT>return ans<EOL><DEDENT>def fill_space(data):<EOL><INDENT>best = <NUM_LIT><EOL>for it in range(iterations):<EOL><INDENT>d = euclid_distance(data)<EOL>if d < best:<EOL><INDENT>d_opt = d<EOL>data_opt = data.copy()<EOL><DEDENT>data = _mix(data)<EOL><DEDENT>print("<STR_LIT>", d_opt)<EOL>return data_opt<EOL><DEDENT>if hasattr(dist, "<STR_LIT>"):  <EOL><INDENT>nvars = len(dist)<EOL>x = np.vstack((np.zeros(nvars), np.ones(nvars)))<EOL>unif_data = fill_space(_lhs(x, samples=size))<EOL>dist_data = np.empty_like(unif_data)<EOL>for i, d in enumerate(dist):<EOL><INDENT>dist_data[:, i] = d.ppf(unif_data[:, i])<EOL><DEDENT><DEDENT>else:  <EOL><INDENT>nvars = dims<EOL>x = np.vstack((np.zeros(nvars), np.ones(nvars)))<EOL>unif_data = fill_space(_lhs(x, samples=size))<EOL>dist_data = np.empty_like(unif_data)<EOL>for i in range(nvars):<EOL><INDENT>dist_data[:, i] = dist.ppf(unif_data[:, i])<EOL><DEDENT><DEDENT><DEDENT>elif form is "<STR_LIT>":<EOL><INDENT>raise NotImplementedError(<EOL>"<STR_LIT>"<EOL>)<EOL><DEDENT>else:<EOL><INDENT>raise ValueError('<STR_LIT>' % (form))<EOL><DEDENT>if dist_data.shape[<NUM_LIT:1>] > <NUM_LIT:1>:<EOL><INDENT>cor_matrix = np.zeros((nvars, nvars))<EOL>for i in range(nvars):<EOL><INDENT>for j in range(nvars):<EOL><INDENT>x_data = dist_data[:, i].copy()<EOL>y_data = dist_data[:, j].copy()<EOL>x_mean = x_data.mean()<EOL>y_mean = y_data.mean()<EOL>num = np.sum((x_data - x_mean) * (y_data - y_mean))<EOL>den = np.sqrt(<EOL>np.sum((x_data - x_mean) ** <NUM_LIT:2>) * np.sum((y_data - y_mean) ** <NUM_LIT:2>)<EOL>)<EOL>cor_matrix[i, j] = num / den<EOL>cor_matrix[j, i] = num / den<EOL><DEDENT><DEDENT>inv_cor_matrix = np.linalg.pinv(cor_matrix)<EOL>VIF = np.max(np.diag(inv_cor_matrix))<EOL>if showcorrelations:<EOL><INDENT>print("<STR_LIT>", cor_matrix)<EOL>print("<STR_LIT>", inv_cor_matrix)<EOL>print("<STR_LIT>", VIF)<EOL><DEDENT><DEDENT>return dist_data<EOL>
Create a Latin-Hypercube sample design based on distributions defined in the
`scipy.stats` module

Parameters
----------
dist: array_like
    frozen scipy.stats.rv_continuous or rv_discrete distribution objects 
    that are defined previous to calling LHD

size: int
    integer value for the number of samples to generate for each 
    distribution object

dims: int, optional
    if dist is a single distribution object, and dims > 1, the one
    distribution will be used to generate a size-by-dims sampled design

form: str, optional (non-functional at the moment)
    determines how the sampling is to occur, with the following optional 
    values:
        - 'randomized' - completely randomized sampling
        - 'spacefilling' - space-filling sampling (generally gives a more 
          accurate sampling of the design when the number of sample points 
          is small)
        - 'orthogonal' - balanced space-filling sampling (experimental)

    The 'spacefilling' and 'orthogonal' forms require some iterations to 
    determine the optimal sampling pattern. 

iterations: int, optional (non-functional at the moment)
    used to control the number of allowable search iterations for generating
    'spacefilling' and 'orthogonal' designs

Returns
-------
out: 2d-array, 
    A 2d-array where each column corresponds to each input distribution and 
    each row is a sample in the design

Examples
--------

Single distribution: 
    - uniform distribution, low = -1, width = 2

>>> import scipy.stats as ss
>>> d0 = ss.uniform(loc=-1,scale=2)
>>> print lhd(dist=d0,size=5)
[[ 0.51031081]
 [-0.28961427]
 [-0.68342107]
 [ 0.69784371]
 [ 0.12248842]]

Single distribution for multiple variables:
    - normal distribution, mean = 0, stdev = 1

>>> d1 = ss.norm(loc=0,scale=1)
>>> print lhd(dist=d1,size=7,dims=5)
[[-0.8612785   0.23034412  0.21808001]
 [ 0.0455778   0.07001606  0.31586419]
 [-0.978553    0.30394663  0.78483995]
 [-0.26415983  0.15235896  0.51462024]
 [ 0.80805686  0.38891031  0.02076505]
 [ 1.63028931  0.52104917  1.48016008]] 

Multiple distributions:
    - beta distribution, alpha = 2, beta = 5
    - exponential distribution, lambda = 1.5

>>> d2 = ss.beta(2,5)
>>> d3 = ss.expon(scale=1/1.5)
>>> print lhd(dist=(d1,d2,d3),size=6)
[[-0.8612785   0.23034412  0.21808001]
 [ 0.0455778   0.07001606  0.31586419]
 [-0.978553    0.30394663  0.78483995]
 [-0.26415983  0.15235896  0.51462024]
 [ 0.80805686  0.38891031  0.02076505]
 [ 1.63028931  0.52104917  1.48016008]]