Coverage for src / sdynpy / modal / sdynpy_smac.py: 10%
825 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Implementation of the Synthesize Modes and Correlate (SMAC) curve fitter for Python
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
25from ..core.sdynpy_data import TransferFunctionArray, GUIPlot
26from ..core.sdynpy_geometry import Geometry
27from ..core.sdynpy_shape import ShapeArray
28from ..signal_processing.sdynpy_correlation import mac, matrix_plot
29from scipy.signal import find_peaks
30from scipy.ndimage import maximum_filter
31from .sdynpy_modeshape import compute_residues as modeshape_compute_residues, compute_shapes as modeshape_compute_shapes, ShapeSelection
32import matplotlib.pyplot as plt
33import matplotlib.cm as cm
34import os
35from enum import Enum
37from qtpy import QtWidgets, uic, QtGui
38from qtpy.QtWidgets import QApplication, QMainWindow
39import pyqtgraph
40pyqtgraph.setConfigOption('background', 'w')
41pyqtgraph.setConfigOption('foreground', 'k')
44def correlation_coefficient(x, y, axis=-1):
45 mx = np.mean(x, axis=axis, keepdims=True)
46 my = np.mean(y, axis=axis, keepdims=True)
47 xn = x - mx
48 yn = y - my
49 return np.sum(xn * yn, axis=axis) / np.sqrt(np.sum(xn**2, axis=axis) * np.sum(yn**2, axis=axis))
52class ConvergenceException(Exception):
53 pass
56class AutoFitTypes(Enum):
57 PARABOLOID = 0
58 ALTERNATE = 1
61class SMAC:
62 def __init__(self, frfs: TransferFunctionArray, min_frequency=None, max_frequency=None, complex_modes=False,
63 displacement_derivative=2):
64 self.frfs = frfs.reshape_to_matrix()
65 self.min_frequency = min_frequency
66 self.max_frequency = max_frequency
67 abscissa_indices = np.ones(self.frequencies.shape, dtype=bool)
68 if min_frequency is not None:
69 abscissa_indices &= (self.frequencies >= min_frequency)
70 if max_frequency is not None:
71 abscissa_indices &= (self.frequencies <= max_frequency)
72 abscissa = self.frequencies[abscissa_indices]
73 freq_range = np.array((np.min(abscissa), np.max(abscissa)))
74 index_range = np.argmin(np.abs(self.frequencies - freq_range[:, np.newaxis]), axis=1)
75 self.frequency_slice = slice(index_range[0], index_range[1] + 1)
76 self.displacement_derivative = displacement_derivative
77 # Other information
78 self.pinv_frfs = None
79 self.complex_modes = complex_modes
80 self.initial_rootlist = None
81 self.rootlist = None
82 self.residuals = None
83 self.initial_correlation = None
85 def save(self, filename):
86 np.savez(filename, frfs=self.frfs,
87 frequency_slice=self.frequency_slice,
88 pinv_frfs=self.pinv_frfs,
89 complex_modes=self.complex_modes,
90 initial_rootlist=self.initial_rootlist,
91 rootlist=self.rootlist)
93 @property
94 def frequencies(self):
95 return self.frfs[0, 0].abscissa
97 @property
98 def angular_frequencies(self):
99 return 2 * np.pi * self.frequencies
101 @property
102 def reference_coordinates(self):
103 return self.frfs[0, :].reference_coordinate
105 @property
106 def frequency_spacing(self):
107 return np.mean(np.diff(self.frequencies))
109 def __repr__(self):
110 if self.pinv_frfs is None:
111 next_step = "Compute pseudoinverse with compute_pseudoinverse"
112 elif self.initial_rootlist is None:
113 next_step = "Compute initial rootlist with compute_initial_rootlist"
114 elif self.rootlist is None:
115 next_step = "Iterate on the initial rootlist with autofit_roots"
116 else:
117 next_step = 'Analysis Finished'
118 return 'SMAC analysis of {:}. Next step: {:}'.format(self.frfs.__repr__(), next_step)
120 def compute_pseudoinverse(self):
121 ordinate_for_pinv = self.frfs.ordinate[..., self.frequency_slice].transpose(1, 2, 0)
122 if self.complex_modes:
123 self.pinv_frfs = np.linalg.pinv(ordinate_for_pinv)
124 else:
125 self.pinv_frfs = np.linalg.pinv(np.concatenate(
126 (ordinate_for_pinv.real, ordinate_for_pinv.imag), axis=1))
128 def compute_correlation_matrix(self, low_frequency=None, high_frequency=None,
129 frequency_samples=None, frequency_resolution=None,
130 low_damping=0.0025, high_damping=0.05, damping_samples=21,
131 frequency_lines_for_correlation=20, plot=False):
132 if self.pinv_frfs is None:
133 raise RuntimeError(
134 'Pseudoinverse must be calculated first by calling compute_pseudoinverse')
135 if low_frequency is None:
136 low_frequency = self.frequencies[self.frequency_slice].min()
137 if high_frequency is None:
138 high_frequency = self.frequencies[self.frequency_slice].max()
139 if frequency_samples is None:
140 if frequency_resolution is None:
141 frequency_samples = int(
142 np.round((high_frequency - low_frequency) / self.frequency_spacing)) + 1
143 else:
144 frequency_samples = int(
145 np.round((high_frequency - low_frequency) / frequency_resolution)) + 1
146 correlation_frequencies = np.linspace(low_frequency, high_frequency, frequency_samples)
147 correlation_dampings = np.linspace(low_damping, high_damping, damping_samples)
148 # Set up for a big broadcast
149 omega_frequency_line = self.angular_frequencies[self.frequency_slice, np.newaxis]
150 omega_correlation_frequencies = 2 * np.pi * \
151 correlation_frequencies[:, np.newaxis, np.newaxis, np.newaxis]
152 zeta = correlation_dampings[:, np.newaxis, np.newaxis]
153 frfs_sdof = ((1j * omega_frequency_line)**self.displacement_derivative /
154 (omega_correlation_frequencies**2
155 + 1j * 2 * zeta * omega_frequency_line * omega_correlation_frequencies
156 - omega_frequency_line**2))
157 if self.complex_modes:
158 psi = self.pinv_frfs[:, np.newaxis, np.newaxis] @ frfs_sdof
159 else:
160 psi = self.pinv_frfs[:, np.newaxis,
161 np.newaxis] @ np.concatenate((frfs_sdof.real, frfs_sdof.imag), axis=2)
162 frfs_sdof_reprojected = self.frfs.ordinate.transpose(
163 1, 2, 0)[:, np.newaxis, np.newaxis, self.frequency_slice, :] @ psi
164 correlation_matrix = np.zeros(frfs_sdof_reprojected.shape[:3])
165 for i, frequency in enumerate(correlation_frequencies):
166 closest_index = np.argmin(np.abs(self.frequencies[self.frequency_slice] - frequency))
167 abscissa_slice = slice(max([0, closest_index - frequency_lines_for_correlation]),
168 closest_index + frequency_lines_for_correlation) # Note that this will automatically stop at the end of the array
169 correlation_matrix[:, i] = correlation_coefficient(np.abs(frfs_sdof[i, :, abscissa_slice, 0]),
170 np.abs(frfs_sdof_reprojected[:, i, :, abscissa_slice, 0]))
171 if plot:
172 X, Y = np.meshgrid(correlation_dampings, correlation_frequencies)
173 for i, (matrix, coordinate) in enumerate(zip(correlation_matrix, self.frfs[0, :].reference_coordinate)):
174 if np.all([shape > 1 for shape in matrix.shape]):
175 fig, ax = plt.subplots(2, 1, num='Reference {:} Correlation'.format(
176 str(coordinate)), gridspec_kw={'height_ratios': [3, 1]})
177 mat2plot = np.log10(1 - matrix)
178 surf = ax[0].contourf(X, Y, -mat2plot, -np.linspace(-1, mat2plot.min(), 20))
179 c = fig.colorbar(surf, ax=ax) # , shrink=0.5, aspect=5)
180 c.ax.set_yticklabels(['{:0.4f}'.format(1 - 10**(-x))
181 for x in c.ax.get_yticks()])
182 ax[1].plot(correlation_frequencies, np.max(matrix, axis=1))
183 ax[1].set_ylim(0, 1)
184 else:
185 fig, ax = plt.subplots(
186 1, 1, num='Reference {:} Correlation'.format(str(coordinate)))
187 ax.plot(correlation_frequencies, np.max(matrix, axis=1))
188 ax.set_ylim(0, 1)
190 return correlation_frequencies, correlation_dampings, correlation_matrix
192 def find_peaks(self, correlation_matrix, size=3, threshold=0.9):
193 return [np.where((maximum_filter(matrix, size=size) == matrix) & (matrix > threshold)) for matrix in correlation_matrix]
195 def compute_initial_rootlist(self,
196 frequency_samples=None, frequency_resolution=None,
197 low_damping=0.0025, high_damping=0.05, damping_samples=21,
198 frequency_lines_for_correlation=20, peak_finder_filter_size=3,
199 correlation_threshold=0.9, num_roots_mif='cmif',
200 num_roots_frequency_threshold=0.005, plot_correlation=False):
201 self.initial_correlation_frequencies, self.initial_correlation_dampings, self.initial_correlation_matrix = self.compute_correlation_matrix(
202 None, None,
203 frequency_samples, frequency_resolution,
204 low_damping, high_damping, damping_samples,
205 frequency_lines_for_correlation, plot=plot_correlation
206 )
207 peaklists = self.find_peaks(self.initial_correlation_matrix,
208 peak_finder_filter_size, correlation_threshold)
209 # Collapse identical roots
210 root_index_list = {}
211 for reference_index, (matrix, peaklist) in enumerate(zip(self.initial_correlation_matrix, peaklists)):
212 for key in zip(*peaklist):
213 value = matrix[key]
214 if key not in root_index_list or root_index_list[key][0] < value:
215 root_index_list[key] = (value, reference_index)
216 root_list = np.ndarray(len(root_index_list), dtype=[('frequency', 'float64'),
217 ('damping', 'float64'),
218 ('correlation', 'float64'),
219 ('reference_index', 'int'),
220 ('num_roots', 'int')])
221 for index, (frequency_index, damping_index) in enumerate(sorted(root_index_list)):
222 frequency = self.initial_correlation_frequencies[frequency_index]
223 damping = self.initial_correlation_dampings[damping_index]
224 root_list[index] = (frequency, damping,) + \
225 root_index_list[frequency_index, damping_index] + (0,)
226 root_list['num_roots'] = self.get_num_roots(
227 root_list['frequency'], num_roots_mif, num_roots_frequency_threshold)
228 root_list = root_list[root_list['num_roots'] > 0]
229 self.initial_rootlist = root_list
231 def get_num_roots(self, frequencies, mif_type, frequency_threshold=0.005, plot=False):
232 # Find the closest index for each frequency in the data
233 abscissa = self.frequencies
234 differences = np.abs(abscissa - frequencies[:, np.newaxis])
235 frequency_indices = np.argmin(differences, axis=-1)
236 nroots = np.zeros(len(frequencies), dtype=int)
237 # MMIF
238 if mif_type.lower() == 'mmif':
239 mmif = self.frfs.compute_mmif()
240 mmif_inv = 1 / mmif
241 peaks = [find_peaks(mmif_inv.ordinate[i])[0]
242 for i in range(self.reference_coordinates.size)]
243 # Now see if we are within the tolerance of a peak
244 for i in range(self.reference_coordinates.size):
245 for j, (freq, index) in enumerate(zip(frequencies, frequency_indices)):
246 argument = np.abs(abscissa[peaks[i]] - freq)
247 peak_index = np.argmin(argument)
248 diff = argument[peak_index]
249 if diff / freq < frequency_threshold or diff / self.frequency_spacing < 1.25: # Gives it the ability to be one frequency line off
250 # Check if it is repeated
251 if i > 0: # check if we are on our second mif line or greater
252 # To be repeated, the higher MIF must drop to at least 60% of the principal mif value
253 # The MMIF value should also be below 0.9
254 # The primary mif must also have a peak -- I don't think this is true
255 if ((1 - mmif[i].ordinate[peaks[i][peak_index]]) / (1 - mmif[i - 1].ordinate[peaks[i][peak_index]]) >= 0.6
256 and mmif[i].ordinate[peaks[i][peak_index]] < 0.9
257 # and nroots[j] > 0
258 ): # noqa: E124
259 nroots[j] += 1
260 else:
261 nroots[j] += 1
262 if plot:
263 ax = mmif.plot()
264 for nroot, freq_ind in zip(nroots, frequency_indices):
265 ax.plot(abscissa[freq_ind], mmif.ordinate[0, freq_ind], 'r*')
266 ax.text(abscissa[freq_ind], mmif.ordinate[0, freq_ind], '{:}'.format(nroot))
267 # CMIF
268 elif mif_type.lower() == 'cmif':
269 cmif = self.frfs.compute_cmif()
270 peaks = [find_peaks(cmif.ordinate[i])[0]
271 for i in range(self.reference_coordinates.size)]
272 # Now see if we are within the tolerance of a peak
273 for i in range(self.reference_coordinates.size):
274 for j, (freq, index) in enumerate(zip(frequencies, frequency_indices)):
275 argument = np.abs(abscissa[peaks[i]] - freq)
276 peak_index = np.argmin(argument)
277 diff = argument[peak_index]
278 if diff / freq < frequency_threshold or diff / self.frequency_spacing < 1.25: # Gives it the ability to be one frequency line off
279 # Check if it is repeated
280 if i > 0: # check if we are on our second mif line or greater
281 # To be repeated, the lower MIF must be at least 10% of the primary
282 # The primary mif must also have a peak - I don't think this is true
283 if ((cmif[i].ordinate[peaks[i][peak_index]]) / (cmif[i - 1].ordinate[peaks[i][peak_index]]) >= 0.1
284 # and nroots[j] > 0
285 ): # noqa: E124
286 nroots[j] += 1
287 else:
288 nroots[j] += 1
289 if plot:
290 ax = cmif.plot()
291 for nroot, freq_ind in zip(nroots, frequency_indices):
292 ax.plot(abscissa[freq_ind], cmif.ordinate[0, freq_ind], 'r*')
293 ax.text(abscissa[freq_ind], cmif.ordinate[0, freq_ind], '{:}'.format(nroot))
294 return nroots
296 def autofit_roots(self, frequency_range=0.01, frequency_points=21, frequency_convergence=0.00025,
297 damping_low=0.0025, damping_high=0.05, damping_points=21, damping_convergence=0.02,
298 frequency_lines_for_correlation=20, max_iter=200,
299 zoom_rate=0.75, plot_convergence=False,
300 autofit_type=AutoFitTypes.ALTERNATE):
301 rootlist_dtype = [('frequency', 'float64'),
302 ('damping', 'float64'),
303 ('correlation', 'float64'),
304 ('initial_frequency', 'float64'),
305 ('initial_damping', 'float64'),
306 ('initial_correlation', 'float64'),
307 ('shape_reference', 'int'),
308 ('drive_point_coefficient', 'int')]
309 root_list = np.ndarray((0,), rootlist_dtype)
310 for i, (initial_frequency,
311 initial_damping,
312 initial_correlation,
313 initial_reference_index,
314 initial_num_roots) in enumerate(self.initial_rootlist):
315 print('Fitting peak {:}: {:0.2f} Hz {:0.2f} % Damping'.format(
316 i, initial_frequency, initial_damping * 100))
317 if initial_num_roots > 1:
318 print(' Multiple roots detected. Splitting up peak into multiple seeds.')
319 initial_frequencies = initial_frequency * \
320 (1 + frequency_range / 2 * np.array((-1, 0, 1)))
321 initial_dampings = np.ones(3) * np.mean([damping_low, damping_high])
322 for initial_frequency, initial_damping in zip(initial_frequencies, initial_dampings):
323 try:
324 if autofit_type == AutoFitTypes.ALTERNATE:
325 converged_frequency, converged_damping, converged_correlation = self.autofit_root_alternate(
326 initial_frequency, initial_damping, frequency_range,
327 frequency_points, frequency_convergence, damping_low,
328 damping_high, damping_points, damping_convergence,
329 frequency_lines_for_correlation, max_iter, zoom_rate,
330 plot_convergence)
331 elif autofit_type == AutoFitTypes.PARABOLOID:
332 converged_frequency, converged_damping, converged_correlation = self.autofit_root_paraboloid(
333 initial_frequency, initial_damping, frequency_range,
334 frequency_points, frequency_convergence, damping_low,
335 damping_high, damping_points, damping_convergence,
336 frequency_lines_for_correlation, max_iter, zoom_rate,
337 plot_convergence)
338 except ConvergenceException:
339 print(' Root diverged. Skipping.')
340 continue
341 # Check if the root is already in there
342 same_roots = np.where(((np.abs(root_list['frequency'] - converged_frequency) / converged_frequency < frequency_convergence)
343 & (np.abs(root_list['damping'] - converged_damping) / converged_damping < damping_convergence)))[0]
344 if same_roots.size > 0:
345 print(' Root already exists.')
346 # Just grab the first one?
347 same_roots = same_roots[0]
348 previous_correlation = root_list['correlation'][same_roots]
349 if converged_correlation > previous_correlation:
350 # Replace the previous one
351 print(' Replacing previous root')
352 root_list[same_roots] = (converged_frequency,
353 converged_damping,
354 converged_correlation,
355 initial_frequency,
356 initial_damping,
357 initial_correlation,
358 0, 0)
359 # Otherwise we don't add it
360 else:
361 print(' Not adding this root.')
362 else:
363 print(' Appending to root list')
364 root_list = np.append(root_list, np.array((converged_frequency,
365 converged_damping,
366 converged_correlation,
367 initial_frequency,
368 initial_damping,
369 initial_correlation,
370 0, 0), rootlist_dtype))
371 else:
372 try:
373 if autofit_type == AutoFitTypes.ALTERNATE:
374 converged_frequency, converged_damping, converged_correlation = self.autofit_root_alternate(
375 initial_frequency, initial_damping, frequency_range,
376 frequency_points, frequency_convergence, damping_low,
377 damping_high, damping_points, damping_convergence,
378 frequency_lines_for_correlation, max_iter, zoom_rate,
379 plot_convergence)
380 elif autofit_type == AutoFitTypes.PARABOLOID:
381 converged_frequency, converged_damping, converged_correlation = self.autofit_root_paraboloid(
382 initial_frequency, initial_damping, frequency_range,
383 frequency_points, frequency_convergence, damping_low,
384 damping_high, damping_points, damping_convergence,
385 frequency_lines_for_correlation, max_iter, zoom_rate,
386 plot_convergence)
387 except ConvergenceException:
388 print(' Root diverged. Skipping.')
389 continue
390 # Check if the root is already in there
391 same_roots = np.where(((np.abs(root_list['frequency'] - converged_frequency) / converged_frequency < frequency_convergence)
392 & (np.abs(root_list['damping'] - converged_damping) / converged_damping < damping_convergence)))[0]
393 if same_roots.size > 0:
394 print(' Root already exists.')
395 # Just grab the first one?
396 same_roots = same_roots[0]
397 previous_correlation = root_list['correlation'][same_roots]
398 if converged_correlation > previous_correlation:
399 # Replace the previous one
400 print(' Replacing previous root')
401 root_list[same_roots] = (converged_frequency,
402 converged_damping,
403 converged_correlation,
404 initial_frequency,
405 initial_damping,
406 initial_correlation,
407 0, 0)
408 # Otherwise we don't add it
409 else:
410 print(' Not adding this root.')
411 else:
412 print(' Appending to root list')
413 root_list = np.append(root_list, np.array((converged_frequency,
414 converged_damping,
415 converged_correlation,
416 initial_frequency,
417 initial_damping,
418 initial_correlation,
419 0, 0), rootlist_dtype))
420 self.rootlist = np.sort(root_list)
422 def autofit_root_paraboloid(self, initial_frequency, initial_damping,
423 frequency_range=0.01, frequency_points=21,
424 frequency_convergence=0.00025,
425 damping_low=0.0025, damping_high=0.05,
426 damping_points=21, damping_convergence=0.02,
427 frequency_lines_for_correlation=20, max_iter=200,
428 zoom_rate=0.75, plot_convergence=False):
429 niter = 0
430 frequency_bounds = initial_frequency * (1 + frequency_range / 2 * np.array([-1, 1]))
431 damping_bounds = np.array([damping_low, damping_high])
432 last_frequency = initial_frequency
433 last_damping = initial_damping
434 frequency_range = frequency_bounds[1] - frequency_bounds[0]
435 damping_range = damping_bounds[1] - damping_bounds[0]
436 if plot_convergence:
437 fig = plt.figure('Root {:.2f} Hz, {:.2f} % Convergence'.format(initial_frequency, 100 * initial_damping),
438 figsize=plt.figaspect(1 / 3))
439 fig.clf()
440 ax_3d = fig.add_subplot(1, 3, 1, projection='3d')
441 ax_3d.set_title('Correlation Coefficient')
442 ax_3d.set_xlabel('Damping')
443 ax_3d.set_ylabel('Frequency')
444 ax_freq = fig.add_subplot(1, 3, 2)
445 ax_freq.set_title('Frequency Convergence')
446 ax_freq.set_ylabel('Frequency')
447 ax_damp = fig.add_subplot(1, 3, 3)
448 ax_damp.set_title('Damping Convergence')
449 ax_damp.set_ylabel('Damping')
450 ax_freq.plot(niter, last_frequency, 'b.')
451 ax_damp.plot(niter, last_damping, 'b.')
452 convergence_count = 0
453 while niter < max_iter:
454 # Compute the correlation matrix
455 freq, damp, matrix = self.compute_correlation_matrix(
456 low_frequency=frequency_bounds[0],
457 high_frequency=frequency_bounds[1],
458 frequency_samples=frequency_points,
459 low_damping=damping_bounds[0],
460 high_damping=damping_bounds[1],
461 damping_samples=damping_points,
462 frequency_lines_for_correlation=frequency_lines_for_correlation
463 )
464 x, y = np.meshgrid(damp, freq)
465 z = matrix
466 peak_matrix, peak_frequency, peak_damping = np.unravel_index(np.argmax(z), z.shape)
467 peak_corr = z[peak_matrix, peak_frequency, peak_damping]
468 peak_damping = damp[peak_damping]
469 peak_frequency = freq[peak_frequency]
470 (next_damping, next_frequency, corr_coef), point_fits = self.fit_paraboloid(x, y, matrix)
471 best_correlation_index = np.argmax(corr_coef)
472 next_damping = next_damping[best_correlation_index]
473 next_frequency = next_frequency[best_correlation_index]
474 next_correlation = corr_coef[best_correlation_index]
475 if plot_convergence:
476 ax_3d.cla()
477 ax_3d.plot_surface(x, y, z[best_correlation_index],
478 color='b', linewidth=0, alpha=0.5)
479 ax_3d.plot_surface(point_fits[0].reshape(x.shape),
480 point_fits[1].reshape(y.shape),
481 point_fits[2].reshape(z.shape)[best_correlation_index],
482 color='r', linewidth=0, alpha=0.5)
483 ax_3d.plot(next_damping, next_frequency, corr_coef[best_correlation_index], 'k.')
484 # Check if it is within the bounds
485 if (next_damping > damping_bounds[1] or next_damping < damping_bounds[0]):
486 next_damping = peak_damping
487 else:
488 damping_range *= zoom_rate
489 if (next_frequency > frequency_bounds[1] or next_frequency < frequency_bounds[0]):
490 next_frequency = peak_frequency
491 else:
492 frequency_range *= zoom_rate
493 frequency_bounds = np.array((next_frequency - frequency_range / 2,
494 next_frequency + frequency_range / 2))
495 damping_bounds = np.array((next_damping - damping_range / 2,
496 next_damping + damping_range / 2))
497 niter += 1
498 if plot_convergence:
499 ax_freq.plot(niter, next_frequency, 'b.')
500 ax_damp.plot(niter, next_damping, 'b.')
501 plt.pause(0.001)
502 if (np.abs(next_damping - last_damping) / last_damping < damping_convergence and
503 np.abs(next_frequency - last_frequency) / last_frequency < frequency_convergence):
504 convergence_count += 1
505 if convergence_count > 4:
506 break
507 else:
508 convergence_count = 0
509 if ( # np.abs(next_damping - initial_damping)/initial_damping > 10 or
510 np.abs(next_frequency - initial_frequency) / initial_frequency > 0.1 or
511 next_damping < 0 or
512 next_frequency < 0):
513 print('Diverged: Damping Change {:}\nFrequency Change {:}'.format(
514 np.abs(next_damping - initial_damping) / initial_damping,
515 np.abs(next_frequency - initial_frequency) / initial_frequency))
516 raise ConvergenceException('Peak Diverged!')
517 last_damping = next_damping
518 last_frequency = next_frequency
519 # break
520 return next_frequency, next_damping, next_correlation
522 def autofit_root_alternate(self, initial_frequency, initial_damping,
523 frequency_range=0.01, frequency_points=21,
524 frequency_convergence=0.00025,
525 damping_low=0.0025, damping_high=0.05,
526 damping_points=21, damping_convergence=0.02,
527 frequency_lines_for_correlation=20, max_iter=200,
528 zoom_rate=0.75, plot_convergence=False):
529 niter = 0
530 frequency_bounds = initial_frequency * (1 + frequency_range / 2 * np.array([-1, 1]))
531 damping_bounds = np.array([damping_low, damping_high])
532 last_frequency = initial_frequency
533 last_damping = initial_damping
534 frequency_range = frequency_bounds[1] - frequency_bounds[0]
535 damping_range = damping_bounds[1] - damping_bounds[0]
536 if plot_convergence:
537 fig = plt.figure('Root {:.2f} Hz, {:.2f} % Convergence'.format(initial_frequency, 100 * initial_damping),
538 figsize=plt.figaspect(1 / 3))
539 fig.clf()
540 ax_freq_fit = fig.add_subplot(2, 2, 1)
541 ax_freq_fit.set_title('Frequency Correlation')
542 ax_freq_fit.set_xlabel('Frequency')
543 ax_freq_fit.set_ylabel('Correlation')
544 ax_damp_fit = fig.add_subplot(2, 2, 2)
545 ax_damp_fit.set_title('Damping Correlation')
546 ax_damp_fit.set_xlabel('Damping')
547 ax_damp_fit.set_ylabel('Correlation')
548 ax_freq = fig.add_subplot(2, 2, 3)
549 ax_freq.set_title('Frequency Convergence')
550 ax_freq.set_ylabel('Frequency')
551 ax_damp = fig.add_subplot(2, 2, 4)
552 ax_damp.set_title('Damping Convergence')
553 ax_damp.set_ylabel('Damping')
554 ax_freq.plot(niter, last_frequency, 'b.')
555 ax_damp.plot(niter, last_damping, 'b.')
556 convergence_count = 0
557 while niter < max_iter:
558 next_frequency, freqs_checked, freq_matrix = self.fit_frequency(*frequency_bounds, last_damping,
559 frequency_points, frequency_lines_for_correlation)
560 next_damping, damps_checked, damp_matrix = self.fit_damping(*damping_bounds, last_frequency,
561 damping_points, frequency_lines_for_correlation)
562 next_correlation = np.max([np.max(freq_matrix), np.max(damp_matrix)])
563 if plot_convergence:
564 ax_freq_fit.cla()
565 ax_freq_fit.plot(freqs_checked, freq_matrix.T)
566 ax_damp_fit.cla()
567 ax_damp_fit.plot(damps_checked, damp_matrix.T)
568 # Check if it is within the bounds
569 if (next_damping <= damping_bounds[1] and next_damping >= damping_bounds[0]):
570 damping_range *= zoom_rate
571 if (next_frequency <= frequency_bounds[1] or next_frequency >= frequency_bounds[0]):
572 frequency_range *= zoom_rate
573 frequency_bounds = np.array((next_frequency - frequency_range / 2,
574 next_frequency + frequency_range / 2))
575 damping_bounds = np.array((next_damping - damping_range / 2,
576 next_damping + damping_range / 2))
577 niter += 1
578 if plot_convergence:
579 ax_freq.plot(niter, next_frequency, 'b.')
580 ax_damp.plot(niter, next_damping, 'b.')
581 plt.pause(0.001)
582 if (np.abs(next_damping - last_damping) / last_damping < damping_convergence and
583 np.abs(next_frequency - last_frequency) / last_frequency < frequency_convergence):
584 convergence_count += 1
585 if convergence_count > 4:
586 break
587 else:
588 convergence_count = 0
589 if ( # np.abs(next_damping - initial_damping)/initial_damping > 10 or
590 np.abs(next_frequency - initial_frequency) / initial_frequency > 0.1 or
591 next_damping < 0 or
592 next_frequency < 0):
593 print('Diverged: Damping Change {:}\nFrequency Change {:}'.format(
594 np.abs(next_damping - initial_damping) / initial_damping,
595 np.abs(next_frequency - initial_frequency) / initial_frequency))
596 raise ConvergenceException('Peak Diverged!')
597 last_damping = next_damping
598 last_frequency = next_frequency
599 return next_frequency, next_damping, next_correlation
601 def fit_frequency(self, min_freq, max_freq, damping, frequency_points=21,
602 frequency_lines_for_correlation=20):
603 freq, damp, matrix = self.compute_correlation_matrix(
604 low_frequency=min_freq, high_frequency=max_freq,
605 frequency_samples=frequency_points,
606 low_damping=damping, high_damping=damping, damping_samples=1,
607 frequency_lines_for_correlation=frequency_lines_for_correlation, plot=False)
608 # Find the maximum frequency
609 max_freq_ind = np.argmax(np.max(matrix[..., 0], 0))
610 max_freq = freq[max_freq_ind]
611 return max_freq, freq, matrix[..., 0]
613 def fit_damping(self, min_damp, max_damp, frequency, damping_points=21,
614 frequency_lines_for_correlation=20):
615 freq, damp, matrix = self.compute_correlation_matrix(
616 low_frequency=frequency, high_frequency=frequency,
617 frequency_samples=1,
618 low_damping=min_damp, high_damping=max_damp, damping_samples=damping_points,
619 frequency_lines_for_correlation=frequency_lines_for_correlation, plot=False)
620 # Find the maximum frequency
621 max_damping_ind = np.argmax(np.max(matrix[..., 0, :], 0))
622 max_damp = damp[max_damping_ind]
623 return max_damp, damp, matrix[..., 0, :]
625 def compute_residues(self, roots,
626 residuals=True, weighting='magnitude'):
627 self.natural_frequencies = roots['frequency']
628 self.damping_ratios = roots['damping']
629 self.residues, self.frfs_synth_residue, self.frfs_synth_residual = (
630 modeshape_compute_residues(
631 self.frfs, self.natural_frequencies, self.damping_ratios,
632 not self.complex_modes, residuals, self.min_frequency, self.max_frequency,
633 weighting, self.displacement_derivative))
635 def compute_shapes(self):
636 self.shapes, negative_drive_points = modeshape_compute_shapes(
637 self.natural_frequencies, self.damping_ratios,
638 self.frfs.coordinate, self.residues, ShapeSelection.DRIVE_POINT_COEFFICIENT)
639 self.negative_drive_points = negative_drive_points
641 def frf_sdof_real(self, frequencies, root_frequencies, root_dampings):
642 omega = 2 * np.pi * frequencies
643 omega_n = 2 * np.pi * root_frequencies
644 zeta = root_dampings
645 return -omega[:, np.newaxis]**2 / (-omega[:, np.newaxis]**2 + omega_n**2 + 2 * 1j * zeta * omega_n * omega[:, np.newaxis])
647 def frf_sdof_complex(self, frequencies, root_frequencies, root_dampings):
648 raise NotImplementedError('Complex Shapes are not implemented yet!')
650 def fit_paraboloid(self, x, y, z):
651 x = np.array(x)
652 original_x_shape = x.shape
653 x = x.flatten()
654 y = np.array(y).flatten()
655 z = np.array(z)
656 new_z_shape = z.shape[:-len(original_x_shape)] + (np.prod(x.shape), 1)
657 z = z.reshape(*new_z_shape)
658 # Now do some normalization
659 x_range = x.max() - x.min()
660 x_mean = (x.max() + x.min()) / 2
661 y_range = y.max() - y.min()
662 y_mean = (y.max() + y.min()) / 2
663 # z_range = z.max() - z.min()
664 # z_mean = (z.max() + z.min())/2
665 # x_range = 1
666 # x_mean = 0
667 # y_range = 1
668 # y_mean = 0
669 z_range = 1
670 z_mean = 0
671 x = (x - x_mean) / x_range
672 y = (y - y_mean) / y_range
673 z = (z - z_mean) / z_range
674 A = np.array((
675 x**2,
676 y**2,
677 x * y,
678 x,
679 y,
680 np.ones(x.shape)
681 )).T
682 a, b, c, d, e, f = (np.linalg.pinv(A) @ z)[..., 0].T
683 x_c = (-2 * b * d + c * e) / (4 * a * b - c**2)
684 y_c = (-2 * a * e + c * d) / (4 * a * b - c**2)
685 z_c = (a * x_c**2 + b * y_c**2 + c * x_c * y_c + d * x_c + e * y_c + f)
686 x_c = x_c * x_range + x_mean
687 y_c = y_c * y_range + y_mean
688 z_c = z_c * z_range + z_mean
689 x_fit = x * x_range + x_mean
690 y_fit = y * y_range + y_mean
691 z_fit = (a[:, np.newaxis] * x**2 + b[:, np.newaxis] * y**2 + c[:, np.newaxis] * x *
692 y + d[:, np.newaxis] * x + e[:, np.newaxis] * y + f[:, np.newaxis]) # *z_range+z_mean
693 return (x_c, y_c, z_c), (x_fit, y_fit, z_fit)
696class AddRootDialog(QtWidgets.QDialog):
697 def __init__(self, parent):
698 super(QtWidgets.QDialog, self).__init__(parent)
699 uic.loadUi(os.path.join(os.path.abspath(os.path.dirname(
700 os.path.abspath(__file__))), 'smac_add_root.ui'), self)
701 self.setWindowTitle('PySMAC: Add Root')
703 self.low_damping.setValue(parent.minimum_damping_autofit_selector.value())
704 self.high_damping.setValue(parent.maximum_damping_autofit_selector.value())
705 self.low_frequency.setValue(parent.low_frequency)
706 self.high_frequency.setValue(parent.high_frequency)
707 for widget in [self.low_frequency, self.high_frequency]:
708 widget.setMinimum(parent.low_frequency)
709 widget.setMaximum(parent.high_frequency)
710 self.frequency_samples.setValue(parent.frequency_autofit_values_in_range.value())
711 self.damping_samples.setValue(parent.damping_autofit_number_samples_selector.value())
713 self.correlation_plot = self.correlation_view.getPlotItem()
714 self.correlation_plot.setLimits(xMin=parent.low_frequency,
715 xMax=parent.high_frequency,
716 yMin=0.001,
717 yMax=99.99)
718 self.correlation_plot.setRange(xRange=(parent.low_frequency, parent.high_frequency),
719 yRange=(parent.minimum_damping_autofit_selector.value(),
720 parent.maximum_damping_autofit_selector.value()), padding=0.0)
721 self.correlation_plot.getAxis('left').setLabel('Damping (%)')
722 self.correlation_plot.getAxis('bottom').setLabel('Frequency (Hz)')
723 self.SMAC_GUI = parent
725 # Compute initial correlation
726 self.image = None
727 self.image_log_scale = False
728 self.recompute_correlation()
730 self.max_correlation = None
731 self.max_freq = None
732 self.max_damp = None
734 self.connect_callbacks()
735 print('Completed init')
737 def connect_callbacks(self):
738 self.recompute_correlation_button.clicked.connect(self.recompute_correlation)
739 self.linear_plot.toggled.connect(self.switch_log_plot)
740 self.correlation_plot.sigRangeChanged.connect(self.update_range_selectors)
741 self.low_frequency.valueChanged.connect(self.update_plot_range)
742 self.high_frequency.valueChanged.connect(self.update_plot_range)
743 self.low_damping.valueChanged.connect(self.update_plot_range)
744 self.high_damping.valueChanged.connect(self.update_plot_range)
746 def update_range_selectors(self):
747 self.low_frequency.blockSignals(True)
748 self.high_frequency.blockSignals(True)
749 self.low_damping.blockSignals(True)
750 self.high_damping.blockSignals(True)
751 ((xmin, xmax), (ymin, ymax)) = self.correlation_plot.vb.viewRange()
752 self.low_frequency.setValue(xmin)
753 self.high_frequency.setValue(xmax)
754 self.low_damping.setValue(ymin)
755 self.high_damping.setValue(ymax)
756 self.low_frequency.blockSignals(False)
757 self.high_frequency.blockSignals(False)
758 self.low_damping.blockSignals(False)
759 self.high_damping.blockSignals(False)
761 def update_plot_range(self):
762 self.correlation_plot.blockSignals(True)
763 lf = self.low_frequency.value()
764 hf = self.high_frequency.value()
765 ld = self.low_damping.value()
766 hd = self.high_damping.value()
767 self.correlation_plot.setRange(xRange=(lf, hf),
768 yRange=(ld, hd), padding=0.0)
769 self.correlation_plot.blockSignals(False)
771 def switch_log_plot(self):
772 if self.image_log_scale and self.linear_plot.isChecked():
773 # Undo log scale back to linear
774 self.image.setImage((-(10**(-self.image.image))) + 1)
775 self.image_log_scale = False
776 if not self.image_log_scale and self.log_plot.isChecked():
777 # put into log scale
778 self.image.setImage(-np.log10(1 - self.image.image))
779 self.image_log_scale = True
781 def recompute_correlation(self):
782 self.correlation_plot.clear()
783 lf = self.low_frequency.value()
784 hf = self.high_frequency.value()
785 ld = self.low_damping.value()
786 hd = self.high_damping.value()
787 fs = self.frequency_samples.value()
788 ds = self.damping_samples.value()
789 df = (hf - lf) / (fs - 1)
790 dd = (hd - ld) / (ds - 1)
791 freq, damp, matrix = self.SMAC_GUI.smac.compute_correlation_matrix(
792 low_frequency=lf, high_frequency=hf,
793 frequency_samples=fs,
794 low_damping=ld / 100,
795 high_damping=hd / 100,
796 damping_samples=ds,
797 frequency_lines_for_correlation=self.SMAC_GUI.frequency_autofit_lines_to_use_in_correlation.value())
798 matrix = np.max(matrix, axis=0)
799 freq_ind, damp_ind = np.unravel_index(matrix.argmax(), matrix.shape)
800 self.max_correlation = matrix[freq_ind, damp_ind]
801 self.max_freq = freq[freq_ind]
802 self.max_damp = damp[damp_ind]
803 self.correlation_description.setText('Correlation:\n{:0.5f}\nFrequency:\n{:0.3f}\nDamping:\n{:0.3f}%'.format(
804 self.max_correlation, self.max_freq, self.max_damp * 100))
805 if self.linear_plot.isChecked():
806 self.image_log_scale = False
807 else:
808 self.image_log_scale = True
809 matrix = -np.log10(1 - matrix)
810 self.image = pyqtgraph.ImageItem(
811 matrix, rect=[lf - (df / 2), ld - (dd / 2), hf - lf + df, hd - ld + dd])
812 self.correlation_plot.addItem(self.image)
814 @staticmethod
815 def add_root(parent):
816 dialog = AddRootDialog(parent)
817 parent.add_root_dialog = dialog
818 result = (dialog.exec_() == QtWidgets.QDialog.Accepted)
819 print('Continued: {:}'.format(result))
820 if result:
821 return dialog.max_freq, dialog.max_damp, dialog.max_correlation
822 else:
823 return None, None, None
824 parent.add_root_dialog = None
827class SMAC_GUI(QMainWindow):
828 def __init__(self, frf_data: TransferFunctionArray):
829 super(SMAC_GUI, self).__init__()
830 uic.loadUi(os.path.join(os.path.abspath(
831 os.path.dirname(os.path.abspath(__file__))), 'smac.ui'), self)
833 # Store FRFs
834 # TODO Make it so we can load date if none is passed
835 self.frfs = frf_data.reshape_to_matrix()
837 # Set up colormaps
838 self.cm = cm.Dark2
839 self.cm_mod = 8
840 self.compare_cm = cm.Paired
841 self.compare_cm_mod = 12
843 # Initialize Plots
844 self.mif_plot = self.mif_plot_widget.getPlotItem()
845 self.correlation_plot = self.corrrelation_coefficient_plot_widget.getPlotItem()
847 self.mif_plot.getAxis('left').setLabel('Mode Indicator Function')
848 self.correlation_plot.getAxis('left').setLabel('Correlation Coefficient')
849 self.mif_plot.setXLink(self.correlation_plot)
850 self.correlation_coefficient_selector = pyqtgraph.LinearRegionItem(values=(0, 0.9),
851 orientation='horizontal',
852 bounds=(0, 1))
854 self.connect_callbacks()
856 # Set defaults
857 self.high_frequency = self.frfs.abscissa.max()
858 self.low_frequency = self.frfs.abscissa.min()
859 self.low_frequency_selector.setValue(self.low_frequency)
860 self.high_frequency_selector.setValue(self.high_frequency)
861 self.frequency_spacing_selector.setValue(np.mean(np.diff(self.frfs[0, 0].abscissa)))
862 self.num_responses_label.setText('{:} Responses'.format(self.frfs.shape[0]))
863 self.num_references_label.setText('{:} References'.format(self.frfs.shape[1]))
864 self.num_frequencies_label.setText('{:} Frequencies'.format(self.frfs.num_elements))
865 self.smac = None
866 self.initial_roots = []
867 self.initial_roots_for_autofit = []
868 self.initial_correlation = None
869 self.root_markers = []
870 self.crosshair_v_mif = pyqtgraph.InfiniteLine(angle=90, movable=False)
871 self.crosshair_v_corr = pyqtgraph.InfiniteLine(angle=90, movable=False)
872 self.mif_plot.addItem(self.crosshair_v_mif)
873 self.correlation_plot.addItem(self.correlation_coefficient_selector)
874 self.correlation_plot.addItem(self.crosshair_v_corr)
875 self.rootlist = []
876 self.guiplots = []
877 self.shapes = None
878 self.geometry = None
879 self.resynthesized_frfs = None
880 self.add_root_dialog = None
881 self.setWindowTitle('PySMAC')
883 self.show()
885 def connect_callbacks(self):
886 self.frequency_spacing_selector.valueChanged.connect(self.update_frequency_line_label)
887 self.compute_pseudoinverse_button.clicked.connect(self.compute_pseudoinverse)
888 self.compute_correlation_matrix_button.clicked.connect(self.compute_correlation_matrix)
889 self.correlation_plot.scene().sigMouseMoved.connect(self.update_crosshairs)
890 self.correlation_plot.scene().sigMouseClicked.connect(self.add_initial_root)
891 self.roots_delete_button.clicked.connect(self.delete_initial_roots)
892 self.rootlist_tablewidget.itemSelectionChanged.connect(self.paint_markers)
893 self.minimum_coefficient_entry.valueChanged.connect(self.update_selector_and_refind)
894 self.correlation_coefficient_selector.sigRegionChanged.connect(
895 self.update_coefficient_and_refind)
896 self.confirm_initial_rootlist_button.clicked.connect(self.confirm_initial_roots_for_autofit)
897 self.autofit_roots_button.clicked.connect(self.autofit_roots)
898 self.delete_button.clicked.connect(self.delete_roots)
899 self.add_button.clicked.connect(self.add_root)
900 self.resynthesize_button.clicked.connect(self.compute_shapes)
901 self.load_geometry_button.clicked.connect(self.load_geometry)
902 self.plot_shapes_button.clicked.connect(self.plot_shapes)
903 self.plot_mac_button.clicked.connect(self.plot_mac)
904 self.save_shapes_button.clicked.connect(self.save_shapes)
905 self.export_mode_table_button.clicked.connect(self.export_mode_table)
906 self.merge_button.clicked.connect(self.load_shapes)
908 def plot_mac(self):
909 mac_matrix = mac(self.shapes.shape_matrix.T)
910 matrix_plot(mac_matrix)
912 def save_shapes(self):
913 if self.shapes is None:
914 QtWidgets.QMessageBox.warning(self, 'PySMAC', 'Shapes not created!')
915 return
916 filename, file_filter = QtWidgets.QFileDialog.getSaveFileName(
917 self, 'Save Shapes', filter='Numpy Files (*.npy)')
918 if filename == '':
919 return
920 self.shapes.save(filename)
922 def load_shapes(self):
923 filename, file_filter = QtWidgets.QFileDialog.getOpenFileName(
924 self, 'Open Shapes', filter='Numpy Files (*.npy);;Universal Files (*.unv *.uff)')
925 if filename == '':
926 return
927 new_shapes = ShapeArray.load(filename).flatten()
928 new_root = np.zeros((new_shapes.size,), dtype=self.rootlist.dtype)
929 new_root['frequency'] = new_shapes.frequency
930 new_root['damping'] = new_shapes.damping
931 new_root['correlation'] = 0
932 new_root['shape_reference'] = -1
933 self.rootlist = np.sort(np.concatenate((self.rootlist, new_root)))
934 self.update_rootlist_table()
936 def export_mode_table(self):
937 QtWidgets.QMessageBox.warning(self, 'PySMAC', 'Not Implemented')
939 def load_geometry(self):
940 filename, file_filter = QtWidgets.QFileDialog.getOpenFileName(
941 self, 'Open Geometry', filter='Numpy Files (*.npz);;Universal Files (*.unv *.uff)')
942 if filename == '':
943 return
944 self.geometry = Geometry.load(filename)
946 def plot_shapes(self):
947 if self.geometry is None:
948 self.load_geometry()
949 if self.geometry is None:
950 QtWidgets.QMessageBox.warning(self, 'PySMAC', 'No Geometry Loaded!')
951 return
952 if self.shapes is None:
953 QtWidgets.QMessageBox.warning(self, 'PySMAC', 'Shapes not created!')
954 return
955 self.geometry.plot_shape(self.shapes)
957 def compute_shapes(self):
958 try:
959 row_indices = [i for i in range(self.root_table.rowCount())
960 if self.root_table.cellWidget(i, 0).isChecked()]
961 logical_index = np.zeros(self.rootlist.shape, dtype=bool)
962 logical_index[row_indices] = True
963 rootlist = self.rootlist[logical_index]
964 if rootlist.size == 0:
965 QtWidgets.QMessageBox.warning(
966 self, 'PySMAC', 'No Roots Selected!\n\nPlease select roots to use in resynthesis!')
967 return
968 self.smac.compute_residues(rootlist, residuals=self.use_residuals_checkbox.isChecked())
969 self.smac.compute_shapes()
970 self.resynthesized_frfs = (self.smac.shapes.compute_frf(self.frfs[0, 0].abscissa[self.smac.frequency_slice],
971 np.unique(
972 self.frfs.response_coordinate),
973 np.unique(
974 self.frfs.reference_coordinate),
975 self.datatype_selector.currentIndex())
976 + self.smac.frfs_synth_residual)
977 if self.collapse_to_real_checkbox.isChecked():
978 self.shapes = self.smac.shapes.to_real()
979 else:
980 self.shapes = self.smac.shapes
981 for row_index in np.arange(self.root_table.rowCount()):
982 item = QtWidgets.QTableWidgetItem('')
983 self.root_table.setItem(row_index, 6, item)
984 for row_index, shape in zip(np.sort(row_indices), self.shapes):
985 item = QtWidgets.QTableWidgetItem(str(shape.comment1))
986 self.root_table.setItem(row_index, 6, item)
987 # Now plot the shapes requested
988 self.guiplots = []
989 if self.nmif_checkbox.isChecked():
990 self.guiplots.append(
991 GUIPlot(self.frfs.compute_nmif(),
992 self.resynthesized_frfs.compute_nmif())
993 )
994 if self.mmif_checkbox.isChecked():
995 self.guiplots.append(
996 GUIPlot(self.frfs.compute_mmif(),
997 self.resynthesized_frfs.compute_mmif())
998 )
999 if self.cmif_checkbox.isChecked():
1000 self.guiplots.append(
1001 GUIPlot(self.frfs.compute_cmif(),
1002 self.resynthesized_frfs.compute_cmif())
1003 )
1004 if self.frf_checkbox.isChecked():
1005 self.guiplots.append(
1006 GUIPlot(self.frfs,
1007 self.resynthesized_frfs)
1008 )
1009 except Exception as e:
1010 print('Error: {:}'.format(str(e)))
1012 def add_root(self):
1013 freq, damp, corr = AddRootDialog.add_root(self)
1014 if freq is not None:
1015 # add root to the table
1016 new_root = np.zeros((1,), dtype=self.rootlist.dtype)
1017 new_root['frequency'] = freq
1018 new_root['damping'] = damp
1019 new_root['correlation'] = corr
1020 new_root['shape_reference'] = -1
1021 self.rootlist = np.sort(np.concatenate((self.rootlist, new_root)))
1022 self.update_rootlist_table()
1024 def delete_roots(self):
1025 select = self.root_table.selectionModel()
1026 row_indices = [val.row() for val in select.selectedRows()]
1027 logical_index = np.ones(self.rootlist.shape, dtype=bool)
1028 logical_index[row_indices] = False
1029 self.rootlist = self.rootlist[logical_index]
1030 self.update_rootlist_table()
1032 def autofit_roots(self):
1033 self.smac.autofit_roots(
1034 self.frequency_autofit_range_selector.value() / 100,
1035 self.frequency_autofit_values_in_range.value(),
1036 self.frequency_autofit_convergence_percentage_selector.value() / 100,
1037 self.minimum_damping_autofit_selector.value() / 100,
1038 self.maximum_damping_autofit_selector.value() / 100,
1039 self.damping_autofit_number_samples_selector.value(),
1040 self.damping_convergence_percentage_selector.value() / 100,
1041 self.frequency_autofit_lines_to_use_in_correlation.value(),
1042 self.maximum_iterations_selector.value(),
1043 self.zoom_rate_selector.value(),
1044 self.display_convergence_graphs_checkbox.isChecked())
1045 self.rootlist = self.smac.rootlist
1046 self.rootlist['shape_reference'] = -1
1047 # Initialize table
1048 self.update_rootlist_table()
1049 self.smac_tabs.setCurrentIndex(4)
1051 def update_rootlist_table(self):
1052 self.shapes = None
1053 self.resynthesized_frfs = None
1054 self.root_table.clearContents()
1055 self.root_table.setRowCount(self.rootlist.size)
1056 for i, root in enumerate(self.rootlist):
1057 checkbox = QtWidgets.QCheckBox()
1058 checkbox.setChecked(True)
1059 self.root_table.setCellWidget(i, 0, checkbox)
1060 item = QtWidgets.QTableWidgetItem('{:0.2f}'.format(root['frequency']))
1061 self.root_table.setItem(i, 1, item)
1062 item = QtWidgets.QTableWidgetItem('{:0.3f}'.format(root['damping'] * 100))
1063 self.root_table.setItem(i, 2, item)
1064 item = QtWidgets.QTableWidgetItem('{:0.4f}'.format(root['correlation']))
1065 self.root_table.setItem(i, 3, item)
1066 item = QtWidgets.QTableWidgetItem('{:0.2f}'.format(root['initial_frequency']))
1067 self.root_table.setItem(i, 4, item)
1068 item = QtWidgets.QTableWidgetItem('{:0.3f}'.format(root['initial_damping'] * 100))
1069 self.root_table.setItem(i, 5, item)
1070 item = QtWidgets.QTableWidgetItem('{:0.4f}'.format(root['initial_correlation']))
1071 self.root_table.setItem(i, 6, item)
1072 if root['shape_reference'] >= 0:
1073 item = QtWidgets.QTableWidgetItem(
1074 str(np.unique(self.frfs.reference_coordinate)[root['shape_reference']]))
1075 self.root_table.setItem(i, 7, item)
1076 item = QtWidgets.QTableWidgetItem('{:0.2f}'.format(root['drive_point_coefficient']))
1077 self.root_table.setItem(i, 8, item)
1079 def confirm_initial_roots_for_autofit(self):
1080 self.initial_roots_for_autofit = self.initial_roots.copy()
1081 self.smac.initial_rootlist = self.initial_roots_for_autofit
1082 self.smac_tabs.setCurrentIndex(3)
1084 def update_coefficient_and_refind(self):
1085 self.minimum_coefficient_entry.blockSignals(True)
1086 self.minimum_coefficient_entry.setValue(
1087 self.correlation_coefficient_selector.getRegion()[1])
1088 self.minimum_coefficient_entry.blockSignals(False)
1089 self.refind_peaks()
1091 def update_selector_and_refind(self):
1092 self.correlation_coefficient_selector.blockSignals(True)
1093 self.correlation_coefficient_selector.setRegion((0, self.minimum_coefficient_entry.value()))
1094 self.correlation_coefficient_selector.blockSignals(False)
1095 self.refind_peaks()
1097 def refind_peaks(self):
1098 num_roots_mif = 'cmif'
1099 num_roots_frequency_threshold = 0.005
1100 peaklists = self.smac.find_peaks(self.initial_correlation[..., np.newaxis], size=3,
1101 threshold=self.minimum_coefficient_entry.value())
1102 # Collapse identical roots
1103 root_index_list = {}
1104 for reference_index, (matrix, peaklist) in enumerate(zip(self.initial_correlation[..., np.newaxis], peaklists)):
1105 for key in zip(*peaklist):
1106 value = matrix[key]
1107 if key not in root_index_list or root_index_list[key][0] < value:
1108 root_index_list[key] = (value, reference_index)
1109 root_list = np.ndarray(len(root_index_list), dtype=[('frequency', 'float64'),
1110 ('damping', 'float64'),
1111 ('correlation', 'float64'),
1112 ('reference_index', 'int'),
1113 ('num_roots', 'int')])
1114 for index, (frequency_index, damping_index) in enumerate(sorted(root_index_list)):
1115 frequency = self.initial_correlation_frequencies[frequency_index]
1116 damping = self.initial_correlation_damping
1117 root_list[index] = (frequency, damping,) + root_index_list[frequency_index, 0] + (0,)
1118 root_list['num_roots'] = self.smac.get_num_roots(
1119 root_list['frequency'], num_roots_mif, num_roots_frequency_threshold)
1120 root_list = root_list[root_list['num_roots'] > 0]
1121 self.initial_roots = root_list
1122 self.update_initial_rootlist_tab(no_reset_axes=True)
1124 def paint_markers(self):
1125 select = self.rootlist_tablewidget.selectionModel()
1126 row_indices = [val.row() for val in select.selectedRows()]
1127 logical_index = np.zeros(self.initial_roots.shape, dtype=bool)
1128 logical_index[row_indices] = True
1129 for markers, selected in zip(self.root_markers, logical_index):
1130 if selected:
1131 for marker in markers:
1132 marker.setBrush((0, 0, 255))
1133 else:
1134 for marker in markers:
1135 marker.setBrush((0, 0, 0, 0))
1137 def delete_initial_roots(self):
1138 select = self.rootlist_tablewidget.selectionModel()
1139 row_indices = [val.row() for val in select.selectedRows()]
1140 logical_index = np.ones(self.initial_roots.shape, dtype=bool)
1141 logical_index[row_indices] = False
1142 self.initial_roots = self.initial_roots[logical_index]
1143 self.update_initial_rootlist_tab(no_reset_axes=True)
1145 def compute_correlation_matrix(self):
1146 self.smac.compute_initial_rootlist(
1147 frequency_resolution=self.frequency_spacing_selector.value(),
1148 low_damping=self.initial_damping_selector.value() / 100,
1149 damping_samples=1,
1150 frequency_lines_for_correlation=self.correlation_lines_selector.value())
1151 self.initial_roots = self.smac.initial_rootlist
1152 self.initial_correlation = self.smac.initial_correlation_matrix[..., 0]
1153 self.initial_correlation_frequencies = self.smac.initial_correlation_frequencies
1154 self.initial_correlation_damping = self.smac.initial_correlation_dampings[0]
1155 self.update_initial_rootlist_tab()
1156 self.smac_tabs.setCurrentIndex(2)
1158 def update_initial_rootlist_tab(self, no_reset_axes=False):
1159 if no_reset_axes:
1160 mif_range = self.mif_plot.vb.viewRange()
1161 corr_range = self.correlation_plot.vb.viewRange()
1162 self.mif_plot.clear()
1163 self.mif_plot.addItem(self.crosshair_v_mif)
1164 self.correlation_plot.clear()
1165 self.correlation_plot.addItem(self.correlation_coefficient_selector)
1166 self.correlation_plot.addItem(self.crosshair_v_corr)
1167 if self.cmif_radio_button.isChecked():
1168 self.mif_plot.setLogMode(False, True)
1169 mif = self.frfs.compute_cmif()
1170 for i, curve in enumerate(mif):
1171 pen = pyqtgraph.mkPen(color=[int(255 * v) for v in self.cm(i % self.cm_mod)])
1172 self.mif_plot.plot(curve.abscissa, curve.ordinate, pen=pen)
1173 elif self.mmif_radio_button.isChecked():
1174 self.mif_plot.setLogMode(False, False)
1175 mif = self.frfs.compute_mmif()
1176 for i, curve in enumerate(mif):
1177 pen = pyqtgraph.mkPen(color=[int(255 * v) for v in self.cm(i % self.cm_mod)])
1178 self.mif_plot.plot(curve.abscissa, curve.ordinate, pen=pen)
1179 elif self.nmif_radio_button.isChecked():
1180 self.mif_plot.setLogMode(False, False)
1181 mif = self.frfs.compute_nmif()
1182 pen = pyqtgraph.mkPen(color=[int(255 * v) for v in self.cm(0 % self.cm_mod)])
1183 self.mif_plot.plot(mif.abscissa, mif.ordinate)
1184 # Now plot the correlation coefficient
1185 if self.initial_correlation is not None:
1186 for i, curve in enumerate(self.initial_correlation):
1187 pen = pyqtgraph.mkPen(color=[int(255 * v) for v in self.cm(i % self.cm_mod)])
1188 self.correlation_plot.plot(self.initial_correlation_frequencies, curve, pen=pen)
1189 self.correlation_plot.setLimits(yMin=0, yMax=1.1)
1190 if no_reset_axes:
1191 self.correlation_plot.setRange(
1192 xRange=corr_range[0], yRange=corr_range[1], padding=0)
1193 self.mif_plot.setRange(xRange=mif_range[0], yRange=mif_range[1], padding=0)
1194 else:
1195 self.correlation_plot.setRange(xRange=(self.low_frequency, self.high_frequency),
1196 yRange=(0, 1))
1197 self.rootlist_tablewidget.setRowCount(len(self.initial_roots))
1198 self.root_markers = []
1199 try:
1200 mif = mif[0]
1201 except IndexError:
1202 pass
1203 root_mif_ordinates = np.interp(
1204 self.initial_roots['frequency'], mif.abscissa, mif.ordinate)
1205 root_corr_ordinates = np.interp(
1206 self.initial_roots['frequency'], self.initial_correlation_frequencies, np.max(self.initial_correlation, axis=0))
1207 for index, (root, mif_ord, corr_ord) in enumerate(zip(self.initial_roots, root_mif_ordinates, root_corr_ordinates)):
1208 frequency = root['frequency']
1209 item = QtWidgets.QTableWidgetItem('{:0.2f}'.format(frequency))
1210 self.rootlist_tablewidget.setItem(index, 0, item)
1211 item = QtWidgets.QTableWidgetItem('{:0.4f}'.format(root['correlation']))
1212 self.rootlist_tablewidget.setItem(index, 1, item)
1213 item = QtWidgets.QTableWidgetItem('{:}'.format(root['num_roots']))
1214 self.rootlist_tablewidget.setItem(index, 2, item)
1215 mif_item = pyqtgraph.ScatterPlotItem([frequency], [np.log10(
1216 mif_ord)], symbol='o', symbolPen=(0, 0, 255), brush=(0, 0, 0, 0))
1217 self.mif_plot.addItem(mif_item)
1218 corr_item = pyqtgraph.ScatterPlotItem(
1219 [frequency], [corr_ord], symbol='o', symbolPen=(0, 0, 255), brush=(0, 0, 0, 0))
1220 self.correlation_plot.addItem(corr_item)
1221 self.root_markers.append((mif_item, corr_item))
1223 def update_crosshairs(self, position):
1224 if self.correlation_plot.vb.sceneBoundingRect().contains(position):
1225 mouse_point = self.correlation_plot.vb.mapSceneToView(position)
1226 self.crosshair_v_corr.setPos(mouse_point.x())
1227 self.crosshair_v_mif.setPos(mouse_point.x())
1229 def compute_pseudoinverse(self):
1230 self.smac = SMAC(self.frfs, self.low_frequency_selector.value(),
1231 self.high_frequency_selector.value(), self.complex_modes_selector.isChecked(),
1232 self.datatype_selector.currentIndex())
1233 self.smac.compute_pseudoinverse()
1234 self.low_frequency = self.low_frequency_selector.value()
1235 self.high_frequency = self.high_frequency_selector.value()
1236 self.update_frequency_line_label()
1237 self.smac_tabs.setCurrentIndex(1)
1239 def update_frequency_line_label(self):
1240 num_samples = int(np.round((self.high_frequency - self.low_frequency) /
1241 self.frequency_spacing_selector.value())) + 1
1242 self.frequency_spacing_label.setText(
1243 ' Hz Frequency Spacing ({:} Frequency Samples)'.format(num_samples))
1245 def add_initial_root(self, event):
1246 position = event.scenePos()
1247 if self.correlation_plot.vb.sceneBoundingRect().contains(position):
1248 mouse_point = self.correlation_plot.vb.mapSceneToView(position)
1249 self.crosshair_v_corr.setPos(mouse_point.x())
1250 self.crosshair_v_mif.setPos(mouse_point.x())
1251 frequency_index = np.argmin(
1252 np.abs(self.initial_correlation_frequencies - mouse_point.x()))
1253 new_root = np.ndarray((1,), dtype=self.initial_roots.dtype)
1254 new_root['frequency'] = self.initial_correlation_frequencies[frequency_index]
1255 new_root['reference_index'] = np.argmax(self.initial_correlation[:, frequency_index])
1256 new_root['correlation'] = self.initial_correlation[new_root['reference_index'], frequency_index]
1257 new_root['damping'] = self.initial_correlation_damping
1258 new_root['num_roots'] = self.smac.get_num_roots(
1259 np.array([new_root['frequency']]), 'cmif')[0]
1260 self.initial_roots = np.sort(np.concatenate((self.initial_roots, new_root)))
1261 self.update_initial_rootlist_tab(no_reset_axes=True)