From 06607da39de2ae4c65780165470fc12d569406dd Mon Sep 17 00:00:00 2001 From: johnjasa Date: Thu, 12 Oct 2023 11:54:57 -0500 Subject: [PATCH 1/4] jac fix for timeseries derivs --- .../pseudospectral_timeseries_output_comp.py | 34 +++++++++++++------ 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py index c56a7ed0f..fb2021bfe 100644 --- a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py +++ b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py @@ -150,35 +150,49 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals self._vars[name] = (input_name, output_name, shape, rate) size = np.prod(shape) - jac = np.zeros((output_num_nodes, size, input_num_nodes, size)) if rate: mat = self.differentiation_matrix else: mat = self.interpolation_matrix - for i in range(size): - if _USE_SPARSE: - jac[:, i, :, i] = mat.toarray() - else: - jac[:, i, :, i] = mat + # Preallocate lists to hold data for jac + jac_data = [] + jac_indices = [] + jac_indptr = [0] + + # Iterate over the dense dimension 'size' + for s in range(size): + # Extend the data and indices using the CSR attributes of mat + jac_data.extend(mat.data) + jac_indices.extend(mat.indices + s * input_num_nodes) + + # For every non-zero row in mat, update jac's indptr + new_indptr = mat.indptr[1:] + s * len(mat.data) + jac_indptr.extend(new_indptr) + + # Correct the last entry of jac_indptr + jac_indptr[-1] = len(jac_data) - jac = jac.reshape((output_num_nodes * size, input_num_nodes * size), order='C') + # Construct the sparse jac matrix in CSR format + jac = sp.csr_matrix((jac_data, jac_indices, jac_indptr), shape=(output_num_nodes * size, input_num_nodes * size)) - jac_rows, jac_cols = np.nonzero(jac) + # Now, if you want to get the row and column indices of the non-zero entries in the jac matrix: + jac_rows, jac_cols = jac.nonzero() # There's a chance that the input for this output was pulled from another variable with # different units, so account for that with a conversion. + val = np.squeeze(np.array(jac[jac_rows, jac_cols])) if input_units is None or units is None: self.declare_partials(of=output_name, wrt=input_name, - rows=jac_rows, cols=jac_cols, val=jac[jac_rows, jac_cols]) + rows=jac_rows, cols=jac_cols, val=val) else: scale, offset = unit_conversion(input_units, units) self._conversion_factors[output_name] = scale, offset self.declare_partials(of=output_name, wrt=input_name, rows=jac_rows, cols=jac_cols, - val=scale * jac[jac_rows, jac_cols]) + val=scale * val) return added_source From 074bde7e02be1756e3c771e2ae671fec5ccabcd9 Mon Sep 17 00:00:00 2001 From: johnjasa Date: Thu, 12 Oct 2023 12:11:52 -0500 Subject: [PATCH 2/4] Also updated analytic phase deriv calc --- .../analytic_timeseries_output_comp.py | 34 +++++++++++++------ 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py index 23dab54d6..f515e7efa 100644 --- a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py +++ b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py @@ -150,35 +150,49 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals self._vars[name] = (input_name, output_name, shape, rate) size = np.prod(shape) - jac = np.zeros((output_num_nodes, size, input_num_nodes, size)) if rate: mat = self.differentiation_matrix else: mat = self.interpolation_matrix - for i in range(size): - if _USE_SPARSE: - jac[:, i, :, i] = mat.toarray() - else: - jac[:, i, :, i] = mat + # Preallocate lists to hold data for jac + jac_data = [] + jac_indices = [] + jac_indptr = [0] + + # Iterate over the dense dimension 'size' + for s in range(size): + # Extend the data and indices using the CSR attributes of mat + jac_data.extend(mat.data) + jac_indices.extend(mat.indices + s * input_num_nodes) + + # For every non-zero row in mat, update jac's indptr + new_indptr = mat.indptr[1:] + s * len(mat.data) + jac_indptr.extend(new_indptr) + + # Correct the last entry of jac_indptr + jac_indptr[-1] = len(jac_data) - jac = jac.reshape((output_num_nodes * size, input_num_nodes * size), order='C') + # Construct the sparse jac matrix in CSR format + jac = sp.csr_matrix((jac_data, jac_indices, jac_indptr), shape=(output_num_nodes * size, input_num_nodes * size)) - jac_rows, jac_cols = np.nonzero(jac) + # Now, if you want to get the row and column indices of the non-zero entries in the jac matrix: + jac_rows, jac_cols = jac.nonzero() # There's a chance that the input for this output was pulled from another variable with # different units, so account for that with a conversion. + val = np.squeeze(np.array(jac[jac_rows, jac_cols])) if input_units is None or units is None: self.declare_partials(of=output_name, wrt=input_name, - rows=jac_rows, cols=jac_cols, val=jac[jac_rows, jac_cols]) + rows=jac_rows, cols=jac_cols, val=val) else: scale, offset = unit_conversion(input_units, units) self._conversion_factors[output_name] = scale, offset self.declare_partials(of=output_name, wrt=input_name, rows=jac_rows, cols=jac_cols, - val=scale * jac[jac_rows, jac_cols]) + val=scale * val) return added_source From 04530a9c9449f3740bcae1a9dd77a524321879d5 Mon Sep 17 00:00:00 2001 From: johnjasa Date: Mon, 16 Oct 2023 11:23:39 -0500 Subject: [PATCH 3/4] Removed use_sparse flag --- .../analytic/analytic_timeseries_output_comp.py | 11 ++--------- .../pseudospectral_timeseries_output_comp.py | 11 ++--------- 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py index f515e7efa..5514edb51 100644 --- a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py +++ b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py @@ -9,9 +9,6 @@ from ..common.timeseries_output_comp import TimeseriesOutputCompBase -_USE_SPARSE = True - - class AnalyticTimeseriesOutputComp(TimeseriesOutputCompBase): """ Class definition of the timeseries output comp for AnalyticPhase. @@ -86,12 +83,8 @@ def setup(self): L_blocks.append(L) D_blocks.append(D) - if _USE_SPARSE: - self.interpolation_matrix = sp.block_diag(L_blocks, format='csr') - self.differentiation_matrix = sp.block_diag(D_blocks, format='csr') - else: - self.interpolation_matrix = block_diag(*L_blocks) - self.differentiation_matrix = block_diag(*D_blocks) + self.interpolation_matrix = sp.block_diag(L_blocks, format='csr') + self.differentiation_matrix = sp.block_diag(D_blocks, format='csr') self.add_input('dt_dstau', shape=(self.input_num_nodes,), units=self.options['time_units']) diff --git a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py index fb2021bfe..d10aeef55 100644 --- a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py +++ b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py @@ -9,9 +9,6 @@ from ...common.timeseries_output_comp import TimeseriesOutputCompBase -_USE_SPARSE = True - - class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase): """ Class definition of the PseudospectralTimeseriesOutputComp. @@ -86,12 +83,8 @@ def setup(self): L_blocks.append(L) D_blocks.append(D) - if _USE_SPARSE: - self.interpolation_matrix = sp.block_diag(L_blocks, format='csr') - self.differentiation_matrix = sp.block_diag(D_blocks, format='csr') - else: - self.interpolation_matrix = block_diag(*L_blocks) - self.differentiation_matrix = block_diag(*D_blocks) + self.interpolation_matrix = sp.block_diag(L_blocks, format='csr') + self.differentiation_matrix = sp.block_diag(D_blocks, format='csr') self.add_input('dt_dstau', shape=(self.input_num_nodes,), units=self.options['time_units']) From fba061dc70e6b9b6c00ec5fba5a143251cee2133 Mon Sep 17 00:00:00 2001 From: johnjasa Date: Mon, 16 Oct 2023 11:46:32 -0500 Subject: [PATCH 4/4] pycodestyle updates --- .../analytic/analytic_timeseries_output_comp.py | 5 +++-- .../components/pseudospectral_timeseries_output_comp.py | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py index 5514edb51..ab4373ad5 100644 --- a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py +++ b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py @@ -159,7 +159,7 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals # Extend the data and indices using the CSR attributes of mat jac_data.extend(mat.data) jac_indices.extend(mat.indices + s * input_num_nodes) - + # For every non-zero row in mat, update jac's indptr new_indptr = mat.indptr[1:] + s * len(mat.data) jac_indptr.extend(new_indptr) @@ -168,7 +168,8 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals jac_indptr[-1] = len(jac_data) # Construct the sparse jac matrix in CSR format - jac = sp.csr_matrix((jac_data, jac_indices, jac_indptr), shape=(output_num_nodes * size, input_num_nodes * size)) + jac = sp.csr_matrix((jac_data, jac_indices, jac_indptr), + shape=(output_num_nodes * size, input_num_nodes * size)) # Now, if you want to get the row and column indices of the non-zero entries in the jac matrix: jac_rows, jac_cols = jac.nonzero() diff --git a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py index d10aeef55..c2e029edb 100644 --- a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py +++ b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py @@ -159,7 +159,7 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals # Extend the data and indices using the CSR attributes of mat jac_data.extend(mat.data) jac_indices.extend(mat.indices + s * input_num_nodes) - + # For every non-zero row in mat, update jac's indptr new_indptr = mat.indptr[1:] + s * len(mat.data) jac_indptr.extend(new_indptr) @@ -168,7 +168,8 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals jac_indptr[-1] = len(jac_data) # Construct the sparse jac matrix in CSR format - jac = sp.csr_matrix((jac_data, jac_indices, jac_indptr), shape=(output_num_nodes * size, input_num_nodes * size)) + jac = sp.csr_matrix((jac_data, jac_indices, jac_indptr), + shape=(output_num_nodes * size, input_num_nodes * size)) # Now, if you want to get the row and column indices of the non-zero entries in the jac matrix: jac_rows, jac_cols = jac.nonzero()