Coverage for src / sdynpy / fem / sdynpy_dof.py: 12%
48 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Functions for optimally down-selecting degrees of freedom from a candidate set.
4"""
5"""
6Copyright 2022 National Technology & Engineering Solutions of Sandia,
7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
8Government retains certain rights in this software.
10This program is free software: you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation, either version 3 of the License, or
13(at your option) any later version.
15This program is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""
24import numpy as np
25import time
27UPDATE_TIME = 5
30def by_condition_number(sensors_to_keep, shape_matrix, return_condition_numbers=False):
31 '''Get the best set of degrees of freedom by mode shape condition number
33 This function accepts a shape matrix and returns the set of degrees of
34 freedom that corresponds to the lowest condition number.
36 Parameters
37 ----------
38 sensors_to_keep : int
39 The number of sensors to keep.
40 shape_matrix : np.ndarray
41 A 2D or 3D numpy array. If it is a 2D array, the row indices should
42 correspond to the degree of freedom and the column indices should
43 correspond to the mode. For a 3D array, the first index corresponds to
44 a "bundle" of channels (e.g. a triaxial accelerometer) that must be
45 kept or removed as one unit.
46 return_condition_numbers : bool (default False)
47 If True, return a second value that is the condition number at each
48 iteration of the technique.
50 Returns
51 -------
52 indices : np.array
53 A 1d array corresponding to the indices to keep in the first dimension
54 of the shape_matrix array (e.g. new_shape_matrix =
55 shape_matrix[incies,...])
56 returned_condition_numbers : list
57 The condition number at each iteration. Returned only if
58 return_condition_numbers is True.
59 '''
60 shape_matrix = shape_matrix.copy()
61 keep_indices = np.arange(shape_matrix.shape[0])
62 if return_condition_numbers:
63 returned_condition_numbers = [np.linalg.cond(
64 shape_matrix.reshape(-1, shape_matrix.shape[-1]))]
65 start_time = time.time()
66 while shape_matrix.shape[0] > sensors_to_keep:
67 condition_numbers = [np.linalg.cond(shape_matrix[np.arange(shape_matrix.shape[0]) != removed_dof_index, ...].reshape(
68 -1, shape_matrix.shape[-1])) for removed_dof_index in range(shape_matrix.shape[0])]
69 dof_to_remove = np.argmin(condition_numbers)
70# print('Condition Numbers {:}'.format(condition_numbers))
71# print('Removing DOF {:}'.format(dof_to_remove))
72 shape_matrix = np.delete(shape_matrix.copy(), dof_to_remove, axis=0)
73 keep_indices = np.delete(keep_indices, dof_to_remove, axis=0)
74 if return_condition_numbers:
75 returned_condition_numbers.append(condition_numbers[dof_to_remove])
76 new_time = time.time()
77 if new_time - start_time > UPDATE_TIME:
78 print('{:} DoFs Remaining'.format(shape_matrix.shape[0]))
79 start_time += UPDATE_TIME
80 if return_condition_numbers:
81 return keep_indices, returned_condition_numbers
82 else:
83 return keep_indices
86def by_effective_independence(sensors_to_keep, shape_matrix, return_efi=False):
87 '''Get the best set of degrees of freedom by mode shape effective independence
89 This function accepts a shape matrix and returns the set of degrees of
90 freedom that corresponds to the maximum effective independence.
92 Parameters
93 ----------
94 sensors_to_keep : int
95 The number of sensors to keep.
96 shape_matrix : np.ndarray
97 A 2D or 3D numpy array. If it is a 2D array, the row indices should
98 correspond to the degree of freedom and the column indices should
99 correspond to the mode. For a 3D array, the first index corresponds to
100 a "bundle" of channels (e.g. a triaxial accelerometer) that must be
101 kept or removed as one unit.
102 return_efi : bool (default False)
103 If True, return a second value that is the effective independence at
104 each iteration of the technique.
106 Returns
107 -------
108 indices : np.array
109 A 1d array corresponding to the indices to keep in the first dimension
110 of the shape_matrix array (e.g. new_shape_matrix =
111 shape_matrix[incies,...])
112 returned_efi : list
113 The effective independence at each iteration. Returned only if
114 return_efi is set to True.
115 '''
116 shape_matrix = shape_matrix.copy()
117 keep_indices = np.arange(shape_matrix.shape[0])
118 if return_efi:
119 returned_efi = []
120 start_time = time.time()
121 while shape_matrix.shape[0] > sensors_to_keep:
122 Q = shape_matrix.reshape(-1, shape_matrix.shape[-1]
123 ).T @ shape_matrix.reshape(-1, shape_matrix.shape[-1])
124 if return_efi:
125 returned_efi.append(np.linalg.det(Q))
126 Qinv = np.linalg.inv(Q)
127 if shape_matrix.ndim == 2:
128 EfIs = np.diag(shape_matrix @ Qinv @ shape_matrix.T)
129 else:
130 EfIs = 1 - np.linalg.det(np.eye(3) - np.einsum('ijk,kl,iml->ijm',
131 shape_matrix, Qinv, shape_matrix))
132 dof_to_remove = np.argmin(EfIs)
133# print('Effective Independences {:}'.format(EfI3s))
134# print('Removing DOF {:}'.format(dof_to_remove))
135 shape_matrix = np.delete(shape_matrix.copy(), dof_to_remove, axis=0)
136 keep_indices = np.delete(keep_indices, dof_to_remove, axis=0)
137 new_time = time.time()
138 if new_time - start_time > UPDATE_TIME:
139 print('{:} DoFs Remaining'.format(shape_matrix.shape[0]))
140 start_time += UPDATE_TIME
141 if return_efi:
142 return keep_indices, returned_efi
143 else:
144 return keep_indices