From c973a787128b29ccaa75ae89f23e1edb21eb8e60 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 23 Oct 2024 13:29:32 -0700 Subject: [PATCH 01/22] Started documenting latent_space.py I also removed the multheadedattention activation type (plus the associated apply_attention function) because they did not make sense in this class. I have now documented the MLP class, but need to add documentation to the autoencoder class. --- src/lasdi/latent_space.py | 387 ++++++++++++++++++++++++++++---------- 1 file changed, 285 insertions(+), 102 deletions(-) diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index fdf8fef..d45b495 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -1,146 +1,323 @@ -import torch -import numpy as np +# ------------------------------------------------------------------------------------------------- +# Imports and setup +# ------------------------------------------------------------------------------------------------- + +import torch +import numpy as np + +from .physics import Physics # activation dict -act_dict = {'ELU': torch.nn.ELU, - 'hardshrink': torch.nn.Hardshrink, - 'hardsigmoid': torch.nn.Hardsigmoid, - 'hardtanh': torch.nn.Hardtanh, - 'hardswish': torch.nn.Hardswish, - 'leakyReLU': torch.nn.LeakyReLU, - 'logsigmoid': torch.nn.LogSigmoid, - 'multihead': torch.nn.MultiheadAttention, - 'PReLU': torch.nn.PReLU, - 'ReLU': torch.nn.ReLU, - 'ReLU6': torch.nn.ReLU6, - 'RReLU': torch.nn.RReLU, - 'SELU': torch.nn.SELU, - 'CELU': torch.nn.CELU, - 'GELU': torch.nn.GELU, - 'sigmoid': torch.nn.Sigmoid, - 'SiLU': torch.nn.SiLU, - 'mish': torch.nn.Mish, - 'softplus': torch.nn.Softplus, - 'softshrink': torch.nn.Softshrink, - 'tanh': torch.nn.Tanh, - 'tanhshrink': torch.nn.Tanhshrink, - 'threshold': torch.nn.Threshold, - } - -def initial_condition_latent(param_grid, physics, autoencoder): - - ''' - - Outputs the initial condition in the latent space: Z0 = encoder(U0) - - ''' - - n_param = param_grid.shape[0] - Z0 = [] - - sol_shape = [1, 1] + physics.qgrid_size +act_dict = {'ELU' : torch.nn.ELU, + 'hardshrink' : torch.nn.Hardshrink, + 'hardsigmoid' : torch.nn.Hardsigmoid, + 'hardtanh' : torch.nn.Hardtanh, + 'hardswish' : torch.nn.Hardswish, + 'leakyReLU' : torch.nn.LeakyReLU, + 'logsigmoid' : torch.nn.LogSigmoid, + 'PReLU' : torch.nn.PReLU, + 'ReLU' : torch.nn.ReLU, + 'ReLU6' : torch.nn.ReLU6, + 'RReLU' : torch.nn.RReLU, + 'SELU' : torch.nn.SELU, + 'CELU' : torch.nn.CELU, + 'GELU' : torch.nn.GELU, + 'sigmoid' : torch.nn.Sigmoid, + 'SiLU' : torch.nn.SiLU, + 'mish' : torch.nn.Mish, + 'softplus' : torch.nn.Softplus, + 'softshrink' : torch.nn.Softshrink, + 'tanh' : torch.nn.Tanh, + 'tanhshrink' : torch.nn.Tanhshrink, + 'threshold' : torch.nn.Threshold,} + + + +# ------------------------------------------------------------------------------------------------- +# initial_conditions_latent function +# ------------------------------------------------------------------------------------------------- + +def initial_condition_latent(param_grid : np.ndarray, + physics : Physics, + autoencoder : torch.nn.Module) -> list[np.ndarray]: + """ + This function maps a set of initial conditions for the fom to initial conditions for the + latent space dynamics. Specifically, we take in a set of possible parameter values. For each + set of parameter values, we recover the fom IC (from physics), then map this fom IC to a + latent space IC (by encoding it using the autoencoder). We do this for each parameter + combination and then return a list housing the latent space ICs. + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + param_grid: A 2d numpy.ndarray object of shape (number of parameter combination) x (number of + parameters). + + physics: A "Physics" object that stores the datasets for each parameter combination. + autoencoder: The actual autoencoder object that we use to map the ICs into the latent space. + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + A list of numpy ndarray objects whose i'th element holds the latent space initial condition + for the i'th set of parameters in the param_grid. That is, if we let U0_i denote the fom IC for + the i'th set of parameters, then the i'th element of the returned list is Z0_i = encoder(U0_i). + """ + + # Figure out how many parameters we are testing at. + n_param : int = param_grid.shape[0] + Z0 : list[np.ndarray] = [] + sol_shape : list[int] = [1, 1] + physics.qgrid_size + + # Cycle through the parameters. for i in range(n_param): - u0 = physics.initial_condition(param_grid[i]) + # TODO(kevin): generalize parameter class. + + # Fetch the IC for the i'th set of parameters. Then map it to a tensor. + u0 : np.ndarray = physics.initial_condition(param_grid[i]) u0 = u0.reshape(sol_shape) u0 = torch.Tensor(u0) + + # Encode the IC, then map the encoding to a numpy array. z0 = autoencoder.encoder(u0) z0 = z0[0, 0, :].detach().numpy() + + # Append the new IC to the list of latent ICs Z0.append(z0) + # Return the list of latent ICs. return Z0 + + +# ------------------------------------------------------------------------------------------------- +# MLP class +# ------------------------------------------------------------------------------------------------- + class MultiLayerPerceptron(torch.nn.Module): + def __init__( self, + layer_sizes : list[int], + act_type : str = 'sigmoid', + reshape_index : int = None, + reshape_shape : tuple[int] = None, + threshold : float = 0.1, + value : float = 0.0) -> None: + """ + This class defines a standard multi-layer network network. - def __init__(self, layer_sizes, - act_type='sigmoid', reshape_index=None, reshape_shape=None, - threshold=0.1, value=0.0, num_heads=1): - super(MultiLayerPerceptron, self).__init__() - # including input, hidden, output layers - self.n_layers = len(layer_sizes) - self.layer_sizes = layer_sizes - # Linear features between layers - self.fcs = [] - for k in range(self.n_layers-1): - self.fcs += [torch.nn.Linear(layer_sizes[k], layer_sizes[k + 1])] - self.fcs = torch.nn.ModuleList(self.fcs) - self.init_weight() + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + layer_sizes: A list of integers specifying the widths of the layers (including the + dimensionality of the domain of each layer, as well as the co-domain of the final layer). + Suppose this list has N elements. Then the network will have N - 1 layers. The i'th layer + maps from \mathbb{R}^{layer_sizes[i]} to \mathbb{R}^{layers_sizes[i]}. Thus, the i'th + element of this list represents the domain of the i'th layer AND the co-domain of the + i-1'th layer. + + act_type: A string specifying which activation function we want to use at the end of each + layer (except the final one). We use the same activation for each layer. + + reshape_index: This argument specifies if we should reshape the network's input or output + (or neither). If the user specifies reshape_index, then it must be either 0 or -1. Further, + in this case, they must also specify reshape_shape (you need to specify both together). If + it is 0, then reshape_shape specifies how we reshape the input before passing it through + the network (the input to the first layer). If reshape_index is -1, then reshape_shape + specifies how we reshape the network output before returning it (the output to the last + layer). + + reshape_shape: These are lists of k integers specifying the final k dimensions of the shape + of the input to the first layer (if reshape_index == 0) or the output of the last layer + (if reshape_index == -1). You must specify this argument if and only if you specify + reshape_index. + + threshold: You should only specify this argument if you are using the "Threshold" + activation type. In this case, it specifies the threshold value for that activation (if x + is the input and x > threshold, then the function returns x. Otherwise, it returns the + value). + + value: You should only specify this argument if you are using the "Threshold" activation + type. In this case, "value" specifies the value that the function returns if the input is + below the threshold. - # Reshape input or output layer + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + # Run checks. + # Make sure the reshape index is either 0 (input to 1st layer) or -1 (output of last + # layer). Also make sure that that product of the dimensions in the reshape_shape match + # the input dimension for the 1st layer (if reshape_index == 0) or the output dimension of + # the last layer (if reshape_index == 1). + # + # Why do we need the later condition? Well, suppose that reshape_shape has a length of k. + # If reshape_index == 0, then we squeeze the final k dimensions of the input and feed that + # into the first layer. Thus, in this case, we need the last dimension of the squeezed + # array to match the input domain of the first layer. On the other hand, reshape_index == -1 + # then we reshape the final dimension of the output to match the reshape_shape. Thus, in + # both cases, we need the product of the components of reshape_shape to match a + # corresponding element of layer_sizes. assert((reshape_index is None) or (reshape_index in [0, -1])) assert((reshape_shape is None) or (np.prod(reshape_shape) == layer_sizes[reshape_index])) - self.reshape_index = reshape_index - self.reshape_shape = reshape_shape - # Initalize activation function - self.act_type = act_type - self.use_multihead = False + super(MultiLayerPerceptron, self).__init__() + + # Note that layer_sizes specifies the dimensionality of the domains and co-domains of each + # layer. Specifically, the i'th element specifies the input dimension of the i'th layer, + # while the final element specifies the dimensionality of the co-domain of the final layer. + # Thus, the number of layers is one less than the length of layer_sizes. + self.n_layers : int = len(layer_sizes) - 1; + self.layer_sizes : list[int] = layer_sizes + + # Set up the affine parts of the layers. + self.layers : list[torch.nn.Module] = []; + for k in range(self.n_layers): + self.layers += [torch.nn.Linear(layer_sizes[k], layer_sizes[k + 1])] + self.layers = torch.nn.ModuleList(self.layers) + + # Now, initialize the weight matrices and bias vectors in the affine portion of each layer. + self.init_weight() + + # Reshape input to the 1st layer or output of the last layer. + self.reshape_index : int = reshape_index + self.reshape_shape : list[int] = reshape_shape + + # Set up the activation function. Note that we need to handle the threshold activation type + # separately because it requires an extra set of parameters to initialize. + self.act_type : str = act_type if act_type == "threshold": - self.act = act_dict[act_type](threshold, value) - - elif act_type == "multihead": - self.use_multihead = True - if (self.n_layers > 3): # if you have more than one hidden layer - self.act = [] - for i in range(self.n_layers-2): - self.act += [act_dict[act_type](layer_sizes[i+1], num_heads)] - else: - self.act = [torch.nn.Identity()] # No additional activation - self.act = torch.nn.ModuleList(self.fcs) - - #all other activation functions initialized here + self.act : torch.nn.Module = act_dict[act_type](threshold, value) else: - self.act = act_dict[act_type]() + self.act : torch.nn.Module = act_dict[act_type]() + + # All done! return - def forward(self, x): + + + def forward(self, x : torch.Tensor) -> torch.Tensor: + """ + This function defines the forward pass through self. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + x: A tensor holding a batch of inputs. We pass this tensor through the network's layers + and then return the result. If self.reshape_index == 0 and self.reshape_shape has k + elements, then the final k elements of x's shape must match self.reshape_shape. + + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + The image of x under the network's layers. If self.reshape_index == -1 and + self.reshape_shape has k elements, then we reshape the output so that the final k elements + of its shape match those of self.reshape_shape. + """ + + # If the reshape_index is 0, we need to reshape x before passing it through the first + # layer. if (self.reshape_index == 0): - # make sure the input has a proper shape + # Make sure the input has a proper shape. There is a lot going on in this line; let's + # break it down. If reshape_index == 0, then we need to reshape the input, x, before + # passing it through the layers. Let's assume that reshape_shape has k elements. Then, + # we need to squeeze the final k dimensions of the input, x, so that the resulting + # tensor has a final dimension size that matches the input dimension size for the first + # layer. The check below makes sure that the final k dimensions of the input, x, match + # the stored reshape_shape. assert(list(x.shape[-len(self.reshape_shape):]) == self.reshape_shape) - # we use torch.Tensor.view instead of torch.Tensor.reshape, - # in order to avoid data copying. + + # Now that we know the final k dimensions of x have the correct shape, let's squeeze + # them into 1 dimension (so that we can pass the squeezed tensor through the first + # layer). To do this, we reshape x by keeping all but the last k dimensions of x, and + # replacing the last k with a single dimension whose size matches the dimensionality of + # the domain of the first layer. Note that we use torch.Tensor.view instead of + # torch.Tensor.reshape in order to avoid data copying. x = x.view(list(x.shape[:-len(self.reshape_shape)]) + [self.layer_sizes[self.reshape_index]]) - for i in range(self.n_layers-2): - x = self.fcs[i](x) # apply linear layer - if (self.use_multihead): - x = self.apply_attention(self, x, i) - else: - x = self.act(x) + # Pass x through the network layers (except for the final one, which has no activation + # function). + for i in range(self.n_layers - 1): + x : torch.Tensor = self.layers[i](x) # apply linear layer + x : torch.Tensor = self.act(x) # apply activation - x = self.fcs[-1](x) + # Apply the final (output) layer. + x = self.layers[-1](x) + # If the reshape_index is -1, then we need to reshape the output before returning. if (self.reshape_index == -1): - # we use torch.Tensor.view instead of torch.Tensor.reshape, - # in order to avoid data copying. + # In this case, we need to split the last dimension of x, the output of the final + # layer, to match the reshape_shape. This is precisely what the line below does. Note + # that we use torch.Tensor.view instead of torch.Tensor.reshape in order to avoid data + # copying. x = x.view(list(x.shape[:-1]) + self.reshape_shape) + # All done! return x - def apply_attention(self, x, act_idx): - x = x.unsqueeze(1) # Add sequence dimension for attention - x, _ = self.act[act_idx](x, x, x) # apply attention - x = x.squeeze(1) # Remove sequence dimension - return x - - def init_weight(self): + + + def init_weight(self) -> None: + """ + This function initializes the weight matrices and bias vectors in self's layers. + + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + None! + + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + # TODO(kevin): support other initializations? - for fc in self.fcs: - torch.nn.init.xavier_uniform_(fc.weight) + for layer in self.layers: + torch.nn.init.xavier_uniform_(layer.weight) + torch.nn.init.zeros_(layer.bias) + + # All done! return -class Autoencoder(torch.nn.Module): - def __init__(self, physics, config): + +# ------------------------------------------------------------------------------------------------- +# Autoencoder class +# ------------------------------------------------------------------------------------------------- + +class Autoencoder(torch.nn.Module): + def __init__(self, physics : Physics, config : dict): super(Autoencoder, self).__init__() - self.qgrid_size = physics.qgrid_size - self.space_dim = np.prod(self.qgrid_size) - hidden_units = config['hidden_units'] - n_z = config['latent_dimension'] - self.n_z = n_z + self.qgrid_size = physics.qgrid_size; + self.space_dim : np.ndarray = np.prod(self.qgrid_size); + hidden_units : int = config['hidden_units']; + n_z : int = config['latent_dimension']; + self.n_z : int = n_z layer_sizes = [self.space_dim] + hidden_units + [n_z] #grab relevant initialization values from config @@ -159,6 +336,8 @@ def __init__(self, physics, config): return + + def forward(self, x): x = self.encoder(x) @@ -166,10 +345,14 @@ def forward(self, x): return x + + def export(self): dict_ = {'autoencoder_param': self.cpu().state_dict()} return dict_ + + def load(self, dict_): self.load_state_dict(dict_['autoencoder_param']) return \ No newline at end of file From 0104f9e5830cb3bb4d2ebd6d9c9a67b7ee842d0f Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 23 Oct 2024 14:56:54 -0700 Subject: [PATCH 02/22] Finished adding documentation to latent_space.py I added comments + doc strings to the Autoencoder class. --- src/lasdi/latent_space.py | 185 +++++++++++++++++++++++++++++++------- 1 file changed, 154 insertions(+), 31 deletions(-) diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index d45b495..fb664fe 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -111,7 +111,6 @@ def __init__( self, This class defines a standard multi-layer network network. - ------------------------------------------------------------------------------------------- Arguments ------------------------------------------------------------------------------------------- @@ -149,7 +148,6 @@ def __init__( self, below the threshold. - ------------------------------------------------------------------------------------------- Returns ------------------------------------------------------------------------------------------- @@ -222,7 +220,6 @@ def forward(self, x : torch.Tensor) -> torch.Tensor: elements, then the final k elements of x's shape must match self.reshape_shape. - ------------------------------------------------------------------------------------------- Returns ------------------------------------------------------------------------------------------- @@ -278,7 +275,6 @@ def init_weight(self) -> None: """ This function initializes the weight matrices and bias vectors in self's layers. - ------------------------------------------------------------------------------------------- Arguments @@ -287,7 +283,6 @@ def init_weight(self) -> None: None! - ------------------------------------------------------------------------------------------- Returns ------------------------------------------------------------------------------------------- @@ -310,49 +305,177 @@ def init_weight(self) -> None: # ------------------------------------------------------------------------------------------------- class Autoencoder(torch.nn.Module): - def __init__(self, physics : Physics, config : dict): + def __init__(self, physics : Physics, config : dict) -> None: + """ + Initializes an Autoencoder object. An Autoencoder consists of two networks, an encoder, + E : \mathbb{R}^F -> \mathbb{R}^L, and a decoder, D : \mathbb{R}^L -> \marthbb{R}^F. We + assume that the dataset consists of samples of a parameterized L-manifold in + \mathbb{R}^F. The idea then is that E and D act like the inverse coordinate patch and + coordinate patch, respectively. In our case, E and D are trainable neural networks. We + try to train E and map data in \mathbb{R}^F to elements of a low dimensional latent + space (\mathbb{R}^L) which D can send back to the original data. (thus, E, and D should + act like inverses of one another). + + The Autoencoder class implements this model as a trainable torch.nn.Module object. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + physics: A "Physics" object that holds the fom solution frames. We use this object to + determine the shape of each fom solution frame. Recall that each Physics object has a + corresponding PDE. We + + config: A dictionary representing the loaded .yml configuration file. We expect it to have + the following keys/: + hidden_units: A list of integers specifying the dimension of the co-domain of each + encoder layer except for the final one. Thus, if the k'th layer maps from + \mathbb{R}^{n(k)} to \mathbb{R}^{n(k + 1)} and there are K layers (indexed 0, 1, ... , + K - 1), then hidden_units should specify n(1), ... , n(K - 1). + + latent_dimension: The dimensionality of the Autoencoder's latent space. Equivalently, + the dimensionality of the co-domain of the encoder (i.e., the dimensionality of the + co-domain of the last layer of the encoder) and the domain of the decoder (i.e., the + dimensionality of the domain of the first layer of the decoder). + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + super(Autoencoder, self).__init__() - self.qgrid_size = physics.qgrid_size; + # A Physics object's qgrid_size is a list of integers specifying the shape of each frame of + # the fom solution. If the solution is scalar valued, then this is just a list whose i'th + # element specifies the number of grid points along the i'th spatial axis. If the solution + # is vector valued, however, we prepend the dimensionality of the vector field to the list + # from the scalar list (so the 0 element represents the dimension of the vector field at + # each point). + self.qgrid_size : list[int] = physics.qgrid_size; + + # The product of the elements of qgrid_size is the number of dimensions in each fom + # solution frame. This number represents the dimensionality of the input to the encoder + # (since we pass a flattened fom frame as input). self.space_dim : np.ndarray = np.prod(self.qgrid_size); - hidden_units : int = config['hidden_units']; + + # Fetch information about the domain/co-domain of each encoder layer. + hidden_units : list[int] = config['hidden_units']; n_z : int = config['latent_dimension']; - self.n_z : int = n_z - - layer_sizes = [self.space_dim] + hidden_units + [n_z] - #grab relevant initialization values from config - act_type = config['activation'] if 'activation' in config else 'sigmoid' - threshold = config["threshold"] if "threshold" in config else 0.1 - value = config["value"] if "value" in config else 0.0 - num_heads = config['num_heads'] if 'num_heads' in config else 1 - - self.encoder = MultiLayerPerceptron(layer_sizes, act_type, - reshape_index=0, reshape_shape=self.qgrid_size, - threshold=threshold, value=value, num_heads=num_heads) - - self.decoder = MultiLayerPerceptron(layer_sizes[::-1], act_type, - reshape_index=-1, reshape_shape=self.qgrid_size, - threshold=threshold, value=value, num_heads=num_heads) + self.n_z : int = n_z; + + # Build the "layer_sizes" argument for the MLP class. This consists of the dimensions of + # each layers' domain + the dimension of the co-domain of the final layer. + layer_sizes = [self.space_dim] + hidden_units + [n_z]; + + # Use the settings to set up the activation information for the encoder. + act_type = config['activation'] if 'activation' in config else 'sigmoid' + threshold = config["threshold"] if "threshold" in config else 0.1 + value = config["value"] if "value" in config else 0.0 + + # Now, build the encoder. + self.encoder = MultiLayerPerceptron( + layer_sizes = layer_sizes, + act_type = act_type, + reshape_index = 0, # We need to flatten the spatial dimensions of each fom frame. + reshape_shape = self.qgrid_size, + threshold = threshold, + value = value); + + self.decoder = MultiLayerPerceptron( + latent_sizes = layer_sizes[::-1], # Reverses the order of the the list. + act_type = act_type, + reshape_index = -1, + reshape_shape = self.qgrid_size, # We need to reshape the network output to a fom frame. + threshold = threshold, + value = value) + # All done! return - def forward(self, x): + def forward(self, x : torch.Tensor) -> torch.Tensor: + """ + This function defines the forward pass through self. - x = self.encoder(x) - x = self.decoder(x) + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- - return x + x: A tensor holding a batch of inputs. We pass this tensor through the encoder + decoder + and then return the result. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + The image of x under the encoder + decoder. + """ + + # Encoder the input + z : torch.Tensor = self.encoder(x) + + # Now decode z. + y : torch.Tensor = self.decoder(z) + + # All done! Hopefully y \approx x. + return y - def export(self): - dict_ = {'autoencoder_param': self.cpu().state_dict()} + def export(self) -> dict: + """ + This function extracts self's parameters and returns them in a dictionary. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + None! + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + The A dictionary housing self's state dictionary. + """ + + # TO DO: deep export which includes all information needed to re-initialize self from + # scratch. This would probably require changing the initializer. + + dict_ = { 'autoencoder_param' : self.cpu().state_dict()} return dict_ - def load(self, dict_): + def load(self, dict_ : dict) -> None: + """ + This function loads self's state dictionary. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + dict_: This should be a dictionary with the key "autoencoder_param" whose corresponding + value is the state dictionary of an autoencoder which has the same architecture (i.e., + layer sizes) as self. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + self.load_state_dict(dict_['autoencoder_param']) return \ No newline at end of file From 1b198ac7bd744adea3449c71a5f9fc4dc8e4e786 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 23 Oct 2024 15:24:03 -0700 Subject: [PATCH 03/22] Fixed bugs gplasdi had an instance of np.Inf (which is depricated). I also fixed a typo in latent_spaces.py --- src/lasdi/gplasdi.py | 2 +- src/lasdi/latent_space.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lasdi/gplasdi.py b/src/lasdi/gplasdi.py index ac3d046..5131f24 100644 --- a/src/lasdi/gplasdi.py +++ b/src/lasdi/gplasdi.py @@ -138,7 +138,7 @@ def __init__(self, physics, autoencoder, latent_dynamics, param_space, config): else: self.device = 'cpu' - self.best_loss = np.Inf + self.best_loss = np.inf self.best_coefs = None self.restart_iter = 0 diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index fb664fe..2e78542 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -386,7 +386,7 @@ def __init__(self, physics : Physics, config : dict) -> None: value = value); self.decoder = MultiLayerPerceptron( - latent_sizes = layer_sizes[::-1], # Reverses the order of the the list. + layer_sizes = layer_sizes[::-1], # Reverses the order of the the list. act_type = act_type, reshape_index = -1, reshape_shape = self.qgrid_size, # We need to reshape the network output to a fom frame. From 7d96482e81fdbe78f3354051f9d4c9c5ac3667b6 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 23 Oct 2024 15:46:57 -0700 Subject: [PATCH 04/22] Added documentation to param.py --- src/lasdi/param.py | 548 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 461 insertions(+), 87 deletions(-) diff --git a/src/lasdi/param.py b/src/lasdi/param.py index a8a2223..208a311 100644 --- a/src/lasdi/param.py +++ b/src/lasdi/param.py @@ -1,83 +1,310 @@ -import numpy as np -from .inputs import InputParser +# ------------------------------------------------------------------------------------------------- +# Imports +# ------------------------------------------------------------------------------------------------- -def get_1dspace_from_list(config): - Nx = len(config['list']) - paramRange = np.array(config['list']) +import numpy as np +from .inputs import InputParser + + + +# ------------------------------------------------------------------------------------------------- +# Helper functions +# ------------------------------------------------------------------------------------------------- + +def get_1dspace_from_list(param_dict : dict) -> tuple[int, np.ndarray]: + """ + This function generates the parameter range (set of possible parameter values) for a parameter + that uses the list type test space. That is, "test_space_type" should be a key for the + parameter dictionary and the corresponding value should be "list". The parameter dictionary + should also have a "list" key whose value is a list of the possible parameter values. + + We parse this list and turn it into a numpy ndarray. + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + param_dict: A dictionary specifying one of the parameters. We should fetch this from the + configuration yaml file. It must have a "list" key whose corresponding value is a list of + floats. + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + Two arguments: Nx and paramRange. paramRange is a 1d numpy ndarray (whose ith value is the + i'th element of param_dict["list"]). Nx is the length of paramRange. + """ + + # In this case, the parameter dictionary should have a "list" attribute which should list the + # parameter values we want to test. Fetch it (it's length is Nx) and use it to generate an + # array of possible values. + Nx : int = len(param_dict['list']) + paramRange : np.ndarray = np.array(param_dict['list']) + + # All done! return Nx, paramRange -def create_uniform_1dspace(config): - Nx = config['sample_size'] - minval = config['min'] - maxval = config['max'] - if (config['log_scale']): - paramRange = np.exp(np.linspace(np.log(minval), np.log(maxval), Nx)) + + +def create_uniform_1dspace(param_dict : dict) -> tuple[int, np.ndarray]: + """ + This function generates the parameter range (set of possible parameter values) for a parameter + that uses the uniform type test space. That is, "test_space_type" should be a key for the + parameter dictionary and the corresponding value should be "uniform". The parameter dictionary + should also have the following keys: + "min" + "max" + "sample_size" + "log_scale" + "min" and "max" specify the minimum and maximum value of the parameter, respectively. + "sample_size" specifies the number of parameter values we generate. Finally, log_scale, if + true, specifies if we should use a uniform or logarithmic spacing between samples of the + parameter. + + The values corresponding to "min" and "max" should be floats while the values corresponding to + "sample_size" and "log_scale" should be an int and a bool, respectively. + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + param_dict: A dictionary specifying one of the parameters. We should fetch this from the + configuration yaml file. It must have a "min", "max", "sample_size", and "log_scale" + keys (see above). + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + Two arguments: Nx and paramRange. paramRange is a 1d numpy ndarray (whose ith value is the + i'th possible value of the parameter. Thus, paramRange[0] = param_dict["min"] and + paramRange[-1] = param_dict["max"]). Nx is the length of paramRange or, equivalently + param_dict["sample_size"]. + """ + + # Fetch the number of samples and the min/max value for this parameter. + Nx : int = param_dict['sample_size'] + minval : float = param_dict['min'] + maxval : float = param_dict['max'] + + # Generate the range of parameter values. Note that we have to generate a uniform grid in the + # log space, then exponentiate it to generate logarithmic spacing. + if (param_dict['log_scale']): + paramRange : np.ndarray = np.exp(np.linspace(np.log(minval), np.log(maxval), Nx)) else: - paramRange = np.linspace(minval, maxval, Nx) + paramRange : np.ndarray = np.linspace(minval, maxval, Nx) + + # All done! return Nx, paramRange -getParam1DSpace = {'list': get_1dspace_from_list, - 'uniform': create_uniform_1dspace} + + +# A macro that allows us to switch function we use to generate generate a parameter's range. +getParam1DSpace : dict[str, callable] = {'list' : get_1dspace_from_list, + 'uniform' : create_uniform_1dspace} + + + +# ------------------------------------------------------------------------------------------------- +# ParameterSpace Class +# ------------------------------------------------------------------------------------------------- class ParameterSpace: - param_list = [] - param_name = [] - n_param = 0 - train_space = None - test_space = None - n_init = 0 - test_grid_sizes = [] - test_meshgrid = None - - def __init__(self, config): - assert('parameter_space' in config) - parser = InputParser(config['parameter_space'], name="param_space_input") - - self.param_list = parser.getInput(['parameters'], datatype=list) - self.n_param = len(self.param_list) - - self.param_name = [] + # Initialize class variables + param_list : list[dict] = [] # A list housing the parameter dictionaries (from the yml file) + param_name : list[str] = [] # A list housing the parameter names. + n_param : int = 0 # The number of parameters. + train_space : np.ndarray = None # A 2D array of shape (n_train, n_param) whose i,j element is the j'th parameter value in the i'th combination of training parameters. + test_space : np.ndarray = None # A 2D array of shape (n_test, n_param) whose i,j element is the j'th parameter value in the i'th combination of testing parameters. + n_init : int = 0 # The number of combinations of parameters in the training set??? + test_grid_sizes : list[int] = [] # A list whose i'th element is the number of different values of the i'th parameter in the test instances. + test_meshgrid : tuple[np.ndarray] = None + + + + def __init__(self, config : dict) -> None: + """ + Initializes a ParameterSpace object using the settings passed in the conf dictionary (which + should have been read from a yaml file). + + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + config: This is a dictionary that houses the settings we want to use to run the code. This + should have been read from a yaml file. We assume it contains the following keys. If one + or more keys are tabbed over relative to one key above them, then the one above is a + dictionary and the ones below should be keys within that dictionary. + - parameter_space + - parameters (this should have at least one parameter defined!) + - test_space + - type (should be "grid") + """ + + # Make sure the configuration dictionary has a "parameter_space" setting. This should house + # information about which variables are present in the code, as well as how we want to test + # the various possible parameter values. + assert('parameter_space' in config); + + # Load the parameter_space settings into an InputParser object (which makes it easy to + # fetch setting values). We then fetch the list of parameters. Each parameters has a name, + # min and max, and information on how many instances we want. + parser = InputParser(config['parameter_space'], name = "param_space_input") + self.param_list : list[dict] = parser.getInput(['parameters'], datatype = list) + + # Fetch the parameter names. + self.param_name : list[str] = [] for param in self.param_list: - self.param_name += [param['name']] + self.param_name += [param['name']]; - self.train_space = self.createInitialTrainSpace(self.param_list) - self.n_init = self.train_space.shape[0] + # First, let's fetch the set of possible parameter values. This yields a 2^k x k matrix, + # where k is the number of parameters. The i,j entry of this matrix gives the value of the + # j'th parameter on the i'th instance. + self.train_space = self.createInitialTrainSpace(self.param_list) + self.n_init = self.train_space.shape[0] - test_space_type = parser.getInput(['test_space', 'type'], datatype=str) + # Next, let's make a set of possible parameter values to test at. + test_space_type = parser.getInput(['test_space', 'type'], datatype = str) if (test_space_type == 'grid'): + # Generate the set possible parameter combinations. See the docstring for + # "createTestGridSpace" for details. self.test_grid_sizes, self.test_meshgrid, self.test_space = self.createTestGridSpace(self.param_list) + # All done! return - def n_train(self): + + + def n_train(self) -> int: + """ + Returns the number of combinations of parameters in the training set. + """ + return self.train_space.shape[0] - def n_test(self): + + + def n_test(self) -> int: + """ + Returns the number of combinations of parameters in the testing set. + """ + return self.test_space.shape[0] - - def createInitialTrainSpace(self, param_list): - paramRanges = [] - for param in param_list: - minval = param['min'] - maxval = param['max'] + + + def createInitialTrainSpace(self, param_list : list[dict]) -> np.ndarray: + """ + Sets up a grid of parameter values to train at. Note that we only use the min and max value + of each parameter when setting up this grid. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param_list: A list of parameter dictionaries. Each entry should be a dictionary with the + following keys: + - name + - min + - max + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A 2d array of shape ((2)^k, k), where k is the number of parameters (k == len(param_list)). + The i'th column is the flattened i'th mesh_grid array we when we create a mesh grid using + the min and max value of each parameter as the argument. See "createHyperMeshGrid" for + details. + + Specifically, we return exactly what "createHyperGridSpace" returns. See the doc-string + for that function for further details. + """ + + # We need to know the min and max value for each parameter to set up the grid of possible + # parameter values. + paramRanges : list[np.ndarray] = [] + + for param in param_list: + # Fetch the min, max value of the current parameter. + minval : float = param['min'] + maxval : float = param['max'] + + # Store these values in an array and add them to the list. paramRanges += [np.array([minval, maxval])] - mesh_grids = self.createHyperMeshGrid(paramRanges) + # Use the ranges to set up a grid of possible parameter values. + mesh_grids : tuple[np.ndarray] = self.createHyperMeshGrid(paramRanges) + + # Now combine these grids into a 2d return self.createHyperGridSpace(mesh_grids) - def createTestGridSpace(self, param_list): - paramRanges = [] - gridSizes = [] + + def createTestGridSpace(self, param_list : list[dict]) -> tuple[list[int], tuple[np.ndarray], np.ndarray]: + """ + This function sets up a grid of parameter values to test at. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param_list: A list of parameter dictionaries. Each dictionary should either use the + "uniform" or "list" format. See create_uniform_1dspace and get_1dspace_from_list, + respectively. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A three element tuple. + + The first is a list whose i'th element specifies the number of distinct values of the i'th + parameter we consider (this is the length of the i'th element of "paramRanges" below). + + The second is a a tuple of k numpy ndarrays (where k = len(param_list)), the i'th one of + which is a k-dimensional array with shape (N0, ... , N{k - 1}), where Ni = + param_list[i].size whose i(0), ... , i(k - 1) element specifies the value of the i'th + parameter in the i(0), ... , i(k - 1)'th unique combination of parameter values. + + The third one is a 2d array of parameter values. It has shape (M, k), where + M = \prod_{i = 0}^{k - 1} param_list[i].size. + """ + + # Set up arrays to hold the parameter values + number of parameter values for each + # parameter. + paramRanges : np.ndarray = [] + gridSizes : list[int] = [] + + # Cycle through the parameters for param in param_list: - Nx, paramRange = getParam1DSpace[param['test_space_type']](param) - gridSizes += [Nx] - paramRanges += [paramRange] + # Fetch the set of possible parameter values (paramRange) + the size of this set (Nx) + Nx, paramRange = getParam1DSpace[param['test_space_type']](param) + + # Add Nx, ParamRange to their corresponding lists + gridSizes += [Nx] + paramRanges += [paramRange] - mesh_grids = self.createHyperMeshGrid(paramRanges) + # Now that we have the set of parameter values for each parameter, turn it into a grid. + mesh_grids : tuple[np.ndarray] = self.createHyperMeshGrid(paramRanges) + + # All done. Return a list specifying the number of possible values for each parameter, the + # grids of parameter values, and the flattened/2d version of it. return gridSizes, mesh_grids, self.createHyperGridSpace(mesh_grids) + + def getParameter(self, param_vector): ''' convert numpy array parameter vector to a dict. @@ -91,63 +318,210 @@ def getParameter(self, param_vector): return param - def createHyperMeshGrid(self, param_ranges): + + + def createHyperMeshGrid(self, param_ranges : list[np.ndarray]) -> tuple[np.ndarray]: ''' - param_ranges: list of numpy 1d arrays, each corresponding to 1d parameter grid space. - The list size is equal to the number of parameters. - - Output: paramSpaces - - tuple of numpy nd arrays, corresponding to each parameter. - Dimension of the array equals to the number of parameters + This function generates arrays of parameter values. Specifically, if there are k + parameters (param_ranges has k elements), then we return k k-d arrays, the i'th one of + which is a k-dimensional array whose i(0), ... , i(k - 1) element specifies the value of + the i'th variable in the i(0), ... , i(k - 1)'th unique combination of parameter values. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param_ranges: list of numpy 1d arrays, each corresponding to 1d parameter grid space. The + i'th element of this list should be a 2-element numpy.ndarray object housing the max and + min value for the i'th parameter. The list size should equal the number of parameters. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + the "paramSpaces" tuple. This is a tuple of numpy ndarray objects, the i'th one of which + gives the grid of parameter values for the i'th parameter. Specifically, if there are + k parameters and if param_range[i].size = Ni, then the j'th return array has shape + (N0, ... , N{k - 1}) and the i(0), ... , i(k - 1) element of this array houses the i(j)'th + value of the j'th parameter. + + Thus, if there are k parameters, the returned tuple has k elements, each one of + which is an array of shape (N0, ... , N{k - 1}). ''' - args = () + + # Fetch the ranges, add them to a tuple (this is what the meshgrid function needs). + args : tuple[np.ndarray] = () for paramRange in param_ranges: args += (paramRange,) - paramSpaces = np.meshgrid(*args, indexing='ij') + # Use numpy's meshgrid function to generate the grids of parameter values. + paramSpaces : tuple[np.ndarray] = np.meshgrid(*args, indexing='ij') + + # All done! return paramSpaces - def createHyperGridSpace(self, mesh_grids): + + + def createHyperGridSpace(self, mesh_grids : tuple[np.ndarray]) -> np.ndarray: ''' - mesh_grids: tuple of numpy nd arrays, corresponding to each parameter. - Dimension of the array equals to the number of parameters - - Output: param_grid - - numpy 2d array of size (grid size x number of parameters). + Flattens the mesh_grid numpy.ndarray objects returned by createHyperMeshGrid and combines + them into a single 2d array of shape (grid size, number of parameters) (see below). + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- - grid size is the size of a numpy nd array. + mesh_grids: tuple of numpy nd arrays, corresponding to each parameter. This should ALWAYS + be the output of the "CreateHyperMeshGrid" function. See the return section of that + function's docstring for details. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + The param_grid. This is a 2d numpy.ndarray object of shape (grid size, number of + parameters). If each element of mesh_grids is a numpy.ndarray object of shape (N(1), ... , + N(k)) (k parameters), then (grid size) = N(1)*N(2)*...*N(k) and (number of parameters) = k. ''' - param_grid = None + + # For each parameter, we flatten its mesh_grid into a 1d array (of length (grid size)). We + # horizontally stack these flattened grids to get the final param_grid array. + param_grid : np.ndarray = None for k, paramSpace in enumerate(mesh_grids): + # Special treatment for the first parameter to initialize param_grid if (k == 0): - param_grid = paramSpace.reshape(-1, 1) + param_grid : np.ndarray = paramSpace.reshape(-1, 1) continue - param_grid = np.hstack((param_grid, paramSpace.reshape(-1, 1))) + # Flatten the new mesh grid and append it to the grid. + param_grid : np.ndarray = np.hstack((param_grid, paramSpace.reshape(-1, 1))) + # All done! return param_grid - def appendTrainSpace(self, param): + + + def appendTrainSpace(self, param : np.ndarray) -> None: + """ + Adds a new parameter to self's train space attribute. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param: A 1d numpy ndarray object. It should have shape (n_param) and should hold a + parameter value that we want to add to the training set. + + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + # Make sure param has n_param components/can be appended to the set of training parameters. assert(self.train_space.shape[1] == param.size) - self.train_space = np.vstack((self.train_space, param)) + # Add the new parameter to the training space by appending it as a new row to + # self.train_space + self.train_space : np.ndarray = np.vstack((self.train_space, param)) + + # All done! return - def export(self): - dict_ = {'train_space': self.train_space, - 'test_space': self.test_space, - 'test_grid_sizes': self.test_grid_sizes, - 'test_meshgrid': self.test_meshgrid, - 'n_init': self.n_init} + + + def export(self) -> dict: + """ + This function packages the testing/training examples into a dictionary, which it returns. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + None! + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A dictionary with 4 keys. Below is a list of the keys with a short description of each + corresponding value. + train_space: self.train_space, a 2d array of shape (n_train, n_param) whose i,j element + holds the value of the j'th parameter in the i'th training case. + + test_space: self.test_space, a 2d array of shape (n_test, n_param) whose i,j element + holds the value of the j'th parameter in the i'th testing case. + + test_grid_sizes: A list whose i'th element specifies how many distinct parameter values + we use for the i'th parameter. + + test_meshgrid: a tuple of n_param numpy.ndarray array objects whose i'th element is a + n_param-dimensional array whose i(1), i(2), ... , i(n_param) element holds the value of + the i'th parameter in the i(1), ... , i(n_param) combination of parameter values in the + testing test. + + n_init: The number of combinations of training parameters in the training set. + """ + + # Build the dictionary + dict_ = {'train_space' : self.train_space, + 'test_space' : self.test_space, + 'test_grid_sizes' : self.test_grid_sizes, + 'test_meshgrid' : self.test_meshgrid, + 'n_init' : self.n_init} + + # All done! return dict_ - def load(self, dict_): - self.train_space = dict_['train_space'] - self.test_space = dict_['test_space'] - self.test_grid_sizes = dict_['test_grid_sizes'] - self.test_meshgrid = dict_['test_meshgrid'] - - assert(self.n_init == dict_['n_init']) - assert(self.train_space.shape[1] == self.n_param) - assert(self.test_space.shape[1] == self.n_param) + + + def load(self, dict_ : dict) -> None: + """ + This function builds a parameter space object from a dictionary. This dictionary should + be one that was returned by th export method, or a loaded copy of a dictionary that was + returned by the export method. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + dict_: This should be a dictionary with the following keys: + - train_space + - test_space + - test_grid_sizes + - test_meshgrid + - n_init + This dictionary should have been returned by the export method. We use the values in this + dictionary to set up self. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + # Extract information from the dictionary. + self.train_space : np.ndarray = dict_['train_space'] + self.test_space : np.ndarray = dict_['test_space'] + self.test_grid_sizes : list[int] = dict_['test_grid_sizes'] + self.test_meshgrid : tuple[np.ndarray] = dict_['test_meshgrid'] + + # Run checks + assert(self.n_init == dict_['n_init']) + assert(self.train_space.shape[1] == self.n_param) + assert(self.test_space.shape[1] == self.n_param) + + # All done! return From 362344ad0546640126b1efcd0e20f536f3fb9ea7 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 23 Oct 2024 15:52:07 -0700 Subject: [PATCH 05/22] Added documentation to gp.py --- src/lasdi/gp.py | 184 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 152 insertions(+), 32 deletions(-) diff --git a/src/lasdi/gp.py b/src/lasdi/gp.py index 3af98ff..19e7763 100644 --- a/src/lasdi/gp.py +++ b/src/lasdi/gp.py @@ -1,80 +1,200 @@ -import numpy as np -from sklearn.gaussian_process.kernels import ConstantKernel, Matern, RBF -from sklearn.gaussian_process import GaussianProcessRegressor +import numpy as np +from sklearn.gaussian_process.kernels import ConstantKernel, Matern, RBF +from sklearn.gaussian_process import GaussianProcessRegressor -def fit_gps(X, Y): - ''' - Trains each GP given the interpolation dataset. - X: (n_train, n_param) numpy 2d array - Y: (n_train, n_coef) numpy 2d array + + +# ------------------------------------------------------------------------------------------------- +# Gaussian Process functions! +# ------------------------------------------------------------------------------------------------- + +def fit_gps(X : np.ndarray, Y : np.ndarray) -> list[GaussianProcessRegressor]: + """ + Trains a GP for each column of Y. If Y has shape N x k, then we train k GP regressors. In this + case, we assume that X has shape N x M. Thus, the Input to the GP is in \mathbb{R}^M. For each + k, we train a GP where the i'th row of X is the input and the i,k component of Y is the + corresponding target. Thus, we return a list of k GP Regressor objects, the k'th one of which + makes predictions for the k'th coefficient in the latent dynamics. + We assume each target coefficient is independent with each other. - gp_dictionnary is a dataset containing the trained GPs (as sklearn objects) - ''' - n_coef = 1 if (Y.ndim == 1) else Y.shape[1] + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + X: A 2d numpy array of shape (n_train, input_dim), where n_train is the number of training + examples and input_dim is the number of components in each input (e.g., the number of + parameters) + + Y: A 2d numpy array of shape (n_train, n_coef), where n_train is the number of training + examples and n_coef is the number of coefficients in the latent dynamics. + + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + A list of trained GP regressor objects. If Y has k columns, then the returned list has k + elements. It's i'th element holds a trained GP regressor object whose training inputs are the + columns of X and whose corresponding target values are the elements of the i'th column of Y. + """ + + # Determine the number of components (columns) of Y. Since this is a regression task, we will + # perform a GP regression fit on each component (column) of Y. + n_coef : int = 1 if (Y.ndim == 1) else Y.shape[1] + + # Transpose Y so that each row corresponds to a particular coefficient. This allows us to + # iterate over the coefficients by iterating through the rows of Y. if (n_coef > 1): Y = Y.T + # Sklearn requires X to be a 2D array... so make sure this holds. if X.ndim == 1: X = X.reshape(-1, 1) - gp_dictionnary = [] + # Initialize a list to hold the trained GP objects. + gp_list : list[GaussianProcessRegressor] = [] + # Cycle through the rows of Y (which we transposed... so this just cycles through the + # coefficients) for yk in Y: + # Make the kernel. # kernel = ConstantKernel() * Matern(length_scale_bounds = (0.01, 1e5), nu = 1.5) - kernel = ConstantKernel() * RBF(length_scale_bounds = (0.1, 1e5)) + kernel = ConstantKernel() * RBF(length_scale_bounds = (0.1, 1e5)) - gp = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = 10, random_state = 1) - gp.fit(X, yk) + # Initialize the GP object. + gp = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = 10, random_state = 1) - gp_dictionnary += [gp] + # Fit it to the data (train), then add it to the list of trained GPs + gp.fit(X, yk) + gp_list += [gp] - return gp_dictionnary + # All done! + return gp_list -def eval_gp(gp_dictionnary, param_grid): - ''' +def eval_gp(gp_list : list[GaussianProcessRegressor], param_grid : np.ndarray) -> tuple: + """ Computes the GPs predictive mean and standard deviation for points of the parameter space grid - ''' - n_coef = len(gp_dictionnary) + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + gp_list: a list of trained GP regressor objects. The number of elements in this list should + match the number of columns in param_grid. The i'th element of this list is a GP regressor + object that predicts the i'th coefficient. + + param_grid: A 2d numpy.ndarray object of shape (number of parameter combination, number of + parameters). The i,j element of this array specifies the value of the j'th parameter in the + i'th combination of parameters. We use this as the testing set for the GP evaluation. + + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + A two element tuple. Both are 2d numpy arrays of shape (number of parameter combinations, + number of coefficients). The two arrays hold the predicted means and std's for each parameter + at each training example, respectively. + + Thus, the i,j element of the first return variable holds the predicted mean of the j'th + coefficient in the latent dynamics at the i'th training example. Likewise, the i,j element of + the second return variable holds the standard deviation in the predicted distribution for the + j'th coefficient in the latent dynamics at the i'th combination of parameter values. + """ + + # Fetch the numbers coefficients. Since there is one GP Regressor per SINDy coefficient, this + # just the length of the gp_list. + n_coef : int = len(gp_list) + + # Fetch the number of parameters, make sure the grid is 2D. if (param_grid.ndim == 1): param_grid = param_grid.reshape(1, -1) n_points = param_grid.shape[0] + # Initialize arrays to hold the mean, STD. pred_mean, pred_std = np.zeros([n_points, n_coef]), np.zeros([n_points, n_coef]) - for k, gp in enumerate(gp_dictionnary): + + # Cycle through the GPs (one for each coefficient in the SINDy coefficients!). + for k, gp in enumerate(gp_list): + # Make predictions using the parameters in the param_grid. pred_mean[:, k], pred_std[:, k] = gp.predict(param_grid, return_std = True) + # All done! return pred_mean, pred_std -def sample_coefs(gp_dictionnary, param, n_samples): - ''' - Generates sample sets of ODEs for one given parameter. - coef_samples is a list of length n_samples, where each terms is a matrix of SINDy coefficients sampled from the GP predictive - distributions +def sample_coefs( gp_list : list[GaussianProcessRegressor], + param : np.ndarray, + n_samples : int): + """ + Generates sets of ODE (SINDy) coefficients sampled from the predictive distribution for those + coefficients at the specified parameter value (parma). Specifically, for the k'th SINDy + coefficient, we draw n_samples samples of the predictive distribution for the k'th coefficient + when param is the parameter. + + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + gp_list: a list of trained GP regressor objects. The number of elements in this list should + match the number of columns in param_grid. The i'th element of this list is a GP regressor + object that predicts the i'th coefficient. + + param: A combination of parameter values. i.e., a single test example. We evaluate each GP in + the gp_list at this parameter value (getting a prediction for each coefficient). + + n_samples: Number of samples of the predicted latent dynamics used to build ensemble of fom + predictions. N_s in the paper. + + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + A 2d numpy ndarray object called coef_samples. It has shape (n_samples, n_coef), where n_coef + is the number of coefficients (length of gp_list). The i,j element of this list is the i'th + sample of the j'th SINDy coefficient. + """ - ''' + # Fetch the number of coefficients (since there is one GP Regressor per coefficient, this is + # just the length of the gp_list). + n_coef : int = len(gp_list) - n_coef = len(gp_dictionnary) - coef_samples = np.zeros([n_samples, n_coef]) + # Initialize an array to hold the coefficient samples. + coef_samples : np.ndarray = np.zeros([n_samples, n_coef]) + # Make sure param is a 2d array with one row, we need this when evaluating the GP Regressor + # object. if param.ndim == 1: param = param.reshape(1, -1) - n_points = param.shape[0] + + # Make sure we only have a single sample. + n_points : int = param.shape[0] assert(n_points == 1) - pred_mean, pred_std = eval_gp(gp_dictionnary, param) + # Evaluate the predicted mean and std at the parameter value. + pred_mean, pred_std = eval_gp(gp_list, param) pred_mean, pred_std = pred_mean[0], pred_std[0] + # Cycle through the samples and coefficients. For each sample of the k'th coefficient, we draw + # a sample from the normal distribution with mean pred_mean[k] and std pred_std[k]. for s in range(n_samples): for k in range(n_coef): coef_samples[s, k] = np.random.normal(pred_mean[k], pred_std[k]) + # All done! return coef_samples \ No newline at end of file From cc1efc2fb1ec4d045f820e1eab844131f571780a Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Thu, 24 Oct 2024 07:39:28 -0700 Subject: [PATCH 06/22] Added documentation to inputs.py I also reformatted enums.py --- src/lasdi/enums.py | 16 ++--- src/lasdi/inputs.py | 160 ++++++++++++++++++++++++++++++++++---------- 2 files changed, 134 insertions(+), 42 deletions(-) diff --git a/src/lasdi/enums.py b/src/lasdi/enums.py index 33b9bd8..ce24936 100644 --- a/src/lasdi/enums.py +++ b/src/lasdi/enums.py @@ -1,13 +1,13 @@ from enum import Enum class NextStep(Enum): - Train = 1 - PickSample = 2 - RunSample = 3 - CollectSample = 4 + Train = 1 + PickSample = 2 + RunSample = 3 + CollectSample = 4 class Result(Enum): - Unexecuted = 1 - Success = 2 - Fail = 3 - Complete = 4 \ No newline at end of file + Unexecuted = 1 + Success = 2 + Fail = 3 + Complete = 4 \ No newline at end of file diff --git a/src/lasdi/inputs.py b/src/lasdi/inputs.py index 0a654b6..fa030bb 100644 --- a/src/lasdi/inputs.py +++ b/src/lasdi/inputs.py @@ -1,67 +1,159 @@ +# ------------------------------------------------------------------------------------------------- +# Imports +# ------------------------------------------------------------------------------------------------- + from warnings import warn -verbose = False +verbose : bool = False + + +# ------------------------------------------------------------------------------------------------- +# Input Parser class +# ------------------------------------------------------------------------------------------------- class InputParser: - dict_ = None - name = "" + """ + A InputParser objects acts as a wrapper around a dictionary of settings. Thus, each setting is + a key and the corresponding value is the setting's value. Because one setting may itself be + a dictionary (we often group settings; each group has a name but several constituent settings), + the underlying dictionary is structured as a sequence of nested dictionaries. This class allows + the user to select a specific setting from that structure by specifying (via a list of strings) + where in that nested structure the desired setting lives. + """ + dict_ : dict = None + name : str = "" + + + + def __init__(self, dict : dict, name : str = "") -> None: + """" + Initializes an InputParser object by setting the underlying dictionary of settings as an + attribute. + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + dict: The dictionary of settings. To avoid the risk of the user accidentally changing one + of the settings after wrapping it, we store a deep copy of dict in self. + """ - def __init__(self, dict, name = ""): + # A shallow copy could cause issues if the user changes dict's keys/values after + # initializing this object. We store a deep copy to avoid this risk. from copy import deepcopy self.dict_ = deepcopy(dict) + self.name = name - return - def getInput(self, keys, fallback=None, datatype=None): + + + def getInput(self, keys : list, fallback = None, datatype = None): ''' - Find the value corresponding to the list of keys. - If the specified keys do not exist, use the fallback value. - If the fallback value does not exist, returns an error. - If the datatype is specified, enforce the output value has the right datatype. + A InputParser object acts as a wrapper around a dictionary of settings. That is, self.dict_ + is structured as a nested family of dictionaries. Each setting corresponds to a key in + self.dict_. The setting's value is the corresponding value in self.dict_. In many cases, + a particular setting may be nested within others. That is, a setting's value may itself be + another dictionary housing various sub-settings. This function allows us to fetch a + specific setting from this nested structure. + + Specifically, we specify a list of strings. keys[0] should be a key in self.dict_ + If so, we set val = self.dict_[keys[0]]. If there are more keys, then val should be a + dictionary and keys[1] should be a key in this dictionary. In this case, we replace val + with val[key[1]] and so on. This continues until we have exhausted all keys. There is one + important exception: + + If at some point in the process, there are more keys but val is not a dictionary, or if + there are more keys and val is a dictionary but the next key is not a key in that + dictionary, then we return the fallback value. If the fallback value does not exist, + returns an error. + + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + keys: A list of keys we want to fetch from self.dict. keys[0] should be a key in self.dict_ + If so, we set val = self.dict_[keys[0]]. If there are more keys, then val should be a + dictionary and keys[1] should be a key in this dictionary. In this case, we replace val + with val[key[1]] and so on. This continues until we have exhausted all keys. + + fallback: A sort of default value. If at some point, val is not a dictionary (and there are + more keys) or val is a dictionary but the next key is not a valid key in that dictionary, + then we return the fallback value. + + datatype: If not None, then we require that the final val has this datatype. If the final + val does not have the desired datatype, we raise an exception. + + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + The final val value as outlined by the process described above. ''' + + # Concatenate the keys together. This is for debugging purposes. keyString = "" for key_ in keys: keyString += key_ + "/" + # Begin by initializing val to self.dict_ val = self.dict_ + + # Cycle through the keys for key in keys: + + # Check if the current key is a key in val (this assumes val is a dictionary). If so, + # update val. Otherwise, return the fallback (if it is present) or raise an exception. if key in val: val = val[key] elif (fallback != None): return fallback else: raise RuntimeError("%s does not exist in the input dictionary %s!" % (keyString, self.name)) - + + # Check if the fallback and final val have the same type. if (fallback != None): if (type(val) != type(fallback)): warn("%s does not match the type with the fallback value %s!" % (str(type(val)), str(type(fallback)))) - + + # Check thast the final val matches the desired datatype if (datatype != None): if (type(val) != datatype): raise RuntimeError("%s does not match the specified datatype %s!" % (str(type(val)), str(datatype))) else: if verbose: warn("InputParser Warning: datatype is not checked.\n key: %s\n value type: %s" % (keys, type(val))) + + # All done! return val -def getDictFromList(list_, inputDict): - ''' - get a dict with {key: val} from a list of dicts - NOTE: it returns only the first item in the list, - even if the list has more than one dict with {key: val}. - ''' - dict_ = None - for item in list_: - isDict = True - for key, val in inputDict.items(): - if key not in item: - isDict = False - break - if (item[key] != val): - isDict = False - break - if (isDict): - dict_ = item - break - if (dict_ == None): - raise RuntimeError('Given list does not have a dict with {%s: %s}!' % (key, val)) - return dict_ + + +# ------------------------------------------------------------------------------------------------- +# Unused: getDictFromList function +# ------------------------------------------------------------------------------------------------- + +#def getDictFromList(list_, inputDict): +# ''' +# get a dict with {key: val} from a list of dicts +# NOTE: it returns only the first item in the list, +# even if the list has more than one dict with {key: val}. +# ''' +# dict_ = None +# for item in list_: +# isDict = True +# for key, val in inputDict.items(): +# if key not in item: +# isDict = False +# break +# if (item[key] != val): +# isDict = False +# break +# if (isDict): +# dict_ = item +# break +# if (dict_ == None): +# raise RuntimeError('Given list does not have a dict with {%s: %s}!' % (key, val)) +# return dict_ \ No newline at end of file From b865d871bce12f6b32a74ca9f34e601e123b2f67 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Thu, 24 Oct 2024 08:24:31 -0700 Subject: [PATCH 07/22] Minor comment tweaks to latent_space.py --- src/lasdi/latent_space.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index 2e78542..d7f547c 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -53,7 +53,8 @@ def initial_condition_latent(param_grid : np.ndarray, ----------------------------------------------------------------------------------------------- param_grid: A 2d numpy.ndarray object of shape (number of parameter combination) x (number of - parameters). + parameters). The i,j element of this array holds the value of the j'th parameter in the i'th + combination of parameters. physics: A "Physics" object that stores the datasets for each parameter combination. From fc14c5e57254047c8c0a90c90410932f8d2027bc Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Thu, 24 Oct 2024 08:31:34 -0700 Subject: [PATCH 08/22] Started adding documentation to gplasdi.py --- src/lasdi/gplasdi.py | 334 ++++++++++++++++++++++++++++++--- src/lasdi/physics/burgers1d.py | 7 +- 2 files changed, 314 insertions(+), 27 deletions(-) diff --git a/src/lasdi/gplasdi.py b/src/lasdi/gplasdi.py index 5131f24..68fc64a 100644 --- a/src/lasdi/gplasdi.py +++ b/src/lasdi/gplasdi.py @@ -1,13 +1,52 @@ -from .gp import * -from .latent_space import * -from .enums import * -from .timing import Timer -import torch -import time -import numpy as np +# ------------------------------------------------------------------------------------------------- +# Imports +# ------------------------------------------------------------------------------------------------- -def average_rom(autoencoder, physics, latent_dynamics, gp_dictionary, param_grid): +import time +import torch +import numpy as np + +from .gp import eval_gp, sample_coefs, fit_gps +from .latent_space import initial_condition_latent, Autoencoder +from .enums import NextStep, Result +from .physics import Physics +from .latent_dynamics import LatentDynamics +from .timing import Timer +from .param import ParameterSpace + + + +def average_rom(autoencoder : Autoencoder, + physics : Physics, + latent_dynamics : LatentDynamics, + gp_dictionary : dict, + param_grid : np.ndarray): + """ + + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + physics: A "Physics" object that stores the datasets for each parameter combination. + + autoencoder: The actual autoencoder object that we use to map the ICs into the latent space. + + param_grid: A 2d numpy.ndarray object of shape (number of parameter combination) x (number of + parameters). The i,j element of this array holds the value of the j'th parameter in the i'th + combination of parameters. + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + A list of numpy ndarray objects whose i'th element holds the latent space initial condition + for the i'th set of parameters in the param_grid. That is, if we let U0_i denote the fom IC for + the i'th set of parameters, then the i'th element of the returned list is Z0_i = encoder(U0_i). + """ if (param_grid.ndim == 1): param_grid = param_grid.reshape(1, -1) n_test = param_grid.shape[0] @@ -22,7 +61,9 @@ def average_rom(autoencoder, physics, latent_dynamics, gp_dictionary, param_grid return Zis -def sample_roms(autoencoder, physics, latent_dynamics, gp_dictionary, param_grid, n_samples): + + +def sample_roms(autoencoder : Autoencoder, physics, latent_dynamics, gp_dictionary, param_grid, n_samples): ''' Collect n_samples of ROM trajectories on param_grid. gp_dictionary: list of Gaussian process regressors (size of n_test) @@ -49,30 +90,77 @@ def sample_roms(autoencoder, physics, latent_dynamics, gp_dictionary, param_grid return Zis -def get_fom_max_std(autoencoder, Zis): - ''' - Computes the maximum standard deviation accross the parameter space grid and finds the corresponding parameter location +def get_fom_max_std(autoencoder : Autoencoder, Zis : np.ndarray) -> int: + """ + Computes the maximum standard deviation across the parameter space grid and finds the + corresponding parameter index. - ''' - # TODO(kevin): currently this evaluate pointwise maximum standard deviation. + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + autoencoder: The autoencoder. We assume the solved dynamics (whose frames are stored in Zis) + take place in the autoencoder's latent space. We use this to decode the solution frames. + + Zis: A 4d numpy array of shape (n_test, n_samples, n_t, n_z) whose i, j, k, l element holds + the l'th component of the k'th frame of the solution to the latent dynamics when we use the + j'th sample of latent coefficients drawn from the posterior distribution for the i'th testing + parameter. + + + + ----------------------------------------------------------------------------------------------- + Returns: + ----------------------------------------------------------------------------------------------- + + An integer. The index of the testing parameter that gives the largest standard deviation. + Specifically, for each testing parameter, we compute the STD of each component of the fom + solution at each frame generated by samples from the posterior coefficient distribution for + that parameter. We compute the maximum of these STDs and pair that number with the parameter. + We then return the index of the parameter whose corresponding maximum std (the number we pair + with it) is greatest. + """ + # TODO(kevin): currently this evaluate point-wise maximum standard deviation. # is this a proper metric? we might want to consider an average, or L2 norm of std. - max_std = 0 - for m, Zi in enumerate(Zis): - Z_m = torch.Tensor(Zi) - X_pred_m = autoencoder.decoder(Z_m).detach().numpy() - X_pred_m_std = X_pred_m.std(0) - max_std_m = X_pred_m_std.max() + max_std : float = 0.0 + # Cycle through the testing parameters. + for m, Zi in enumerate(Zis): + # Zi is a 3d tensor of shape (n_samples, n_t, n_z), where n_samples is the number of + # samples of the posterior distribution per parameter, n_t is the number of time steps in + # the latent dynamics solution, and n_z is the dimension of the latent space. The i,j,k + # element of Zi is the k'th component of the j'th frame of the solution to the latent + # dynamics when the latent dynamics uses the i'th set of sampled parameter values. + Z_m : torch.Tensor = torch.Tensor(Zi) + + # Now decode the frames. + X_pred_m : np.ndarray = autoencoder.decoder(Z_m).detach().numpy() + + # Compute the standard deviation across the sample axis. This gives us an array of shape + # (n_t, n_fom) whose i,j element holds the (sample) standard deviation of the j'th component + # of the i'th frame of the fom solution. In this case, the sample distribution consists of + # the set of j'th components of i'th frames of fom solutions (one for each sample of the + # posterior coefficients). + X_pred_m_std : np.ndarray = X_pred_m.std(0) + + # Now compute the maximum standard deviation across frames/fom components. + max_std_m : np.float32 = X_pred_m_std.max() + + # If this is bigger than the biggest std we have seen so far, update the maximum. if max_std_m > max_std: - m_index = m - max_std = max_std_m + m_index : int = m + max_std : float = max_std_m + # Report the index of the testing parameter that gave the largest maximum std. return m_index + # move optimizer parameters to device def optimizer_to(optim, device): for param in optim.state.values(): @@ -92,7 +180,7 @@ class BayesianGLaSDI: X_train = torch.Tensor([]) X_test = torch.Tensor([]) - def __init__(self, physics, autoencoder, latent_dynamics, param_space, config): + def __init__(self, physics, autoencoder : Autoencoder, latent_dynamics, param_space, config): ''' @@ -280,4 +368,202 @@ def load(self, dict_): self.optimizer.load_state_dict(dict_['optimizer']) if (self.device != 'cpu'): optimizer_to(self.optimizer, self.device) - return \ No newline at end of file + return + + + +class BayesianGLaSDI: + X_train = None + + def __init__(self, physics, autoencoder, latent_dynamics, model_parameters): + + ''' + + This class runs a full GPLaSDI training. It takes into input the autoencoder defined as a PyTorch object and the + dictionnary containing all the training parameters. + The "train" method with run the active learning training loop, compute the reconstruction and SINDy loss, train the GPs, + and sample a new FOM data point. + + ''' + + self.autoencoder = autoencoder + self.latent_dynamics = latent_dynamics + self.physics = physics + self.param_space = physics.param_space + self.timer = Timer() + + self.n_samples = model_parameters['n_samples'] + self.lr = model_parameters['lr'] + self.n_iter = model_parameters['n_iter'] + self.n_greedy = model_parameters['n_greedy'] + self.max_greedy_iter = model_parameters['max_greedy_iter'] + self.sindy_weight = model_parameters['sindy_weight'] + self.coef_weight = model_parameters['coef_weight'] + + self.optimizer = torch.optim.Adam(autoencoder.parameters(), lr = self.lr) + self.MSE = torch.nn.MSELoss() + + self.path_checkpoint = model_parameters['path_checkpoint'] + self.path_results = model_parameters['path_results'] + + from os.path import dirname + from pathlib import Path + Path(dirname(self.path_checkpoint)).mkdir(parents=True, exist_ok=True) + Path(dirname(self.path_results)).mkdir(parents=True, exist_ok=True) + + device = model_parameters['device'] if 'device' in model_parameters else 'cpu' + if (device == 'cuda'): + assert(torch.cuda.is_available()) + self.device = device + elif (device == 'mps'): + assert(torch.backends.mps.is_available()) + self.device = device + else: + self.device = 'cpu' + + self.best_loss = np.Inf + self.restart_iter = 0 + + return + + def train(self): + assert(self.X_train is not None) + assert(self.X_train.size(0) == self.param_space.n_train) + + device = self.device + autoencoder_device = self.autoencoder.to(device) + X_train_device = self.X_train.to(device) + + from pathlib import Path + Path(self.path_checkpoint).mkdir(parents=True, exist_ok=True) + Path(self.path_results).mkdir(parents=True, exist_ok=True) + + ps = self.param_space + ld = self.latent_dynamics + + for iter in range(self.restart_iter, self.n_iter): + self.timer.start("train_step") + + self.optimizer.zero_grad() + Z = autoencoder_device.encoder(X_train_device) + X_pred = autoencoder_device.decoder(Z) + Z = Z.cpu() + + loss_ae = self.MSE(X_train_device, X_pred) + coefs, loss_sindy, loss_coef = ld.calibrate(Z, self.physics.dt, compute_loss=True, numpy=True) + + max_coef = np.abs(coefs).max() + + loss = loss_ae + self.sindy_weight * loss_sindy / ps.n_train + self.coef_weight * loss_coef / ps.n_train + + loss.backward() + self.optimizer.step() + + if loss.item() < self.best_loss: + torch.save(autoencoder_device.state_dict(), self.path_checkpoint + '/' + 'checkpoint.pt') + self.best_coefs = coefs + self.best_loss = loss.item() + + print("Iter: %05d/%d, Loss: %3.10f, Loss AE: %3.10f, Loss SI: %3.10f, Loss COEF: %3.10f, max|c|: %04.1f, " + % (iter + 1, self.n_iter, loss.item(), loss_ae.item(), loss_sindy.item(), loss_coef.item(), max_coef), + end = '') + + if ps.n_train < 6: + print('Param: ' + str(np.round(ps.train_space[0, :], 4)), end = '') + + for i in range(1, ps.n_train - 1): + print(', ' + str(np.round(ps.train_space[i, :], 4)), end = '') + print(', ' + str(np.round(ps.train_space[-1, :], 4))) + + else: + print('Param: ...', end = '') + for i in range(5): + print(', ' + str(np.round(ps.train_space[-6 + i, :], 4)), end = '') + print(', ' + str(np.round(ps.train_space[-1, :], 4))) + + self.timer.end("train_step") + + if ((iter > self.restart_iter) and (iter < self.max_greedy_iter) and (iter % self.n_greedy == 0)): + + if (self.best_coefs.shape[0] == ps.n_train): + coefs = self.best_coefs + + new_sample = self.get_new_sample_point(coefs) + + # TODO(kevin): implement save/load the new parameter + ps.appendTrainSpace(new_sample) + self.restart_iter = iter + next_step, result = NextStep.RunSample, Result.Success + print('New param: ' + str(np.round(new_sample, 4)) + '\n') + # self.timer.end("new_sample") + return result, next_step + + self.timer.start("finalize") + + if (self.best_coefs.shape[0] == ps.n_train): + coefs = self.best_coefs + + gp_dictionnary = fit_gps(ps.train_space, coefs) + + bglasdi_results = {'autoencoder_param': self.autoencoder.state_dict(), 'final_X_train': self.X_train, + 'coefs': coefs, 'gp_dictionnary': gp_dictionnary, 'lr': self.lr, 'n_iter': self.n_iter, + 'n_greedy': self.n_greedy, 'sindy_weight': self.sindy_weight, 'coef_weight': self.coef_weight, + 'n_samples' : self.n_samples, + } + bglasdi_results['physics'] = self.physics.export() + bglasdi_results['parameters'] = self.param_space.export() + # TODO(kevin): restart capability for timer. + bglasdi_results['timer'] = self.timer.export() + bglasdi_results['latent_dynamics'] = self.latent_dynamics.export() + + date = time.localtime() + date_str = "{month:02d}_{day:02d}_{year:04d}_{hour:02d}_{minute:02d}" + date_str = date_str.format(month = date.tm_mon, day = date.tm_mday, year = date.tm_year, hour = date.tm_hour + 3, minute = date.tm_min) + np.save(self.path_results + '/' + 'bglasdi_' + date_str + '.npy', bglasdi_results) + + self.timer.end("finalize") + self.timer.print() + + next_step, result = None, Result.Complete + return result, next_step + + def get_new_sample_point(self, coefs): + self.timer.start("new_sample") + + print('\n~~~~~~~ Finding New Point ~~~~~~~') + # TODO(kevin): william, this might be the place for new sampling routine. + + ae = self.autoencoder.cpu() + ps = self.param_space + ae.load_state_dict(torch.load(self.path_checkpoint + '/' + 'checkpoint.pt')) + + Z0 = initial_condition_latent(ps.test_space, self.physics, ae) + + gp_dictionnary = fit_gps(ps.train_space, coefs) + + coef_samples = [sample_coefs(gp_dictionnary, ps.test_space[i], self.n_samples) for i in range(ps.n_test)] + + Zis = np.zeros([ps.n_test, self.n_samples, self.physics.nt, ae.n_z]) + for i, Zi in enumerate(Zis): + z_ic = Z0[i] + for j, coef_sample in enumerate(coef_samples[i]): + Zi[j] = self.latent_dynamics.simulate(coef_sample, z_ic, self.physics.t_grid) + + m_index = get_fom_max_std(ae, Zis) + + self.timer.end("new_sample") + return ps.test_space[m_index, :] + + def sample_fom(self): + + new_foms = self.param_space.n_train - self.X_train.size(0) + assert(new_foms > 0) + new_params = self.param_space.train_space[-new_foms:, :] + + if not self.physics.offline: + new_X = self.physics.generate_solutions(new_params) + + self.X_train = torch.cat([self.X_train, new_X], dim = 0) + else: + # TODO(kevin): interface for offline FOM simulation + raise RuntimeError("Offline FOM simulation is not supported yet!") \ No newline at end of file diff --git a/src/lasdi/physics/burgers1d.py b/src/lasdi/physics/burgers1d.py index ad2dd32..f7c3bdc 100644 --- a/src/lasdi/physics/burgers1d.py +++ b/src/lasdi/physics/burgers1d.py @@ -21,9 +21,10 @@ def __init__(self, cfg, param_name=None): self.offline = parser.getInput(['offline_driver'], fallback=False) - self.nt = parser.getInput(['number_of_timesteps'], datatype=int) - self.grid_size = parser.getInput(['grid_size'], datatype=list) - self.qgrid_size = self.grid_size + self.nt : int = parser.getInput(['number_of_timesteps'], datatype = int) + self.grid_size : list[int] = parser.getInput(['grid_size'], datatype = list) + self.qgrid_size : list[int] = self.grid_size + assert(self.dim == len(self.grid_size)) self.xmin = parser.getInput(['xmin'], datatype=float) From 370002a1082aa1c9ab4ceb6520d43c1f5002dd43 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Thu, 24 Oct 2024 09:14:33 -0700 Subject: [PATCH 09/22] Added documentation to timing.py --- src/lasdi/param.py | 2 +- src/lasdi/timing.py | 217 +++++++++++++++++++++++++++++--------------- 2 files changed, 147 insertions(+), 72 deletions(-) diff --git a/src/lasdi/param.py b/src/lasdi/param.py index 208a311..ec75e5c 100644 --- a/src/lasdi/param.py +++ b/src/lasdi/param.py @@ -487,7 +487,7 @@ def export(self) -> dict: def load(self, dict_ : dict) -> None: """ This function builds a parameter space object from a dictionary. This dictionary should - be one that was returned by th export method, or a loaded copy of a dictionary that was + be one that was returned by the export method, or a loaded copy of a dictionary that was returned by the export method. diff --git a/src/lasdi/timing.py b/src/lasdi/timing.py index 978045a..f36388b 100644 --- a/src/lasdi/timing.py +++ b/src/lasdi/timing.py @@ -1,135 +1,210 @@ -""" - -.. _NumPy docstring standard: - https://numpydoc.readthedocs.io/en/latest/format.html#docstring-standard - -""" +# ------------------------------------------------------------------------------------------------- +# Imports +# ------------------------------------------------------------------------------------------------- from time import perf_counter -class Timer: - """A light-weight timer class. - """ - def __init__(self): - - self.names = {} - """:obj:`dict(str:int)`: Dictionary that maps job names to job indices.""" - self.calls = [] - """:obj:`list(int)`: List that stores the number of calls for each job.""" +# ------------------------------------------------------------------------------------------------- +# Timer class +# ------------------------------------------------------------------------------------------------- - self.times = [] - """:obj:`list(float)`: List that stores the total time for each job.""" - - self.starts = [] - """ - :obj:`list(float)`: List that stores the start time for each job. - If the job is not running, :obj:`None` is stored instead. - """ +class Timer: + def __init__(self): + # Set instance variables + self.names : dict = {} # Dictionary of named timers. key = name, value = index + self.calls : list = [] # k'th element = # of times we have called the k'th timer + self.times : list = [] # k'th element = total time recorded by the k'th timer + self.starts : list = [] # k'th element = start time for the k'th timer (if running) return + - def start(self, name): - """Start a job named :obj:`name`. - If the job is not listed, register the job in the job list. - Args: - name (:obj:`str`): Name of the job to be started. + def start(self, name : str) -> None: + """ + Starts a specific timer. The user must specify the name of the timer they want to start. + The specified timer can not already be running! + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + name: A string specifying the name of the timer you want to start. - Note: - The job must not have started before calling this method. - Returns: - Does not return a value. + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + Nothing! """ + # If the timer does not already exist, initialize it if name not in self.names: + # If this is the k'th timer that we have added to self (which would be the case if + # self.names has k keys), then we store that information in the "names" dictionary. + # We then store other information about that timer (such as its current value) in the + # k'th element of the calls/times/starts lists. self.names[name] = len(self.names) - self.calls += [0] - self.times += [0.0] - self.starts += [None] + self.calls += [0] + self.times += [0.0] + self.starts += [None] - idx = self.names[name] - # check if the job is already being measured + # Fetch the current timer's number. + idx : int = self.names[name] + + # Make sure the current timer has not already started. If so, it is already running and we + # need to raise an exception. if (self.starts[idx] is not None): raise RuntimeError("Timer.start: %s timer is already ticking!" % name) + + # Set the current timer's start element to the time when we started this timer. self.starts[idx] = perf_counter() + + # All done! return - def end(self, name): - """End a job named :obj:`name`. - Increase the number of calls and the runtime for the job. - Args: - name (:obj:`str`): Name of the job to be ended. - Note: - The job must have started before calling this method. + def end(self, name : str) -> None: + """ + Stops a specific timer. The user must specify the name of the timer they want to stop. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- - Returns: - Does not return a value. + name: A string specifying the name of the timer you want to stop. + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! """ + + # Make sure the requested timer actually exists. assert(name in self.names) + # Fetch the requested timers' index. idx = self.names[name] - # check if the job has started. + + # Make sure the requested timer is actually running. if (self.starts[idx] is None): raise RuntimeError("Timer.end: %s start time is not measured yet!" % name) + # Record the time since the current timer started, add it to the running total for this + # timer. Also increment the number of calls to this timer + reset it's start value to + # None (so we can start it again). self.times[idx] += perf_counter() - self.starts[idx] self.calls[idx] += 1 self.starts[idx] = None + + # All done! return - def print(self): - """Print the list of jobs and their number of calls, total time and time per each call. - Returns: - Does not return a value. - """ + def print(self) -> None: + """ + This function reports information on every timer in self. It has no arguments and returns + nothing. + """ + + # Header print("Function name\tCalls\tTotal time\tTime/call\n") + + # Cycle through timers. for name, idx in self.names.items(): print("%s\t%d\t%.3e\t%.3e\n" % (name, self.calls[idx], self.times[idx], self.times[idx] / self.calls[idx])) + + # All done! return - def export(self): - """Export the list of jobs and their number of calls and total time - into a dictionary. - Note: - All jobs must be ended before calling this method. - Returns: - :obj:`dict` that contains "names", "calls", and "times" as keys + def export(self) -> dict: """ + This function extracts the names, calls, and times attributes of self, stores them in a + dictionary, and then returns that dictionary. If you have another dictionary object, + you can passed the returned dictionary to that object's load method to make that object + into an identical copy of self. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + None! + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A dictionary housing the names, calls, and times attributes of self. The returned + dictionary has three keys: + - names + - calls + - times + + names is a dictionary with string keys whose corresponding values are integer indexes. If a + particular timer was the k'th one added to self, then it's value in names will be k. + + calls is a list whose k'th element specifies how many times the k'th timer was stopped. + times is a list whose k'th element specifies the total time recorded on the k'th timer. + """ + + # Make sure that no timers are currently running. for start in self.starts: if (start is not None): raise RuntimeError('Timer.export: cannot export while Timer is still ticking!') - param_dict = {} + # Set up a dictionary to house the timer information. + param_dict : dict = {} + + # Store names, calls, and timers (but not starts... we don't need that) in the dictionary. param_dict["names"] = self.names param_dict["calls"] = self.calls param_dict["times"] = self.times - return param_dict - - def load(self, dict_): - """Load the list of jobs and their number of calls and total time - from a dictionary. - Args: - `dict_` (:obj:`dict`): Dictionary that contains the list of jobs and their calls and times. + # All done! + return param_dict + - Note: - :obj:`dict_['names']`, :obj:`dict_['calls']` and :obj:`dict_['times']` must have the same size. - Returns: - Does not return a value + def load(self, dict_ : dict) -> None: """ + This function de-serializes a timer object, making self into an identical copy of a + previously serialized timer object. Specifically, we replace self's names, calls, and + times attributes using those in the passed dict_. We use this function to restore a + timer object's state after loading from a checkpoint. + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + dict_: This should be a dictionary with three keys: + - names + - calls + - times + the corresponding values should be the names, calls, and times attributes of another timer + object, respectively. We replace self's attributes with those the values in dict_. dict_ + should be the dictionary returned by calling export on a timer object. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ self.names = dict_['names'] self.calls = dict_['calls'] self.times = dict_['times'] From 541ebfa4f5dfdf7c56f6bcae60ff3fc17679b4b6 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Fri, 25 Oct 2024 07:54:29 -0700 Subject: [PATCH 10/22] Added more documentation to gplasdi.py I still need to document the actual gplasdi class, but have done everything else in the file. I also made small tweaks to comments in other files. --- src/lasdi/gp.py | 8 +- src/lasdi/gplasdi.py | 395 +++++++++++++++----------------------- src/lasdi/latent_space.py | 6 +- 3 files changed, 161 insertions(+), 248 deletions(-) diff --git a/src/lasdi/gp.py b/src/lasdi/gp.py index 19e7763..ddc3fc5 100644 --- a/src/lasdi/gp.py +++ b/src/lasdi/gp.py @@ -20,7 +20,6 @@ def fit_gps(X : np.ndarray, Y : np.ndarray) -> list[GaussianProcessRegressor]: We assume each target coefficient is independent with each other. - ----------------------------------------------------------------------------------------------- Arguments ----------------------------------------------------------------------------------------------- @@ -33,7 +32,6 @@ def fit_gps(X : np.ndarray, Y : np.ndarray) -> list[GaussianProcessRegressor]: examples and n_coef is the number of coefficients in the latent dynamics. - ----------------------------------------------------------------------------------------------- Returns ----------------------------------------------------------------------------------------------- @@ -78,12 +76,11 @@ def fit_gps(X : np.ndarray, Y : np.ndarray) -> list[GaussianProcessRegressor]: -def eval_gp(gp_list : list[GaussianProcessRegressor], param_grid : np.ndarray) -> tuple: +def eval_gp(gp_list : list[GaussianProcessRegressor], param_grid : np.ndarray) -> tuple[np.ndarray, np.ndarray]: """ Computes the GPs predictive mean and standard deviation for points of the parameter space grid - ----------------------------------------------------------------------------------------------- Arguments ----------------------------------------------------------------------------------------------- @@ -97,7 +94,6 @@ def eval_gp(gp_list : list[GaussianProcessRegressor], param_grid : np.ndarray) - i'th combination of parameters. We use this as the testing set for the GP evaluation. - ----------------------------------------------------------------------------------------------- Returns ----------------------------------------------------------------------------------------------- @@ -143,7 +139,6 @@ def sample_coefs( gp_list : list[GaussianProcessRegressor], coefficient, we draw n_samples samples of the predictive distribution for the k'th coefficient when param is the parameter. - ----------------------------------------------------------------------------------------------- Arguments @@ -158,7 +153,6 @@ def sample_coefs( gp_list : list[GaussianProcessRegressor], n_samples: Number of samples of the predicted latent dynamics used to build ensemble of fom predictions. N_s in the paper. - ----------------------------------------------------------------------------------------------- diff --git a/src/lasdi/gplasdi.py b/src/lasdi/gplasdi.py index 68fc64a..b51f7b4 100644 --- a/src/lasdi/gplasdi.py +++ b/src/lasdi/gplasdi.py @@ -5,98 +5,209 @@ import time import torch -import numpy as np +import numpy as np +from sklearn.gaussian_process import GaussianProcessRegressor -from .gp import eval_gp, sample_coefs, fit_gps -from .latent_space import initial_condition_latent, Autoencoder -from .enums import NextStep, Result -from .physics import Physics -from .latent_dynamics import LatentDynamics -from .timing import Timer -from .param import ParameterSpace +from .gp import eval_gp, sample_coefs, fit_gps +from .latent_space import initial_condition_latent, Autoencoder +from .enums import NextStep, Result +from .physics import Physics +from .latent_dynamics import LatentDynamics +from .timing import Timer +from .param import ParameterSpace +# ------------------------------------------------------------------------------------------------- +# Simulate latent dynamics +# ------------------------------------------------------------------------------------------------- + def average_rom(autoencoder : Autoencoder, physics : Physics, latent_dynamics : LatentDynamics, - gp_dictionary : dict, + gp_list : list[GaussianProcessRegressor], param_grid : np.ndarray): """ - + This function simulates the latent dynamics for a collection of testing parameters by using + the mean of the posterior distribution for each coefficient's posterior distribution. + Specifically, for each parameter combination, we determine the mean of the posterior + distribution for each coefficient. We then use this mean to simulate the latent dynamics + forward in time (starting from the latent encoding of the fom initial condition for that + combination of coefficients). ----------------------------------------------------------------------------------------------- Arguments ----------------------------------------------------------------------------------------------- + autoencoder: The actual autoencoder object that we use to map the ICs into the latent space. + physics: A "Physics" object that stores the datasets for each parameter combination. + + latent_dynamics: A LatentDynamics object which describes how we specify the dynamics in the + Autoencoder's latent space. - autoencoder: The actual autoencoder object that we use to map the ICs into the latent space. + gp_list: a list of trained GP regressor objects. The number of elements in this list should + match the number of columns in param_grid. The i'th element of this list is a GP regressor + object that predicts the i'th coefficient. param_grid: A 2d numpy.ndarray object of shape (number of parameter combination) x (number of parameters). The i,j element of this array holds the value of the j'th parameter in the i'th - combination of parameters. + combination of parameters. ----------------------------------------------------------------------------------------------- Returns ----------------------------------------------------------------------------------------------- - A list of numpy ndarray objects whose i'th element holds the latent space initial condition - for the i'th set of parameters in the param_grid. That is, if we let U0_i denote the fom IC for - the i'th set of parameters, then the i'th element of the returned list is Z0_i = encoder(U0_i). + A 3d numpy ndarray whose i, j, k element holds the k'th component of the j'th time step of + the solution to the latent dynamics when we use the latent encoding of the initial condition + from the i'th combination of parameter values """ + + # The param grid needs to be two dimensional, with the first axis corresponding to which + # instance of the parameter values we are using. If there is only one parameter, it may be 1d. + # We can fix that by adding on an axis with size 1. if (param_grid.ndim == 1): param_grid = param_grid.reshape(1, -1) - n_test = param_grid.shape[0] - - Z0 = initial_condition_latent(param_grid, physics, autoencoder) - - pred_mean, _ = eval_gp(gp_dictionary, param_grid) - Zis = np.zeros([n_test, physics.nt, autoencoder.n_z]) - for i in range(n_test): + # Now fetch the number of combinations of parameter values. + n_param : int = param_grid.shape[0] + + # For each parameter in param_grid, fetch the corresponding initial condition and then encode + # it. This gives us a list whose i'th element holds the encoding of the i'th initial condition. + Z0 : list[np.ndarray] = initial_condition_latent(param_grid, physics, autoencoder) + + # Evaluate each GP at each combination of parameter values. This returns two arrays, the + # first of which is a 2d array whose i,j element specifies the mean of the posterior + # distribution for the j'th coefficient at the i'th combination of parameter values. + pred_mean, _ = eval_gp(gp_list, param_grid) + + # For each testing parameter, cycle through the mean value of each coefficient from each + # posterior distribution. For each set of coefficients (combination of parameter values), solve + # the latent dynamics forward in time (starting from the corresponding IC value) and store the + # resulting solution frames in Zis, a 3d array whose i, j, k element holds the k'th component + # of the j'th time step fo the latent solution when we use the coefficients from the posterior + # distribution for the i'th combination of parameter values. + Zis : np.ndarray = np.zeros([n_param, physics.nt, autoencoder.n_z]) + for i in range(n_param): Zis[i] = latent_dynamics.simulate(pred_mean[i], Z0[i], physics.t_grid) + # All done! return Zis -def sample_roms(autoencoder : Autoencoder, physics, latent_dynamics, gp_dictionary, param_grid, n_samples): +def sample_roms(autoencoder : Autoencoder, + physics : Physics, + latent_dynamics : LatentDynamics, + gp_list : list[GaussianProcessRegressor], + param_grid : np.ndarray, + n_samples : int) -> np.ndarray: ''' - Collect n_samples of ROM trajectories on param_grid. - gp_dictionary: list of Gaussian process regressors (size of n_test) - param_grid: numpy 2d array - n_samples: integer - assert(len(gp_dictionnary) == param_grid.shape[0]) + This function samples the latent coefficients, solves the corresponding latent dynamics, and + then returns the resulting latent solutions. + + Specifically, for each combination of parameter values in the param_grid, we draw n_samples + samples of the latent coefficients (from the coefficient posterior distributions evaluated at + that parameter value). This gives us a set of n_samples latent dynamics coefficients. For each + set of coefficients, we solve the corresponding latent dynamics forward in time and store the + resulting solution frames. We do this for each sample and each combination of parameter values, + resulting in an (n_param, n_sample, n_t, n_z) array of solution frames, which is what we + return. - output: np.array of size [n_test, n_samples, physics.nt, autoencoder.n_z] - ''' + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- - if (param_grid.ndim == 1): - param_grid = param_grid.reshape(1, -1) - n_test = param_grid.shape[0] + autoencoder: An autoencoder. We use this to map the fom IC's (stored in Physics) to the + latent space using the autoencoder's encoder. + + physics: A "Physics" object that stores the ICs for each parameter combination. + + latent_dynamics: A LatentDynamics object which describes how we specify the dynamics in the + Autoencoder's latent space. We use this to simulate the latent dynamics forward in time. + + gp_list: a list of trained GP regressor objects. The number of elements in this list should + match the number of columns in param_grid. The i'th element of this list is a GP regressor + object that predicts the i'th coefficient. + + param_grid: A 2d numpy.ndarray object of shape (number of parameter combination) x (number of + parameters). The i,j element of this array holds the value of the j'th parameter in the i'th + combination of parameters. - Z0 = initial_condition_latent(param_grid, physics, autoencoder) + n_samples: The number of samples we want to draw from each posterior distribution for each + coefficient evaluated at each combination of parameter values. + - coef_samples = [sample_coefs(gp_dictionary, param_grid[i], n_samples) for i in range(n_test)] + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + A np.array of size [n_test, n_samples, physics.nt, autoencoder.n_z]. The i, j, k, l element + holds the l'th component of the k'th frame of the solution to the latent dynamics when we use + the j'th sample of latent coefficients drawn from the posterior distribution for the i'th + combination of parameter values (i'th row of param_grid). + ''' - Zis = np.zeros([n_test, n_samples, physics.nt, autoencoder.n_z]) + # The param grid needs to be two dimensional, with the first axis corresponding to which + # instance of the parameter values we are using. If there is only one parameter, it may be 1d. + # We can fix that by adding on an axis with size 1. + if (param_grid.ndim == 1): + param_grid = param_grid.reshape(1, -1) + + # Now fetch the number of combinations of parameter values (rows of param_grid). + n_param : int = param_grid.shape[0] + + # For each parameter in param_grid, fetch the corresponding initial condition and then encode + # it. This gives us a list whose i'th element holds the encoding of the i'th initial condition. + Z0 : list[np.ndarray] = initial_condition_latent(param_grid, physics, autoencoder) + + # Now, for each combination of parameters, draw n_samples samples from the posterior + # distributions for each coefficient at that combination of parameters. We store these samples + # in a list of numpy arrays. The k'th list element is a (n_sample, n_coef) array whose i, j + # element stores the i'th sample from the posterior distribution for the j'th coefficient at + # the k'th combination of parameter values. + coef_samples : list[np.ndarray] = [sample_coefs(gp_list, param_grid[i], n_samples) for i in range(n_param)] + + # For each testing parameter, cycle through the samples of the coefficients for that + # combination of parameter values. For each set of coefficients, solve the corresponding latent + # dynamics forward in time and store the resulting frames in Zis. This is a 4d array whose i, + # j, k, l element holds the l'th component of the k'th frame of the solution to the latent + # dynamics when we use the j'th sample of latent coefficients drawn from the posterior + # distribution for the i'th combination of parameter values. + Zis = np.zeros([n_param, n_samples, physics.nt, autoencoder.n_z]) for i, Zi in enumerate(Zis): z_ic = Z0[i] for j, coef_sample in enumerate(coef_samples[i]): Zi[j] = latent_dynamics.simulate(coef_sample, z_ic, physics.t_grid) + # All done! return Zis def get_fom_max_std(autoencoder : Autoencoder, Zis : np.ndarray) -> int: """ - Computes the maximum standard deviation across the parameter space grid and finds the - corresponding parameter index. - + Computes the maximum standard deviation across the trajectories in Zis and returns the + corresponding parameter index. Specifically, Zis is a 4d tensor of shape (n_test, n_samples, + n_t, n_z). The first axis specifies which parameter combination we're using. For each + combination of parameters, we assume that we drew n_samples of the posterior distribution of + the coefficients at that parameter value, simulated the corresponding dynamics for n_t time + steps, and then recorded the results in Zis[i]. Thus, Zis[i, j, k, :] represents the k'th + time step of the solution to the latent dynamics when we use the coefficients from the j'th + sample of the posterior distribution for the i'th set of parameters. + + Let i \in {1, 2, ... , n_test} and k \in {1, 2, ... , n_t}. For each j, we map the k'th frame + of the j'th solution trajectory for the i'th parameter combination (Zi[i, j, k, :]) to a fom + frame. We do this for each j (the set of samples), which gives us a collection of n_sample + fom frames, representing samples of the distribution of fom frames at the k'th time step + when we use the posterior distribution for the i'th set of parameters. For each l \in {1, 2, + ... , n_fom}, we compute the STD of the set of l'th components of these n_sample fom frames. + We do this for each i and k and then figure out which i, k, l combination gives the largest + STD. We return the corresponding i index. ----------------------------------------------------------------------------------------------- @@ -146,7 +257,7 @@ def get_fom_max_std(autoencoder : Autoencoder, Zis : np.ndarray) -> int: # (n_t, n_fom) whose i,j element holds the (sample) standard deviation of the j'th component # of the i'th frame of the fom solution. In this case, the sample distribution consists of # the set of j'th components of i'th frames of fom solutions (one for each sample of the - # posterior coefficients). + # coefficient posterior distributions). X_pred_m_std : np.ndarray = X_pred_m.std(0) # Now compute the maximum standard deviation across frames/fom components. @@ -161,6 +272,11 @@ def get_fom_max_std(autoencoder : Autoencoder, Zis : np.ndarray) -> int: return m_index + +# ------------------------------------------------------------------------------------------------- +# BayesianGLaSDI class +# ------------------------------------------------------------------------------------------------- + # move optimizer parameters to device def optimizer_to(optim, device): for param in optim.state.values(): @@ -369,201 +485,4 @@ def load(self, dict_): if (self.device != 'cpu'): optimizer_to(self.optimizer, self.device) return - - - -class BayesianGLaSDI: - X_train = None - - def __init__(self, physics, autoencoder, latent_dynamics, model_parameters): - - ''' - - This class runs a full GPLaSDI training. It takes into input the autoencoder defined as a PyTorch object and the - dictionnary containing all the training parameters. - The "train" method with run the active learning training loop, compute the reconstruction and SINDy loss, train the GPs, - and sample a new FOM data point. - - ''' - - self.autoencoder = autoencoder - self.latent_dynamics = latent_dynamics - self.physics = physics - self.param_space = physics.param_space - self.timer = Timer() - - self.n_samples = model_parameters['n_samples'] - self.lr = model_parameters['lr'] - self.n_iter = model_parameters['n_iter'] - self.n_greedy = model_parameters['n_greedy'] - self.max_greedy_iter = model_parameters['max_greedy_iter'] - self.sindy_weight = model_parameters['sindy_weight'] - self.coef_weight = model_parameters['coef_weight'] - - self.optimizer = torch.optim.Adam(autoencoder.parameters(), lr = self.lr) - self.MSE = torch.nn.MSELoss() - - self.path_checkpoint = model_parameters['path_checkpoint'] - self.path_results = model_parameters['path_results'] - - from os.path import dirname - from pathlib import Path - Path(dirname(self.path_checkpoint)).mkdir(parents=True, exist_ok=True) - Path(dirname(self.path_results)).mkdir(parents=True, exist_ok=True) - - device = model_parameters['device'] if 'device' in model_parameters else 'cpu' - if (device == 'cuda'): - assert(torch.cuda.is_available()) - self.device = device - elif (device == 'mps'): - assert(torch.backends.mps.is_available()) - self.device = device - else: - self.device = 'cpu' - - self.best_loss = np.Inf - self.restart_iter = 0 - - return - - def train(self): - assert(self.X_train is not None) - assert(self.X_train.size(0) == self.param_space.n_train) - - device = self.device - autoencoder_device = self.autoencoder.to(device) - X_train_device = self.X_train.to(device) - - from pathlib import Path - Path(self.path_checkpoint).mkdir(parents=True, exist_ok=True) - Path(self.path_results).mkdir(parents=True, exist_ok=True) - - ps = self.param_space - ld = self.latent_dynamics - - for iter in range(self.restart_iter, self.n_iter): - self.timer.start("train_step") - - self.optimizer.zero_grad() - Z = autoencoder_device.encoder(X_train_device) - X_pred = autoencoder_device.decoder(Z) - Z = Z.cpu() - - loss_ae = self.MSE(X_train_device, X_pred) - coefs, loss_sindy, loss_coef = ld.calibrate(Z, self.physics.dt, compute_loss=True, numpy=True) - - max_coef = np.abs(coefs).max() - - loss = loss_ae + self.sindy_weight * loss_sindy / ps.n_train + self.coef_weight * loss_coef / ps.n_train - - loss.backward() - self.optimizer.step() - - if loss.item() < self.best_loss: - torch.save(autoencoder_device.state_dict(), self.path_checkpoint + '/' + 'checkpoint.pt') - self.best_coefs = coefs - self.best_loss = loss.item() - - print("Iter: %05d/%d, Loss: %3.10f, Loss AE: %3.10f, Loss SI: %3.10f, Loss COEF: %3.10f, max|c|: %04.1f, " - % (iter + 1, self.n_iter, loss.item(), loss_ae.item(), loss_sindy.item(), loss_coef.item(), max_coef), - end = '') - - if ps.n_train < 6: - print('Param: ' + str(np.round(ps.train_space[0, :], 4)), end = '') - - for i in range(1, ps.n_train - 1): - print(', ' + str(np.round(ps.train_space[i, :], 4)), end = '') - print(', ' + str(np.round(ps.train_space[-1, :], 4))) - - else: - print('Param: ...', end = '') - for i in range(5): - print(', ' + str(np.round(ps.train_space[-6 + i, :], 4)), end = '') - print(', ' + str(np.round(ps.train_space[-1, :], 4))) - - self.timer.end("train_step") - - if ((iter > self.restart_iter) and (iter < self.max_greedy_iter) and (iter % self.n_greedy == 0)): - - if (self.best_coefs.shape[0] == ps.n_train): - coefs = self.best_coefs - - new_sample = self.get_new_sample_point(coefs) - - # TODO(kevin): implement save/load the new parameter - ps.appendTrainSpace(new_sample) - self.restart_iter = iter - next_step, result = NextStep.RunSample, Result.Success - print('New param: ' + str(np.round(new_sample, 4)) + '\n') - # self.timer.end("new_sample") - return result, next_step - - self.timer.start("finalize") - - if (self.best_coefs.shape[0] == ps.n_train): - coefs = self.best_coefs - - gp_dictionnary = fit_gps(ps.train_space, coefs) - - bglasdi_results = {'autoencoder_param': self.autoencoder.state_dict(), 'final_X_train': self.X_train, - 'coefs': coefs, 'gp_dictionnary': gp_dictionnary, 'lr': self.lr, 'n_iter': self.n_iter, - 'n_greedy': self.n_greedy, 'sindy_weight': self.sindy_weight, 'coef_weight': self.coef_weight, - 'n_samples' : self.n_samples, - } - bglasdi_results['physics'] = self.physics.export() - bglasdi_results['parameters'] = self.param_space.export() - # TODO(kevin): restart capability for timer. - bglasdi_results['timer'] = self.timer.export() - bglasdi_results['latent_dynamics'] = self.latent_dynamics.export() - - date = time.localtime() - date_str = "{month:02d}_{day:02d}_{year:04d}_{hour:02d}_{minute:02d}" - date_str = date_str.format(month = date.tm_mon, day = date.tm_mday, year = date.tm_year, hour = date.tm_hour + 3, minute = date.tm_min) - np.save(self.path_results + '/' + 'bglasdi_' + date_str + '.npy', bglasdi_results) - - self.timer.end("finalize") - self.timer.print() - - next_step, result = None, Result.Complete - return result, next_step - - def get_new_sample_point(self, coefs): - self.timer.start("new_sample") - - print('\n~~~~~~~ Finding New Point ~~~~~~~') - # TODO(kevin): william, this might be the place for new sampling routine. - - ae = self.autoencoder.cpu() - ps = self.param_space - ae.load_state_dict(torch.load(self.path_checkpoint + '/' + 'checkpoint.pt')) - - Z0 = initial_condition_latent(ps.test_space, self.physics, ae) - - gp_dictionnary = fit_gps(ps.train_space, coefs) - - coef_samples = [sample_coefs(gp_dictionnary, ps.test_space[i], self.n_samples) for i in range(ps.n_test)] - - Zis = np.zeros([ps.n_test, self.n_samples, self.physics.nt, ae.n_z]) - for i, Zi in enumerate(Zis): - z_ic = Z0[i] - for j, coef_sample in enumerate(coef_samples[i]): - Zi[j] = self.latent_dynamics.simulate(coef_sample, z_ic, self.physics.t_grid) - - m_index = get_fom_max_std(ae, Zis) - - self.timer.end("new_sample") - return ps.test_space[m_index, :] - - def sample_fom(self): - - new_foms = self.param_space.n_train - self.X_train.size(0) - assert(new_foms > 0) - new_params = self.param_space.train_space[-new_foms:, :] - - if not self.physics.offline: - new_X = self.physics.generate_solutions(new_params) - - self.X_train = torch.cat([self.X_train, new_X], dim = 0) - else: - # TODO(kevin): interface for offline FOM simulation - raise RuntimeError("Offline FOM simulation is not supported yet!") \ No newline at end of file + \ No newline at end of file diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index d7f547c..b84b527 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -81,12 +81,12 @@ def initial_condition_latent(param_grid : np.ndarray, # Fetch the IC for the i'th set of parameters. Then map it to a tensor. u0 : np.ndarray = physics.initial_condition(param_grid[i]) - u0 = u0.reshape(sol_shape) + u0 = u0.reshape(sol_shape) u0 = torch.Tensor(u0) # Encode the IC, then map the encoding to a numpy array. - z0 = autoencoder.encoder(u0) - z0 = z0[0, 0, :].detach().numpy() + z0 : np.ndarray = autoencoder.encoder(u0) + z0 = z0[0, 0, :].detach().numpy() # Append the new IC to the list of latent ICs Z0.append(z0) From ad259136ea6b17e72a3a4d5a62f945af31143c26 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Fri, 25 Oct 2024 10:59:10 -0700 Subject: [PATCH 11/22] Made some minor commenting tweaks to gplasdi.py --- src/lasdi/gplasdi.py | 40 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/src/lasdi/gplasdi.py b/src/lasdi/gplasdi.py index b51f7b4..629fcd7 100644 --- a/src/lasdi/gplasdi.py +++ b/src/lasdi/gplasdi.py @@ -6,6 +6,7 @@ import torch import numpy as np +from torch.optim import Optimizer from sklearn.gaussian_process import GaussianProcessRegressor from .gp import eval_gp, sample_coefs, fit_gps @@ -278,7 +279,28 @@ def get_fom_max_std(autoencoder : Autoencoder, Zis : np.ndarray) -> int: # ------------------------------------------------------------------------------------------------- # move optimizer parameters to device -def optimizer_to(optim, device): +def optimizer_to(optim : Optimizer, device : str) -> None: + """ + This function moves an optimizer object to a specific device. + + + ----------------------------------------------------------------------------------------------- + Arguments + ----------------------------------------------------------------------------------------------- + + optim: The optimizer whose device we want to change. + + device: The device we want to move optim onto. + + + ----------------------------------------------------------------------------------------------- + Returns + ----------------------------------------------------------------------------------------------- + + Nothing. + """ + + # Cycle through the optimizer's parameters. for param in optim.state.values(): # Not sure there are any global tensors in the state dict if isinstance(param, torch.Tensor): @@ -292,6 +314,8 @@ def optimizer_to(optim, device): if subparam._grad is not None: subparam._grad.data = subparam._grad.data.to(device) + + class BayesianGLaSDI: X_train = torch.Tensor([]) X_test = torch.Tensor([]) @@ -351,6 +375,8 @@ def __init__(self, physics, autoencoder : Autoencoder, latent_dynamics, param_sp return + + def train(self): assert(self.X_train.size(0) > 0) assert(self.X_train.size(0) == self.param_space.n_train()) @@ -431,7 +457,9 @@ def train(self): self.timer.print() return - + + + def get_new_sample_point(self): self.timer.start("new_sample") assert(self.X_test.size(0) > 0) @@ -466,7 +494,9 @@ def get_new_sample_point(self): self.timer.end("new_sample") return new_sample - + + + def export(self): dict_ = {'X_train': self.X_train, 'X_test': self.X_test, 'lr': self.lr, 'n_iter': self.n_iter, 'n_samples' : self.n_samples, 'best_coefs': self.best_coefs, 'max_iter': self.max_iter, @@ -474,7 +504,9 @@ def export(self): 'restart_iter': self.restart_iter, 'timer': self.timer.export(), 'optimizer': self.optimizer.state_dict() } return dict_ - + + + def load(self, dict_): self.X_train = dict_['X_train'] self.X_test = dict_['X_test'] From ce8e6c8c00890a6bc123e63a83a496fefce21e9a Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Fri, 25 Oct 2024 16:21:28 -0700 Subject: [PATCH 12/22] Finished documenting gplasdi.py I think I will move onto the LatentDynamics/Physics files before moving onto workflow.py --- src/lasdi/gplasdi.py | 390 ++++++++++++++++++++++++++++---------- src/lasdi/latent_space.py | 35 +--- 2 files changed, 303 insertions(+), 122 deletions(-) diff --git a/src/lasdi/gplasdi.py b/src/lasdi/gplasdi.py index 629fcd7..02cbe1d 100644 --- a/src/lasdi/gplasdi.py +++ b/src/lasdi/gplasdi.py @@ -317,45 +317,86 @@ def optimizer_to(optim : Optimizer, device : str) -> None: class BayesianGLaSDI: - X_train = torch.Tensor([]) - X_test = torch.Tensor([]) - - def __init__(self, physics, autoencoder : Autoencoder, latent_dynamics, param_space, config): - - ''' - - This class runs a full GPLaSDI training. It takes into input the autoencoder defined as a PyTorch object and the - dictionnary containing all the training parameters. - The "train" method with run the active learning training loop, compute the reconstruction and SINDy loss, train the GPs, - and sample a new FOM data point. - - ''' - - self.autoencoder = autoencoder - self.latent_dynamics = latent_dynamics - self.physics = physics - self.param_space = param_space - self.timer = Timer() - - self.n_samples = config['n_samples'] - self.lr = config['lr'] - self.n_iter = config['n_iter'] # number of iterations for one train and greedy sampling - self.max_iter = config['max_iter'] # maximum iterations for overall training - self.max_greedy_iter = config['max_greedy_iter'] # maximum iterations for greedy sampling - self.ld_weight = config['ld_weight'] - self.coef_weight = config['coef_weight'] - - self.optimizer = torch.optim.Adam(autoencoder.parameters(), lr = self.lr) - self.MSE = torch.nn.MSELoss() - - self.path_checkpoint = config['path_checkpoint'] - self.path_results = config['path_results'] + X_train : torch.Tensor = torch.Tensor([]) + X_test : torch.Tensor = torch.Tensor([]) + + def __init__(self, + physics : Physics, + autoencoder : Autoencoder, + latent_dynamics : LatentDynamics, + param_space : ParameterSpace, + config : dict): + """ + This class runs a full GPLaSDI training. As input, it takes the autoencoder defined as a + torch.nn.Module object, a Physics object to recover fom ICs + information on the time + discretization, a + + The "train" method runs the active learning training loop, computes the reconstruction and + SINDy loss, trains the GPs, and samples a new FOM data point. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + physics: A "Physics" object that we use to fetch the fom initial conditions (which we + encode into latent ICs). Each Physics object has + a corresponding PDE with parameters, and a way to generate a solution to that equation + given a particular set of parameter values (and an IC, BCs). We use this object to generate + fom solutions which we then use to train the autoencoder/latent dynamics. + + Autoencoder: An autoencoder object that we use to compress the fom state to a reduced, + latent state. + + latent_dynamics: A LatentDynamics object which describes how we specify the dynamics in the + Autoencoder's latent space. + + param_space: A Parameter space object which holds the set of testing and training + parameters. + + config: A dictionary housing the settings we wna to use to train the model on + the data generated by physics. + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + self.autoencoder = autoencoder + self.latent_dynamics = latent_dynamics + self.physics = physics + self.param_space = param_space + + # Initialize a timer object. We will use this while training. + self.timer = Timer() + + # Extract training/loss hyperparameters from the configuration file. + self.n_samples : int = config['n_samples'] # Number of samples to draw per coefficient per combination of parameters + self.lr : float = config['lr'] # Learning rate for the optimizer. + self.n_iter : int = config['n_iter'] # Number of iterations for one train and greedy sampling + self.max_iter : int = config['max_iter'] # Maximum iterations for overall training + self.max_greedy_iter : int = config['max_greedy_iter'] # Maximum iterations for greedy sampling + self.ld_weight : float = config['ld_weight'] # Weight of the SINDy loss in the loss function. \beta_2 in the paper. + self.coef_weight : float = config['coef_weight'] # Weight of the norm of matrix of latent dynamics coefficients. \beta_3 in the paper. + + # Set up the optimizer and loss function. + self.optimizer : Optimizer = torch.optim.Adam(autoencoder.parameters(), lr = self.lr) + self.MSE = torch.nn.MSELoss() + + # Set paths for checkpointing. + self.path_checkpoint : str = config['path_checkpoint'] + self.path_results : str = config['path_results'] + + # Make sure the checkpoints and results directories exist. from os.path import dirname from pathlib import Path - Path(dirname(self.path_checkpoint)).mkdir(parents=True, exist_ok=True) - Path(dirname(self.path_results)).mkdir(parents=True, exist_ok=True) + Path(dirname(self.path_checkpoint)).mkdir( parents = True, exist_ok = True) + Path(dirname(self.path_results)).mkdir( parents = True, exist_ok = True) + # Set the device to train on. We default to cpu. device = config['device'] if 'device' in config else 'cpu' if (device == 'cuda'): assert(torch.cuda.is_available()) @@ -366,67 +407,121 @@ def __init__(self, physics, autoencoder : Autoencoder, latent_dynamics, param_sp else: self.device = 'cpu' - self.best_loss = np.inf - self.best_coefs = None - self.restart_iter = 0 + # Set up variables to aide checkpointing. + self.best_loss : float = np.inf # The lowest testing loss we have found so far + self.best_coefs : np.ndarray = None # The best coefficients from the iteration with lowest testing loss + self.restart_iter : int = 0 # Iteration number at the end of the last training period + + # Set placeholder tensors to hold the testing and training data. We expect to set up + # X_train to be a tensor of shape (Np, Nt, Nx[0], ... , Nx[Nd - 1]), where Np is the number + # of parameter combinations in the training set, Nt is the number of time steps per fom + # solution, and Nx[0], ... , Nx[Nd - 1] represent the number of steps along the spatial + # axes. X_test has an analogous shape, but it's leading dimension has a size matching the + # number of combinations of parameters in the testing set. + self.X_train : torch.Tensor = torch.Tensor([]) + self.X_test : torch.Tensor = torch.Tensor([]) + + # All done! + return - self.X_train = torch.Tensor([]) - self.X_test = torch.Tensor([]) - return + def train(self) -> None: + """ + Runs a round of training on the autoencoder. + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- - def train(self): + None! + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + # Make sure we have at least one training data point (the 0 axis of X_Train corresponds + # which combination of training parameters we use). assert(self.X_train.size(0) > 0) assert(self.X_train.size(0) == self.param_space.n_train()) - device = self.device - autoencoder_device = self.autoencoder.to(device) - X_train_device = self.X_train.to(device) + # Map everything to self's device. + device : str = self.device + autoencoder_device : Autoencoder = self.autoencoder.to(device) + X_train_device : torch.Tensor = self.X_train.to(device) + # Make sure the checkpoints and results directories exist. from pathlib import Path - Path(self.path_checkpoint).mkdir(parents=True, exist_ok=True) - Path(self.path_results).mkdir(parents=True, exist_ok=True) - - ps = self.param_space - n_train = ps.n_train() - ld = self.latent_dynamics + Path(self.path_checkpoint).mkdir( parents = True, exist_ok = True) + Path(self.path_results).mkdir( parents = True, exist_ok = True) - ''' - determine number of iterations. - Perform n_iter iterations until overall iterations hit max_iter. - ''' - next_iter = min(self.restart_iter + self.n_iter, self.max_iter) + ps : ParameterSpace = self.param_space + n_train : int = ps.n_train() + ld : LatentDynamics = self.latent_dynamics + # Determine number of iterations we should run in this round of training. + next_iter : int = min(self.restart_iter + self.n_iter, self.max_iter) + + # Run the iterations! for iter in range(self.restart_iter, next_iter): + # Begin timing the current training step. self.timer.start("train_step") + # Zero out the gradients. self.optimizer.zero_grad() - Z = autoencoder_device.encoder(X_train_device) - X_pred = autoencoder_device.decoder(Z) - Z = Z.cpu() - loss_ae = self.MSE(X_train_device, X_pred) - coefs, loss_ld, loss_coef = ld.calibrate(Z, self.physics.dt, compute_loss=True, numpy=True) - max_coef = np.abs(coefs).max() + # ------------------------------------------------------------------------------------- + # Forward pass + # Run the forward pass + Z : torch.Tensor = autoencoder_device.encoder(X_train_device) + X_pred : torch.Tensor = autoencoder_device.decoder(Z) + Z : torch.Tensor = Z.cpu() + + # Compute the autoencoder loss. + loss_ae : torch.Tensor = self.MSE(X_train_device, X_pred) + + # Compute the latent dynamics and coefficient losses. Also fetch the current latent + # dynamics coefficients for each training point. The latter is stored in a 3d array + # called "coefs" of shape (n_train, N_z, N_l), where N_{\mu} = n_train = number of + # training parameter combinations, N_z = latent space dimension, and N_l = number of + # terms in the SINDy library. + coefs, loss_ld, loss_coef = ld.calibrate(Z, self.physics.dt, compute_loss = True, numpy = True) + max_coef : np.float32 = np.abs(coefs).max() + + # Compute the final loss. loss = loss_ae + self.ld_weight * loss_ld / n_train + self.coef_weight * loss_coef / n_train + + # ------------------------------------------------------------------------------------- + # Backward Pass + + # Run back propagation and update the model parameters. loss.backward() self.optimizer.step() + # Check if we hit a new minimum loss. If so, make a checkpoint, record the loss and + # the iteration number. if loss.item() < self.best_loss: torch.save(autoencoder_device.cpu().state_dict(), self.path_checkpoint + '/' + 'checkpoint.pt') - autoencoder_device = self.autoencoder.to(device) - self.best_coefs = coefs - self.best_loss = loss.item() + autoencoder_device : Autoencoder = self.autoencoder.to(device) + self.best_coefs : np.ndarray = coefs + self.best_loss : float = loss.item() + + # ------------------------------------------------------------------------------------- + # Report Results from this iteration + # Report the current iteration number and losses print("Iter: %05d/%d, Loss: %3.10f, Loss AE: %3.10f, Loss LD: %3.10f, Loss COEF: %3.10f, max|c|: %04.1f, " % (iter + 1, self.max_iter, loss.item(), loss_ae.item(), loss_ld.item(), loss_coef.item(), max_coef), end = '') + # If there are fewer than 6 training examples, report the set of parameter combinations. if n_train < 6: print('Param: ' + str(np.round(ps.train_space[0, :], 4)), end = '') @@ -434,87 +529,192 @@ def train(self): print(', ' + str(np.round(ps.train_space[i, :], 4)), end = '') print(', ' + str(np.round(ps.train_space[-1, :], 4))) + # Otherwise, report the final 6 parameter combinations. else: print('Param: ...', end = '') for i in range(5): print(', ' + str(np.round(ps.train_space[-6 + i, :], 4)), end = '') print(', ' + str(np.round(ps.train_space[-1, :], 4))) + # We have finished a training step, stop the timer. self.timer.end("train_step") + # We are ready to wrap up the training procedure. self.timer.start("finalize") + # Now that we have completed another round of training, update the restart iteration. self.restart_iter += self.n_iter + # Recover the autoencoder + coefficients which attained the lowest loss. If we recorded + # our bess loss in this round of training, then we replace the autoencoder's parameters + # with those from the iteration that got the best loss. Otherwise, we use the current + # set of coefficients and serialize the current model. if ((self.best_coefs is not None) and (self.best_coefs.shape[0] == n_train)): - state_dict = torch.load(self.path_checkpoint + '/' + 'checkpoint.pt') + state_dict = torch.load(self.path_checkpoint + '/' + 'checkpoint.pt') self.autoencoder.load_state_dict(state_dict) else: - self.best_coefs = coefs + self.best_coefs : np.ndarray = coefs torch.save(autoencoder_device.cpu().state_dict(), self.path_checkpoint + '/' + 'checkpoint.pt') + # Report timing information. self.timer.end("finalize") self.timer.print() + # All done! return - def get_new_sample_point(self): + def get_new_sample_point(self) -> np.ndarray: + """ + This function uses a greedy process to sample a new parameter value. Specifically, it runs + through each combination of parameters in in self.param_space. For the i'th combination of + parameters, we generate a collection of samples of the coefficients in the latent dynamics. + We draw the k'th sample of the j'th coefficient from the posterior distribution for the + j'th coefficient at the i'th combination of parameters. We map the resulting solution back + into the real space and evaluate the standard deviation of the fom frames. We return the + combination of parameters which engenders the largest standard deviation (see the function + get_fom_max_std). + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + None! + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + a 2d numpy ndarray object of shape (1, n_param) whose (0, j) element holds the value of + the j'th parameter in the new sample. + """ + + self.timer.start("new_sample") - assert(self.X_test.size(0) > 0) - assert(self.X_test.size(0) == self.param_space.n_test()) + assert(self.X_test.size(0) > 0) + assert(self.X_test.size(0) == self.param_space.n_test()) assert(self.best_coefs.shape[0] == self.param_space.n_train()) - coefs = self.best_coefs + coefs : np.ndarray = self.best_coefs print('\n~~~~~~~ Finding New Point ~~~~~~~') # TODO(kevin): william, this might be the place for new sampling routine. - ae = self.autoencoder.cpu() - ps = self.param_space - n_test = ps.n_test() + # Move the model to the cpu (this is where all the GP stuff happens) and load the model + # from the last checkpoint. This should be the one that obtained the best loss so far. + # Remember that coefs should specify the coefficients from that iteration. + ae : Autoencoder = self.autoencoder.cpu() + ps : ParameterSpace = self.param_space + n_test : int = ps.n_test() ae.load_state_dict(torch.load(self.path_checkpoint + '/' + 'checkpoint.pt')) - Z0 = initial_condition_latent(ps.test_space, self.physics, ae) - - gp_dictionnary = fit_gps(ps.train_space, coefs) - - coef_samples = [sample_coefs(gp_dictionnary, ps.test_space[i], self.n_samples) for i in range(n_test)] - - Zis = np.zeros([n_test, self.n_samples, self.physics.nt, ae.n_z]) + # Map the initial conditions for the fom to initial conditions in the latent space. + Z0 : list[np.ndarray] = initial_condition_latent(ps.test_space, self.physics, ae) + + # Train the GPs on the training data, get one GP per latent space coefficient. + gp_list : list[GaussianProcessRegressor] = fit_gps(ps.train_space, coefs) + + # For each combination of parameter values in the testing set, for each coefficient, + # draw a set of samples from the posterior distribution for that coefficient evaluated at + # the testing parameters. We store the samples for a particular combination of parameter + # values in a 2d numpy.ndarray of shape (n_sample, n_coef), whose i, j element holds the + # i'th sample of the j'th coefficient. We store the arrays for different parameter values + # in a list of length (number of combinations of parameters in the testing set). + coef_samples : list[np.ndarray] = [sample_coefs(gp_list, ps.test_space[i], self.n_samples) for i in range(n_test)] + + # Now, solve the latent dynamics forward in time for each set of coefficients in + # coef_samples. There are n_test combinations of parameter values, and we have n_samples + # sets of coefficients for each combination of parameter values. For each of those, we want + # to solve the corresponding latent dynamics for n_t time steps. Each one of the frames + # in that solution live in \mathbb{R}^{n_z}. Thus, we need to store the results in a 4d + # array of shape (n_test, n_samples, n_t, n_z) whose i, j, k, l element holds the l'th + # component of the k'th frame of the solution to the latent dynamics when we use the + # j'th sample of the coefficients for the i'th testing parameter value and when the latent + # dynamics uses the encoding of the i'th fom IC as its IC. + Zis : np.ndarray = np.zeros([n_test, self.n_samples, self.physics.nt, ae.n_z]) for i, Zi in enumerate(Zis): z_ic = Z0[i] for j, coef_sample in enumerate(coef_samples[i]): Zi[j] = self.latent_dynamics.simulate(coef_sample, z_ic, self.physics.t_grid) - m_index = get_fom_max_std(ae, Zis) + # Find the index of the parameter with the largest std. + m_index : int = get_fom_max_std(ae, Zis) - new_sample = ps.test_space[m_index, :].reshape(1, -1) + # We have found the testing parameter we want to add to the training set. Fetch it, then + # stop the timer and return the parameter. + new_sample : np.ndarray = ps.test_space[m_index, :].reshape(1, -1) print('New param: ' + str(np.round(new_sample, 4)) + '\n') - self.timer.end("new_sample") + + # All done! return new_sample - def export(self): - dict_ = {'X_train': self.X_train, 'X_test': self.X_test, 'lr': self.lr, 'n_iter': self.n_iter, - 'n_samples' : self.n_samples, 'best_coefs': self.best_coefs, 'max_iter': self.max_iter, - 'max_iter': self.max_iter, 'ld_weight': self.ld_weight, 'coef_weight': self.coef_weight, - 'restart_iter': self.restart_iter, 'timer': self.timer.export(), 'optimizer': self.optimizer.state_dict() - } + def export(self) -> dict: + """ + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A dictionary housing most of the internal variables in self. You can pass this dictionary + to self (after initializing it using ParameterSpace, Autoencoder, and LatentDynamics + objects) to make a GLaSDI object whose internal state matches that of self. + """ + + dict_ = {'X_train' : self.X_train, + 'X_test' : self.X_test, + 'lr' : self.lr, + 'n_iter' : self.n_iter, + 'n_samples' : self.n_samples, + 'best_coefs' : self.best_coefs, + 'max_iter' : self.max_iter, + 'max_iter' : self.max_iter, + 'ld_weight' : self.ld_weight, + 'coef_weight' : self.coef_weight, + 'restart_iter' : self.restart_iter, + 'timer' : self.timer.export(), + 'optimizer' : self.optimizer.state_dict()} return dict_ - def load(self, dict_): - self.X_train = dict_['X_train'] - self.X_test = dict_['X_test'] - self.best_coefs = dict_['best_coefs'] - self.restart_iter = dict_['restart_iter'] + def load(self, dict_ : dict) -> None: + """ + Modifies self's internal state to match the one whose export method generated the dict_ + dictionary. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + dict_: This should be a dictionary returned by calling the export method on another + GLaSDI object. We use this to make self hav the same internal state as the object that + generated dict_. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + # Extract instance variables from dict_. + self.X_train : torch.Tensor = dict_['X_train'] + self.X_test : torch.Tensor = dict_['X_test'] + self.best_coefs : np.ndarray = dict_['best_coefs'] + self.restart_iter : int = dict_['restart_iter'] + + # Load the timer / optimizer. self.timer.load(dict_['timer']) self.optimizer.load_state_dict(dict_['optimizer']) if (self.device != 'cpu'): optimizer_to(self.optimizer, self.device) + + # All done! return \ No newline at end of file diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index b84b527..7f64f65 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -56,7 +56,8 @@ def initial_condition_latent(param_grid : np.ndarray, parameters). The i,j element of this array holds the value of the j'th parameter in the i'th combination of parameters. - physics: A "Physics" object that stores the datasets for each parameter combination. + physics: A "Physics" object that, among other things, stores the IC for each combination of + parameters. autoencoder: The actual autoencoder object that we use to map the ICs into the latent space. @@ -274,21 +275,8 @@ def forward(self, x : torch.Tensor) -> torch.Tensor: def init_weight(self) -> None: """ - This function initializes the weight matrices and bias vectors in self's layers. - - - ------------------------------------------------------------------------------------------- - Arguments - ------------------------------------------------------------------------------------------- - - None! - - - ------------------------------------------------------------------------------------------- - Returns - ------------------------------------------------------------------------------------------- - - Nothing! + This function initializes the weight matrices and bias vectors in self's layers. It takes + no arguments and returns nothing! """ # TODO(kevin): support other initializations? @@ -432,21 +420,14 @@ def forward(self, x : torch.Tensor) -> torch.Tensor: def export(self) -> dict: """ - This function extracts self's parameters and returns them in a dictionary. - - - ------------------------------------------------------------------------------------------- - Arguments - ------------------------------------------------------------------------------------------- - - None! - - ------------------------------------------------------------------------------------------- Returns ------------------------------------------------------------------------------------------- - The A dictionary housing self's state dictionary. + This function extracts self's parameters and returns them in a dictionary. You can pass + the dictionary returned by this function to the load method of another Autoencoder object + (that you initialized to have the same architecture as self) to make the other autoencoder + identical to self. """ # TO DO: deep export which includes all information needed to re-initialize self from From 9f8ba91ac91c7691be131fc07b7bef662f59b57e Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 30 Oct 2024 09:21:31 -0700 Subject: [PATCH 13/22] Made commenting corrections to input.py I forgot to add the "name" argument to the arguments list in the doc-string for the "InputParser" class. I added a breif description of what "name" does. --- src/lasdi/inputs.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/lasdi/inputs.py b/src/lasdi/inputs.py index fa030bb..8f24ac9 100644 --- a/src/lasdi/inputs.py +++ b/src/lasdi/inputs.py @@ -25,25 +25,29 @@ class InputParser: - def __init__(self, dict : dict, name : str = "") -> None: + def __init__(self, dict_ : dict, name : str = "") -> None: """" Initializes an InputParser object by setting the underlying dictionary of settings as an attribute. - ------------------------------------------------------------------------------------------- + + -------------------gi------------------------------------------------------------------------ Arguments ------------------------------------------------------------------------------------------- - dict: The dictionary of settings. To avoid the risk of the user accidentally changing one + dict_: The dictionary of settings. To avoid the risk of the user accidentally changing one of the settings after wrapping it, we store a deep copy of dict in self. + + name: A string that we use toe name the InputParser object. We use this name when + reporting errors (it is purely for debugging purposes). """ # A shallow copy could cause issues if the user changes dict's keys/values after # initializing this object. We store a deep copy to avoid this risk. from copy import deepcopy - self.dict_ = deepcopy(dict) + self.dict_ : dict = deepcopy(dict_) - self.name = name + self.name : str = name From 69096e0aed0aeab9d75b0dbab79c660c7d5237a4 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 30 Oct 2024 12:51:48 -0700 Subject: [PATCH 14/22] Added documentation to fd.py This one was... hard. Wow this code was intracate and difficult to parse. I think I did it, however. I think I documented the Stencil class. I did not document specific sub-classes.. I think the docstring in Stencil should make it clear what the subclasses do. I probably need to give these comments another look... but this took me like 5 hours. I am done for now. --- src/lasdi/fd.py | 253 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 218 insertions(+), 35 deletions(-) diff --git a/src/lasdi/fd.py b/src/lasdi/fd.py index 3a00684..b767edf 100644 --- a/src/lasdi/fd.py +++ b/src/lasdi/fd.py @@ -1,61 +1,236 @@ -import numpy as np -import scipy.sparse as sps -import torch +# ------------------------------------------------------------------------------------------------- +# Imports and Setup +# ------------------------------------------------------------------------------------------------- + +import numpy as np +import scipy.sparse as sps +from scipy.sparse import spmatrix; +import torch + + + +# ------------------------------------------------------------------------------------------------- +# Stencil class +# ------------------------------------------------------------------------------------------------- class Stencil: - leftBdrDepth = 0 - leftBdrWidth = [] - leftBdrStencils = [[]] - leftBdrNorm = [] - - interiorStencils = np.array([]) - interiorIndexes = [] - - def getOperators(self, Nx, periodic=False): - norm = np.ones(Nx,) - periodicOffset = np.zeros(Nx,) - Dxi = sps.diags(self.interiorStencils, - self.interiorIndexes, - shape=(Nx, Nx), format='lil') + # Class variables + leftBdrDepth : int = 0 # Number of time steps (at the start/stop of the time series) we use a special stencil on. + leftBdrWidth : list[int] = [] # i'th element specifies the number of terms in the stencil for the i'th time step + leftBdrStencils : list[list[float]] = [[]] # i'th element is a list wit ki = len(leftBdrWidth[i]) specifying the stencil weights of the first ki time series terms in the stencil for x(t_i) + leftBdrNorm : list[float] = [] # ??? + + """ + Suppose that InteriorStencils is an array of length Ns and interiorIndexes is a list of + length Ns. Then the stencil contains Ns terms representing x(t_0), ... , x(t_{Nx -1}). Further + suppose the underlying time series contains Nx points. In this case, assuming index i is an + "interior index" (not too close to 0 or Nx - 1), then we approximate the time derivative of z + at time t_i as follows: + z'(t_i) \approx c_0 z(t_{i + i(0)}) + ... + c_{Ns - 1} z(t_{i + i(Ns - 1)}) + where c_k = interiorStencils[k] and i(k) = interiorIndexes[k]. Note that the indices may be + negative or positive. Thus, interiorStencils and interiorIndexes tell us how to construct the + finite difference scheme away from the boundary. + + For instance, the central difference scheme corresponds to interiorStencils = [-1/2, 1/2], + interiorIndexes = [-1, 1] and + z'(t_i) \approx (1/2)(-z(t_{i - 1}) + z_{t_{i + 1}}) + + Note: We assume that interiorIndexes is in ASCENDING order. + """ + interiorStencils : np.ndarray = np.array([]) # Specify the weights of the terms in the stencil + interiorIndexes : list[int] = [] # Specify the relative indices of the terms in the stencil + + + + def getOperators(self, Nx : int, periodic : bool = False) -> tuple(spmatrix, torch.Tensor, torch.Tensor): + """ + The stencil class acts as an abstract base class for finite difference schemes. We assume + that the user has a time series, x(t_0), ... , x(t_{Nx - 1}) \in \mathbb{R}^d. We will + further assume that the time stamps are evenly spaced, t_0, ... , t_{Nx - 1}. That is, + there is some \Delta t such that t_k = t_0 + \Delta_t*k. We also assume that the user wants + to approximate the time derivative of x at t_0, ... , t_{Nx - 1} using a finite difference + scheme. + + The getOperators method builds an a sparse Tensor housing the "operator matrix". This is + an Nx x Nx matrix that we use to apply the finite difference scheme to one component of + the time series. How does this work? Suppose the time series is x(t_0), ... , + x(t_{Nx - 1}) \in \mathbb{R}^d. Then, for each i \in \{ 1, 2, ... , d}, we construct Dxi + such that Dxi * xi is the vector whose i'th entry holds the approximation (using the + selected finite difference scheme) of x_i'(t_i). Here, xi = [x_i(t_0), ... , + x_i(t_{Nx - 1})]. The line of code below sets up Dxi such that it gives the correct + approximation to x_i'(t_j) whenever j is an interior index. The rest of this function + adjusts the first/last few rows of Dxi so that it also gives the correct approximation + at the boundaries (j close to to 0 or Nx - 1). + + The Stencil class is not a standalone class; you shouldn't use it directly. Rather, you + should use one of the sub-classes defined below. Each one implements a specific finite + difference scheme to approximate the time derivatives at t_0, ... , t_{Nx - 1} (they may + use different rules on the left and right boundary). Each sub-class should set the + interiorStencils and interiorIndexes attributes. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + Nx: An integer specifying the number of points in the time series we want to apply a + finite difference scheme to. + + periodic: A boolean specifying if we should treat the time series as periodic or not. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Three elements. The first Dxi, the "Operator Matrix" described above. The second holds the + "norm" tensor and the third holds the "PeriodicOffset" tensor + (TODO: what are the last two of these used for?) + """ + + norm = np.ones(Nx,) + periodicOffset = np.zeros(Nx,) + + # Start building the "Operator Matirx" (See docstring) + Dxi : spmatrix = sps.diags(self.interiorStencils, + self.interiorIndexes, + shape = (Nx, Nx), + format = 'lil') + + """ + If we assume the time series is periodic, then we need to adjust the first/last few + rows to "wrap around". In this case, every index is an interior index. For this to work, + however, we may need to be able to fetch x(t_{-k}), or x(t_{Nx + k}). Periodicity + means that we assume x(t_{-k}) = x(t_{Nx - k}) and x(t_{Nx + k}) = x(t_k). This means + that we can use the last few components of x in place of terms like x(t_{-k}). For + instance, if we wanted to use a central difference scheme, then we have + x'(t_0) \approx (1/2)(-x(t_{-1}) + x(t_{1})) = (1/2)(-x(t_{Nx - 1}) + x(t_1)) + In other words, by modifying the first/last few rows of Dxi, we can build an operator + matrix that can recover our approximation of x' near the boundaries (think about it). + """ if (periodic): - bdrIdxes = ([k for k in range(-self.interiorIndexes[0])] + - [k for k in range(-self.interiorIndexes[-1], 0)]) + """ + Get the indices less than zero and greater than Nx that we will need to update the + first/last few rows ofDxi. Remember that the elements of interiorIndexes are in + ascending order. Thus, if interiorIndexes[0] = -k, then we need to update the first + k rows of Dxi to account for the periodic BCs (think about it). Likewise, if + interiorIndexes[-1] = j, then we also need to update the final j rows of Dxi (think + about it). + + If the smallest and largest elements of interiorIndexes are -k and j, respectively, + then the line below produces the list [0, 1, ... , k - 1, -j, -j + 1, ... , -1]. + """ + bdrIdxes : list[int] = ([k for k in range(-self.interiorIndexes[0])] + + [k for k in range(-self.interiorIndexes[-1], 0)]) + + # Cycle through the rows of Dxi that we need to correct. for idx in bdrIdxes: - colIdx = [k + idx for k in self.interiorIndexes] - Dxi[idx, colIdx] = self.interiorStencils + # Determine which columns of Dxi we need to access to apply the stencil for the + # idx'th row of Dxi (corresponding to time t_i). + colIdx : list[int] = [k + idx for k in self.interiorIndexes] + + # Now set the corresponding elements with appropriate wrap around (note that if + # one index is negative, then numpy's indexing rules will automatically wrap it + # since x[-k] = x[x.len - k]). + Dxi[idx, colIdx] = self.interiorStencils + + # Set up the periodicOffset. If idx is negative, we are fixing one of the last + # few rows of Dxi. In this case, periodicOffset holds the sum of the weights + # of the indices at and to the right of the current index in the interiorStencil. + # Otherwise, we set it to the sum of the weights of the indices to the left. if (idx < 0): - temp = [k>=0 for k in colIdx] + temp = [k >= 0 for k in colIdx] periodicOffset[idx] = np.sum(self.interiorStencils[temp]) else: - temp = [k<0 for k in colIdx] + temp = [k < 0 for k in colIdx] periodicOffset[idx] = -np.sum(self.interiorStencils[temp]) + + # Otherwise, we need use special stencils for the left and right boundary terms. else: - Dxi[:self.leftBdrDepth,:] = 0. - Dxi[-self.leftBdrDepth:,:] = 0. + # If leftBrdDepth = k, then we use a special stencil to approximate x'(t_0), ... , + # x'(t_{k - 1}) and for x'(t_{Nx - k}), ... , x'(t_{Nx - 1}). We begin by zeroing out + # the first/last k rows of Dxi. + Dxi[:self.leftBdrDepth,:] = 0. + Dxi[-self.leftBdrDepth:,:] = 0. + + # Now, populate the first/last few rows of Dxi to implement the boundary stencil. for depth in range(self.leftBdrDepth): - width = self.leftBdrWidth[depth] + # leftBdrWidth[depth] = k means that the stencil for x(t_{depth}) depends on + # x(t_0), ... , x(t_{k - 1}). + width : int = self.leftBdrWidth[depth] + + # Update the first width elements of the depth'th row of Dxi using the depth'th + # elements of leftBdrStencils (which holds the stencil weights for this row's + # stencil) Dxi[depth,:width] = self.leftBdrStencils[depth] + + # Update the (depth - 1)'th to last row of Dxi using the reversed stencil (which + # gives us an equivalent stencil for right boundaries). # NOTE: Currently only consider skew-symmetric operators. Dxi[-1-depth,-width:] = -Dxi[depth,(width-1)::-1] - norm[:self.leftBdrDepth] = self.leftBdrNorm - norm[-self.leftBdrDepth:] = norm[(self.leftBdrDepth-1)::-1] + + # ??? + norm[:self.leftBdrDepth] = self.leftBdrNorm + norm[-self.leftBdrDepth:] = norm[(self.leftBdrDepth-1)::-1] + + + # Convert Dxi to a tensor. + Dxf : torch.Tensor = self.convert(sps.coo_matrix(Dxi)) - Dxi = self.convert(sps.coo_matrix(Dxi)) + # All done! return Dxi, torch.Tensor(norm), torch.Tensor(periodicOffset) - def convert(self, scipy_coo): + + + def convert(self, scipy_coo : spmatrix) -> torch.Tensor: + """ + Converts scipy_coo, a sparse numpy array, to a sparse torch Tensor. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + scipy_coo: A sparse numpy array. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A tensor version of scipy_coo. + """ + if (type(scipy_coo) is not sps._coo.coo_matrix): raise RuntimeError("The input sparse matrix must be in scipy COO format!") - values = scipy_coo.data - indices = np.vstack((scipy_coo.row, scipy_coo.col)) - i = torch.LongTensor(indices) - v = torch.FloatTensor(values) + # Fetch the values of the non-zero entries of scipy_coo. + values : np.ndarray = scipy_coo.data + + # Get the row, column numbers of the non-zero entries of scipy_coo. vertically stack them + # into a 2 x N matrix (N being the number of non-zero entries in scipy_coo) whose j'th + # column holds the row and column number of the j'th non-zero cell in scipy_coo. + indices : np.ndarray = np.vstack((scipy_coo.row, scipy_coo.col)) + + # Convert indices to a tensor + i : torch.Tensor = torch.LongTensor(indices) + + # Convert the list of non-zero values in scipy_coo to a tensor. + v : torch.Tensor = torch.FloatTensor(values) + + # Get the shape of scipy.coo shape = scipy_coo.shape + # Convert scipy_coo to a tensor and return! return torch.sparse_coo_tensor(i, v, torch.Size(shape)) + + +# ------------------------------------------------------------------------------------------------- +# Stencil sub-classes +# ------------------------------------------------------------------------------------------------- + class SBP12(Stencil): def __init__(self): self.leftBdrDepth = 1 @@ -66,7 +241,9 @@ def __init__(self): self.interiorStencils = np.array([-0.5, 0.5]) self.interiorIndexes = [-1, 1] return - + + + class SBP24(Stencil): def __init__(self): self.leftBdrDepth = 4 @@ -81,6 +258,8 @@ def __init__(self): self.interiorIndexes = [-2, -1, 1, 2] return + + class SBP36(Stencil): def __init__(self): self.leftBdrDepth = 6 @@ -98,6 +277,8 @@ def __init__(self): self.interiorIndexes = [-3, -2, -1, 1, 2, 3] return + + class SBP48(Stencil): def __init__(self): self.leftBdrDepth = 8 @@ -216,6 +397,8 @@ def __init__(self): self.interiorIndexes = [-4, -3, -2, -1, 1, 2, 3, 4] return + + FDdict = {'sbp12': SBP12(), 'sbp24': SBP24(), 'sbp36': SBP36(), From d853d64b3d741718629d29a8fdc5b69c6b55bb50 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Wed, 30 Oct 2024 12:53:08 -0700 Subject: [PATCH 15/22] Started documenting LatentDynamics I qucikly realized, however, that I could not document these files without understanding fd.py... so I did that one first (see the previous commit). I will hopefully finish the LatentDynamics stuff in the next commit or two. --- src/lasdi/latent_dynamics/__init__.py | 70 ++++++++++++-- src/lasdi/latent_dynamics/sindy.py | 133 ++++++++++++++++++++++---- 2 files changed, 172 insertions(+), 31 deletions(-) diff --git a/src/lasdi/latent_dynamics/__init__.py b/src/lasdi/latent_dynamics/__init__.py index 6efa6d3..de8eb0f 100644 --- a/src/lasdi/latent_dynamics/__init__.py +++ b/src/lasdi/latent_dynamics/__init__.py @@ -1,21 +1,63 @@ +# ------------------------------------------------------------------------------------------------- +# Imports and Setup +# ------------------------------------------------------------------------------------------------- + import numpy as np import torch + + +# ------------------------------------------------------------------------------------------------- +# LatentDynamics base class +# ------------------------------------------------------------------------------------------------- + class LatentDynamics: - dim = -1 - nt = -1 - ncoefs = -1 - # TODO(kevin): do we want to store coefficients as a member variable? - coefs = torch.Tensor([]) - - def __init__(self, dim_, nt_): - self.dim = dim_ - self.nt = nt_ + # Class variables + dim : int = -1 # Dimensionality of the latent space + nt : int = -1 # Number of time steps when solving the latent dynamics + ncoefs : int = -1 # Number of coefficients in the latent space dynamics + + # TODO(kevin): do we want to store coefficients as an instance variable? + coefs : torch.Tensor = torch.Tensor([]) + + + + def __init__(self, dim_ : int, nt_ : int) -> None: + """ + Initializes a LatentDynamics object. Each LatentDynamics object needs to have a + dimensionality (dim), a number of time steps, a model for the latent space dynamics, and + set of coefficients for that model. The model should describe a set of ODEs in + \mathbb{R}^{dim}. We + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + dim_: The number of dimensions in the latent space, where the latent dynamics takes place. + + nt_: The number of time steps we want to generate when solving (numerically) the latent + space dynamics. + """ + + # Set class variables. + self.dim : int = dim_ + self.nt : int = nt_ + + # There must be at least one latent dimension and there must be at least 1 time step. assert(self.dim > 0) assert(self.nt > 0) + + # All done! return - def calibrate(self, Z, dt, compute_loss=True, numpy=False): + + + def calibrate(self, + Z : torch.Tensor, + dt : int, + compute_loss : bool = True, + numpy : bool = False) -> np.ndarray: ''' calibrate coefficients of the latent dynamics and compute loss for one given time series Z. @@ -28,6 +70,8 @@ def calibrate(self, Z, dt, compute_loss=True, numpy=False): else: return coefs + + def simulate(self, coefs, z0, t_grid): ''' time-integrate with one initial condition z0 at time points t_grid, @@ -36,6 +80,8 @@ def simulate(self, coefs, z0, t_grid): raise RuntimeError('Abstract function LatentDynamics.simulate!') return zhist + + def sample(self, coefs_sample, z0_sample, t_grid): ''' Sample time series for given sample initial conditions and coefficients. @@ -53,10 +99,14 @@ def sample(self, coefs_sample, z0_sample, t_grid): return Z_simulated + + def export(self): param_dict = {'dim': self.dim, 'ncoefs': self.ncoefs} return param_dict + + # SINDy does not need to load parameters. # Other latent dynamics might need to. def load(self, dict_): diff --git a/src/lasdi/latent_dynamics/sindy.py b/src/lasdi/latent_dynamics/sindy.py index 25780a4..70c380b 100644 --- a/src/lasdi/latent_dynamics/sindy.py +++ b/src/lasdi/latent_dynamics/sindy.py @@ -1,34 +1,116 @@ -import numpy as np -import torch -from scipy.integrate import odeint -from . import LatentDynamics -from ..inputs import InputParser -from ..fd import FDdict +# ------------------------------------------------------------------------------------------------- +# Imports and Setup +# ------------------------------------------------------------------------------------------------- + +import numpy as np +import torch +from scipy.integrate import odeint + +from . import LatentDynamics +from ..inputs import InputParser +from ..fd import FDdict + + + +# ------------------------------------------------------------------------------------------------- +# SINDy class +# ------------------------------------------------------------------------------------------------- class SINDy(LatentDynamics): - fd_type = '' - fd = None - fd_oper = None + fd_type = '' + fd = None + fd_oper = None + + + + def __init__(self, + dim : int, + nt : int, + config : dict) -> None: + """ + Initializes a SINDy object. This is a subclass of the LatentDynamics class which uses the + SINDy algorithm as its model for the ODE governing the latent state. Specifically, we + assume there is a library of functions, f_1(z), ... , f_N(z), each one of which is a + monomial of the components of the latent space, z, and a set of coefficients c_{i,j}, + i = 1, 2, ... , dim and j = 1, 2, ... , N such that + z_i'(t) = \sum_{j = 1}^{N} c_{i,j} f_j(z) + In this case, we assume that f_1, ... , f_N consists of the set of order <= 1 monomials. + That is, f_1(z), ... , f_N(z) = 1, z_1, ... , z_{dim}. + - def __init__(self, dim, nt, config): - super().__init__(dim, nt) + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- - #TODO(kevin): generalize for high-order dynamics - self.ncoefs = self.dim * (self.dim + 1) + dim: The number of dimensions in the latent space, where the latent dynamics takes place. - assert('sindy' in config) - parser = InputParser(config['sindy'], name='sindy_input') + nt: The number of time steps we want to generate when solving (numerically) the latent + space dynamics. - ''' - fd_type is the string that specifies finite-difference scheme for time derivative: + config: A dictionary housing the settings we need to set up a SINDy object. Specifically, + this dictionary should have a key called "sindy" whose corresponding value is another + dictionary with the following two keys: + - fd_type: A string specifying which finite-difference scheme we should use when + approximating the time derivative of the solution to the latent dynamics at a + specific time. Currently, the following options are allowed: - 'sbp12': summation-by-parts 1st/2nd (boundary/interior) order operator - 'sbp24': summation-by-parts 2nd/4th order operator - 'sbp36': summation-by-parts 3rd/6th order operator - 'sbp48': summation-by-parts 4th/8th order operator - ''' - self.fd_type = parser.getInput(['fd_type'], fallback='sbp12') - self.fd = FDdict[self.fd_type] - self.fd_oper, _, _ = self.fd.getOperators(self.nt) + - coef_norm_order: A string specifying which norm we want to use when computing + the coefficient loss. + """ + + # Run the base class initializer. The only thing this does is set the dim and nt + # attributes. + super().__init__(dim, nt) + + # We only allow library terms of order <= 1. If we let z(t) \in \mathbb{R}^{dim} denote the + # latent state at some time, t, then the possible library terms are 1, z_1(t), ... , + # z_{dim}(t). Since each component function gets its own set of coefficients, there must + # be dim*(dim + 1) total coefficients. + #TODO(kevin): generalize for high-order dynamics + self.ncoefs = self.dim * (self.dim + 1) + + # Now, set up an Input parser to process the contents of the config['sindy'] dictionary. + assert('sindy' in config) + parser = InputParser(config['sindy'], name = 'sindy_input') + + """ + Determine which finite difference scheme we should use to approximate the time derivative + of the latent space dynamics. Currently, we allow the following values for "fd_type": + - 'sbp12': summation-by-parts 1st/2nd (boundary/interior) order operator + - 'sbp24': summation-by-parts 2nd/4th order operator + - 'sbp36': summation-by-parts 3rd/6th order operator + - 'sbp48': summation-by-parts 4th/8th order operator + """ + self.fd_type : str = parser.getInput(['fd_type'], fallback = 'sbp12') + self.fd : callable = FDdict[self.fd_type] + + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + # RESUME HERE + self.fd_oper, _, _ = self.fd.getOperators(self.nt) # NOTE(kevin): by default, this will be L1 norm. self.coef_norm_order = parser.getInput(['coef_norm_order'], fallback=1) @@ -36,8 +118,11 @@ def __init__(self, dim, nt, config): # TODO(kevin): other loss functions self.MSE = torch.nn.MSELoss() + # All done! return + + def calibrate(self, Z, dt, compute_loss=True, numpy=False): ''' loop over all train cases, if Z dimension is 3 ''' if (Z.dim() == 3): @@ -86,6 +171,8 @@ def calibrate(self, Z, dt, compute_loss=True, numpy=False): else: return coefs + + def compute_time_derivative(self, Z, Dt): ''' @@ -101,6 +188,8 @@ def compute_time_derivative(self, Z, Dt): ''' return 1. / Dt * torch.sparse.mm(self.fd_oper, Z) + + def simulate(self, coefs, z0, t_grid): ''' @@ -116,6 +205,8 @@ def simulate(self, coefs, z0, t_grid): return Z_i + + def export(self): param_dict = super().export() param_dict['fd_type'] = self.fd_type From 4748e03e888aa86487ef0aeef6180acca800666a Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Mon, 4 Nov 2024 18:07:15 -0800 Subject: [PATCH 16/22] Added documentation to most of latent_dynamics. --- src/lasdi/fd.py | 49 ++++---- src/lasdi/gplasdi.py | 5 +- src/lasdi/latent_dynamics/sindy.py | 174 +++++++++++++++++++++-------- 3 files changed, 155 insertions(+), 73 deletions(-) diff --git a/src/lasdi/fd.py b/src/lasdi/fd.py index b767edf..a72c6a4 100644 --- a/src/lasdi/fd.py +++ b/src/lasdi/fd.py @@ -18,19 +18,18 @@ class Stencil: leftBdrDepth : int = 0 # Number of time steps (at the start/stop of the time series) we use a special stencil on. leftBdrWidth : list[int] = [] # i'th element specifies the number of terms in the stencil for the i'th time step leftBdrStencils : list[list[float]] = [[]] # i'th element is a list wit ki = len(leftBdrWidth[i]) specifying the stencil weights of the first ki time series terms in the stencil for x(t_i) - leftBdrNorm : list[float] = [] # ??? + leftBdrNorm : list[float] = [] # ??? TODO: what is this? """ Suppose that InteriorStencils is an array of length Ns and interiorIndexes is a list of - length Ns. Then the stencil contains Ns terms representing x(t_0), ... , x(t_{Nx -1}). Further - suppose the underlying time series contains Nx points. In this case, assuming index i is an - "interior index" (not too close to 0 or Nx - 1), then we approximate the time derivative of z - at time t_i as follows: + length Ns. Further suppose the underlying time series contains Nx points, x(t_0), ... , + x(t_{Nx -1}). Assuming index i is an "interior index" (not too close to 0 or Nx - 1), then we + approximate the time derivative of z at time t_i as follows: z'(t_i) \approx c_0 z(t_{i + i(0)}) + ... + c_{Ns - 1} z(t_{i + i(Ns - 1)}) where c_k = interiorStencils[k] and i(k) = interiorIndexes[k]. Note that the indices may be negative or positive. Thus, interiorStencils and interiorIndexes tell us how to construct the finite difference scheme away from the boundary. - + x For instance, the central difference scheme corresponds to interiorStencils = [-1/2, 1/2], interiorIndexes = [-1, 1] and z'(t_i) \approx (1/2)(-z(t_{i - 1}) + z_{t_{i + 1}}) @@ -38,18 +37,18 @@ class Stencil: Note: We assume that interiorIndexes is in ASCENDING order. """ interiorStencils : np.ndarray = np.array([]) # Specify the weights of the terms in the stencil - interiorIndexes : list[int] = [] # Specify the relative indices of the terms in the stencil + interiorIndexes : list[int] = [] # Specify the relative indices of the terms in the stencil. Must be in ascending order. - def getOperators(self, Nx : int, periodic : bool = False) -> tuple(spmatrix, torch.Tensor, torch.Tensor): + def getOperators(self, Nx : int, periodic : bool = False) -> tuple[spmatrix, torch.Tensor, torch.Tensor]: """ The stencil class acts as an abstract base class for finite difference schemes. We assume that the user has a time series, x(t_0), ... , x(t_{Nx - 1}) \in \mathbb{R}^d. We will further assume that the time stamps are evenly spaced, t_0, ... , t_{Nx - 1}. That is, - there is some \Delta t such that t_k = t_0 + \Delta_t*k. We also assume that the user wants - to approximate the time derivative of x at t_0, ... , t_{Nx - 1} using a finite difference - scheme. + there is some \Delta t such that t_k = t_0 + \Delta_t*k for each k. We also assume the user + wants to approximate the time derivative of x at t_0, ... , t_{Nx - 1} using a finite + difference scheme. The getOperators method builds an a sparse Tensor housing the "operator matrix". This is an Nx x Nx matrix that we use to apply the finite difference scheme to one component of @@ -57,10 +56,10 @@ def getOperators(self, Nx : int, periodic : bool = False) -> tuple(spmatrix, tor x(t_{Nx - 1}) \in \mathbb{R}^d. Then, for each i \in \{ 1, 2, ... , d}, we construct Dxi such that Dxi * xi is the vector whose i'th entry holds the approximation (using the selected finite difference scheme) of x_i'(t_i). Here, xi = [x_i(t_0), ... , - x_i(t_{Nx - 1})]. The line of code below sets up Dxi such that it gives the correct - approximation to x_i'(t_j) whenever j is an interior index. The rest of this function - adjusts the first/last few rows of Dxi so that it also gives the correct approximation - at the boundaries (j close to to 0 or Nx - 1). + x_i(t_{Nx - 1})]. To build Dxi, we first set up a matrix to give the correct approximation + to x_i'(t_j) whenever j is an interior index. The rest of this function adjusts the + first/last few rows of Dxi so that it also gives the correct approximation at the + boundaries (j close to to 0 or Nx - 1). The Stencil class is not a standalone class; you shouldn't use it directly. Rather, you should use one of the sub-classes defined below. Each one implements a specific finite @@ -84,14 +83,14 @@ def getOperators(self, Nx : int, periodic : bool = False) -> tuple(spmatrix, tor ------------------------------------------------------------------------------------------- Three elements. The first Dxi, the "Operator Matrix" described above. The second holds the - "norm" tensor and the third holds the "PeriodicOffset" tensor - (TODO: what are the last two of these used for?) + "norm" tensor and the third holds the "PeriodicOffset" tensor. + TODO: what are the last two of these used for? And what are they? """ norm = np.ones(Nx,) periodicOffset = np.zeros(Nx,) - # Start building the "Operator Matirx" (See docstring) + # Start building the "Operator Matrix" (See docstring) Dxi : spmatrix = sps.diags(self.interiorStencils, self.interiorIndexes, shape = (Nx, Nx), @@ -99,10 +98,11 @@ def getOperators(self, Nx : int, periodic : bool = False) -> tuple(spmatrix, tor """ If we assume the time series is periodic, then we need to adjust the first/last few - rows to "wrap around". In this case, every index is an interior index. For this to work, - however, we may need to be able to fetch x(t_{-k}), or x(t_{Nx + k}). Periodicity - means that we assume x(t_{-k}) = x(t_{Nx - k}) and x(t_{Nx + k}) = x(t_k). This means - that we can use the last few components of x in place of terms like x(t_{-k}). For + rows to "wrap around". In this case, every index is an interior index. However, this means + that we will need x(t_{-k}) when approximating x'(t_0), for instance, or x(t_{Nx + k}) + when approximating x'(t_{Nx - 1}). Periodicity means we assume x(t_{-k}) = x(t_{Nx - k}) + and x(t_{Nx + k}) = x(t_k). Thus, we can use the last few components of x in place of terms + like x(t_{-k}). For instance, if we wanted to use a central difference scheme, then we have x'(t_0) \approx (1/2)(-x(t_{-1}) + x(t_{1})) = (1/2)(-x(t_{Nx - 1}) + x(t_1)) In other words, by modifying the first/last few rows of Dxi, we can build an operator @@ -169,13 +169,12 @@ def getOperators(self, Nx : int, periodic : bool = False) -> tuple(spmatrix, tor # NOTE: Currently only consider skew-symmetric operators. Dxi[-1-depth,-width:] = -Dxi[depth,(width-1)::-1] - # ??? + # TODO: what is going on here? norm[:self.leftBdrDepth] = self.leftBdrNorm norm[-self.leftBdrDepth:] = norm[(self.leftBdrDepth-1)::-1] - # Convert Dxi to a tensor. - Dxf : torch.Tensor = self.convert(sps.coo_matrix(Dxi)) + Dxi : torch.Tensor = self.convert(sps.coo_matrix(Dxi)) # All done! return Dxi, torch.Tensor(norm), torch.Tensor(periodicOffset) diff --git a/src/lasdi/gplasdi.py b/src/lasdi/gplasdi.py index 02cbe1d..e5776e7 100644 --- a/src/lasdi/gplasdi.py +++ b/src/lasdi/gplasdi.py @@ -478,7 +478,10 @@ def train(self) -> None: # ------------------------------------------------------------------------------------- # Forward pass - # Run the forward pass + # Run the forward pass. This results in a tensor of shape (Np, Nt, Nz), where Np is the + # number of parameters, Nt is the number of time steps in the time series, and Nz is + # the latent space dimension. X_Pred, should have the same shape as X_Train, (Np, Nt, + # Nx[0], .... , Nx[Nd - 1]). Z : torch.Tensor = autoencoder_device.encoder(X_train_device) X_pred : torch.Tensor = autoencoder_device.decoder(Z) Z : torch.Tensor = Z.cpu() diff --git a/src/lasdi/latent_dynamics/sindy.py b/src/lasdi/latent_dynamics/sindy.py index 70c380b..b4c5e2d 100644 --- a/src/lasdi/latent_dynamics/sindy.py +++ b/src/lasdi/latent_dynamics/sindy.py @@ -22,7 +22,6 @@ class SINDy(LatentDynamics): fd_oper = None - def __init__(self, dim : int, nt : int, @@ -87,33 +86,22 @@ def __init__(self, self.fd_type : str = parser.getInput(['fd_type'], fallback = 'sbp12') self.fd : callable = FDdict[self.fd_type] - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE - # RESUME HERE + """ + Fetch the operator matrix. What does this do? Suppose we have a time series with nt points, + x(t_0), ... , x(t_{nt - 1}) \in \mathbb{R}^d. Further assume that for each j, + t_j = t_0 + j \delta t, where \delta t is some positive constant. Let j \in {0, 1, ... , + nt - 1}. Let xj be j'th vector whose k'th element is x_j(t_k). Then, the i'th element of + M xj holds the approximation to x_j'(t_k) using the stencil we selected above. + + For instance, if we selected sdp12, corresponding to the central difference scheme, then + we have (for j != 0, nt - 1) + [M xj]_i = (x_j(t_{i + 1}) - x(t_{j - 1}))/(2 \delta t). + """ self.fd_oper, _, _ = self.fd.getOperators(self.nt) + # Fetch the norm we are going to use on the sindy coefficients. # NOTE(kevin): by default, this will be L1 norm. - self.coef_norm_order = parser.getInput(['coef_norm_order'], fallback=1) + self.coef_norm_order = parser.getInput(['coef_norm_order'], fallback = 1) # TODO(kevin): other loss functions self.MSE = torch.nn.MSELoss() @@ -123,45 +111,126 @@ def __init__(self, - def calibrate(self, Z, dt, compute_loss=True, numpy=False): - ''' loop over all train cases, if Z dimension is 3 ''' + def calibrate(self, + Z : torch.Tensor, + dt : float, + compute_loss : bool = True, + numpy : bool = False) -> (np.ndarray | torch.Tensor) | tuple[(np.ndarray | torch.Tensor), torch.Tensor, torch.Tensor]: + """ + This function computes the optimal SINDy coefficients using the current latent time + series'. Specifically, let us consider the case when Z has two dimensions (the case when + it has three is identical, just with different coefficients for each instance of the + leading dimension of Z). In this case, we assume that the rows of Z correspond to a + trajectory of latent states. Specifically, we assume the i'th row holds the latent state, + z, at time t_0 + i*dt. We use SINDy to find the coefficients in the dynamical system + z'(t) = C \Phi(z(t)), where C is a matrix of coefficients and \Phi(z(t)) represents a + library of terms. We find the matrix C corresponding to the dynamical system that best + agrees with the data in the rows of Z. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + Z: A 2 or 3d tensor. If Z is a 2d tensor, then it has shape (Nt, Nz) and it's i,j entry + holds the j'th component of the latent state at the time t_0 + i*dt. If it is a 3d tensor, + then it has shape (Np, Nt, Nz). In this case, we assume there at Np different combinations + of parameter values. The i, j, k entry of Z in this case holds the k'th component of the + latent encoding at time t_0 + j*dt when we use the i'th combination of parameter values. + + dt: The time step between time steps. See the description of the "Z" argument. + + compute_loss: A boolean which, if true, will prompt us to calculate the sindy and + coefficient losses, which we will then return with the optimal SINDy coefficients. If not, + we will only return the optimal SINDy coefficients. + + numpy: A boolean. If True, we return the coefficient matrix as a numpy.ndarray object. If + False, we return it as a torch.Tensor object. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + IF compute_loss is True, then we return three variables. The first is a matrix of shape + () + """ + + # ----------------------------------------------------------------------------------------- + # If Z has three dimensions, loop over all train cases. if (Z.dim() == 3): - n_train = Z.size(0) + # Fetch the number of training cases. + n_train : int = Z.size(0) + # Prepare an array to house the flattened coefficient matrices for each combination of + # parameter values. if (numpy): coefs = np.zeros([n_train, self.ncoefs]) else: coefs = torch.Tensor([n_train, self.ncoefs]) + + # Initialize the losses. Note that these are floats which we will replace with + # tensors. loss_sindy, loss_coef = 0.0, 0.0 + # Cycle through the combinations of parameter values. for i in range(n_train): + """" + Get the optimal SINDy coefficients for the i'th combination of parameter values. + Remember that Z is 3d tensor of shape (Np, Nt, Nz) whose (i, j, k) entry holds + the k'th component of the j'th frame of the latent trajectory for the i'th + combination of parameter values. Note that Result is either a torch.Tensor + (if compute_loss = False and numpy = False), a numpy.ndarray (if numpy = True and + compute_loss = True), or a 3 element tuple (if compute_loss = True). + """ result = self.calibrate(Z[i], dt, compute_loss, numpy) + + # If we are computing losses, the 1 and 2 elements of return hold the sindy and + # coefficient losses. Otherwise, the only return variable is the flattened + # coefficient matrix from the i'th combination of parameter values. if (compute_loss): - coefs[i] = result[0] + coefs[i] = result[0] loss_sindy += result[1] - loss_coef += result[2] + loss_coef += result[2] else: coefs[i] = result + # Package everything to return! if (compute_loss): return coefs, loss_sindy, loss_coef else: return coefs - ''' evaluate for one train case ''' + + # ----------------------------------------------------------------------------------------- + # evaluate for one training case. assert(Z.dim() == 2) + + # First, compute the time derivatives. This yields a torch.Tensor object whose i,j entry + # holds an approximation of (d/dt) Z_j(t_0 + i*dt) dZdt = self.compute_time_derivative(Z, dt) time_dim, space_dim = dZdt.shape - Z_i = torch.cat([torch.ones(time_dim, 1), Z], dim = 1) - coefs = torch.linalg.lstsq(Z_i, dZdt).solution - + # Concatenate a column of zeros. + Z_i : torch.Tensor = torch.cat([torch.ones(time_dim, 1), Z], dim = 1) + + # For each j, solve the least squares problem + # min{ || dZdt[:, j] - Z_i c_j|| : C_j \in \mathbb{R}ˆNl } + # where Nl is the number of library terms (in this case, just Nz + 1, since we only allow + # constant and linear terms). We store the resulting solutions in a matrix, coefs, whose + # j'th column holds the results for the j'th column of dZdt. Thus, coefs is a 2d tensor + # with shape (Nl, Nz). + coefs : torch.Tensor = torch.linalg.lstsq(Z_i, dZdt).solution + + # If we need to compute the loss, do so now. if (compute_loss): loss_sindy = self.MSE(dZdt, Z_i @ coefs) # NOTE(kevin): by default, this will be L1 norm. loss_coef = torch.norm(coefs, self.coef_norm_order) - # output of lstsq is not contiguous in memory. + # All done. Prepare coefs and the losses to return. Note that we flatten the coefficient + # matrix. + # Note: output of lstsq is not contiguous in memory. coefs = coefs.detach().flatten() if (numpy): coefs = coefs.numpy() @@ -173,30 +242,41 @@ def calibrate(self, Z, dt, compute_loss=True, numpy=False): - def compute_time_derivative(self, Z, Dt): + def compute_time_derivative(self, Z : torch.Tensor, Dt : float) -> torch.Tensor: + """ + This function builds the SINDy dataset, assuming only linear terms in the SINDy dataset. + The time derivatives are computed through finite difference. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + Z: A 2d tensor of shape (Nt, Nz) whose i, j entry holds the j'th component of the i'th + time step in the latent time series. We assume that Z[i, :] represents the latent state + at time t_0 + i*Dt - ''' + Dt: The time step between latent frames (the time between Z[:, i] and Z[:, i + 1]) - Builds the SINDy dataset, assuming only linear terms in the SINDy dataset. The time derivatives are computed through - finite difference. - Z is the encoder output (2D tensor), with shape [time_dim, space_dim] - Dt is the size of timestep (assumed to be a uniform scalar) + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- - The output dZdt is a 2D tensor with the same shape of Z. + The output dZdt is a 2D tensor with the same shape as Z. It's i, j entry holds an + approximation to (d/dt)Z_j(t_0 + j Dt). We compute this approximation using self's stencil. + """ - ''' - return 1. / Dt * torch.sparse.mm(self.fd_oper, Z) + return (1. / Dt) * torch.sparse.mm(self.fd_oper, Z) def simulate(self, coefs, z0, t_grid): - - ''' + """ Integrates each system of ODEs corresponding to each training points, given the initial condition Z0 = encoder(U0) - ''' + """ # copy is inevitable for numpy==1.26. removed copy=False temporarily. c_i = coefs.reshape([self.dim+1, self.dim]).T dzdt = lambda z, t : c_i[:, 1:] @ z + c_i[:, 0] From 9a5acbdf44956d2e0cf1e6e86bb4d2561350cbd9 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Tue, 5 Nov 2024 12:29:55 -0800 Subject: [PATCH 17/22] Finished documenting latent_dynamics/sindy.py --- src/lasdi/latent_dynamics/sindy.py | 74 ++++++++++++++++++++++++++---- 1 file changed, 66 insertions(+), 8 deletions(-) diff --git a/src/lasdi/latent_dynamics/sindy.py b/src/lasdi/latent_dynamics/sindy.py index b4c5e2d..4a22502 100644 --- a/src/lasdi/latent_dynamics/sindy.py +++ b/src/lasdi/latent_dynamics/sindy.py @@ -271,24 +271,82 @@ def compute_time_derivative(self, Z : torch.Tensor, Dt : float) -> torch.Tensor: - def simulate(self, coefs, z0, t_grid): + def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray): """ + For each combination of training parameters, this function integrates the corresponding + system of ODEs given the initial condition Z0 = encoder(U0) - Integrates each system of ODEs corresponding to each training points, given the initial condition Z0 = encoder(U0) + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + coefs: A one dimensional numpy.ndarray object representing the flattened copy of the array + of latent dynamics coefficients that calibrate returns. + + z0: A numpy ndarray object of shape nz representing the initial condition for the latent + dynamics. Thus, the i'th component of this array should hold the i'th component of the + latent dynamics initial condition. + + t_grid: A 1d numpy ndarray object whose i'th entry holds the value of the i'th time value + where we want to compute the latent solution. The elements of this array should be in + ascending order. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A 2d numpy.ndarray object holding the solution to the latent dynamics at the time values + specified in t_grid. Specifically, this is a 2d array of shape (nt, nz), where nt is the + number of time steps (size of t_grid) and nz is the latent space dimension (self.dim). + Thus, the i,j element of this matrix holds the j'th component of the latent solution at + the time stored in the i'th element of t_grid. """ - # copy is inevitable for numpy==1.26. removed copy=False temporarily. - c_i = coefs.reshape([self.dim+1, self.dim]).T - dzdt = lambda z, t : c_i[:, 1:] @ z + c_i[:, 0] + # First, reshape coefs as a matrix. Since we only allow for linear terms, there are nz + 1 + # library terms and nz equations, where nz = self.dim. + # Note: copy is inevitable for numpy==1.26. removed copy=False temporarily. + c_i = coefs.reshape([self.dim + 1, self.dim]).T + + # Set up a lambda function to approximate dzdt. In SINDy, we learn a coefficient matrix + # C such that the latent state evolves according to the dynamical system + # z'(t) = C \Phi(z(t)), + # where \Phi(z(t)) is the library of terms. Note that the zero column of C corresponds + # to the constant library term, 1. + dzdt = lambda z, t : c_i[:, 1:] @ z + c_i[:, 0] + + # Solve the ODE forward in time. Z_i = odeint(dzdt, z0, t_grid) + # All done! return Z_i def export(self): - param_dict = super().export() - param_dict['fd_type'] = self.fd_type - param_dict['coef_norm_order'] = self.coef_norm_order + """ + This function packages self's contents into a dictionary which it then returns. We can use + this dictionary to create a new SINDy object which has the same internal state as self. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + None! + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A dictionary with two keys: fd_type and coef_norm_order. The former specifies which finite + different scheme we use while the latter specifies which norm we want to use when computing + the coefficient loss. + """ + + param_dict = super().export() + param_dict['fd_type'] = self.fd_type + param_dict['coef_norm_order'] = self.coef_norm_order return param_dict \ No newline at end of file From 7fa7fb132f73521af98f36022f97a28b6c4905c6 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Tue, 5 Nov 2024 13:42:07 -0800 Subject: [PATCH 18/22] Finished documenting latent_dynamics --- src/lasdi/latent_dynamics/__init__.py | 131 +++++++++++++++++++++++--- src/lasdi/latent_dynamics/sindy.py | 19 ++-- 2 files changed, 128 insertions(+), 22 deletions(-) diff --git a/src/lasdi/latent_dynamics/__init__.py b/src/lasdi/latent_dynamics/__init__.py index de8eb0f..bbe56b2 100644 --- a/src/lasdi/latent_dynamics/__init__.py +++ b/src/lasdi/latent_dynamics/__init__.py @@ -27,7 +27,9 @@ def __init__(self, dim_ : int, nt_ : int) -> None: Initializes a LatentDynamics object. Each LatentDynamics object needs to have a dimensionality (dim), a number of time steps, a model for the latent space dynamics, and set of coefficients for that model. The model should describe a set of ODEs in - \mathbb{R}^{dim}. We + \mathbb{R}^{dim}. These ODEs should contain a set of unknown coefficients. We learn those + coefficients using the calibrate function. Once we have learned the coefficients, we can + solve the corresponding set of ODEs forward in time using the simulate function. ------------------------------------------------------------------------------------------- @@ -58,12 +60,54 @@ def calibrate(self, dt : int, compute_loss : bool = True, numpy : bool = False) -> np.ndarray: - ''' - calibrate coefficients of the latent dynamics and compute loss - for one given time series Z. + """ + The user must implement this class on any latent dynamics sub-class. Each latent dynamics + object should parameterize a model for the dynamics in the latent space. Thus, to specify + the latent dynamics, we need to find a set of coefficients that parameterize the equations + in the latent dynamics model. The calibrate function is supposed to do that. + + Specifically, this function should take in a sequence (or sequences) of latent states, a + time step, dt, which specifies the time step between successive terms in the sequence(s) of + latent states, and some optional booleans which control what information we return. + + The function should always return at least one argument: a numpy.ndarray object holding the + optimal coefficients in the latent dynamics model using the data contained in Z. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + Z: A 2d or 3d tensor. If Z is a 2d tensor, then it has shape (Nt, Nz), where Nt specifies + the number of time steps in each sequence of latent states and Nz is the dimension of the + latent space. In this case, the i,j entry of Z holds the j'th component of the latent state + at the time t_0 + i*dt. If it is a 3d tensor, then it has shape (Np, Nt, Nz). In this case, + we assume there at Np different combinations of parameter values. The i, j, k entry of Z in + this case holds the k'th component of the latent encoding at time t_0 + j*dt when we use + he i'th combination of parameter values. + + dt: The time step between time steps. See the description of the "Z" argument. + + compute_loss: A boolean which, if True, this function should return the coefficients and + additional losses based on the set of coefficients we learn. If False, this function should + only return the optimal coefficients for the latent dynamics model using the data in Z. + + numpy: A boolean. If True, this function should return the coefficient matrix as a + numpy.ndarray object. If False, this function should return it as a torch.Tensor object. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A tensor or ndarray (depending on the value of the "numpy" argument) holding the optimal + coefficients for the latent space dynamics given the data stored in Z. If Z is 2d, then + the returned tensor will only contain one set of coefficients. If Z is 3d, with a leading + dimension size of Np (number of combinations of parameter values) then we will return + an array/tensor with a leading dimension of size Np whose i'th entry holds the coefficients + for the sequence of latent states stored in Z[:, ...]. + """ - Z is the encoder output (2D/3D tensor), with shape [time_dim, space_dim] or [n_train, time_dim, space_dim] - ''' raise RuntimeError('Abstract function LatentDynamics.calibrate!') if (compute_loss): return coefs, loss @@ -72,20 +116,79 @@ def calibrate(self, - def simulate(self, coefs, z0, t_grid): - ''' - time-integrate with one initial condition z0 at time points t_grid, - for one given set of coefficients - ''' + def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray) -> np.ndarray: + """ + Time integrates the latent dynamics when it uses the coefficients specified in coefs and + starts from the (single) initial condition in z0. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + coefs: A one dimensional numpy.ndarray object holding the coefficients we want to use + to solve the latent dynamics forward in time. + + z0: A numpy ndarray object of shape nz representing the initial condition for the latent + dynamics. Thus, the i'th component of this array should hold the i'th component of the + latent dynamics initial condition. + + t_grid: A 1d numpy ndarray object whose i'th entry holds the value of the i'th time value + where we want to compute the latent solution. The elements of this array should be in + ascending order. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A 2d numpy.ndarray object holding the solution to the latent dynamics at the time values + specified in t_grid when we use the coefficients in coefs to characterize the latent + dynamics model. Specifically, this is a 2d array of shape (nt, nz), where nt is the + number of time steps (size of t_grid) and nz is the latent space dimension (self.dim). + Thus, the i,j element of this matrix holds the j'th component of the latent solution at + the time stored in the i'th element of t_grid. + """ + raise RuntimeError('Abstract function LatentDynamics.simulate!') return zhist def sample(self, coefs_sample, z0_sample, t_grid): - ''' - Sample time series for given sample initial conditions and coefficients. - ''' + """ + simulate's the latent dynamics for a set of coefficients/initial conditions. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + coefs_sample: A numpy.ndarray object whose leading dimension has size ns (the number of + sets of coefficients/initial conditions/simulations we run). + + z0_sample: A 2d numpy.ndarray object of shape (ns, nz) (where ns is the number of samples + and nz is the dimensionality of the latent space). The i,j entry of z0_sample should hold + the j'th component of the i'th initial condition. + + t_grid: A 1d numpy ndarray object whose i'th entry holds the value of the i'th time value + where we want to compute each latent solution. The elements of this array should be in + ascending order. We use the same array for each set of coefficients. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + + A 3d numpy ndarray object of shape (ns, nt, nz), where ns = the number of samples (the + leading dimension of z0_sample and coefs_sample), nt = the number of time steps (size of + t_grid) and nz is the dimension of the latent space. The i, j, k element of this array + holds the k'th component of the solution of the latent dynamics at the j'th time step (j'th + element of t_grid) when we use the i'th set of coefficients/initial conditions. + """ + + # There needs to be as many initial conditions as sets of coefficients. assert(len(coefs_sample) == len(z0_sample)) for i in range(len(coefs_sample)): diff --git a/src/lasdi/latent_dynamics/sindy.py b/src/lasdi/latent_dynamics/sindy.py index 4a22502..e8c5dba 100644 --- a/src/lasdi/latent_dynamics/sindy.py +++ b/src/lasdi/latent_dynamics/sindy.py @@ -132,11 +132,13 @@ def calibrate(self, Arguments ------------------------------------------------------------------------------------------- - Z: A 2 or 3d tensor. If Z is a 2d tensor, then it has shape (Nt, Nz) and it's i,j entry - holds the j'th component of the latent state at the time t_0 + i*dt. If it is a 3d tensor, - then it has shape (Np, Nt, Nz). In this case, we assume there at Np different combinations - of parameter values. The i, j, k entry of Z in this case holds the k'th component of the - latent encoding at time t_0 + j*dt when we use the i'th combination of parameter values. + Z: A 2d or 3d tensor. If Z is a 2d tensor, then it has shape (Nt, Nz), where Nt specifies + the number of time steps in each sequence of latent states and Nz is the dimension of the + latent space. In this case, the i,j entry of Z holds the j'th component of the latent state + at the time t_0 + i*dt. If it is a 3d tensor, then it has shape (Np, Nt, Nz). In this case, + we assume there at Np different combinations of parameter values. The i, j, k entry of Z in + this case holds the k'th component of the latent encoding at time t_0 + j*dt when we use + he i'th combination of parameter values. dt: The time step between time steps. See the description of the "Z" argument. @@ -273,8 +275,8 @@ def compute_time_derivative(self, Z : torch.Tensor, Dt : float) -> torch.Tensor: def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray): """ - For each combination of training parameters, this function integrates the corresponding - system of ODEs given the initial condition Z0 = encoder(U0) + Time integrates the latent dynamics when it uses the coefficients specified in coefs and + starts from the (single) initial condition in z0. ------------------------------------------------------------------------------------------- @@ -298,7 +300,8 @@ def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray): ------------------------------------------------------------------------------------------- A 2d numpy.ndarray object holding the solution to the latent dynamics at the time values - specified in t_grid. Specifically, this is a 2d array of shape (nt, nz), where nt is the + specified in t_grid when we use the coefficients in coefs to characterize the latent + dynamics model. Specifically, this is a 2d array of shape (nt, nz), where nt is the number of time steps (size of t_grid) and nz is the latent space dimension (self.dim). Thus, the i,j element of this matrix holds the j'th component of the latent solution at the time stored in the i'th element of t_grid. From 44b2da9c548c4f32f69732f4df0b1986ee82b68d Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Tue, 5 Nov 2024 13:42:51 -0800 Subject: [PATCH 19/22] Added a few comments I meant to add these into the previous commit but forgot to save the file. Oops. Anyway, with these additions, the latent_dynamics files are now fully documented. --- src/lasdi/latent_dynamics/__init__.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/lasdi/latent_dynamics/__init__.py b/src/lasdi/latent_dynamics/__init__.py index bbe56b2..fe9f7ad 100644 --- a/src/lasdi/latent_dynamics/__init__.py +++ b/src/lasdi/latent_dynamics/__init__.py @@ -155,7 +155,7 @@ def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray) -> - def sample(self, coefs_sample, z0_sample, t_grid): + def sample(self, coefs_sample : np.ndarray, z0_sample : np.ndarray, t_grid : np.ndarray) -> np.ndarray: """ simulate's the latent dynamics for a set of coefficients/initial conditions. @@ -191,15 +191,21 @@ def sample(self, coefs_sample, z0_sample, t_grid): # There needs to be as many initial conditions as sets of coefficients. assert(len(coefs_sample) == len(z0_sample)) + # Cycle through the set of coefficients for i in range(len(coefs_sample)): - Z_i = self.simulate(coefs_sample[i], z0_sample[i], t_grid) + # Simulate the latent dynamics when we use the i'th set of coefficients + ICs + Z_i : np.ndarray = self.simulate(coefs_sample[i], z0_sample[i], t_grid) + + # Append a leading dimension of size 1. Z_i = Z_i.reshape(1, Z_i.shape[0], Z_i.shape[1]) + # Append the latest trajectory onto the Z_simulated array. if (i == 0): Z_simulated = Z_i else: Z_simulated = np.concatenate((Z_simulated, Z_i), axis = 0) + # All done! return Z_simulated From 8e43fac4aec87a83e8f793fa488b1b279372e1e8 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Tue, 5 Nov 2024 14:53:58 -0800 Subject: [PATCH 20/22] Made formatting tweaks to physics files. --- src/lasdi/physics/__init__.py | 49 ++++++++++++++++++++++++++---- src/lasdi/physics/burgers1d.py | 54 +++++++++++++++++++++++++++++----- 2 files changed, 89 insertions(+), 14 deletions(-) diff --git a/src/lasdi/physics/__init__.py b/src/lasdi/physics/__init__.py index 1aa7e84..b9968b6 100644 --- a/src/lasdi/physics/__init__.py +++ b/src/lasdi/physics/__init__.py @@ -1,16 +1,31 @@ -import numpy as np -import torch +# ------------------------------------------------------------------------------------------------- +# Imports +# ------------------------------------------------------------------------------------------------- + +import numpy as np +import torch + + + +# ------------------------------------------------------------------------------------------------- +# Physics class +# ------------------------------------------------------------------------------------------------- class Physics: # Physical space dimension dim = -1 - # solution (denoted as q) dimension + + # The fom solution can be vector valued. If it is, then qdim specifies the dimensionality of + # the fom solution at each point. If the solution is scalar valued, then qdim = -1. qdim = -1 + # grid_size is the shape of the grid nd-array. grid_size = [] - # the shape of the solution nd-array. + + # the shape of the solution nd-array. This is just the qgrid_size with the qdim prepended onto + # it. qgrid_size = [] - + ''' numpy nd-array, assuming the shape of: - 1d: (space_dim[0],) @@ -34,22 +49,32 @@ class Physics: # list of parameter names to parse parameters. param_name = None + + def __init__(self, cfg, param_name=None): self.param_name = param_name return + + def initial_condition(self, param): raise RuntimeError("Abstract method Physics.initial_condition!") return np.array + + def solve(self, param): raise RuntimeError("Abstract method Physics.solve!") return torch.Tensor + + def export(self): raise RuntimeError("Abstract method Physics.export!") return dict + + def generate_solutions(self, params): ''' Given 2d-array of params, @@ -74,10 +99,18 @@ def generate_solutions(self, params): return X_train + + def residual(self, Xhist): raise RuntimeError("Abstract method Physics.residual!") return res, res_norm - + + + +# ------------------------------------------------------------------------------------------------- +# Offline full order model class +# ------------------------------------------------------------------------------------------------- + class OfflineFOM(Physics): def __init__(self, cfg, param_name=None): super().__init__(cfg, param_name) @@ -107,10 +140,14 @@ def __init__(self, cfg, param_name=None): return + + def generate_solutions(self, params): raise RuntimeError("OfflineFOM does not support generate_solutions!!") return + + def export(self): dict_ = {'t_grid' : self.t_grid, 'x_grid' : self.x_grid, 'dt' : self.dt} return dict_ diff --git a/src/lasdi/physics/burgers1d.py b/src/lasdi/physics/burgers1d.py index f7c3bdc..3dbbf0b 100644 --- a/src/lasdi/physics/burgers1d.py +++ b/src/lasdi/physics/burgers1d.py @@ -1,16 +1,30 @@ -import numpy as np -from scipy.sparse.linalg import spsolve -from scipy.sparse import spdiags -import torch -from ..inputs import InputParser -from . import Physics -from ..fd import FDdict +# ------------------------------------------------------------------------------------------------- +# Inputs +# ------------------------------------------------------------------------------------------------- + +import numpy as np +from scipy.sparse.linalg import spsolve +from scipy.sparse import spdiags +import torch + +from ..inputs import InputParser +from . import Physics +from ..fd import FDdict + + + +# ------------------------------------------------------------------------------------------------- +# Burgers 1D class +# ------------------------------------------------------------------------------------------------- class Burgers1D(Physics): + # Class variables a_idx = None # parameter index for a w_idx = None # parameter index for w - def __init__(self, cfg, param_name=None): + + + def __init__(self, cfg, param_name = None): super().__init__(cfg, param_name) self.qdim = 1 @@ -48,6 +62,8 @@ def __init__(self, cfg, param_name=None): self.w_idx = self.param_name.index('w') return + + def initial_condition(self, param): a, w = 1.0, 1.0 if 'a' in self.param_name: @@ -57,6 +73,8 @@ def initial_condition(self, param): return a * np.exp(- self.x_grid ** 2 / 2 / w / w) + + def solve(self, param): u0 = self.initial_condition(param) @@ -64,10 +82,14 @@ def solve(self, param): new_X = new_X.reshape(1, self.nt, self.grid_size[0]) return torch.Tensor(new_X) + + def export(self): dict_ = {'t_grid' : self.t_grid, 'x_grid' : self.x_grid, 'dt' : self.dt, 'dx' : self.dx} return dict_ + + def residual(self, Xhist): # first axis is time index, and second index is spatial index. dUdx = (Xhist[:, 1:] - Xhist[:, :-1]) / self.dx @@ -78,6 +100,12 @@ def residual(self, Xhist): return r, e + + +# ------------------------------------------------------------------------------------------------- +# Helper functions +# ------------------------------------------------------------------------------------------------- + def residual_burgers(un, uw, c, idxn1): ''' @@ -92,6 +120,8 @@ def residual_burgers(un, uw, c, idxn1): return r + + def jacobian(u, c, idxn1, nx): ''' @@ -110,6 +140,8 @@ def jacobian(u, c, idxn1, nx): return J + + def solver(u0, maxk, convergence_threshold, nt, nx, Dt, Dx): ''' @@ -145,6 +177,12 @@ def solver(u0, maxk, convergence_threshold, nt, nx, Dt, Dx): return u + + +# ------------------------------------------------------------------------------------------------- +# Main function (if running this file as a script). +# ------------------------------------------------------------------------------------------------- + def main(): import argparse import yaml From 08004dc6cfc394372fd9a255c24b05914c4eb776 Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Tue, 5 Nov 2024 14:55:13 -0800 Subject: [PATCH 21/22] Added a few typehints to latent_dynamics files This should have gone in a few commits ago but did not. I only changed a few lines, nothing major. --- src/lasdi/latent_dynamics/__init__.py | 4 ++-- src/lasdi/latent_dynamics/sindy.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/lasdi/latent_dynamics/__init__.py b/src/lasdi/latent_dynamics/__init__.py index fe9f7ad..e415751 100644 --- a/src/lasdi/latent_dynamics/__init__.py +++ b/src/lasdi/latent_dynamics/__init__.py @@ -210,7 +210,7 @@ def sample(self, coefs_sample : np.ndarray, z0_sample : np.ndarray, t_grid : np. - def export(self): + def export(self) -> dict: param_dict = {'dim': self.dim, 'ncoefs': self.ncoefs} return param_dict @@ -218,7 +218,7 @@ def export(self): # SINDy does not need to load parameters. # Other latent dynamics might need to. - def load(self, dict_): + def load(self, dict_ : dict) -> None: assert(self.dim == dict_['dim']) assert(self.ncoefs == dict_['ncoefs']) return diff --git a/src/lasdi/latent_dynamics/sindy.py b/src/lasdi/latent_dynamics/sindy.py index e8c5dba..98d45a0 100644 --- a/src/lasdi/latent_dynamics/sindy.py +++ b/src/lasdi/latent_dynamics/sindy.py @@ -273,7 +273,7 @@ def compute_time_derivative(self, Z : torch.Tensor, Dt : float) -> torch.Tensor: - def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray): + def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray) -> np.ndarray: """ Time integrates the latent dynamics when it uses the coefficients specified in coefs and starts from the (single) initial condition in z0. @@ -327,7 +327,7 @@ def simulate(self, coefs : np.ndarray, z0 : np.ndarray, t_grid : np.ndarray): - def export(self): + def export(self) -> dict: """ This function packages self's contents into a dictionary which it then returns. We can use this dictionary to create a new SINDy object which has the same internal state as self. From 9a3a0eda99448d655807b35bdb044f53cc1802df Mon Sep 17 00:00:00 2001 From: Robert Stephany Date: Tue, 5 Nov 2024 18:06:43 -0800 Subject: [PATCH 22/22] Added documentation to most of Physics I haven't totaly finished, but I added documentation to most of the base Physics class + the Burgers1D subclass. I still do not know what the Offline subclass thing does, and haven't looked too much into how we compute the Burgers' residual... but I think I got most of the documentation in this commit. I also made a minor tweak to latent_space.py... the bulk of the changes in this commit are in the physics sub-dictionary, however. --- src/lasdi/latent_space.py | 2 +- src/lasdi/physics/__init__.py | 194 +++++++++++++++++++++++++++------ src/lasdi/physics/burgers1d.py | 194 ++++++++++++++++++++++++++------- 3 files changed, 318 insertions(+), 72 deletions(-) diff --git a/src/lasdi/latent_space.py b/src/lasdi/latent_space.py index 7f64f65..551c6ff 100644 --- a/src/lasdi/latent_space.py +++ b/src/lasdi/latent_space.py @@ -83,7 +83,7 @@ def initial_condition_latent(param_grid : np.ndarray, # Fetch the IC for the i'th set of parameters. Then map it to a tensor. u0 : np.ndarray = physics.initial_condition(param_grid[i]) u0 = u0.reshape(sol_shape) - u0 = torch.Tensor(u0) + u0 = torch.Tensor(u0) # Encode the IC, then map the encoding to a numpy array. z0 : np.ndarray = autoencoder.encoder(u0) diff --git a/src/lasdi/physics/__init__.py b/src/lasdi/physics/__init__.py index b9968b6..8e31b49 100644 --- a/src/lasdi/physics/__init__.py +++ b/src/lasdi/physics/__init__.py @@ -13,18 +13,18 @@ class Physics: # Physical space dimension - dim = -1 + dim : int= -1 # The fom solution can be vector valued. If it is, then qdim specifies the dimensionality of # the fom solution at each point. If the solution is scalar valued, then qdim = -1. - qdim = -1 + qdim : int = -1 # grid_size is the shape of the grid nd-array. - grid_size = [] + grid_size : list[int] = [] # the shape of the solution nd-array. This is just the qgrid_size with the qdim prepended onto # it. - qgrid_size = [] + qgrid_size : list[int] = [] ''' numpy nd-array, assuming the shape of: @@ -33,75 +33,201 @@ class Physics: - 3d: (3, space_dim[0], space_dim[1], space_dim[2]) - higher dimension... ''' - x_grid = np.array([]) + x_grid : np.ndarray = np.array([]) # the number of time steps, as a positive integer. - nt = -1 - # time step size. assume constant for now. - dt = -1. - # time grid in numpy 1d array - t_grid = np.array([]) + nt : int = -1 + + # time step size. assume constant for now. + dt : float = -1. + + # time grid in numpy 1d array. + t_grid : np.ndarray = np.array([]) # Need an offline FOM simulation or not. # Set at the initialization. - offline = False + offline : bool = False # list of parameter names to parse parameters. - param_name = None + param_name : list[str] = None + + def __init__(self, cfg : dict, param_name : list[str] = None) -> None: + """ + A Physics object acts as a wrapper around a solver for a particular equation. The initial + condition in that function can have named parameters. Each physics object should have a + solve method to solve the underlying equation for a given set of parameters, and an + initial condition function to recover the equation's initial condition for a specific set + of parameters. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + cfg: A dictionary housing the settings for the Physics object. This should be the "physics" + sub-dictionary of the main configuration file. - def __init__(self, cfg, param_name=None): + param_name: A list of strings. There should be one list item for each parameter. The i'th + element of this list should be a string housing the name of the i'th parameter. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + self.param_name = param_name return - def initial_condition(self, param): + def initial_condition(self, param : np.ndarray) -> np.ndarray: + """ + The user should write an instance of this method for their specific Physics sub-class. + It should evaluate and return the initial condition along the spatial grid. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param: A 1d numpy.ndarray object holding the value of self's parameters (necessary to + specify the IC). + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A d-dimensional numpy.ndarray object of shape self.grid_size, where + d = len(self.grid_size). This should hold the IC evaluated on self's spatial grid + (self.x_grid) + """ + raise RuntimeError("Abstract method Physics.initial_condition!") return np.array - def solve(self, param): + def solve(self, param : np.ndarray) -> torch.Tensor: + """ + The user should write an instance of this method for their specific Physics sub-class. + This function should solve the underlying equation when the IC uses the parameters in + param. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param: A 1d numpy.ndarray object with two elements corresponding to the values of the + initial condition parameters. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A (ns + 2)-dimensional torch.Tensor object of shape (1, nt, nx[0], .. , nx[ns - 1]), + where nt is the number of points along the temporal grid and nx = self.grid_size specifies + the number of grid points along the axes in the spatial grid. + """ + raise RuntimeError("Abstract method Physics.solve!") return torch.Tensor - def export(self): + def export(self) -> dict: + """ + This function should return a dictionary that houses self's state. I + """ raise RuntimeError("Abstract method Physics.export!") return dict - def generate_solutions(self, params): - ''' - Given 2d-array of params, - generate solutions of size params.shape[0]. - params.shape[1] must match the required size of - parameters for the specific physics. - ''' + def generate_solutions(self, params : np.ndarray) -> torch.Tensor: + """ + Given 2d-array of params, generate solutions of size params.shape[0]. params.shape[1] must + match the required size of parameters for the specific physics. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param: a 2d numpy.ndarray object of shape (np, n), where np is the number of combinations + of parameters we want to test and n denotes the number of parameters in self's initial + condition function. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A torch.Tensor object of shape (np, nt, nx[0], .. , nx[ns - 1]), where nt is the number of + points along the temporal grid and nx = self.grid_size specifies the number of grid points + along the axes in the spatial grid. + """ + + # Make sure we have a 2d grid of parameter values. assert(params.ndim == 2) - n_param = len(params) + n_param : int = len(params) + + # Report print("Generating %d samples" % n_param) + # Cycle through the parameters. X_train = None for k, param in enumerate(params): - new_X = self.solve(param) + # Solve the underlying equation using the current set of parameter values. + new_X : torch.Tensor = self.solve(param) + + # Now, add this solution to the set of solutions. assert(new_X.size(0) == 1) # should contain one parameter case. if (X_train is None): X_train = new_X else: X_train = torch.cat([X_train, new_X], dim = 0) - print("%d/%d complete" % (k+1, n_param)) + print("%d/%d complete" % (k+1, n_param)) + # All done! return X_train - def residual(self, Xhist): + def residual(self, Xhist : np.ndarray) -> tuple[np.ndarray, float]: + """ + The user should write an instance of this method for their specific Physics sub-class. + This function should compute the PDE residual (difference between the left and right hand + side of of the underlying physics equation when we substitute in the solution in Xhist). + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + Xhist: A (ns + 1)-dimensional numpy.ndarray object of shape self.grid_size = (nt, nx[0], + ... , nx[ns - 1]), where nt is the number of points along the temporal grid and nx = + self.grid_size specifies the number of grid points along the axes in the spatial grid. + The i,j(0), ... , j(ns - 1) element of this array should hold the value of the solution at + the i'th time step and the spatial grid point with index (j(0), ... , j(ns - 1)). + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A two element tuple. The first is a numpy.ndarray object holding the residual on the + spatial and temporal grid. The second should be a float holding the norm of the residual. + """ + raise RuntimeError("Abstract method Physics.residual!") return res, res_norm @@ -112,18 +238,18 @@ def residual(self, Xhist): # ------------------------------------------------------------------------------------------------- class OfflineFOM(Physics): - def __init__(self, cfg, param_name=None): + def __init__(self, cfg, param_name = None): super().__init__(cfg, param_name) self.offline = True assert('offline_fom' in cfg) from ..inputs import InputParser - parser = InputParser(cfg['offline_fom'], name="offline_fom_input") + parser = InputParser(cfg['offline_fom'], name = "offline_fom_input") - self.dim = parser.getInput(['space_dimension'], datatype=int) - self.qdim = parser.getInput(['solution_dimension'], datatype=int) + self.dim = parser.getInput(['space_dimension'], datatype = int) + self.qdim = parser.getInput(['solution_dimension'], datatype = int) - self.grid_size = parser.getInput(['grid_size'], datatype=list) + self.grid_size = parser.getInput(['grid_size'], datatype = list) self.qgrid_size = self.grid_size if (self.qdim > 1): self.qgrid_size = [self.qdim] + self.qgrid_size @@ -134,8 +260,8 @@ def __init__(self, cfg, param_name=None): self.x_grid = None # Assume uniform time stepping for now. - self.nt = parser.getInput(['number_of_timesteps'], datatype=int) - self.dt = parser.getInput(['timestep_size'], datatype=float) + self.nt = parser.getInput(['number_of_timesteps'], datatype = int) + self.dt = parser.getInput(['timestep_size'], datatype = float) self.t_grid = np.linspace(0.0, (self.nt-1) * self.dt, self.nt) return diff --git a/src/lasdi/physics/burgers1d.py b/src/lasdi/physics/burgers1d.py index 3dbbf0b..7a8d162 100644 --- a/src/lasdi/physics/burgers1d.py +++ b/src/lasdi/physics/burgers1d.py @@ -24,80 +24,200 @@ class Burgers1D(Physics): - def __init__(self, cfg, param_name = None): - super().__init__(cfg, param_name) + def __init__(self, cfg : dict, param_name : list[str] = None) -> None: + """ + This is the initializer for the Burgers Physics class. This class essentially acts as a + wrapper around a 1D Burgers solver. - self.qdim = 1 - self.dim = 1 + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- - assert('burgers1d' in cfg) - parser = InputParser(cfg['burgers1d'], name="burgers1d_input") + cfg: A dictionary housing the settings for the Burgers object. This should be the "physics" + sub-dictionary of the configuration file. + + param_name: A list of strings. There should be one list item for each parameter. The i'th + element of this list should be a string housing the name of the i'th parameter. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + Nothing! + """ + + # Call the super class initializer. + super().__init__(cfg, param_name) - self.offline = parser.getInput(['offline_driver'], fallback=False) + # The solution to Burgers' equation is scalar valued, so the qdim is 1. Likewise, since + # there is only one spatial dimension in the 1D burgers example, dim is also 1. + self.qdim : int = 1 + self.dim : int = 1 - self.nt : int = parser.getInput(['number_of_timesteps'], datatype = int) - self.grid_size : list[int] = parser.getInput(['grid_size'], datatype = list) - self.qgrid_size : list[int] = self.grid_size + # Make sure the configuration dictionary is actually for Burgers' equation. + assert('burgers1d' in cfg) + + # Now, get a parser for cfg. + parser : InputParser = InputParser(cfg['burgers1d'], name = "burgers1d_input") + + # Fetch variables from the configuration. + self.offline : bool = parser.getInput(['offline_driver'], fallback = False) # TODO: ??? What does this do ??? + self.nt : int = parser.getInput(['number_of_timesteps'], datatype = int) # number of time steps when solving + self.grid_size : list[int] = parser.getInput(['grid_size'], datatype = list) # number of grid points along each spatial axis + self.qgrid_size : list[int] = self.grid_size + # If there are n spatial dimensions, then the grid needs to have n axes (one for each + # dimension). Make sure this is the case. assert(self.dim == len(self.grid_size)) - self.xmin = parser.getInput(['xmin'], datatype=float) - self.xmax = parser.getInput(['xmax'], datatype=float) - self.dx = (self.xmax - self.xmin) / (self.grid_size[0] - 1) + # Fetch more variables from the + self.xmin = parser.getInput(['xmin'], datatype = float) # Minimum value of the spatial variable in the problem domain + self.xmax = parser.getInput(['xmax'], datatype = float) # Maximum value of the spatial variable in the problem domain + self.dx = (self.xmax - self.xmin) / (self.grid_size[0] - 1) # Spacing between grid points along the spatial axis. assert(self.dx > 0.) - self.tmax = parser.getInput(['simulation_time']) - self.dt = self.tmax / (self.nt - 1) + self.tmax : float = parser.getInput(['simulation_time']) # Final simulation time. We solve form t = 0 to t = tmax + self.dt : float = self.tmax / (self.nt - 1) # step size between successive time steps/the time step we use when solving. - self.x_grid = np.linspace(self.xmin, self.xmax, self.grid_size[0]) - self.t_grid = np.linspace(0, self.tmax, self.nt) + # Set up the spatial, temporal grid. + self.x_grid : np.ndarray = np.linspace(self.xmin, self.xmax, self.grid_size[0]) + self.t_grid : np.ndarray = np.linspace(0, self.tmax, self.nt) - self.maxk = parser.getInput(['maxk'], fallback=10) - self.convergence_threshold = parser.getInput(['convergence_threshold'], fallback=1.e-8) + self.maxk : int = parser.getInput(['maxk'], fallback = 10) # TODO: ??? What is this ??? + self.convergence_threshold : float = parser.getInput(['convergence_threshold'], fallback = 1.e-8) + # Determine which index corresponds to 'a' and 'w' (we pass an array of parameter values, + # we need this information to figure out which element corresponds to which variable). if (self.param_name is not None): if 'a' in self.param_name: self.a_idx = self.param_name.index('a') if 'w' in self.param_name: self.w_idx = self.param_name.index('w') + + # All done! return - def initial_condition(self, param): + def initial_condition(self, param : np.ndarray) -> np.ndarray: + """ + Evaluates the initial condition along the spatial grid. For this class, we use the + following initial condition: + u(x, 0) = a*exp(-x^2 / (2*w^2)) + where a and w are the corresponding parameter values. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param: A 1d numpy.ndarray object with two elements corresponding to the values of the w + and a parameters. self.a_idx and self.w_idx tell us which index corresponds to which + variable. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A 1d numpy.ndarray object of length self.grid_size[0] (the number of grid points along the + spatial axis). + """ + + # Fetch the parameter values. a, w = 1.0, 1.0 if 'a' in self.param_name: a = param[self.a_idx] if 'w' in self.param_name: - w = param[self.w_idx] + w = param[self.w_idx] + # Compute the initial condition and return! return a * np.exp(- self.x_grid ** 2 / 2 / w / w) - def solve(self, param): - u0 = self.initial_condition(param) + def solve(self, param : np.ndarray) -> torch.Tensor: + """ + Solves the 1d burgers equation when the IC uses the parameters in the param array. + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + param: A 1d numpy.ndarray object with two elements corresponding to the values of the w + and a parameters. self.a_idx and self.w_idx tell us which index corresponds to which + variable. + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A 3d torch.Tensor object of shape (1, nt, nx), where nt is the number of points along the + temporal grid and nx is the number along the spatial grid. + """ + + # Fetch the initial condition. + u0 : np.ndarray = self.initial_condition(param) + + # Solve the PDE and then reshape the result to be a 3d tensor with a leading dimension of + # size 1. new_X = solver(u0, self.maxk, self.convergence_threshold, self.nt - 1, self.grid_size[0], self.dt, self.dx) new_X = new_X.reshape(1, self.nt, self.grid_size[0]) + + # All done! return torch.Tensor(new_X) - def export(self): - dict_ = {'t_grid' : self.t_grid, 'x_grid' : self.x_grid, 'dt' : self.dt, 'dx' : self.dx} + def export(self) -> dict: + """ + Returns a dictionary housing self's internal state. You can use this dictionary to + effectively serialize self. + """ + + dict_ : dict = {'t_grid' : self.t_grid, 'x_grid' : self.x_grid, 'dt' : self.dt, 'dx' : self.dx} return dict_ - def residual(self, Xhist): + def residual(self, Xhist : np.ndarray) -> tuple[np.ndarray, float]: + """ + This function computes the PDE residual (difference between the left and right hand side + of Burgers' equation when we substitute in the solution in Xhist). + + + ------------------------------------------------------------------------------------------- + Arguments + ------------------------------------------------------------------------------------------- + + Xhist: A 2d numpy.ndarray object of shape (nt, nx), where nt is the number of points along + the temporal axis and nx is the number of points along the spatial axis. The i,j element of + this array should have the j'th component of the solution at the i'th time step. + + + ------------------------------------------------------------------------------------------- + Returns + ------------------------------------------------------------------------------------------- + + A two element tuple. The first is a numpy.ndarray object of shape (nt - 2, nx - 2) whose + i, j element holds the residual at the i + 1'th temporal grid point and the j + 1'th + spatial grid point. + """ + + # First, approximate the spatial and teporal derivatives. # first axis is time index, and second index is spatial index. dUdx = (Xhist[:, 1:] - Xhist[:, :-1]) / self.dx dUdt = (Xhist[1:, :] - Xhist[:-1, :]) / self.dt - r = dUdt[:, :-1] - Xhist[:-1, :-1] * dUdx[:-1, :] - e = np.linalg.norm(r) + # compute the residual + the norm of the residual. + r : np.ndarray = dUdt[:, :-1] - Xhist[:-1, :-1] * dUdx[:-1, :] + e : float = np.linalg.norm(r) + # All done! return r, e @@ -135,7 +255,7 @@ def jacobian(u, c, idxn1, nx): subdiag_comp = np.ones(nx - 1) subdiag_comp[:-1] = -c * u[1:] data = np.array([diag_comp, subdiag_comp]) - J = spdiags(data, [0, -1], nx - 1, nx - 1, format='csr') + J = spdiags(data, [0, -1], nx - 1, nx - 1, format = 'csr') J[0, -1] = -c * u[0] return J @@ -190,8 +310,8 @@ def main(): import sys parser = argparse.ArgumentParser(description = "", formatter_class=argparse.RawTextHelpFormatter) - parser.add_argument('config_file', metavar='string', type=str, - help='config file to run LasDI workflow.\n') + parser.add_argument('config_file', metavar = 'string', type = str, + help = 'config file to run LasDI workflow.\n') args = parser.parse_args(sys.argv[1:]) print("config file: %s" % args.config_file) @@ -206,28 +326,28 @@ def main(): physics = Burgers1D(config['physics'], param_space.param_name) # read training parameter points - train_param_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'train_param_file'], datatype=str) - train_sol_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'train_sol_file'], datatype=str) + train_param_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'train_param_file'], datatype = str) + train_sol_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'train_sol_file'], datatype = str) with h5py.File(train_param_file, 'r') as f: new_train_params = f['train_params'][...] # generate and write FOM solution new_X = physics.generate_solutions(new_train_params) with h5py.File(train_sol_file, 'w') as f: - f.create_dataset("train_sol", new_X.shape, data=new_X) + f.create_dataset("train_sol", new_X.shape, data = new_X) # check if test parameter points exist - test_param_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'test_param_file'], datatype=str) + test_param_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'test_param_file'], datatype = str) import os.path if (os.path.isfile(test_param_file)): # read test parameter points - test_sol_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'test_sol_file'], datatype=str) + test_sol_file = cfg_parser.getInput(['workflow', 'offline_greedy_sampling', 'test_sol_file'], datatype = str) with h5py.File(test_param_file, 'r') as f: new_test_params = f['test_params'][...] # generate and write FOM solution new_X = physics.generate_solutions(new_test_params) with h5py.File(test_sol_file, 'w') as f: - f.create_dataset("test_sol", new_X.shape, data=new_X) + f.create_dataset("test_sol", new_X.shape, data = new_X) return \ No newline at end of file