romtools.hyper_reduction.deim

  1#
  2# ************************************************************************
  3#
  4#                         ROM Tools and Workflows
  5# Copyright 2019 National Technology & Engineering Solutions of Sandia,LLC
  6#                              (NTESS)
  7#
  8# Under the terms of Contract DE-NA0003525 with NTESS, the
  9# U.S. Government retains certain rights in this software.
 10#
 11# ROM Tools and Workflows is licensed under BSD-3-Clause terms of use:
 12#
 13# Redistribution and use in source and binary forms, with or without
 14# modification, are permitted provided that the following conditions
 15# are met:
 16#
 17# 1. Redistributions of source code must retain the above copyright
 18# notice, this list of conditions and the following disclaimer.
 19#
 20# 2. Redistributions in binary form must reproduce the above copyright
 21# notice, this list of conditions and the following disclaimer in the
 22# documentation and/or other materials provided with the distribution.
 23#
 24# 3. Neither the name of the copyright holder nor the names of its
 25# contributors may be used to endorse or promote products derived
 26# from this software without specific prior written permission.
 27#
 28# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 29# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 30# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 31# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 32# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 33# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 34# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 35# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 36# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 37# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
 38# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 39# POSSIBILITY OF SUCH DAMAGE.
 40#
 41# Questions? Contact Eric Parish (ejparis@sandia.gov)
 42#
 43# ************************************************************************
 44#
 45
 46'''Implementation of DEIM technique for hyper-reduction'''
 47
 48import numpy as np
 49import romtools.linalg.linalg as la
 50
 51
 52
 53def _deim_get_indices_sharedmem(U):
 54    '''
 55    Implementation of the discrete empirical method as described in Algorithm 1 of
 56    S. Chaturantabut and D. C. Sorensen, "Discrete Empirical Interpolation for
 57    nonlinear model reduction," doi: 10.1109/CDC.2009.5400045.
 58
 59    Args:
 60        $\\mathbf{U} \\in \\mathbb{R}^{m \\times n}$, where
 61            $m$ is the number of DOFs and
 62            $n$ is the number of samples.
 63            Function basis in matrix format
 64
 65    Returns:
 66        $\\mathrm{indices} \\in \\mathbb{I}^{n}$: sample mesh indices
 67    '''
 68
 69    m = np.shape(U)[1]
 70    first_index = la.argmax(np.abs(U[:, 0]))
 71    indices = first_index
 72    for ell in range(1, m):
 73        LHS = U[indices, 0:ell]
 74        RHS = U[indices, ell]
 75        if ell == 1:
 76            LHS = np.ones((1, 1))*LHS
 77            RHS = np.ones(1)*RHS
 78        C = np.linalg.solve(LHS, RHS)
 79
 80        residual = U[:, ell] - U[:, 0:ell] @ C
 81        index_to_add = la.argmax(np.abs(residual))
 82        indices = np.append(indices, index_to_add)
 83    return indices
 84
 85
 86class _dist_deim_data:
 87    def __init__(self, i, r):
 88        self.local_indices = np.array([int(i)])
 89        self.owning_ranks = np.array([int(r)])
 90
 91    def append(self, i, r):
 92        '''Adds the local index and rank to self.local_indices and self.owning_ranks.'''
 93        self.local_indices = np.append(self.local_indices , int(i))
 94        self.owning_ranks = np.append(self.owning_ranks, int(r))
 95
 96
 97def _deim_get_indices_distributed(U, comm):
 98    m = np.shape(U)[1]
 99    local_index, foundRank = la.argmax(np.abs(U[:, 0]), comm)
100    result = _dist_deim_data(local_index, foundRank)
101    if m == 1:
102        return result.local_indices, result.owning_ranks
103
104    myRank = comm.Get_rank()
105    LHS, RHS, C = np.array([]), np.array([]), np.array([])
106    for ell in range(1, m):
107        indices = result.local_indices[result.owning_ranks==myRank]
108        LHS = np.array([]) if indices.size == 0 else U[indices, 0:ell]
109        RHS = np.array([]) if indices.size == 0 else U[indices, ell]
110
111        A, b = la.move_distributed_linear_system_to_rank_zero(LHS, RHS, comm)
112        if myRank == 0:
113            C = np.linalg.solve(A, b)
114        C = comm.bcast(C, root=0)
115
116        residual = U[:, ell] - U[:, 0:ell] @ C
117        local_index, foundRank = la.argmax(np.abs(residual), comm)
118        result.append(local_index, foundRank)
119
120    return result.local_indices, result.owning_ranks
121
122
123def deim_get_indices(U, comm=None):
124    '''
125    Implementation of the discrete empirical method as described in Algorithm 1 of
126    S. Chaturantabut and D. C. Sorensen, "Discrete Empirical Interpolation for
127    nonlinear model reduction," doi: 10.1109/CDC.2009.5400045.
128
129    Args:
130        $\\mathbf{U} \\in \\mathbb{R}^{m \\times n}$, where
131            $m$ is the number of DOFs and
132            $n$ is the number of samples.
133            Function basis in matrix format.
134        comm: Optional communicator object. If none, algorithm assumes shared-memory data.
135
136    Returns:
137         if comm==None:
138            $\\mathrm{indices} \\in \\mathbb{I}^{n}$: sample mesh indices
139
140         else:
141            sample mesh local indices and the corresponding owing ranks
142    '''
143
144    assert np.shape(U)[1] >= 1, "deim requires a basis matrix with at least one basis vector (one column)"
145
146    if comm:
147        return _deim_get_indices_distributed(U, comm)
148
149    return _deim_get_indices_sharedmem(U)
150
151
152def _deim_multi_state_get_indices_sharedmem(U):
153    all_indices = np.zeros(0, dtype=int)
154    n_var = U.shape[0]
155    for i in range(0, n_var):
156        data_matrix = U[i]
157        indices = deim_get_indices(data_matrix)
158        all_indices = np.unique(np.append(all_indices, indices))
159    return all_indices
160
161
162def _deim_multi_state_get_indices_distributed(U, comm):
163    all_local_indices = np.zeros(0, dtype=int)
164    all_ranks = np.zeros(0, dtype=int)
165    n_var = U.shape[0]
166    for i in range(0, n_var):
167        data_matrix = U[i]
168        indices, ranks = deim_get_indices(data_matrix, comm)
169        # print(my_rank, ":::", indices, ranks)
170        all_local_indices = np.append(all_local_indices, indices)
171        all_ranks = np.append(all_ranks, ranks)
172
173    # from here we need to operate on ranks and indices together
174    # so stack them for convenience
175    M = np.vstack((all_local_indices, all_ranks))
176
177    # first, sort based on rankID
178    M = np.unique(M, axis=1)
179    inds = M[1,:].argsort()
180    M = M[:, inds]
181
182    # once we are here, M is sorted based on the ranks, but could look like this:
183
184    # M = [ 0 14 13  5 15  4  2  6  8 10  4  4 14  7  8  3  6]
185    #     [ 0  0  0  0  0  0  0  1  1  1  1  2  2  2  2  2  2]
186
187    # we do an additional step where we sort based on the local index within each rank
188
189    rank_ids = np.unique(M[1,:])
190    for rank in rank_ids:
191        locs = np.where(M[1,:] == rank)
192        # only need to sort the 0-th row since we are sorting locally for each rank with same rank id
193        M[0, locs] = np.sort(M[0, locs])
194
195    # return local indices and ranks separately
196    return M[0,:], M[1,:]
197
198
199def multi_state_deim_get_indices(U, comm=None):
200    '''
201    Version of DEIM for multi-state systems.
202
203    We perform DEIM on each state variable, and
204    then return the union of all indices.
205    Repeated indices are removed.
206
207
208    Args:
209         $\\mathbf{U} \\in \\mathbb{R}^{l \\times m \\times n}$, where
210            $l$ is the number of variables,
211            $m$ is the number of DOFs and
212            $n$ the number of samples.
213            Multi-dimensional function basis in tensor format.
214
215    Returns:
216         if comm==None:
217            $\\mathrm{indices} \\in \\mathbb{I}^{n}$: sample mesh indices
218
219         else:
220            sample mesh local indices and the corresponding owing ranks
221    '''
222    shape = np.shape(U)
223    assert len(shape) == 3
224    assert shape[2] >= 1, "deim requires a basis matrix with at least one basis vector (one column)"
225
226    if comm:
227        return _deim_multi_state_get_indices_distributed(U, comm)
228
229    return _deim_multi_state_get_indices_sharedmem(U)
230
231
232
233def deim_get_approximation_matrix(function_basis, sample_indices):
234    '''
235    Given a function basis $\\mathbf{U}$ and sample indices defining $\\mathbf{P}$, we compute
236    $$ \\mathbf{U} \\mathrm{pinv}( \\mathbf{P}^T \\mathbf{U})$$
237    which comprises the matrix needed for the DEIM approximation to $\\mathbf{f}$
238
239    Args:
240        function_basis: ($m$, $n$) array, where $m$ is the number of DOFs and $n$ the number of basis functions.
241            Basis for function to be approximated.
242        sample_indices: ($n_s$, ) array, where $n_s$ is the number of sample points. Sampling points.
243
244    Returns:
245        deim_matrix: ($n$,$n_s$) array. DEIM approximation basis
246    '''
247    sampled_function_basis = function_basis[sample_indices]
248    PU_pinv = np.linalg.pinv(sampled_function_basis)
249    deim_matrix = function_basis @ PU_pinv
250    return deim_matrix
251
252
253def multi_state_deim_get_test_basis(test_basis, function_basis, sample_indices):
254    '''
255    For multistate systems. Constructs an independent DEIM basis for each state variable using uniform sample indices
256    Args:
257        test_basis: ($n_{var}$, $m$, $k$) array, where
258            $n_{var}$ is the number of state variables,
259            $m$ is the number of DOFs and
260            $k$ is the number of basis functions. Test basis in projection scheme
261        function_basis: ($n_{var}$, m, n) array, where
262            $n_{var}$ is the number of state variables,
263            $m$ is the number of DOFs and
264            $n$ is the number of basis functions.
265            Basis for function to be approximated.
266        sample_indices: ($n_s$, ) array, where $n_s$ is the number of sample points. Sampling points.
267
268    Returns:
269        deim_test_basis: ($n_{var}$, $n_s$, $k$) array, where
270            $n_{var}$ is the number of state variables,
271            $n_s$ is the number of sample points, and
272            $k$ the number of basis functions.
273            DEIM test basis matrix.
274
275    '''
276    n_var = function_basis.shape[0]
277    deim_test_basis = deim_get_test_basis(
278        test_basis[0], function_basis[0], sample_indices
279    )
280    deim_test_basis = deim_test_basis[None]
281    for i in range(1, n_var):
282        deim_test_basis_i = deim_get_test_basis(
283            test_basis[i], function_basis[i], sample_indices
284        )
285        deim_test_basis = np.append(deim_test_basis, deim_test_basis_i[None], axis=0)
286    return deim_test_basis
287
288
289def deim_get_test_basis(test_basis, function_basis, sample_indices, comm=None):
290    '''
291    Given a test basis $\\mathbf{\\Phi}$, a function basis $\\mathbf{U}$, and
292    sample indices defining $\\mathbf{P}$, we compute
293    $$[ \\mathbf{\\Phi}^T \\mathbf{U} \\mathrm{pinv}( \\mathbf{P}^T \\mathbf{U}) ]^T$$
294    which comprises the "test basis" for the DEIM approximation for
295    $\\mathbf{\\Phi}^T \\mathbf{f}$
296
297    Args:
298        test_basis: ($m$, $k$) array, where
299            $m$ is the number of DOFs and
300            $k$ is the number of basis functions. Test basis in projection scheme
301        function_basis: ($m$, $n$) array, where
302            $m$ is the number of DOFs and
303            $n$ is the number of basis functions. Basis for function to be approximated.
304        sample_indices: ($n_s$, ) array, where $n_s$ is the number of sample points. Sampling points.
305        comm: Optional communicator object. If none, algorithm assumes shared-memory data.
306
307    Returns:
308        deim_test_basis: ($n_s$, $k$) array, where
309            $n_s$ is the number of sample points and
310            $k$ the number of basis functions. DEIM test basis matrix.
311
312    '''
313    # Determine sampled_function_basis, allowing for empty ranks
314    sampled_function_basis = np.empty((0, function_basis.shape[1]))
315    if len(sample_indices) > 0:
316        sampled_function_basis = function_basis[sample_indices]
317
318    # pinv(P^T U) ^T
319    PU_pinv_transpose = la.pinv(sampled_function_basis, comm=comm)
320
321    # Phi^T U (distributed operation)
322    phi_U = np.empty((test_basis.shape[1], function_basis.shape[1]))
323    la.product("T", "N", 1, test_basis, function_basis, 0, phi_U, comm=comm)
324
325    # Gather pinvs if distributed
326    if comm is not None:
327        local_pinvs = comm.allgather(PU_pinv_transpose)
328        PU_pinv_transpose = np.vstack(local_pinvs)
329
330    # phi_U pinv (local operation)
331    deim_test_basis = np.empty((phi_U.shape[0], PU_pinv_transpose.shape[0]))
332    la.product("N", "T", 1, phi_U, PU_pinv_transpose, 0, deim_test_basis)
333
334    return deim_test_basis.transpose()