Skip to content
Navigation Menu
{{ message }}
forked from nasa/SROMPy
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathSROM.py
More file actions
362 lines (276 loc) · 13.4 KB
/
Copy pathSROM.py
File metadata and controls
362 lines (276 loc) · 13.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
# Copyright 2018 United States Government as represented by the Administrator of
# the National Aeronautics and Space Administration. No copyright is claimed in
# the United States under Title 17, U.S. Code. All Other Rights Reserved.
# The Stochastic Reduced Order Models with Python (SROMPy) platform is licensed
# under the Apache License, Version 2.0 (the "License"); you may not use this
# file except in compliance with the License. You may obtain a copy of the
# License at http://www.apache.org/licenses/LICENSE-2.0.
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
# License for the specific language governing permissions and limitations
# under the License.
"""
Define stochastic reduced order model (SROM) class
"""
import copy
import os
import numpy as np
from SROMPy.optimize import Optimizer
from SROMPy.target.RandomEntity import RandomEntity
class SROM(object):
"""
This is the primary SROMPy class for defining and utilizing a
stochastic reduced order model (SROM). Main capability is optimizing
for the defining SROM parameters to model a given target random quantity.
Other functions provided to calculate SROM statistics, set/get defining
parameters directly, and store/load SROM to/from file.
:param _size: SROM size
:type size: int
:param _dim: dimension of random quantity being modeled
:type dim: int
"""
def __init__(self, size, dim):
"""
Initialize SROM w/ specified size for random vector of dimension dim
m = SROM size, d = dimension;
"""
if size <= 0:
raise(ValueError("SROM size must be greater than zero."))
if dim <= 0:
raise(ValueError("SROM dimension must be greater than 0."))
self._size = int(size)
self._dim = int(dim)
self._samples = None
self._probabilities = None
def set_params(self, samples, probabilities):
"""
Set defining SROM parameters - samples & corresponding probabilities.
:param samples: Array of SROM samples
:type samples: 2d Numpy array, size - (SROM size) x (dim)
:param probabilities: Array of SROM probabilities
:type probabilities: 1d Numpy array, size - (SROM size) x 1
The sample/probability arrays have the following convention (srom sample
index as rows, components of sample as columns):
Samples:
| [[ x_1^(1), x_2^(1), ..., x_d^(1)],
| [x_1^(2), x_2^(2), ..., x_d^(2)],
| ... ... ... ....
| [x_1^(m), x_2^(m), ... x_d^(m)]]
Probabilities:
| [p^(1), p^(2), ..., p^(m)]^T
"""
# Handle 1 dimension case, adjust shape:
if len(samples.shape) == 1:
samples.shape = (len(samples), 1)
# Verify dimensions of samples/probabilities.
(size, dim) = samples.shape
if size != self._size and dim != self._dim:
msg = "SROM samples have wrong dimension, must be (srom_size x dim)"
raise ValueError(msg)
if len(probabilities) != self._size:
raise ValueError("SROM probabilities must have dim. equal to srom "
"size")
self._samples = copy.deepcopy(samples)
self._probabilities = \
copy.deepcopy(probabilities.reshape((self._size, 1)))
def get_params(self):
"""
Returns: tuple of SROM sample & probability arrays. Samples array
has size (SROM size x dim) and probability array has length (SROM size)
The sample/probability arrays have the following convention (srom sample
index as rows, components of sample as columns):
Samples:
| [[ x_1^(1), x_2^(1), ..., x_d^(1)],
| [x_1^(2), x_2^(2), ..., x_d^(2)],
| ... ... ... ....
| [x_1^(m), x_2^(m), ... x_d^(m)]]
Probabilities:
| [p^(1), p^(2), ..., p^(m)]^T
"""
return self._samples, self._probabilities
def compute_moments(self, max_order):
"""
Calculates and returns SROM moments.
:param max_order: Maximum order of moments to return
:type max_order: int
Returns (max_order x dim) size Numpy array with SROM moments for
each dimension.
"""
# Make sure SROM has been properly initialized.
if self._samples is None or self._probabilities is None:
raise ValueError("Must initialize SROM before computing moments")
max_order = int(max_order)
moments = np.zeros((max_order, self._dim))
for q in range(max_order):
# moment_q = sum_{k=1}^m p(k) * x(k)^q.
moment_q = np.zeros((1, self._dim))
for k, sample in enumerate(self._samples):
moment_q = moment_q + self._probabilities[k] * pow(sample, q + 1)
moments[q, :] = moment_q
return moments
def compute_cdf(self, x_grid):
"""
Computes the SROM marginal CDF values in each dimension.
:param x_grid: Grid of points to compute CDF values on. If 1d array is
provided, the same points are used to evaluate CDF in each
dimension. If 2d array is provided, calculates CDF values on
different points, but must have same # points for each dimension.
Size is (# grid pts) x (dim) or (# grid pts) x (1).
:type x_grid: Numpy array.
Returns: Numpy array of CDF values at x_grid points. Size is (# grid
pts) x (dim).
Note:
* Increasing the number of grid points can significantly slow
down the SROM optimization problem.
* Providing a 2d array for x_grid can specify a different range
of values for each dimension, but must use the same number of pts.
"""
# Make sure SROM has been properly initialized
if self._samples is None or self._probabilities is None:
raise ValueError("Must initialize SROM before computing CDF")
if len(x_grid.shape) == 1:
x_grid = x_grid.reshape((len(x_grid), 1))
(num_pts, dim) = x_grid.shape
# If only one grid was provided for multiple dims, repeat to generalize.
if (dim == 1) and (self._dim > 1):
x_grid = np.repeat(x_grid, self._dim, axis=1)
cdf_values = np.zeros((num_pts, self._dim))
# Vectorized indicator implementation for CDF.
# CDF(x) = sum_{k=1}^m 1( sample^(k) < x) prob^(k).
for i, grid in enumerate(x_grid.T):
for k, sample in enumerate(self._samples):
indices = grid >= sample[i]
cdf_values[indices, i] += self._probabilities[k]
return cdf_values
def compute_corr_mat(self):
"""
Returns the SROM correlation matrix as (dim x dim) numpy array
srom_corr = sum_{k=1}^m [ x^(k) * (x^(k))^T ] * p^(k)
"""
# Make sure SROM has been properly initialized
if self._samples is None or self._probabilities is None:
raise ValueError("Must initialize SROM before computing moments")
corr = np.zeros((self._dim, self._dim))
for k, sample in enumerate(self._samples):
corr = corr + np.outer(sample, sample) * self._probabilities[k]
return corr
def optimize(self, target_random_variable,
weights=None,
num_test_samples=50,
error='SSE',
max_moment=5,
cdf_grid_pts=100,
tolerance=None,
options=None,
method=None,
joint_opt=False):
"""
Optimize for the SROM samples & probabilities to best match the
target random vector statistics. The main functionality provided
by the SROM class. Solves SROM the optimization problem and sets
the samples and probabilities for the SROM object to the optimized
values.
:param target_random_variable: the target random quantity
(variable/vector) being modeled by the SROM.
:type target_random_variable: SROMPy target object
(AnalyticRandomVector, SampleRandomVector, or random variable class)
:param weights: relative weights specifying importance of matching
CDFs, moments, and correlation of the target during optimization.
Default is equal weights [1,1,1].
:type weights: 1d Numpy array (length = 3)
:param num_test_samples: Number of sample sets (iterations) to run
optimization.
:type num_test_samples: int
:param error: Type of error metric to use in objective ("SSE", "MAX",
"MEAN").
:type error: string
:param max_moment: Max. number of target moments to consider matching
:type max_moment: int
:param cdf_grid_pts: Number of points to evaluate CDF error on
:type cdf_grid_pts: int
:param tolerance: tolerance for scipy optimization algorithm (TODO)
:type tolerance: float
:param options: scipy optimization algorithm options see scipy
documentation. (TODO)
:type options: dict
:param method: method used for scipy optimization (TODO)
:type method: string
:param joint_opt: Flag to optimize jointly for samples & probabilities.
:type joint_opt: bool
Returns: None. Sets samples/probabilities member variables.
Assumes the targetRV object has been properly initialized beforehand.
The optimization for SROM samples & probabilities is currently
performed sequentially - a random set of samples are first drawn and the
probabilities are then optimization for those samples. The input
"num_test_samples" is the number of random sample sets this is
performed for before terminating. The random sample set and optimal
probabilities found that produce the lowest objective function value
are used as the optimal parameters. The joint_opt input flag can
specify to do the optimization over samples and probabilities
simultaneously.
"""
if not isinstance(target_random_variable, RandomEntity):
raise TypeError("target_random_variable must inherit from "
"RandomEntity.")
# Use optimizer to form SROM objective func & gradient and minimize:
opt = Optimizer(target_random_variable,
self,
weights,
error,
max_moment,
cdf_grid_pts)
(samples, probabilities) = opt.get_optimal_params(num_test_samples,
tolerance,
options,
method,
joint_opt)
self.set_params(samples, probabilities)
def save_params(self, outfile="srom_params.txt", delimiter=' '):
"""
Write the SROM parameters to file.
:param outfile: output file name
:type outfile: string
:param delimiter: delimiter used in output file (default - whitespace)
:type delimiter: string
Returns: None. Produces output file.
Writes output file with the following format (samples in each row with
prob after):
| x_1^(1), x_2^(1), ..., x_d^(1), p^(1)
| x_1^(2), x_2^(2), ..., x_d^(2), p^(2)
| ... ... ... .... ...
| x_1^(m), x_2^(m), ... x_d^(m), p^(m)
"""
# Make sure SROM has been properly initialized
if self._samples is None or self._probabilities is None:
raise ValueError("Must initialize SROM before saving to disk")
srom_params = np.hstack((self._samples, self._probabilities))
np.savetxt(outfile, srom_params, delimiter=delimiter)
def load_params(self, infile="srom_params.txt", delimiter=' '):
"""
Load SROM parameters from file.
:param infile: input file name containing SROM parameters
:type infile: string
:param delimiter: delimiter used in input file (default - whitespace)
:type delimiter: string
Returns: None. Sets sample/probability member variables.
Assumes input file has the following format (samples in each row with
prob after):
| x_1^(1), x_2^(1), ..., x_d^(1), p^(1)
| x_1^(2), x_2^(2), ..., x_d^(2), p^(2)
| ... ... ... .... ...
| x_1^(m), x_2^(m), ... x_d^(m), p^(m)
The dimension of the samples and probabilities arrays must be
compatible with the SROM size and dimension that was used to initialize
the SROM class.
"""
if not os.path.isfile(infile):
raise IOError("SROM parameter input file does not exist: " + infile)
srom_params = np.genfromtxt(infile, delimiter=delimiter)
(size, dim) = srom_params.shape
dim -= 1 # Account for probabilities in last col.
if size != self._size and dim != self._dim:
msg = "Dimension mismatch when loading SROM params from file"
raise ValueError(msg)
self._samples = srom_params[:, :-1]
self._probabilities = srom_params[:, -1]
You can’t perform that action at this time.