Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
3 changed files
with
322 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,82 @@ | ||
''' | ||
Created on 30 Jun 2015 | ||
@author: will2 | ||
''' | ||
|
||
import numpy as np | ||
from . import common_args | ||
from ..util import read_param_file, unscale_samples | ||
from ..sample.ff import generate_contrast, sample, extend_bounds | ||
|
||
def analyze(problem, X, Y, second_order=False, print_to_console=False): | ||
|
||
problem = extend_bounds(problem) | ||
num_vars = problem['num_vars'] | ||
|
||
X = generate_contrast(problem) | ||
|
||
main_effect = (1. / (2 * num_vars)) * np.dot(Y, X) | ||
|
||
Si = dict((k, [None] * num_vars) | ||
for k in ['names', 'ME']) | ||
Si['ME'] = main_effect | ||
Si['names'] = problem['names'] | ||
|
||
if print_to_console: | ||
print("Parameter ME") | ||
for j in range(num_vars): | ||
print("%s %f" % (problem['names'][j], Si['ME'][j])) | ||
|
||
if second_order == True: | ||
interaction_names, interaction_effects = interactions(problem, | ||
Y, | ||
print_to_console) | ||
|
||
Si['names'].append(interaction_names) | ||
Si['IE'] = interaction_effects | ||
|
||
return Si | ||
|
||
|
||
def interactions(problem, Y, print_to_console=False): | ||
''' | ||
Computes the second order effects (interactions) between | ||
all combinations of pairs of input factors | ||
''' | ||
|
||
names = problem['names'] | ||
num_vars = problem['num_vars'] | ||
|
||
X = generate_contrast(problem) | ||
|
||
ie_names = [] | ||
IE = [] | ||
|
||
for col in range(X.shape[1]): | ||
for col_2 in range(col): | ||
x = X[:, col] * X[:, col_2] | ||
var_names = names[col_2] + names[col] | ||
ie_names.append(var_names) | ||
IE.append((1. / (2 * num_vars)) * np.dot(Y, x)) | ||
if print_to_console: | ||
print[('%s %f \n' % (n, i) ) for (n, i) in zip(ie_names, IE) ] | ||
|
||
return ie_names, IE | ||
|
||
if __name__ == "__main__": | ||
|
||
parser = common_args.create() | ||
parser.add_argument('-X', '--model-input-file', type=str, | ||
required=True, default=None, help='Model input file') | ||
args = parser.parse_args() | ||
|
||
problem = read_param_file(args.paramfile) | ||
|
||
Y = np.loadtxt(args.model_output_file, delimiter=args.delimiter, usecols=(args.column,)) | ||
X = np.loadtxt(args.model_input_file, delimiter=args.delimiter, ndmin=2) | ||
if len(X.shape) == 1: | ||
X = X.reshape((len(X), 1)) | ||
|
||
analyze(problem, X, Y, num_resamples=args.resamples, print_to_console=True, | ||
num_levels=args.levels, grid_jump=args.grid_jump) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,70 @@ | ||
''' | ||
Created on 30 Jun 2015 | ||
@author: will2 | ||
''' | ||
|
||
from scipy.linalg import hadamard | ||
import numpy as np | ||
from . import common_args | ||
from ..util import scale_samples, read_param_file | ||
|
||
|
||
def find_smallest(num_vars): | ||
''' | ||
Find the smallest exponent of two that is greater than the number | ||
of variables | ||
''' | ||
for x in range(10): | ||
if num_vars <= 2 ** x: | ||
return x | ||
|
||
|
||
def extend_bounds(problem): | ||
|
||
num_vars = problem['num_vars'] | ||
num_ff_vars = 2 ** find_smallest(num_vars) | ||
num_dummy_variables = num_ff_vars - num_vars | ||
|
||
bounds = list(problem['bounds']) | ||
names = problem['names'] | ||
if num_dummy_variables > 0: | ||
bounds.extend([[0, 1] for x in range(num_dummy_variables)]) | ||
names.extend(["dummy_" + str(var) for var in range(num_dummy_variables)]) | ||
problem['bounds'] = bounds | ||
problem['names'] = names | ||
problem['num_vars'] = num_ff_vars | ||
|
||
return problem | ||
|
||
def generate_contrast(problem): | ||
|
||
num_vars = problem['num_vars'] | ||
|
||
# Find the smallest n, such that num_vars < k | ||
k = [2 ** n for n in range(16)] | ||
k_chosen = 2 ** find_smallest(num_vars) | ||
|
||
# Generate the fractional factorial contrast | ||
contrast = np.vstack([hadamard(k_chosen), -hadamard(k_chosen)]) | ||
|
||
return contrast | ||
|
||
def sample(problem): | ||
|
||
contrast = generate_contrast(problem) | ||
sample = np.array((contrast + 1.) / 2, dtype=np.float) | ||
problem = extend_bounds(problem) | ||
scale_samples(sample, problem['bounds']) | ||
return sample | ||
|
||
if __name__ == "__main__": | ||
|
||
parser = common_args.create() | ||
args = parser.parse_args() | ||
|
||
problem = read_param_file(args.paramfile) | ||
param_values = sample(problem) | ||
|
||
np.savetxt(args.output, param_values, delimiter=args.delimiter, | ||
fmt='%.' + str(args.precision) + 'e') |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,170 @@ | ||
''' | ||
Created on 30 Jun 2015 | ||
@author: will2 | ||
''' | ||
import numpy as np | ||
from numpy.testing import assert_equal, assert_allclose | ||
from ..sample.ff import sample, find_smallest, extend_bounds | ||
from ..analyze.ff import analyze, interactions | ||
|
||
|
||
def test_find_smallest(): | ||
''' | ||
''' | ||
num_vars = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16, 17, 31, 32, 33] | ||
expected = [0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6] | ||
for x, y in zip(num_vars, expected): | ||
actual = find_smallest(x) | ||
assert_equal(actual, y) | ||
|
||
|
||
def test_extend_bounds(): | ||
problem = {'bounds': np.repeat([-1, 1], 12).reshape(2, 12).T, | ||
'num_vars': 12, | ||
'names': ["x" + str(x + 1) for x in range(12)] | ||
} | ||
actual = extend_bounds(problem) | ||
expected = {'names': ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x12', 'dummy_0', 'dummy_1', 'dummy_2', 'dummy_3'], | ||
'bounds': [np.array([-1, 1]), np.array([-1, 1]), | ||
np.array([-1, 1]), np.array([-1, 1]), | ||
np.array([-1, 1]), np.array([-1, 1]), | ||
np.array([-1, 1]), np.array([-1, 1]), | ||
np.array([-1, 1]), np.array([-1, 1]), | ||
np.array([-1, 1]), np.array([-1, 1]), | ||
np.array([0, 1]), np.array([0, 1]), | ||
np.array([0, 1]), np.array([0, 1])], | ||
'num_vars': 16} | ||
|
||
assert_equal(actual, expected) | ||
|
||
def test_ff_sample(): | ||
problem = {'bounds': [[0., 1.], [0., 1.], [0., 1.], [0., 1.]], | ||
'num_vars': 4, | ||
'names': ['x1', 'x2', 'x3', 'x4']} | ||
actual = sample(problem) | ||
expected = np.array([[ 1, 1, 1, 1], | ||
[ 1, 0, 1, 0], | ||
[ 1, 1, 0, 0], | ||
[ 1, 0, 0, 1], | ||
[0, 0, 0, 0], | ||
[0, 1, 0, 1], | ||
[0, 0, 1, 1], | ||
[0, 1, 1, 0]], dtype=np.float) | ||
assert_equal(actual, expected) | ||
|
||
|
||
def test_ff_sample_scaled(): | ||
''' | ||
''' | ||
problem = {'bounds': [[0., 2.5], [0., 1.], [0., 1.], [0., 1.]], | ||
'num_vars': 4, | ||
'names': ['x1', 'x2', 'x3', 'x4']} | ||
actual = sample(problem) | ||
expected = np.array([[ 2.5, 1, 1, 1], | ||
[ 2.5, 0, 1, 0], | ||
[ 2.5, 1, 0, 0], | ||
[ 2.5, 0, 0, 1], | ||
[0, 0, 0, 0], | ||
[0, 1, 0, 1], | ||
[0, 0, 1, 1], | ||
[0, 1, 1, 0]], dtype=np.float) | ||
assert_equal(actual, expected) | ||
|
||
|
||
def test_ff_analyze(): | ||
''' | ||
''' | ||
|
||
problem = {'bounds': [[0., 2.5], [0., 1.], [0., 1.], [0., 1.]], | ||
'num_vars': 4, | ||
'names': ['x1', 'x2', 'x3', 'x4']} | ||
X = np.array([[ 1, 1, 1, 1], | ||
[ 1, 0, 1, 0], | ||
[ 1, 1, 0, 0], | ||
[ 1, 0, 0, 1], | ||
[0, 0, 0, 0], | ||
[0, 1, 0, 1], | ||
[0, 0, 1, 1], | ||
[0, 1, 1, 0]], dtype=np.float) | ||
Y = np.array([1.5, 1, 1.5, 1, 2, 2.5, 2, 2.5], dtype=np.float) | ||
actual = analyze(problem, X, Y) | ||
expected = {'ME': np.array([ -0.5 , 0.25, 0. , 0. ]), 'names': ['x1', 'x2', 'x3', 'x4']} | ||
assert_equal(actual, expected) | ||
|
||
|
||
def test_ff_example(): | ||
''' | ||
''' | ||
|
||
problem = {'bounds': np.repeat([-1, 1], 12).reshape(2, 12).T, | ||
'num_vars': 12, | ||
'names': ["x" + str(x + 1) for x in range(12)] | ||
} | ||
|
||
X = sample(problem) | ||
Y = X[:, 0] + 2 * X[:, 1] + 3 * X[:, 2] + 4 * X[:, 6] * X[:, 11] | ||
|
||
expected = np.array([10, -2, 4, -8, 2, 6, -4, | ||
0, 2, 6, -4, 0, 10, -2, 4, -8, | ||
- 2, -6, 4, 0, -10, 2, -4, 8, | ||
- 10, 2, -4, 8, -2, -6, 4, 0]) | ||
|
||
assert_equal(Y, expected) | ||
|
||
Si = analyze(problem, X, Y) | ||
|
||
expected = np.array([1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=np.float) | ||
assert_equal(expected, Si['ME']) | ||
|
||
|
||
def test_interactions_from_saltelli(): | ||
''' | ||
''' | ||
problem = {'bounds': np.repeat([-1, 1], 12).reshape(2, 12).T, | ||
'num_vars': 12, | ||
'names': ["x" + str(x + 1) for x in range(12)] | ||
} | ||
|
||
X = sample(problem) | ||
|
||
Y = np.array([10, -2, 4, -8, 2, 6, -4, 0, | ||
2, 6, -4, 0, 10, -2, 4, -8, | ||
- 2, -6, 4, 0, -10, 2, -4, 8, | ||
- 10, 2, -4, 8, -2, -6, 4, 0]) | ||
|
||
Si = analyze(problem, X, Y, second_order=True) | ||
actual = Si['IE'] | ||
expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | ||
0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, | ||
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, | ||
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] | ||
|
||
assert_equal(actual, expected) | ||
|
||
|
||
def test_interactions(): | ||
''' | ||
''' | ||
problem = {'bounds': [[0., 2.5], [0., 1.], [0., 1.], [0., 1.]], | ||
'num_vars': 4, | ||
'names': ['x1', 'x2', 'x3', 'x4']} | ||
X = np.array([[ 2.5, 1.0, 1.0, 1.0], | ||
[ 2.5, 0, 1.0, 0], | ||
[ 2.5, 1.0, 0, 0], | ||
[ 2.5, 0, 0, 1.0], | ||
[0, 0, 0, 0], | ||
[0, 1.0, 0, 1.0], | ||
[0, 0, 1.0, 1.0], | ||
[0, 1.0, 1.0, 0]], dtype=np.float) | ||
Y = X[:, 0] + (0.1 * X[:, 1]) + ((1.2 * X[:, 2]) * (1.3 + X[:, 3])) | ||
# Y = np.array([1.5, 1, 1.5, 1, 2, 2.5, 2, 2.5], dtype=np.float) | ||
ie_names, ie = interactions(problem, Y, print_to_console=True) | ||
actual = ie | ||
assert_allclose(actual, [0.3, 0, 0, 0, 0, 0.3], rtol=1e-4, atol=1e-4) |