first version

data/file1.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:90713e007e25e2b4467711981274dec5f15548666bf7867a30cdf5e189b94b80
+size 7200262

data/file10.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:123c044fc18c7e1b0671261fc081aa7e9a60eab726f35b2758b75f3e70ebe76e
+size 7200262

data/file2.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ebdbdcb3465e49f16073eb828dc18163a12e8b5968cb7f4661e3931c13ac0cea
+size 7200262

data/file3.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6b14cf26774c287dd0fb2351a6fc6ce3a2eae135a57c265a93c03f017d479d3a
+size 7200262

data/file4.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7810f50b9c30a19fbfcb3db77bf838b02df4fc54f6caa63e8dd3ca0f2abba6c1
+size 7200262

data/file5.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34f3485194b57325b0f32554c4f14218389a6d47d5b6edb7e37626fc48d6aae8
+size 7200262

data/file6.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d59084be5aff0eb578c7ff5ee62027ef853d1f5d8d2794d6230d5b706aa0f6aa
+size 7200262

data/file7.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9035c00bb34d17940b033e3bae40097296c493a4852630fbcf963c6f80391f5a
+size 7200262

data/file8.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:955406a5ec5af34373c4ae8573195ca580bfcfbbef13efdcb410fc05d58d66f8
+size 7200262

data/file9.dat.npz

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3ea42398e0062a3280981795ea8fd3a49bd5e104de178decfad9e69e41c0de8c
+size 7200262

model/lstm.ckpt

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1080dbdc37acdb1e9e6a29c140711908d2426b39735617eafbadd49fb5772ef4
+size 7286190

model/lstm.py

@@ -0,0 +1,22 @@
+import pytorch_lightning as pl
+import torch
+import torch.nn as nn
+
+class LSTMModel(pl.LightningModule):
+    def __init__(self, **config):
+        super(LSTMModel, self).__init__()
+        self.save_hyperparameters(config)
+
+        self.lstm = nn.LSTM(input_size=21,hidden_size=512,num_layers=3,proj_size=21,batch_first=True)
+        self.linear = nn.Linear(in_features=21, out_features=7)
+
+    def forward(self, x):
+        outputs = []
+        hidden, cell = None, None
+        for i in range(20):
+            if i == 0:
+                output, (hidden, cell) = self.lstm(x)
+            else:
+                output, (hidden, cell) = self.lstm(output[:, -1, :].unsqueeze(1), (hidden, cell))
+            outputs.append(self.linear(output[:, -1, :]))
+        return torch.stack(outputs, dim=1)

model/tcn.ckpt

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fc975af786670412a2e419978636038452467676e1a9dac59d2ed77033f9f67b
+size 43742454

model/tcn.py

@@ -0,0 +1,40 @@
+import os
+os.environ["KERAS_BACKEND"] = "torch"
+import pytorch_lightning as pl
+import torch
+import torch.nn as nn
+
+from keras.layers import Input, Dense
+from keras.models import Model
+from model.tcn_module import TCN
+
+
+class TCNModel(pl.LightningModule):
+    def __init__(self, **config):
+        super(TCNModel, self).__init__()
+        self.save_hyperparameters(config)
+
+        input_layer = Input(shape=(self.hparams.windows_size, self.hparams.input_size))
+        self.tcn = TCN(input_shape=(self.hparams.windows_size, self.hparams.input_size))(input_layer)
+        self.linear = Dense(7)(self.tcn)
+        self.model = Model(inputs=input_layer, outputs=self.linear)
+
+    def forward(self, x):
+        output = self.model(x)
+        return torch.stack([output], dim=1)
+
+def move_custom_layers_to_device(model, device):
+    for name, module in model.named_children():
+        # 如果是标准层，named_children已经处理了
+        if isinstance(module, nn.Module):
+            continue
+
+        # 对于非标准层，例如包含在列表或字典中的层
+        if isinstance(module, list):
+            for sub_module in module:
+                if isinstance(sub_module, nn.Module):
+                    sub_module.to(device)
+        elif isinstance(module, dict):
+            for sub_module in module.values():
+                if isinstance(sub_module, nn.Module):
+                    sub_module.to(device)

model/tcn_module.py

@@ -0,0 +1,511 @@
+import inspect
+from typing import List
+
+import os
+os.environ["KERAS_BACKEND"] = "torch"
+import keras
+
+# from keras_core import backend as K, Model, Input, optimizers
+# from keras_core import backend as Model, Input, optimizers
+# from keras_core import backend as K
+
+from keras import Model
+from keras import optimizers
+from keras import ops as K
+from keras import config as KK
+
+from keras import layers
+from keras.layers import Input, Layer, Conv1D, Dense, BatchNormalization, LayerNormalization, Activation, SpatialDropout1D, Lambda
+
+
+def is_power_of_two(num: int):
+    return num != 0 and ((num & (num - 1)) == 0)
+
+
+def adjust_dilations(dilations: list):
+    if all([is_power_of_two(i) for i in dilations]):
+        return dilations
+    else:
+        new_dilations = [2 ** i for i in dilations]
+        return new_dilations
+
+
+class ResidualBlock(Layer):
+
+    def __init__(self,
+                 dilation_rate: int,
+                 nb_filters: int,
+                 kernel_size: int,
+                 padding: str,
+                 activation: str = 'relu',
+                 dropout_rate: float = 0,
+                 kernel_initializer: str = 'he_normal',
+                 use_batch_norm: bool = False,
+                 use_layer_norm: bool = False,
+                 use_weight_norm: bool = False,
+                 **kwargs):
+        """Defines the residual block for the WaveNet TCN
+        Args:
+            x: The previous layer in the model
+            training: boolean indicating whether the layer should behave in training mode or in inference mode
+            dilation_rate: The dilation power of 2 we are using for this residual block
+            nb_filters: The number of convolutional filters to use in this block
+            kernel_size: The size of the convolutional kernel
+            padding: The padding used in the convolutional layers, 'same' or 'causal'.
+            activation: The final activation used in o = Activation(x + F(x))
+            dropout_rate: Float between 0 and 1. Fraction of the input units to drop.
+            kernel_initializer: Initializer for the kernel weights matrix (Conv1D).
+            use_batch_norm: Whether to use batch normalization in the residual layers or not.
+            use_layer_norm: Whether to use layer normalization in the residual layers or not.
+            use_weight_norm: Whether to use weight normalization in the residual layers or not.
+            kwargs: Any initializers for Layer class.
+        """
+
+        self.dilation_rate = dilation_rate
+        self.nb_filters = nb_filters
+        self.kernel_size = kernel_size
+        self.padding = padding
+        self.activation = activation
+        self.dropout_rate = dropout_rate
+        self.use_batch_norm = use_batch_norm
+        self.use_layer_norm = use_layer_norm
+        self.use_weight_norm = use_weight_norm
+        self.kernel_initializer = kernel_initializer
+        self.layers = []
+        self.shape_match_conv = None
+        self.res_output_shape = None
+        self.final_activation = None
+
+        super(ResidualBlock, self).__init__(**kwargs)
+
+    def _build_layer(self, layer):
+        """Helper function for building layer
+        Args:
+            layer: Appends layer to internal layer list and builds it based on the current output
+                   shape of ResidualBlocK. Updates current output shape.
+        """
+        self.layers.append(layer)
+        self.layers[-1].build(self.res_output_shape)
+        self.res_output_shape = self.layers[-1].compute_output_shape(self.res_output_shape)
+
+    def build(self, input_shape):
+
+        #with K.name_scope(self.name):  # name scope used to make sure weights get unique names
+        self.layers = []
+        self.res_output_shape = input_shape
+
+        for k in range(2):  # dilated conv block.
+            name = 'conv1D_{}'.format(k)
+            # with K.name_scope(name):  # name scope used to make sure weights get unique names
+            conv = Conv1D(
+                filters=self.nb_filters,
+                kernel_size=self.kernel_size,
+                dilation_rate=self.dilation_rate,
+                padding=self.padding,
+                name=name,
+                kernel_initializer=self.kernel_initializer
+            )
+            if self.use_weight_norm:
+                from tensorflow_addons.layers import WeightNormalization
+                # wrap it. WeightNormalization API is different than BatchNormalization or LayerNormalization.
+                #with K.name_scope('norm_{}'.format(k)):
+                conv = WeightNormalization(conv)
+            self._build_layer(conv)
+
+            #with K.name_scope('norm_{}'.format(k)):
+            if self.use_batch_norm:
+                self._build_layer(BatchNormalization())
+            elif self.use_layer_norm:
+                self._build_layer(LayerNormalization())
+            elif self.use_weight_norm:
+                pass  # done above.
+
+            # with K.name_scope('act_and_dropout_{}'.format(k)):
+            self._build_layer(Activation(self.activation, name='Act_Conv1D_{}'.format(k)))
+            self._build_layer(SpatialDropout1D(rate=self.dropout_rate, name='SDropout_{}'.format(k)))
+
+        if self.nb_filters != input_shape[-1]:
+            # 1x1 conv to match the shapes (channel dimension).
+            name = 'matching_conv1D'
+            #with K.name_scope(name):
+                # make and build this layer separately because it directly uses input_shape.
+                # 1x1 conv.
+            self.shape_match_conv = Conv1D(
+                filters=self.nb_filters,
+                kernel_size=1,
+                padding='same',
+                name=name,
+                kernel_initializer=self.kernel_initializer
+            )
+        else:
+            name = 'matching_identity'
+            self.shape_match_conv = Lambda(lambda x: x, name=name)
+
+        #with K.name_scope(name):
+        self.shape_match_conv.build(input_shape)
+        self.res_output_shape = self.shape_match_conv.compute_output_shape(input_shape)
+
+        self._build_layer(Activation(self.activation, name='Act_Conv_Blocks'))
+        self.final_activation = Activation(self.activation, name='Act_Res_Block')
+        self.final_activation.build(self.res_output_shape)  # probably isn't necessary
+
+        # this is done to force Keras to add the layers in the list to self._layers
+        for layer in self.layers:
+            self.__setattr__(layer.name, layer)
+        self.__setattr__(self.shape_match_conv.name, self.shape_match_conv)
+        self.__setattr__(self.final_activation.name, self.final_activation)
+
+        super(ResidualBlock, self).build(input_shape)  # done to make sure self.built is set True
+
+    def call(self, inputs, training=None, **kwargs):
+        """
+        Returns: A tuple where the first element is the residual model tensor, and the second
+                 is the skip connection tensor.
+        """
+        # https://arxiv.org/pdf/1803.01271.pdf  page 4, Figure 1 (b).
+        # x1: Dilated Conv -> Norm -> Dropout (x2).
+        # x2: Residual (1x1 matching conv - optional).
+        # Output: x1 + x2.
+        # x1 -> connected to skip connections.
+        # x1 + x2 -> connected to the next block.
+        #       input
+        #     x1      x2
+        #   conv1D    1x1 Conv1D (optional)
+        #    ...
+        #   conv1D
+        #    ...
+        #       x1 + x2
+        x1 = inputs
+        for layer in self.layers:
+            training_flag = 'training' in dict(inspect.signature(layer.call).parameters)
+            x1 = layer(x1, training=training) if training_flag else layer(x1)
+        x2 = self.shape_match_conv(inputs)
+        x1_x2 = self.final_activation(layers.add([x2, x1], name='Add_Res'))
+        return [x1_x2, x1]
+
+    def compute_output_shape(self, input_shape):
+        return [self.res_output_shape, self.res_output_shape]
+
+
+class TCN(Layer):
+    """Creates a TCN layer.
+        Input shape:
+            A tensor of shape (batch_size, timesteps, input_dim).
+        Args:
+            nb_filters: The number of filters to use in the convolutional layers. Can be a list.
+            kernel_size: The size of the kernel to use in each convolutional layer.
+            dilations: The list of the dilations. Example is: [1, 2, 4, 8, 16, 32, 64].
+            nb_stacks : The number of stacks of residual blocks to use.
+            padding: The padding to use in the convolutional layers, 'causal' or 'same'.
+            use_skip_connections: Boolean. If we want to add skip connections from input to each residual blocK.
+            return_sequences: Boolean. Whether to return the last output in the output sequence, or the full sequence.
+            activation: The activation used in the residual blocks o = Activation(x + F(x)).
+            dropout_rate: Float between 0 and 1. Fraction of the input units to drop.
+            kernel_initializer: Initializer for the kernel weights matrix (Conv1D).
+            use_batch_norm: Whether to use batch normalization in the residual layers or not.
+            use_layer_norm: Whether to use layer normalization in the residual layers or not.
+            use_weight_norm: Whether to use weight normalization in the residual layers or not.
+            kwargs: Any other arguments for configuring parent class Layer. For example "name=str", Name of the model.
+                    Use unique names when using multiple TCN.
+        Returns:
+            A TCN layer.
+        """
+
+    def __init__(self,
+                 nb_filters=256,
+                 kernel_size=5,
+                 nb_stacks=1,
+                 dilations=(1, 2, 4, 8, 16, 32),
+                 padding='causal',
+                 use_skip_connections=True,
+                 dropout_rate=0.0,
+                 return_sequences=False,
+                 activation='relu',
+                 kernel_initializer='he_normal',
+                 use_batch_norm=False,
+                 use_layer_norm=False,
+                 use_weight_norm=False,
+                 **kwargs):
+        print("nb_filters:", nb_filters, "kernel_size", kernel_size)
+        self.return_sequences = return_sequences
+        self.dropout_rate = dropout_rate
+        self.use_skip_connections = use_skip_connections
+        self.dilations = dilations
+        self.nb_stacks = nb_stacks
+        self.kernel_size = kernel_size
+        self.nb_filters = nb_filters
+        self.activation_name = activation
+        self.padding = padding
+        self.kernel_initializer = kernel_initializer
+        self.use_batch_norm = use_batch_norm
+        self.use_layer_norm = use_layer_norm
+        self.use_weight_norm = use_weight_norm
+        self.skip_connections = []
+        self.residual_blocks = []
+        self.layers_outputs = []
+        self.build_output_shape = None
+        self.slicer_layer = None  # in case return_sequence=False
+        self.output_slice_index = None  # in case return_sequence=False
+        self.padding_same_and_time_dim_unknown = False  # edge case if padding='same' and time_dim = None
+
+        if self.use_batch_norm + self.use_layer_norm + self.use_weight_norm > 1:
+            raise ValueError('Only one normalization can be specified at once.')
+
+        if isinstance(self.nb_filters, list):
+            assert len(self.nb_filters) == len(self.dilations)
+            if len(set(self.nb_filters)) > 1 and self.use_skip_connections:
+                raise ValueError('Skip connections are not compatible '
+                                 'with a list of filters, unless they are all equal.')
+
+        if padding != 'causal' and padding != 'same':
+            raise ValueError("Only 'causal' or 'same' padding are compatible for this layer.")
+
+        # initialize parent class
+        super(TCN, self).__init__(**kwargs)
+
+    @property
+    def receptive_field(self):
+        return 1 + 2 * (self.kernel_size - 1) * self.nb_stacks * sum(self.dilations)
+
+    def build(self, input_shape):
+
+        # member to hold current output shape of the layer for building purposes
+        self.build_output_shape = input_shape
+
+        # list to hold all the member ResidualBlocks
+        self.residual_blocks = []
+        total_num_blocks = self.nb_stacks * len(self.dilations)
+        if not self.use_skip_connections:
+            total_num_blocks += 1  # cheap way to do a false case for below
+
+        for s in range(self.nb_stacks):
+            for i, d in enumerate(self.dilations):
+                res_block_filters = self.nb_filters[i] if isinstance(self.nb_filters, list) else self.nb_filters
+                self.residual_blocks.append(ResidualBlock(dilation_rate=d,
+                                                          nb_filters=res_block_filters,
+                                                          kernel_size=self.kernel_size,
+                                                          padding=self.padding,
+                                                          activation=self.activation_name,
+                                                          dropout_rate=self.dropout_rate,
+                                                          use_batch_norm=self.use_batch_norm,
+                                                          use_layer_norm=self.use_layer_norm,
+                                                          use_weight_norm=self.use_weight_norm,
+                                                          kernel_initializer=self.kernel_initializer,
+                                                          name='residual_block_{}'.format(len(self.residual_blocks))))
+                # build newest residual block
+                self.residual_blocks[-1].build(self.build_output_shape)
+                self.build_output_shape = self.residual_blocks[-1].res_output_shape
+
+        # this is done to force keras to add the layers in the list to self._layers
+        for layer in self.residual_blocks:
+            self.__setattr__(layer.name, layer)
+
+        self.output_slice_index = None
+        if self.padding == 'same':
+            time = self.build_output_shape.as_list()[1]
+            if time is not None:  # if time dimension is defined. e.g. shape = (bs, 500, input_dim).
+                self.output_slice_index = int(self.build_output_shape.as_list()[1] / 2)
+            else:
+                # It will known at call time. c.f. self.call.
+                self.padding_same_and_time_dim_unknown = True
+
+        else:
+            self.output_slice_index = -1  # causal case.
+        self.slicer_layer = Lambda(lambda tt: tt[:, self.output_slice_index, :], name='Slice_Output')
+
+        if type(self.build_output_shape) == tuple:
+            static = list(self.build_output_shape)
+        else:
+            static = self.build_output_shape.as_list()
+        self.slicer_layer.build(static)
+
+    def compute_output_shape(self, input_shape):
+        """
+        Overridden in case keras uses it somewhere... no idea. Just trying to avoid future errors.
+        """
+        if not self.built:
+            self.build(input_shape)
+        if not self.return_sequences:
+            batch_size = self.build_output_shape[0]
+            batch_size = batch_size.value if hasattr(batch_size, 'value') else batch_size
+            nb_filters = self.build_output_shape[-1]
+            return [batch_size, nb_filters]
+        else:
+            # Compatibility tensorflow 1.x
+            return [v.value if hasattr(v, 'value') else v for v in self.build_output_shape]
+
+    def call(self, inputs, training=None, **kwargs):
+        x = inputs
+        self.layers_outputs = [x]
+        self.skip_connections = []
+        for res_block in self.residual_blocks:
+            # try:
+            #     x, skip_out = res_block(x, training=training)
+            # except TypeError:  # compatibility with tensorflow 1.x
+            #     x, skip_out = res_block(K.cast(x, 'float32'), training=training)
+            x, skip_out = res_block(x, training=training)
+
+            self.skip_connections.append(skip_out)
+            self.layers_outputs.append(x)
+
+        if self.use_skip_connections:
+            x = layers.add(self.skip_connections, name='Add_Skip_Connections')
+            self.layers_outputs.append(x)
+
+        if not self.return_sequences:
+            # case: time dimension is unknown. e.g. (bs, None, input_dim).
+            if self.padding_same_and_time_dim_unknown:
+                self.output_slice_index = K.shape(self.layers_outputs[-1])[1] // 2
+            x = self.slicer_layer(x)
+            self.layers_outputs.append(x)
+        return x
+
+    def get_config(self):
+        """
+        Returns the config of a the layer. This is used for saving and loading from a model
+        :return: python dictionary with specs to rebuild layer
+        """
+        config = super(TCN, self).get_config()
+        config['nb_filters'] = self.nb_filters
+        config['kernel_size'] = self.kernel_size
+        config['nb_stacks'] = self.nb_stacks
+        config['dilations'] = self.dilations
+        config['padding'] = self.padding
+        config['use_skip_connections'] = self.use_skip_connections
+        config['dropout_rate'] = self.dropout_rate
+        config['return_sequences'] = self.return_sequences
+        config['activation'] = self.activation_name
+        config['use_batch_norm'] = self.use_batch_norm
+        config['use_layer_norm'] = self.use_layer_norm
+        config['use_weight_norm'] = self.use_weight_norm
+        config['kernel_initializer'] = self.kernel_initializer
+        return config
+
+
+def compiled_tcn(num_feat,  # type: int
+                 num_classes,  # type: int
+                 nb_filters,  # type: int
+                 kernel_size,  # type: int
+                 dilations,  # type: List[int]
+                 nb_stacks,  # type: int
+                 max_len,  # type: int
+                 output_len=1,  # type: int
+                 padding='causal',  # type: str
+                 use_skip_connections=False,  # type: bool
+                 return_sequences=True,
+                 regression=False,  # type: bool
+                 dropout_rate=0.05,  # type: float
+                 name='tcn',  # type: str,
+                 kernel_initializer='he_normal',  # type: str,
+                 activation='relu',  # type:str,
+                 opt='adam',
+                 lr=0.002,
+                 use_batch_norm=False,
+                 use_layer_norm=False,
+                 use_weight_norm=False):
+    # type: (...) -> Model
+    """Creates a compiled TCN model for a given task (i.e. regression or classification).
+    Classification uses a sparse categorical loss. Please input class ids and not one-hot encodings.
+    Args:
+        num_feat: The number of features of your input, i.e. the last dimension of: (batch_size, timesteps, input_dim).
+        num_classes: The size of the final dense layer, how many classes we are predicting.
+        nb_filters: The number of filters to use in the convolutional layers.
+        kernel_size: The size of the kernel to use in each convolutional layer.
+        dilations: The list of the dilations. Example is: [1, 2, 4, 8, 16, 32, 64].
+        nb_stacks : The number of stacks of residual blocks to use.
+        max_len: The maximum sequence length, use None if the sequence length is dynamic.
+        padding: The padding to use in the convolutional layers.
+        use_skip_connections: Boolean. If we want to add skip connections from input to each residual blocK.
+        return_sequences: Boolean. Whether to return the last output in the output sequence, or the full sequence.
+        regression: Whether the output should be continuous or discrete.
+        dropout_rate: Float between 0 and 1. Fraction of the input units to drop.
+        activation: The activation used in the residual blocks o = Activation(x + F(x)).
+        name: Name of the model. Useful when having multiple TCN.
+        kernel_initializer: Initializer for the kernel weights matrix (Conv1D).
+        opt: Optimizer name.
+        lr: Learning rate.
+        use_batch_norm: Whether to use batch normalization in the residual layers or not.
+        use_layer_norm: Whether to use layer normalization in the residual layers or not.
+        use_weight_norm: Whether to use weight normalization in the residual layers or not.
+    Returns:
+        A compiled keras TCN.
+    """
+
+    dilations = adjust_dilations(dilations)
+
+    input_layer = Input(shape=(max_len, num_feat))
+
+    x = TCN(nb_filters, kernel_size, nb_stacks, dilations, padding,
+            use_skip_connections, dropout_rate, return_sequences,
+            activation, kernel_initializer, use_batch_norm, use_layer_norm,
+            use_weight_norm, name=name)(input_layer)
+
+    print('x.shape=', x.shape)
+
+    def get_opt():
+        if opt == 'adam':
+            return optimizers.Adam(lr=lr, clipnorm=1.)
+        elif opt == 'rmsprop':
+            return optimizers.RMSprop(lr=lr, clipnorm=1.)
+        else:
+            raise Exception('Only Adam and RMSProp are available here')
+
+    if not regression:
+        # classification
+        print('asdasfdasfa')
+        x = Dense(num_classes)(x)
+        x = Activation('softmax')(x)
+        output_layer = x
+        model = Model(input_layer, output_layer)
+
+        # https://github.com/keras-team/keras/pull/11373
+        # It's now in Keras@master but still not available with pip.
+        # TODO remove later.
+        def accuracy(y_true, y_pred):
+            # reshape in case it's in shape (num_samples, 1) instead of (num_samples,)
+            if K.ndim(y_true) == K.ndim(y_pred):
+                y_true = K.squeeze(y_true, -1)
+            # convert dense predictions to labels
+            y_pred_labels = K.argmax(y_pred, axis=-1)
+            y_pred_labels = K.cast(y_pred_labels, KK.floatx())
+            return K.cast(K.equal(y_true, y_pred_labels), KK.floatx())
+
+        model.compile(get_opt(), loss='sparse_categorical_crossentropy', metrics=[accuracy])
+    else:
+        # regression
+        x = Dense(output_len)(x)
+        x = Activation('linear')(x)
+        output_layer = x
+        model = Model(input_layer, output_layer)
+        model.compile(get_opt(), loss='mean_squared_error')
+    print('model.x = {}'.format(input_layer.shape))
+    print('model.y = {}'.format(output_layer.shape))
+    return model
+
+
+def tcn_full_summary(model: Model, expand_residual_blocks=True):
+
+    layers = model._layers.copy()  # store existing layers
+    model._layers.clear()  # clear layers
+
+    for i in range(len(layers)):
+        if isinstance(layers[i], TCN):
+            for layer in layers[i]._layers:
+                if not isinstance(layer, ResidualBlock):
+                    if not hasattr(layer, '__iter__'):
+                        model._layers.append(layer)
+                else:
+                    if expand_residual_blocks:
+                        for lyr in layer._layers:
+                            if not hasattr(lyr, '__iter__'):
+                                model._layers.append(lyr)
+                    else:
+                        model._layers.append(layer)
+        else:
+            model._layers.append(layers[i])
+
+    model.summary()  # print summary
+
+    # restore original layers
+    model._layers.clear()
+    [model._layers.append(lyr) for lyr in layers]

pages/inference.py

@@ -0,0 +1,547 @@
+import os
+import time
+
+import streamlit as st
+import torch
+import torch.nn.functional as F
+import random
+from sklearn.preprocessing import MinMaxScaler
+import numpy as np
+import pandas as pd
+
+from model.lstm import LSTMModel
+from model.tcn import TCNModel
+from model.tcn import move_custom_layers_to_device
+from utils.lowlevel import LowLevel
+from utils.highlevel import HighLevel
+from utils.midpoint import MidPoint
+
+from utils.transform import compute_gradient
+
+st.set_page_config(page_title="Inference", page_icon=":chart_with_upwards_trend:", layout="wide", initial_sidebar_state="auto")
+
+def uniform_sampling(data, n_sample):
+    k = len(data) // n_sample
+    return data[::k]
+
+def low_level(option_time, slider_sample_orbit, progress_bar):
+    time.sleep(0.1)
+    low_level_total_start_time = time.time()
+    low_level_30000_start_time = time.time()
+
+    lowlevelhelper = LowLevel(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = lowlevelhelper.initial()
+
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(option_time) - t) / h
+    current_iteration = 0
+
+    original_low_level_data = []
+
+    while t < float(option_time):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+
+        t = t + h
+
+        if jn % 10 == 0:
+            original_low_level_data.append([b, x, y, z, px, py, pz])
+            # Update progress bar
+            progress_percentage = int((current_iteration / total_iterations) * 100)
+            progress_bar.progress(progress_percentage)
+
+        if jn == 300000:
+            low_level_30000_end_time = time.time()
+            low_level_30000_execute_time = low_level_30000_end_time - low_level_30000_start_time
+            low_level_2000_start_time = time.time()
+        jn = jn + 1
+        current_iteration += 1
+
+    progress_bar.progress(100)
+
+    low_level_2000_end_time = time.time()
+    low_level_2000_execute_time = low_level_2000_end_time - low_level_2000_start_time
+    low_level_total_end_time = time.time()
+    low_level_total_execute_time = low_level_total_end_time - low_level_total_start_time
+
+    result = uniform_sampling(np.array(original_low_level_data), n_sample=int(option_time/100))
+
+    return low_level_30000_execute_time, low_level_2000_execute_time, low_level_total_execute_time, result
+
+def high_level(option_time, slider_sample_orbit, progress_bar):
+    time.sleep(0.1)
+    high_level_total_start_time = time.time()
+    high_level_30000_start_time = time.time()
+
+    highlevelhelper = HighLevel(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = highlevelhelper.initial()
+
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(option_time) - t) / h
+    current_iteration = 0
+
+    original_high_level_data = []
+
+    while t < float(option_time):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = highlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = highlevelhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = highlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = highlevelhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+
+        t = t + h
+        vx, vy, vz, vpx, vpy, vpz, e = highlevelhelper.f(x, y, z, px, py, pz, b)
+        en = np.asarray(e).astype(np.float64)
+
+        if jn % 10 == 0:
+            original_high_level_data.append([b, x, y, z, px, py, pz])
+        if jn == 300000:
+            high_level_30000_end_time = time.time()
+            high_level_30000_execute_time = high_level_30000_end_time - high_level_30000_start_time
+            high_level_2000_start_time = time.time()
+        jn = jn + 1
+
+        # Update progress bar
+        progress_percentage = int((current_iteration / total_iterations) * 100)
+        progress_bar.progress(progress_percentage)
+        current_iteration += 1
+
+    progress_bar.progress(100)
+    high_level_2000_end_time = time.time()
+    high_level_2000_execute_time = high_level_2000_end_time - high_level_2000_start_time
+    high_level_total_end_time = time.time()
+    high_level_total_execute_time = high_level_total_end_time - high_level_total_start_time
+
+    result = uniform_sampling(np.array(original_high_level_data), n_sample=int(option_time / 100))
+
+    return high_level_30000_execute_time, high_level_2000_execute_time, high_level_total_execute_time, result
+
+def midpoint(option_time, slider_sample_orbit, progress_bar):
+    time.sleep(0.1)
+    mid_point_total_start_time = time.time()
+    mid_point_30000_start_time = time.time()
+
+    midpointhelper = MidPoint(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = midpointhelper.initial()
+
+    #en0 = np.asarray(e0).astype(np.float64)
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(option_time) - t) / h
+    current_iteration = 0
+
+    original_mid_point_data = []
+
+    while t < float(option_time):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+
+        t = t + h
+
+        if jn % 10 == 0:
+            original_mid_point_data.append([b, x, y, z, px, py, pz])
+        if jn == 300000:
+            mid_point_30000_end_time = time.time()
+            mid_point_30000_execute_time = mid_point_30000_end_time - mid_point_30000_start_time
+            mid_point_2000_start_time = time.time()
+        jn = jn + 1
+
+        # Update progress bar
+        progress_percentage = int((current_iteration / total_iterations) * 100)
+        progress_bar.progress(progress_percentage)
+        current_iteration += 1
+
+    #mid_point_df.to_excel('mid_point_df_output.xlsx', index=False)
+    progress_bar.progress(100)
+    mid_point_2000_end_time = time.time()
+    mid_point_2000_execute_time = mid_point_2000_end_time - mid_point_2000_start_time
+    mid_point_total_end_time = time.time()
+    mid_point_total_execute_time = mid_point_total_end_time - mid_point_total_start_time
+
+    result = uniform_sampling(np.array(original_mid_point_data), n_sample=int(option_time / 100))
+
+    return mid_point_30000_execute_time, mid_point_2000_execute_time, mid_point_total_execute_time, result
+
+def low_level_lstm(slider_sample_orbit, lstm_progress_bar):
+    time.sleep(0.1)
+    total_start_time = time.time()
+
+    lstm_ckpt_file = os.path.join("model", "lstm.ckpt")
+    lstm_model = LSTMModel.load_from_checkpoint(lstm_ckpt_file)
+    lstm_model.to("cpu")
+    lstm_model.eval()
+
+    # Initialize variables for the classical method
+    lowlevelhelper = LowLevel(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = lowlevelhelper.initial()
+
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(30000) - t) / h
+    current_iteration = 0
+
+    original_low_level_data = []
+
+    low_level_start_time = time.time()
+
+    # Perform classical method prediction for the initial segment
+    while t < float(30000):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+        t = t + h
+
+        if jn % 10 == 0:
+            original_low_level_data.append([b, x, y, z, px, py, pz])
+            # Update progress bar
+            progress_percentage = int((current_iteration / total_iterations) * 100)
+            lstm_progress_bar.progress(progress_percentage)
+
+        jn = jn + 1
+        current_iteration += 1
+
+    original_low_level_data = np.array(original_low_level_data)
+    low_level_end_time = time.time()
+    low_level_data = original_low_level_data.copy()
+    low_level_data = uniform_sampling(low_level_data, n_sample=300)
+    scaler = MinMaxScaler()
+    low_level_data = scaler.fit_transform(low_level_data)
+    low_level_data = torch.tensor(np.stack(low_level_data)).float()
+    low_level_data = torch.stack([compute_gradient(i, degree=2) for i in low_level_data]).unsqueeze(0)
+
+    lstm_start_time = time.time()
+    with torch.no_grad():
+        lstm_preds = lstm_model(low_level_data[:, 100:300, :])
+        lstm_innv_preds = scaler.inverse_transform(lstm_preds.squeeze().cpu().numpy())
+
+    original_low_level_data = uniform_sampling(original_low_level_data, n_sample=300)
+
+    lstm_end_time = time.time()
+    lstm_progress_bar.progress(100)
+
+    combined_preds = np.concatenate([original_low_level_data, lstm_innv_preds], axis=0)
+
+    lstm_total_time = lstm_end_time - lstm_start_time
+    low_level_total_time = low_level_end_time - low_level_start_time
+
+    total_end_time = time.time()
+    total_time = total_end_time - total_start_time
+
+    return low_level_total_time, lstm_total_time, total_time, combined_preds
+
+def mid_point_lstm(slider_sample_orbit, lstm_progress_bar):
+    time.sleep(0.1)
+    total_start_time = time.time()
+
+    lstm_ckpt_file = os.path.join("model", "lstm.ckpt")
+    lstm_model = LSTMModel.load_from_checkpoint(lstm_ckpt_file)
+    lstm_model.to("cpu")
+    lstm_model.eval()
+
+    midpointhelper = MidPoint(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = midpointhelper.initial()
+
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(30000) - t) / h
+    current_iteration = 0
+
+    original_mid_point_data = []
+
+    mid_point_start_time = time.time()
+
+    while t < float(30000):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+
+        t = t + h
+
+        if jn % 10 == 0:
+            original_mid_point_data.append([b, x, y, z, px, py, pz])
+            # Update progress bar
+            progress_percentage = int((current_iteration / total_iterations) * 100)
+            lstm_progress_bar.progress(progress_percentage)
+        jn = jn + 1
+        current_iteration += 1
+
+    original_mid_point_data = np.array(original_mid_point_data)
+    mid_point_end_time = time.time()
+    mid_point_data = original_mid_point_data.copy()
+    mid_point_data = uniform_sampling(mid_point_data, n_sample=300)
+    scaler = MinMaxScaler()
+    mid_point_data = scaler.fit_transform(mid_point_data)
+    mid_point_data = torch.tensor(np.stack(mid_point_data)).float()
+    mid_point_data = torch.stack([compute_gradient(i, degree=2) for i in mid_point_data]).unsqueeze(0)
+
+    lstm_start_time = time.time()
+    with torch.no_grad():
+        lstm_preds = lstm_model(mid_point_data[:, 100:300, :])
+        lstm_innv_preds = scaler.inverse_transform(lstm_preds.squeeze().cpu().numpy())
+
+    original_mid_point_data = uniform_sampling(original_mid_point_data, n_sample=300)
+
+    lstm_end_time = time.time()
+    lstm_progress_bar.progress(100)
+
+    combined_preds = np.concatenate([original_mid_point_data, lstm_innv_preds], axis=0)
+
+    lstm_total_time = lstm_end_time - lstm_start_time
+    mid_point_total_time = mid_point_end_time - mid_point_start_time
+
+    total_end_time = time.time()
+    total_time = total_end_time - total_start_time
+
+    return mid_point_total_time, lstm_total_time, total_time, combined_preds
+
+def low_level_tcn(slider_sample_orbit, tcn_progress_bar):
+    time.sleep(0.1)
+    total_start_time = time.time()
+
+    tcn_ckpt_file = os.path.join("model", "tcn.ckpt")
+    tcn_model = TCNModel.load_from_checkpoint(tcn_ckpt_file)
+    move_custom_layers_to_device(tcn_model, "cpu")
+    tcn_model.eval()
+
+    # Initialize variables for the classical method
+    lowlevelhelper = LowLevel(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = lowlevelhelper.initial()
+
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(30000) - t) / h
+    current_iteration = 0
+
+    original_low_level_data = []
+
+    low_level_start_time = time.time()
+
+    # Perform classical method prediction for the initial segment
+    while t < float(30000):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = lowlevelhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza)
+
+        t = t + h
+
+        if jn % 10 == 0:
+            original_low_level_data.append([b, x, y, z, px, py, pz])
+            # Update progress bar
+            progress_percentage = int((current_iteration / total_iterations) * 100)
+            tcn_progress_bar.progress(progress_percentage)
+
+        jn = jn + 1
+        current_iteration += 1
+
+    original_low_level_data = np.array(original_low_level_data)
+    low_level_end_time = time.time()
+    low_level_data = original_low_level_data.copy()
+    low_level_data = uniform_sampling(low_level_data, n_sample=300)
+    scaler = MinMaxScaler()
+    low_level_data = scaler.fit_transform(low_level_data)
+    low_level_data = torch.tensor(np.stack(low_level_data)).float()
+    low_level_data = torch.stack([compute_gradient(i, degree=2) for i in low_level_data]).unsqueeze(0)
+
+    tcn_start_time = time.time()
+    with torch.no_grad():
+        tcn_preds = None
+        for i in range(20):
+            if i == 0:
+                tcn_preds = tcn_model(low_level_data[:, :300, :])
+            else:
+                gd_y_hat = compute_gradient(tcn_preds[:, :i, :], degree=2).to('cpu')
+                output = tcn_model(torch.cat([low_level_data[:, i:300, :], gd_y_hat], dim=1).to('cpu'))
+                tcn_preds = torch.cat([tcn_preds, output], dim=1)
+        tcn_innv_preds = scaler.inverse_transform(tcn_preds.squeeze().cpu().numpy())
+
+    original_low_level_data = uniform_sampling(original_low_level_data, n_sample=300)
+
+    tcn_end_time = time.time()
+    tcn_progress_bar.progress(100)
+
+    combined_preds = np.concatenate([original_low_level_data, tcn_innv_preds], axis=0)
+
+    tcn_total_time = tcn_end_time - tcn_start_time
+    low_level_total_time = low_level_end_time - low_level_start_time
+
+    total_end_time = time.time()
+    total_time = total_end_time - total_start_time
+
+    return low_level_total_time, tcn_total_time, total_time, combined_preds
+
+def mid_point_tcn(slider_sample_orbit, tcn_progress_bar):
+    time.sleep(0.1)
+    total_start_time = time.time()
+
+    tcn_ckpt_file = os.path.join("model", "tcn.ckpt")
+    tcn_model = TCNModel.load_from_checkpoint(tcn_ckpt_file)
+    move_custom_layers_to_device(tcn_model, "cpu")
+    tcn_model.eval()
+
+    # Initialize variables for the classical method
+    midpointhelper = MidPoint(j=slider_sample_orbit)
+    j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = midpointhelper.initial()
+
+    a1 = 1 / (2 - 2 ** (1 / 3))
+    a2 = 1 - 2 * a1
+    jn = 0
+    t = 0.1
+
+    # Calculate the total number of iterations for the progress bar update
+    total_iterations = (float(30000) - t) / h
+    current_iteration = 0
+
+    original_mid_point_data = []
+
+    mid_point_start_time = time.time()
+
+    # Perform classical method prediction for the initial segment
+    while t < float(30000):
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
+        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza)
+
+        t = t + h
+
+        if jn % 10 == 0:
+            original_mid_point_data.append([b, x, y, z, px, py, pz])
+            # Update progress bar
+            progress_percentage = int((current_iteration / total_iterations) * 100)
+            tcn_progress_bar.progress(progress_percentage)
+        jn = jn + 1
+        current_iteration += 1
+
+    original_mid_point_data = np.array(original_mid_point_data)
+    mid_point_end_time = time.time()
+    mid_point_data = original_mid_point_data.copy()
+    mid_point_data = uniform_sampling(mid_point_data, n_sample=300)
+    scaler = MinMaxScaler()
+    mid_point_data = scaler.fit_transform(mid_point_data)
+    mid_point_data = torch.tensor(np.stack(mid_point_data)).float()
+    mid_point_data = torch.stack([compute_gradient(i, degree=2) for i in mid_point_data]).unsqueeze(0)
+
+    tcn_start_time = time.time()
+    with torch.no_grad():
+        tcn_preds = None
+        for i in range(20):
+            if i == 0:
+                tcn_preds = tcn_model(mid_point_data[:, :300, :])
+            else:
+                gd_y_hat = compute_gradient(tcn_preds[:, :i, :], degree=2).to('cpu')
+                output = tcn_model(torch.cat([mid_point_data[:, i:300, :], gd_y_hat], dim=1).to('cpu'))
+                tcn_preds = torch.cat([tcn_preds, output], dim=1)
+        tcn_innv_preds = scaler.inverse_transform(tcn_preds.squeeze().cpu().numpy())
+
+    original_mid_point_data = uniform_sampling(original_mid_point_data, n_sample=300)
+
+    tcn_end_time = time.time()
+    tcn_progress_bar.progress(100)
+
+    combined_preds = np.concatenate([original_mid_point_data, tcn_innv_preds], axis=0)
+
+    tcn_total_time = tcn_end_time - tcn_start_time
+    mid_point_total_time = mid_point_end_time - mid_point_start_time
+
+    total_end_time = time.time()
+    total_time = total_end_time - total_start_time
+
+    return mid_point_total_time, tcn_total_time, total_time, combined_preds
+
+container = st.container()
+container1, container2 = st.columns(2)
+plot_container = st.container()
+
+with st.sidebar:
+    slider_sample_orbit = st.slider('Orbit Sample ID', 1, 10, 1)
+    option_time = 32000
+    st.write(f'Total Time Step: {option_time}')
+    options_method = st.multiselect(
+        'Compared Methods',
+        ['Low-Level', 'High-Level', 'Midpoint', 'Low-Level with LSTM', 'Low-Level with TCN', 'Midpoint with LSTM', 'Midpoint with TCN'],
+        ['Low-Level'])
+    btn_go = st.button("Go", type="primary", use_container_width=True)
+
+if btn_go:
+    if 'Low-Level' in options_method:
+        with container1:
+            st.write('Low Level Progress Bar')
+            low_level_progress_bar = st.progress(0)
+        low_level_30000_time, low_level_2000_time, low_level_total_time, low_level_result = low_level(option_time, slider_sample_orbit, low_level_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"Low Level", '30000 Time Steps (s)': [low_level_30000_time], '2000 Time Steps (s)': [low_level_2000_time], 'Total Time (s)': [low_level_total_time]}))
+    if 'High-Level' in options_method:
+        with container1:
+            st.write('High Level Progress Bar')
+            high_level_progress_bar = st.progress(0)
+        high_level_30000_time, high_level_2000_time, high_level_total_time, high_level_result = high_level(option_time, slider_sample_orbit, high_level_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"High Level", '30000 Time Steps (s)': [high_level_30000_time], '2000 Time Steps (s)': [high_level_2000_time], 'Total Time (s)': [high_level_total_time]}))
+    if 'Midpoint' in options_method:
+        with container1:
+            st.write('Midpoint Progress Bar')
+            mid_point_progress_bar = st.progress(0)
+        mid_point_30000_time, mid_point_2000_time, mid_point_total_time, mid_point_result = midpoint(option_time, slider_sample_orbit, mid_point_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"Midpoint", '30000 Time Steps (s)': [mid_point_30000_time], '2000 Time Steps (s)': [mid_point_2000_time], 'Total Time (s)': [mid_point_total_time]}))
+    if 'Low-Level with LSTM' in options_method:
+        with container1:
+            st.write('Low Level LSTM Progress Bar')
+            low_level_lstm_progress_bar = st.progress(0)
+        lstm_30000_time, lstm_2000_time, lstm_total_time, lstm_result = low_level_lstm(slider_sample_orbit, low_level_lstm_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"Low Level + LSTM", '30000 Time Steps (s)': [lstm_30000_time], '2000 Time Steps (s)': [lstm_2000_time], 'Total Time (s)': [lstm_total_time]}))
+    if 'Low-Level with TCN' in options_method:
+        with container1:
+            st.write('Low Level TCN Progress Bar')
+            low_level_tcn_progress_bar = st.progress(0)
+        tcn_30000_time, tcn_2000_time, tcn_total_time, tcn_result = low_level_tcn(slider_sample_orbit, low_level_tcn_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"Low Level + TCN", '30000 Time Steps (s)': [tcn_30000_time], '2000 Time Steps (s)': [tcn_2000_time], 'Total Time (s)': [tcn_total_time]}))
+    if 'Midpoint with LSTM' in options_method:
+        with container1:
+            st.write('Midpoint LSTM Progress Bar')
+            mid_point_lstm_progress_bar = st.progress(0)
+        md_lstm_30000_time, md_lstm_2000_time, md_lstm_total_time, md_lstm_result = mid_point_lstm(slider_sample_orbit, mid_point_lstm_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"Midpoint + LSTM", '30000 Time Steps (s)': [md_lstm_30000_time], '2000 Time Steps (s)': [md_lstm_2000_time], 'Total Time (s)': [md_lstm_total_time]}))
+    if 'Midpoint with TCN' in options_method:
+        with container1:
+            st.write('Midpoint TCN Progress Bar')
+            mid_point_tcn_progress_bar = st.progress(0)
+        md_tcn_30000_time, md_tcn_2000_time, md_tcn_total_time, md_tcn_result = mid_point_tcn(slider_sample_orbit, mid_point_tcn_progress_bar)
+        with container2:
+            st.table(pd.DataFrame({'Model':"Midpoint + TCN", '30000 Time Steps (s)': [md_tcn_30000_time], '2000 Time Steps (s)': [md_tcn_2000_time], 'Total Time (s)': [md_tcn_total_time]}))
+            

prediction.py

@@ -0,0 +1,196 @@
+import streamlit as st
+import os
+import numpy as np
+import pandas as pd
+from sklearn.preprocessing import MinMaxScaler
+import torch
+import time
+
+from utils.transform import compute_gradient
+from model.lstm import LSTMModel
+from model.tcn import TCNModel
+from model.tcn import move_custom_layers_to_device
+from utils.metrics import calculate_metrics
+
+each_feature_name = ["q_1","q_2","q_3","p_1","p_2","p_3"]
+
+def uniform_sampling(data, n_sample):
+    k = len(data) // n_sample
+    return data[::k]
+
+st.set_page_config(page_title="Prediction", page_icon=":chart_with_upwards_trend:", layout="wide", initial_sidebar_state="auto")
+
+#st.title("Prediction")
+
+with st.sidebar:
+    slider_predict_step = st.slider('Predicted Step', 0, 20, 20)
+
+    number_input_sample_id = st.number_input("Select Sample ID 1~10", value=1, placeholder="Type a number...", min_value=1, max_value=10, step=1)
+
+    squences_start_idx = st.slider('Squences Start Index', 0, 700 - slider_predict_step, 0)
+
+    st.subheader("Model Configuration")
+    st.write("LSTM Window Size: ", 200)
+    st.write("TCN Window Size: ", 300)
+    st.write("Predicted Step: ", slider_predict_step)
+    st.write("Feature Augmentation: ", "Second-order derivative")
+
+file_path = os.path.join("data", "file"+str(number_input_sample_id)+".dat.npz")
+data = pd.DataFrame(np.load(file_path)['data'])
+scaler = MinMaxScaler()
+uniform_data = uniform_sampling(data, n_sample=1000).sort_index().values[:, 1:8]
+normal_uniform_data = scaler.fit_transform(uniform_data)
+data_sequences = torch.tensor(np.stack(normal_uniform_data)).float()
+original_data_sequences = torch.tensor(np.stack(uniform_data)).float()
+selected_data = data_sequences[squences_start_idx:squences_start_idx+300+slider_predict_step]
+original_selected_data = original_data_sequences[squences_start_idx:squences_start_idx+300+slider_predict_step]
+input_data = torch.stack([compute_gradient(i, degree=2) for i in selected_data]).unsqueeze(0)
+
+with st.sidebar:
+    st.subheader("Data Configuration")
+    st.write("Sample ID: ", number_input_sample_id)
+    #st.write("Origianl Shape: ", data_sequences.shape)
+    st.write("Squences Start Index: ", squences_start_idx)
+    #st.write("Selected Shape: ", selected_data.shape)
+    #st.write("Input Shape: ", input_data.shape)
+
+# st.write(selected_data[0])
+# st.write(input_data[0][0])
+
+#################################################
+## LSTM GPU Inference
+#################################################
+lstm_ckpt_file = os.path.join("model", "lstm.ckpt")
+lstm_model = LSTMModel.load_from_checkpoint(lstm_ckpt_file)
+lstm_model.eval()
+lstm_start_time = time.time()
+with torch.no_grad():
+    lstm_preds = lstm_model(input_data[:, 100:300, :].cuda())
+    lstm_end_time = time.time()
+    lstm_innv_preds = scaler.inverse_transform(lstm_preds.squeeze().cpu().numpy())
+    lstm_normal_preds = lstm_preds.squeeze().cpu().numpy()
+
+    #lstm_model.to_onnx("model/lstm.onnx", torch.randn((1, 200, 21)), export_params=True)
+
+del lstm_model
+
+#################################################
+## LSTM CPU Inference
+#################################################
+lstm_cpu_ckpt_file = os.path.join("model", "lstm.ckpt")
+lstm_cpu_model = LSTMModel.load_from_checkpoint(lstm_ckpt_file)
+lstm_cpu_model.to("cpu")
+lstm_cpu_model.eval()
+lstm_cpu_start_time = time.time()
+with torch.no_grad():
+    lstm_cpu_preds = lstm_cpu_model(input_data[:, 100:300, :])
+    lstm_cpu_end_time = time.time()
+
+del lstm_cpu_model
+
+#################################################
+## TCN GPU Inference
+#################################################
+tcn_ckpt_file = os.path.join("model", "tcn.ckpt")
+tcn_model = TCNModel.load_from_checkpoint(tcn_ckpt_file)
+tcn_model.eval()
+tcn_start_time = time.time()
+with torch.no_grad():
+    input_data_cuda = input_data[:,:300,:].cuda()
+    y_hat = tcn_model(input_data_cuda)
+    for i in range(1, slider_predict_step):
+        gd_y_hat = compute_gradient(y_hat[:, :i, :], degree=2).cuda()
+        output = tcn_model(torch.cat([input_data[:, i:300, :].cuda(), gd_y_hat], dim=1)).cuda()
+        y_hat = torch.cat([y_hat, output], dim=1)
+    tcn_end_time = time.time()
+    tcn_preds = y_hat
+    tcn_innv_preds = scaler.inverse_transform(tcn_preds.squeeze().cpu().numpy())
+    tcn_normal_preds = tcn_preds.squeeze().cpu().numpy()
+
+    #tcn_model.to_onnx("model/tcn.onnx", torch.randn((1, 300, 21)), export_params=True)
+
+del tcn_model
+del y_hat, gd_y_hat, output
+
+#################################################
+## TCN CPU Inference
+#################################################
+input_data_cpu = input_data.to("cpu")
+tcn_cpu_ckpt_file = os.path.join("model", "tcn.ckpt")
+tcn_cpu_model = TCNModel.load_from_checkpoint(tcn_cpu_ckpt_file)
+move_custom_layers_to_device(tcn_cpu_model, "cpu")
+tcn_cpu_model.eval()
+tcn_cpu_start_time = time.time()
+with torch.no_grad():
+    y_hat = None
+    for i in range(slider_predict_step):
+        if i == 0:
+            y_hat = tcn_cpu_model(input_data_cpu[:,:300,:])
+        else:
+            gd_y_hat = compute_gradient(y_hat[:, :i, :], degree=2).to('cpu')
+            output = tcn_cpu_model(torch.concatenate([input_data_cpu[:, i:300, :], gd_y_hat], dim=1).to('cpu'))
+            y_hat = torch.concatenate([y_hat, output], dim=1)
+    tcn_cpu_preds = y_hat
+    tcn_cpu_end_time = time.time()
+
+del tcn_cpu_model
+
+st.subheader("Normalized Prediction")
+
+i = 1
+for each_col in st.columns(6):
+    with each_col:
+        raw_data = selected_data[:, i]
+        lstm_data = [np.nan] * 300 + lstm_normal_preds[:slider_predict_step, :][:, i].tolist()
+        tcn_data = [np.nan] * 300 + tcn_normal_preds[:, i].tolist()
+        st.markdown(f"<div style='text-align: center'>{each_feature_name[i-1]}</div>", unsafe_allow_html=True)
+        #st.write(np.array(raw_data).shape, np.array(lstm_data).shape, np.array(tcn_data).shape)
+        st.line_chart(pd.DataFrame({"Original": raw_data, "LSTM": lstm_data, "TCN": tcn_data}),
+                      color=["#EE4035", "#0077BB", "#7BC043"])
+    i += 1
+
+# with st.sidebar:
+#     st.write("Predicted Shape: ", lstm_preds.shape)
+
+st.subheader("Inverse Normalized Prediction")
+
+i = 1
+for each_col in st.columns(6):
+    with each_col:
+        raw_data = original_selected_data[:, i]
+        lstm_data = [np.nan] * 300 + lstm_innv_preds[:slider_predict_step, :][:, i].tolist()
+        tcn_data = [np.nan] * 300 + tcn_innv_preds[:, i].tolist()
+        st.markdown(f"<div style='text-align: center'>{each_feature_name[i - 1]}</div>", unsafe_allow_html=True)
+        st.line_chart(pd.DataFrame({"Original": raw_data, "LSTM": lstm_data, "TCN": tcn_data}),
+                      color=["#EE4035", "#0077BB", "#7BC043"])
+    i += 1
+
+LSTM_SMAPE, LSTM_MSE, LSTM_RMSE, LSTM_MAE, LSTM_R2, LSTM_PSD = calculate_metrics(selected_data[300:300+slider_predict_step, :].cpu().numpy(), lstm_normal_preds[:slider_predict_step, :])
+TCN_SMAPE, TCN_MSE, TCN_RMSE, TCN_MAE, TCN_R2, TCN_PSD = calculate_metrics(selected_data[300:300+slider_predict_step, :].cpu().numpy(), tcn_normal_preds)
+
+results_df = pd.DataFrame({
+    "Model": ["LSTM", "TCN"],
+    "SMAPE": [LSTM_SMAPE, TCN_SMAPE],
+    "MSE": [LSTM_MSE, TCN_MSE],
+    "RMSE": [LSTM_RMSE, TCN_RMSE],
+    "MAE": [LSTM_MAE, TCN_MAE],
+    "R2": [LSTM_R2, TCN_R2],
+    "PSD": [LSTM_PSD, TCN_PSD]
+})
+
+time_df = pd.DataFrame({
+    "Model": ["LSTM-GPU", "TCN-GPU", "LSTM-CPU", "TCN-CPU"],
+    "Time(ms)": [(lstm_end_time - lstm_start_time)*1000,
+                 (tcn_end_time - tcn_start_time)*1000,
+                 (lstm_cpu_end_time - lstm_cpu_start_time)*1000,
+                 (tcn_cpu_end_time - tcn_cpu_start_time)*1000]
+})
+
+col1, col2 = st.columns(2)
+with col1:
+    st.subheader("Evaluation Metrics")
+    st.write(results_df)
+
+with col2:
+    st.subheader("Prediction Time")
+    st.write(time_df)

requirements.txt

@@ -0,0 +1,7 @@
+pandas==1.5.3
+scikit-learn==1.3.2
+torch==2.1.1
+streamlit==1.32.2
+keras==3.0.0
+torchvision==0.16.1
+pytorch-lightning==2.1.2

utils/highlevel.py

@@ -0,0 +1,160 @@
+import sympy as sp
+
+class HighLevel():
+
+    def __init__(self, j):
+        self.j = j
+
+    def initial(self):
+        j = self.j
+
+        #init parameters
+        h, b, n, u = 0.1, None, None, None
+        x, y, z, px, py, pz = None, 0.1, 0.001, 0.01, None, 0.0001
+        xa, ya, za, pxa, pya, pza = None, None, None, None, None, None
+
+        if j == 1:
+            b = 5.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 10.0
+            py = 0.5
+        elif j == 2:
+            b = 3.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 8.3
+            py = 0.6
+        elif j == 3:
+            b = 3.0 / 2
+            x = 12.0
+            py = 0.4
+        elif j == 4:
+            b = 7.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 15.0
+            py = 0.35
+        elif j == 5:
+            b = 1.0
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 18.0
+            py = 0.3
+        elif j == 6:
+            b = 3.0 / 5
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 20.0
+            py = 0.25
+        elif j == 7:
+            b = 5.0 / 7
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 22.0
+            py = 0.22
+        elif j == 8:
+            b = 2.0
+            x = 26.0
+            py = 0.2
+        elif j == 9:
+            b = 0.5
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 30.0
+            y = 0.5
+            z = 0.1
+            pz = 0.01
+        elif j == 10:
+            b = 5.0
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 35.0
+            y = 2.0
+            z = 0.1
+            pz = 0.03
+            py = 0.15
+
+        xa, ya, za, pxa, pya, pza = x, y, z, px, py, pz
+        return j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza
+
+    def f(self, x, y, z, px, py, pz, b):
+        x_val, y_val, z_val, px_val, py_val, pz_val, b_val = x, y, z, px, py, pz, b
+        x, y, z, px, py, pz, b = sp.symbols('x y z px py pz b')
+
+        c = 1.0
+
+        u = 1 / (1 / b + b + 2)
+        ht = px ** 2 / 2 + py ** 2 / 2 + pz ** 2 / 2
+        hv = -1 / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        h1pn = 1 / (2 * x ** 2 + 2 * y ** 2 + 2 * z ** 2) - (((u + 3) * (px ** 2 + py ** 2 + pz ** 2)) / 2 + (u * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) ** 2) / 2) / (
+                           x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + ((3 * u) / 8 - 1 / 8) * (
+                           px ** 2 + py ** 2 + pz ** 2) ** 2
+
+        e = ht + hv + h1pn
+
+        de_dx = sp.diff(e, x)
+        de_dy = sp.diff(e, y)
+        de_dz = sp.diff(e, z)
+        de_dpx = sp.diff(e, px)
+        de_dpy = sp.diff(e, py)
+        de_dpz = sp.diff(e, pz)
+
+        de_dx_val = de_dx.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+        de_dy_val = de_dy.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+        de_dz_val = de_dz.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+        de_dpx_val = de_dpx.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+        de_dpy_val = de_dpy.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+        de_dpz_val = de_dpz.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+
+        e_val = e.subs({x: x_val, y: y_val, z: z_val, px: px_val, py: py_val, pz: pz_val, b: b_val})
+
+        return de_dx_val, de_dy_val, de_dz_val, de_dpx_val, de_dpy_val, de_dpz_val, e_val
+
+    def rejust(self, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza):
+
+        x = (x + xa) / 2
+        y = (y + ya) / 2
+        z = (z + za) / 2
+
+        px = (px + pxa) / 2
+        py = (py + pya) / 2
+        pz = (pz + pza) / 2
+        xa = x
+        ya = y
+        za = z
+        pxa = px
+        pya = py
+        pza = pz
+
+        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza
+
+    def symplectic(self, h, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza, b):
+
+        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
+        x = x + h / 2 * vpx
+        y = y + h / 2 * vpy
+        z = z + h / 2 * vpz
+        pxa = pxa - h / 2 * vxa
+        pya = pya - h / 2 * vya
+        pza = pza - h / 2 * vza
+
+        vx, vy, vz, vpxa, vpya, vpza, e = self.f(x, y, z, pxa, pya, pza, b)
+        xa = xa + h * vpxa
+        ya = ya + h * vpya
+        za = za + h * vpza
+        px = px - h * vx
+        py = py - h * vy
+        pz = pz - h * vz
+
+        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
+        x = x + h / 2 * vpx
+        y = y + h / 2 * vpy
+        z = z + h / 2 * vpz
+        pxa = pxa - h / 2 * vxa
+        pya = pya - h / 2 * vya
+        pza = pza - h / 2 * vza
+
+        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza

utils/lowlevel.py

@@ -0,0 +1,158 @@
+import pandas as pd
+import numpy as np
+
+class LowLevel():
+
+    def __init__(self, j):
+        self.j = j
+
+    def initial(self):
+        j = self.j
+
+        #init parameters
+        h, b, n, u = 0.1, None, None, None
+        x, y, z, px, py, pz = None, 0.1, 0.001, 0.01, None, 0.0001
+        xa, ya, za, pxa, pya, pza = None, None, None, None, None, None
+
+        if j == 1:
+            b = 5.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 10.0
+            py = 0.5
+        elif j == 2:
+            b = 3.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 8.3
+            py = 0.6
+        elif j == 3:
+            b = 3.0 / 2
+            x = 12.0
+            py = 0.4
+        elif j == 4:
+            b = 7.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 15.0
+            py = 0.35
+        elif j == 5:
+            b = 1.0
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 18.0
+            py = 0.3
+        elif j == 6:
+            b = 3.0 / 5
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 20.0
+            py = 0.25
+        elif j == 7:
+            b = 5.0 / 7
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 22.0
+            py = 0.22
+        elif j == 8:
+            b = 2.0
+            x = 26.0
+            py = 0.2
+        elif j == 9:
+            b = 0.5
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 30.0
+            y = 0.5
+            z = 0.1
+            pz = 0.01
+        elif j == 10:
+            b = 5.0
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 35.0
+            y = 2.0
+            z = 0.1
+            pz = 0.03
+            py = 0.15
+
+        xa, ya, za, pxa, pya, pza = x, y, z, px, py, pz
+        return j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza
+
+    def f(self, x, y, z, px, py, pz, b):
+        n = b / (1 + b)**2
+        u = 1 / (1 / b + b + 2)
+        ht = px**2 / 2 + py**2 / 2 + pz**2 / 2
+        hv = -1 / (x**2 + y**2 + z**2)**(1/2)
+        h1pn = 1/(2*x**2 + 2*y**2 + 2*z**2) - (((u + 3)*(px**2 + py**2 +pz**2))/2 + (u*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/ (x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))**2)/2)/(x**2 + y**2 + z**2)**(1/2) + ((3*u)/8 -1/8)*(px**2 + py**2 + pz**2)**2
+
+        e = ht + hv + h1pn
+
+        vnpx=px
+        v1pnpx=4*px*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (px*(n + 3) + (n*x*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
+        vpx=vnpx+v1pnpx
+
+        vnpy=py
+        v1pnpy=4*py*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (py*(n + 3) + (n*y*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
+        vpy=vnpy+v1pnpy
+
+        vnpz=pz
+        v1pnpz=4*pz*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (pz*(n + 3) + (n*z*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
+        vpz=vnpz+v1pnpz
+
+        vnx=x/(x**2 + y**2 + z**2)**(3/2)
+        v1pnx=(x*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*x)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((px*x**2)/(x**2 + y**2 + z**2)**(3/2) -px/(x**2 + y**2 + z**2)**(1/2) + (py*x*y)/(x**2 + y**2 + z**2)**(3/2) + (pz*x*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
+        vx=vnx+v1pnx
+
+        vny=y/(x**2 + y**2 + z**2)**(3/2)
+        v1pny=(y*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*y)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((py*y**2)/(x**2 + y**2 + z**2)**(3/2) -py/(x**2 + y**2 + z**2)**(1/2) + (px*x*y)/(x**2 + y**2 + z**2)**(3/2) + (pz*y*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
+        vy=vny+v1pny
+
+        vnz=z/(x**2 + y**2 + z**2)**(3/2)
+        v1pnz=(z*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*z)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((pz*z**2)/(x**2 + y**2 + z**2)**(3/2) -pz/(x**2 + y**2 + z**2)**(1/2) + (px*x*z)/(x**2 + y**2 + z**2)**(3/2) + (py*y*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
+        vz=vnz+v1pnz
+
+        return vx,vy,vz,vpx,vpy,vpz,e
+
+    def rejust(self, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza):
+        x = (x + xa) / 2
+        y = (y + ya) / 2
+        z = (z + za) / 2
+
+        px = (px + pxa) / 2
+        py = (py + pya) / 2
+        pz = (pz + pza) / 2
+        xa = x
+        ya = y
+        za = z
+        pxa = px
+        pya = py
+        pza = pz
+        return x,y,z,px,py,pz,xa,ya,za,pxa,pya,pza
+
+    def symplectic(self, h, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza, b):
+        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
+        x = x + h / 2 * vpx
+        y = y + h / 2 * vpy
+        z = z + h / 2 * vpz
+        pxa = pxa - h / 2 * vxa
+        pya = pya - h / 2 * vya
+        pza = pza - h / 2 * vza
+
+        vx, vy, vz, vpxa, vpya, vpza, e = self.f(x, y, z, pxa, pya, pza, b)
+        xa = xa + h * vpxa
+        ya = ya + h * vpya
+        za = za + h * vpza
+        px = px - h * vx
+        py = py - h * vy
+        pz = pz - h * vz
+
+        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
+        x = x + h / 2 * vpx
+        y = y + h / 2 * vpy
+        z = z + h / 2 * vpz
+        pxa = pxa - h / 2 * vxa
+        pya = pya - h / 2 * vya
+        pza = pza - h / 2 * vza
+
+        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza

utils/metrics.py

@@ -0,0 +1,28 @@
+import numpy as np
+from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
+
+def calculate_metrics(y, y_hat, y_train=None):
+    def smape(a, f):
+        return 1/len(a) * np.sum(2 * np.abs(f - a) / (np.abs(a) + np.abs(f) + np.finfo(float).eps))
+
+    def mase(y_actual, y_pred, y_train):
+        n = y_train.shape[1]
+        d = np.abs(np.diff(y_train)).sum() / (n - 1)
+        errors = np.abs(y_actual - y_pred)
+        return errors.mean() / d
+
+    def phase_space_distance(y_actual, y_pred):
+        return np.sqrt(np.sum(np.square(y_actual - y_pred)))
+
+    SMAPE = np.mean([smape(yi.reshape(-1), y_hati.reshape(-1)) for yi, y_hati in zip(y, y_hat)])
+    MSE = np.mean([mean_squared_error(yi.reshape(-1), y_hati.reshape(-1)) for yi, y_hati in zip(y, y_hat)])
+    RMSE = np.mean([np.sqrt(mean_squared_error(yi.reshape(-1), y_hati.reshape(-1))) for yi, y_hati in zip(y, y_hat)])
+    MAE = np.mean([mean_absolute_error(yi.reshape(-1), y_hati.reshape(-1)) for yi, y_hati in zip(y, y_hat)])
+    R2 = np.mean([r2_score(yi.reshape(-1), y_hati.reshape(-1)) for yi, y_hati in zip(y, y_hat)])
+    PSD = np.mean([phase_space_distance(yi.reshape(-1), y_hati.reshape(-1)) for yi, y_hati in zip(y, y_hat)])
+
+    if y_train is None:
+        return SMAPE, MSE, RMSE, MAE, R2, PSD
+    else:
+        MASE = np.mean([mase(yi, y_hati, yt) for yi, y_hati, yt in zip(y, y_hat, y_train)])
+        return SMAPE, MSE, RMSE, MAE, R2, MASE, PSD

utils/midpoint.py

@@ -0,0 +1,164 @@
+
+class MidPoint():
+
+    def __init__(self, j):
+        self.j = j
+
+    def initial(self):
+        j = self.j
+
+        #init parameters
+        h, b, n, u = 0.1, None, None, None
+        x, y, z, px, py, pz = None, 0.1, 0.001, 0.01, None, 0.0001
+        xa, ya, za, pxa, pya, pza = None, None, None, None, None, None
+
+        if j == 1:
+            b = 5.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 10.0
+            py = 0.5
+        elif j == 2:
+            b = 3.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 8.3
+            py = 0.6
+        elif j == 3:
+            b = 3.0 / 2
+            x = 12.0
+            py = 0.4
+        elif j == 4:
+            b = 7.0 / 4
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 15.0
+            py = 0.35
+        elif j == 5:
+            b = 1.0
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 18.0
+            py = 0.3
+        elif j == 6:
+            b = 3.0 / 5
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 20.0
+            py = 0.25
+        elif j == 7:
+            b = 5.0 / 7
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 22.0
+            py = 0.22
+        elif j == 8:
+            b = 2.0
+            x = 26.0
+            py = 0.2
+        elif j == 9:
+            b = 0.5
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 30.0
+            y = 0.5
+            z = 0.1
+            pz = 0.01
+        elif j == 10:
+            b = 5.0
+            n = b / (1 + b) ** 2
+            u = 1.0 / (1.0 / b + b + 2.0)
+            x = 35.0
+            y = 2.0
+            z = 0.1
+            pz = 0.03
+            py = 0.15
+
+        xa, ya, za, pxa, pya, pza = x, y, z, px, py, pz
+        return j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza
+
+    def f(self, x, y, z, px, py, pz, b):
+
+        n = b / (1 + b) ** 2
+        u = 1 / (1 / b + b + 2)
+        ht = px ** 2 / 2 + py ** 2 / 2 + pz ** 2 / 2
+        hv = -1 / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        h1pn = 1 / (2 * x ** 2 + 2 * y ** 2 + 2 * z ** 2) - (((u + 3) * (px ** 2 + py ** 2 + pz ** 2)) / 2 + (u * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) ** 2) / 2) / (
+                           x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + ((3 * u) / 8 - 1 / 8) * (
+                           px ** 2 + py ** 2 + pz ** 2) ** 2
+
+        e = ht + hv + h1pn
+
+        vnpx=px
+        v1pnpx=4*px*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (px*(n + 3) + (n*x*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
+        vpx=vnpx+v1pnpx
+
+        vnpy=py
+        v1pnpy=4*py*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (py*(n + 3) + (n*y*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
+        vpy=vnpy+v1pnpy
+
+        vnpz=pz
+        v1pnpz=4*pz*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (pz*(n + 3) + (n*z*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
+        vpz=vnpz+v1pnpz
+
+        vnx=x/(x**2 + y**2 + z**2)**(3/2)
+        v1pnx=(x*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*x)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((px*x**2)/(x**2 + y**2 + z**2)**(3/2) -px/(x**2 + y**2 + z**2)**(1/2) + (py*x*y)/(x**2 + y**2 + z**2)**(3/2) + (pz*x*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
+        vx=vnx+v1pnx
+
+        vny=y/(x**2 + y**2 + z**2)**(3/2)
+        v1pny=(y*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*y)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((py*y**2)/(x**2 + y**2 + z**2)**(3/2) -py/(x**2 + y**2 + z**2)**(1/2) + (px*x*y)/(x**2 + y**2 + z**2)**(3/2) + (pz*y*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
+        vy=vny+v1pny
+
+        vnz=z/(x**2 + y**2 + z**2)**(3/2)
+        v1pnz=(z*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*z)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((pz*z**2)/(x**2 + y**2 + z**2)**(3/2) -pz/(x**2 + y**2 + z**2)**(1/2) + (px*x*z)/(x**2 + y**2 + z**2)**(3/2) + (py*y*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
+        vz=vnz+v1pnz
+
+        return vx,vy,vz,vpx,vpy,vpz,e
+
+    def rejust(self, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza):
+
+        x = (x + xa) / 2
+        y = (y + ya) / 2
+        z = (z + za) / 2
+
+        px = (px + pxa) / 2
+        py = (py + pya) / 2
+        pz = (pz + pza) / 2
+        xa = x
+        ya = y
+        za = z
+        pxa = px
+        pya = py
+        pza = pz
+
+        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza
+
+    def symplectic(self, h, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza, b):
+
+        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
+        x = x + h / 2 * vpx
+        y = y + h / 2 * vpy
+        z = z + h / 2 * vpz
+        pxa = pxa - h / 2 * vxa
+        pya = pya - h / 2 * vya
+        pza = pza - h / 2 * vza
+
+        vx, vy, vz, vpxa, vpya, vpza, e = self.f(x, y, z, pxa, pya, pza, b)
+        xa = xa + h * vpxa
+        ya = ya + h * vpya
+        za = za + h * vpza
+        px = px - h * vx
+        py = py - h * vy
+        pz = pz - h * vz
+
+        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
+        x = x + h / 2 * vpx
+        y = y + h / 2 * vpy
+        z = z + h / 2 * vpz
+        pxa = pxa - h / 2 * vxa
+        pya = pya - h / 2 * vya
+        pza = pza - h / 2 * vza
+
+        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza

utils/transform.py

@@ -0,0 +1,8 @@
+import torch
+
+def compute_gradient(x, degree):
+    gradients = [x]
+    for i in range(degree):
+        x = torch.diff(x, dim=-1, prepend=x[..., 0:1])
+        gradients.append(x)
+    return torch.concatenate(gradients, dim=-1)