Merge pull request #4 from Steriva/main

Addition of Test and Automatic Setup

.github/workflows/testing.yml

@@ -0,0 +1,52 @@
+# This workflow will install Python dependencies, run tests and lint with a single version of Python
+# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python
+
+name: Testing pyforce
+
+on:
+  workflow_dispatch:
+  pull_request:
+    branches: [ "main" ]
+
+permissions:
+  contents: read
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v4
+
+    - name: Set up Miniconda
+      uses: conda-incubator/setup-miniconda@v2
+      with:
+        auto-update-conda: true
+        python-version: "3.10"
+        channels: conda-forge
+        use-mamba: true  # Optional: to speed up the environment solving
+
+    - name: Create Conda environment
+      run: |
+        conda env create -f pyforce/environment.yml
+
+    - name: Activate Conda environment
+      run: |
+        conda activate pyforce-env
+        python -m pip install pytest
+        python -m pip install pyforce/
+      shell: bash -l {0}
+
+    # - name: Lint with flake8
+    #   run: |
+    #     # stop the build if there are Python syntax errors or undefined names
+    #     flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+    #     # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+    #     flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Test with pytest
+      run: |
+        conda init
+        conda activate pyforce-env
+        pytest
+      shell: bash -l {0}

.gitignore

@@ -75,6 +75,9 @@ instance/
 # Sphinx documentation
 docs/_build/
 
+# Testing Notebooks
+tests/ongoing/
+
 # PyBuilder
 target/
 

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2023 ROSE
+Copyright (c) 2024 ROSE-pyforce
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -4,9 +4,9 @@
   </a>
 </p>
 
-[![Reference Paper 1](https://img.shields.io/badge/Reference%20Paper%201-arXiv:%202401.07300-gray?labelColor=blue&style=flat&link=https://arxiv.org/abs/2401.07300)](https://arxiv.org/abs/2401.07300) [![Reference Paper 2](https://img.shields.io/badge/Reference%20Paper%202-10.1016/j.nucengdes.2024.113105-gray?labelColor=blue&style=flat&link=https://www.sciencedirect.com/science/article/pii/S002954932400205X)](https://www.sciencedirect.com/science/article/pii/S002954932400205X) [![Docs](https://img.shields.io/badge/Docs-green?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/intro.html)](https://ermete-lab.github.io/ROSE-pyforce/intro.html) [![Tutorials](https://img.shields.io/badge/Tutorials-red?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/tutorials.html)](https://ermete-lab.github.io/ROSE-pyforce/tutorials.html)
+[![Reference Paper 1](https://img.shields.io/badge/Reference%20Paper%201-arXiv:%202401.07300-gray?labelColor=blue&style=flat&link=https://arxiv.org/abs/2401.07300)](https://arxiv.org/abs/2401.07300) [![Reference Paper 2](https://img.shields.io/badge/Reference%20Paper%202-10.1016/j.nucengdes.2024.113105-gray?labelColor=blue&style=flat&link=https://www.sciencedirect.com/science/article/pii/S002954932400205X)](https://www.sciencedirect.com/science/article/pii/S002954932400205X) [![Testing pyforce](https://github.com/Steriva/ROSE-pyforce/actions/workflows/testing.yml/badge.svg)](https://github.com/Steriva/ROSE-pyforce/actions/workflows/testing.yml) [![Docs](https://img.shields.io/badge/Docs-green?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/intro.html)](https://ermete-lab.github.io/ROSE-pyforce/intro.html) [![Tutorials](https://img.shields.io/badge/Tutorials-red?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/tutorials.html)](https://ermete-lab.github.io/ROSE-pyforce/tutorials.html)
 
-**pyforce: Python Framework data-driven model Order Reduction for multi-physiCs problEms** 
+**pyforce: Python Framework data-driven model Order Reduction for multi-physiCs problEms**
 
 - [Description](#description)
 - [How to cite *pyforce*](#how-to-cite-pyforce)
@@ -17,7 +17,7 @@
 
 ## Description
 
-*pyforce* is a Python package implementing some Data-Driven Reduced Order Modelling (DDROM) techniques for applications to multi-physics problems, mainly set in the **Nuclear Engineering** world. These techniques have been implemented upon the [dolfinx](https://github.com/FEniCS/dolfinx) package (currently v0.6.0), part of the [FEniCSx](https://fenicsproject.org/) project, to handle mesh generation, integral calculation and functions storage. The package is part of the **ROSE (Reduced Order modelling with data-driven techniques for multi-phySics problEms)**: mathematical algorithms aimed at reducing the complexity of multi-physics models (for nuclear reactors applications), at searching for optimal sensor positions and at integrating real measures to improve the knowledge on the physical systems.
+*pyforce* is a Python package implementing Data-Driven Reduced Order Modelling (DDROM) techniques for applications to multi-physics problems, mainly set in the **Nuclear Engineering** world. These techniques have been implemented upon the [dolfinx](https://github.com/FEniCS/dolfinx) package (currently v0.6.0), part of the [FEniCSx](https://fenicsproject.org/) project, to handle mesh generation, integral calculation and functions storage. The package is part of the **ROSE (Reduced Order modelling with data-driven techniques for multi-phySics problEms)**: mathematical algorithms aimed at reducing the complexity of multi-physics models (for nuclear reactors applications), at searching for optimal sensor positions and at integrating real measures to improve the knowledge on the physical systems.
 
 The techniques implemented here follow the same underlying idea expressed in the following figure: in the offline (training) phase, a dimensionality reduction process retrieves a reduced coordinate system onto which encodes the information of the mathematical model; the sensor positioning algorithm then uses this set to select the optimal location of sensors according to some optimality criterion, which depends on the adopted algorithm. In the online phase, the DA process begins, retrieving a novel set of reduced variables and then computing the reconstructed state through a decoding step.
 
@@ -33,11 +33,11 @@ At the moment, the following techniques have been implemented:
 - **Parameterised-Background Data-Weak formulation**
 - an **Indirect Reconstruction** algorithm to reconstruct non-observable fields
 
-This package is aimed to be a valuable tool for other researchers, engineers, and data scientists working in various fields, not only restricted in the Nuclear Engineering world. 
+This package is aimed to be a valuable tool for other researchers, engineers, and data scientists working in various fields, not only restricted in the Nuclear Engineering world.
 
 ## How to cite *pyforce*
 
-If you are going to use *pyforce* in your research work, please cite the following articles. 
+If you are going to use *pyforce* in your research work, please cite the following articles.
 The authors would be pleased if you could cite the relevant papers:
 
 1. Stefano Riva, Carolina Introini, and Antonio Cammi. Multi-physics model bias correction with data-driven reduced order modelling techniques: Application to nuclear case studies, 2024. [arXiv:2401.07300](http://arxiv.org/abs/2401.07300).
@@ -48,7 +48,7 @@ For LaTeX users:
 ```bibtex
 
 @misc{RMP_2024,
-      title={Multi-Physics Model Bias Correction with Data-Driven Reduced Order Modelling Techniques: Application to Nuclear Case Studies}, 
+      title={Multi-Physics Model Bias Correction with Data-Driven Reduced Order Modelling Techniques: Application to Nuclear Case Studies},
       author={Stefano Riva and Carolina Introini and Antonio Cammi},
       year={2024},
       eprint={2401.07300},
@@ -78,15 +78,20 @@ keywords = {Hybrid Data-Assimilation, Generalized Empirical Interpolation Method
 - Stefano Riva, Sophie Deanesi, Carolina Introini, Stefano Lorenzi, and Antonio Cammi. Neutron Flux Re- construction from Out-Core Sparse Measurements using Data-Driven Reduced Order Modelling. In International Conference on Physics of Reactors (PHYSOR24), San Francisco, USA, April 2024.
 
 ## Installation
-The package can be installed using `pip`, make sure all the dependencies are installed (following these [steps](https://ermete-lab.github.io/ROSE-pyforce/installation.html#set-up-a-conda-environment-for-pyforce)).
-The requirements are listed [here](https://github.com/ERMETE-Lab/ROSE-pyforce/blob/main/pyforce/requirements.txt).
+The package can be installed using `pip`, make sure all the dependencies are installed (following these [steps](https://ermete-lab.github.io/ROSE-pyforce/installation.html#set-up-a-conda-environment-for-pyforce)). The requirements are listed [here](https://github.com/ERMETE-Lab/ROSE-pyforce/blob/main/pyforce/requirements.txt).
 
-Clone the repository 
+It is suggested to create a conda environment: at first, clone the repository
 ```bash
-git clone https://github.com/ROSE-Polimi/pyforce.git
+git clone https://github.com/ERMETE-Lab/ROSE-pyforce.git
 ```
-Change directory to *pyforce* and install using `pip`
-```bash 
+create a conda environment using `environment.yml`
+```bash
+cd ROSE-pyforce
+conda env create -f pyforce/environment.yml
+```
+activate the environment and then install the package using `pip`
+```bash
+conda activate pyforce-env
 python -m pip install pyforce/
 ```
 

docs/Tutorials/01_LaminarNS/01_generateFlowData.ipynb

@@ -19,7 +19,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+    "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -49,7 +49,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+    "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -106,7 +106,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -396,7 +396,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.10.14"
   }
  },
  "nbformat": 4,

pyforce/environment.yml

@@ -0,0 +1,21 @@
+name: pyforce-env
+channels:
+  - conda-forge
+dependencies:
+  - python=3.10
+  - numpy=1.23.5
+  - scipy=1.13.0
+  - matplotlib=3.8.0
+  - pyvista=0.42.2
+  - fenics-dolfinx=0.6.0
+  - mpi4py=3.1.4
+  - petsc4py=3.18.5
+  - h5py=3.9.0
+  - pip
+  - pip:
+    - tqdm
+    - gmsh
+    - gmsh-api
+    - fluidfoam
+    - scikit-learn
+    - setuptools==62.0.0

pyforce/pyforce/offline/pod.py

@@ -1,13 +1,14 @@
 # Offline Phase: Proper Orthogonal Decomposition
 # Author: Stefano Riva, PhD Student, NRG, Politecnico di Milano
-# Latest Code Update: 05 March 2024
-# Latest Doc  Update: 05 March 2024
+# Latest Code Update: 24 May 2024
+# Latest Doc  Update: 24 May 2024
 
 import numpy as np
 import scipy
 import warnings
 
 from dolfinx.fem import (Function, FunctionSpace)
+from sklearn.utils.extmath import randomized_svd
 
 from pyforce.tools.backends import norms, LoopProgress
 from pyforce.tools.functions_list import *
@@ -317,23 +318,22 @@ class DiscretePOD():
       Name of the field.
     Nmax : int, optional (default=None)
       If `None` the full matrices are stored, else only the first `Nmax`.
+    random : bool, optional (default = False)
+      If True and if `Nmax` is provided, the randomised SVD is used.
     """
-  def __init__(self, train_snap: FunctionsMatrix, name: str, Nmax = None) -> None:
+  def __init__(self, train_snap: FunctionsList, name: str, Nmax = None, random = False) -> None:
 
     self.Ns = len(train_snap)
     self.Nh = len(train_snap(0))
     self.name = name
-
-    # Creating the snapshots matrix - check the variable type
-    if isinstance(train_snap, FunctionsList):
-      snapMatrix = fun_list_2_fun_matrix(train_snap)
-    else:
-      snapMatrix = train_snap
-
-    self.modes = FunctionsMatrix(snapMatrix.dofs)
+
+    self.modes = FunctionsList(dofs = train_snap.fun_shape)
 
     # Performing SVD of the snapshot matrix
-    U, Sigma, V_T = np.linalg.svd(snapMatrix.return_matrix(), full_matrices=False)
+    if random and Nmax is not None:
+      U, Sigma, V_T  = randomized_svd(train_snap.return_matrix(), n_components=Nmax, n_iter='auto')
+    else:
+      U, Sigma, V_T = np.linalg.svd(train_snap.return_matrix(), full_matrices=False)
     if sum(Sigma < 0) > 0: 
         warnings.warn("Check singular values: some of them are negative!")   
 
@@ -404,7 +404,7 @@ def reconstruct(self, Vh_star: np.ndarray):
 
     return np.dot(self.modes.return_matrix()[:, :N] * self.sing_vals[:N], Vh_star)
 
-  def train_error(self, train_snap: FunctionsMatrix, maxBasis: int, verbose = False):
+  def train_error(self, train_snap: FunctionsList, maxBasis: int, verbose = False):
     r"""
     The maximum absolute :math:`E_N` and relative :math:`\varepsilon_N` error on the train set is computed, by projecting it onto the reduced space in :math:`l^2`-sense
 

pyforce/pyforce/offline/sensors.py

@@ -369,7 +369,9 @@ def generate(self, N: int, Mmax: int, tol: float = 0.2,
     if m < Mmax:
         if verbose:
             print(' ')
-            print('Starting approximation loop', end = "\r")
+            print('-----------------------------------------')
+            print('Starting approximation loop')
+            print(' ')
         self.approx_loop(Mmax, is_H1=is_H1)
         m = Mmax
 

pyforce/pyforce/online/pod_interpolation.py

@@ -1,7 +1,7 @@
 # Synthetic Online Phase: Proper Orthogonal Decomposition with Interpolation
 # Author: Stefano Riva, PhD Student, NRG, Politecnico di Milano
-# Latest Code Update: 06 December 2023
-# Latest Doc  Update: 06 December 2023
+# Latest Code Update: 30 April 2024
+# Latest Doc  Update: 30 April 2024
 
 import numpy as np
 import scipy.linalg as la
@@ -20,7 +20,7 @@ class PODI():
     modes : FunctionsList
         List of POD modes computed during the offline phase.
     maps : list
-        List of maps for the POD modal coefficients, they must be callable.
+        List of maps for the POD modal coefficients, they must be callable. If `None`, the reduced coefficient must be provided as input later!
     name : str
         Name of the snapshots (e.g., temperature T)
     """
@@ -56,12 +56,11 @@ def synt_test_error(self, test_snap: FunctionsList, mu_estimated: np.ndarray, ma
         test_snap : FunctionsList
             List of snapshots onto which the test error of the POD basis is performed.
         mu_estimated : np.ndarray
-            Arrays with the estimated parameters from the Parameter Estimation phase, it must have dimension `[Ns, p]` in which `Ns` is the number of test snapshots and `p` the number of parameters.
+            Arrays with the estimated parameters from the Parameter Estimation phase, it must have dimension `[Ns, p]` in which `Ns` is the number of test snapshots and `p` the number of parameters. If `None`, the reduced coefficients `alpha_coeffs` must be given.
         maxBasis : int
             Integer input indicating the maximum number of modes to use.
         alpha_coeff : np.ndarray (optional, Default: None)
             Matrix with the estimated coefficients :math:`\alpha_n`, they will be used if the input `alpha_coeffs` is not `None`.
-        
         verbose : boolean, optional (Default = False) 
             If `True`, print of the progress is enabled.
 
@@ -80,6 +79,8 @@ def synt_test_error(self, test_snap: FunctionsList, mu_estimated: np.ndarray, ma
         if mu_estimated is not None:
             assert (mu_estimated.shape[0] == Ns_test)
             n_feature = mu_estimated.shape[1]
+        else:
+            assert alpha_coeffs is not None, 'Both inputs mu_estimated and alpha_coeffs are None'
 
         absErr = np.zeros((Ns_test, maxBasis))
         relErr = np.zeros_like(absErr)
@@ -90,7 +91,6 @@ def synt_test_error(self, test_snap: FunctionsList, mu_estimated: np.ndarray, ma
         # Variables to store the computational times
         computational_time = dict()
         computational_time['CoeffEstimation'] = np.zeros((Ns_test, maxBasis))
-        computational_time['LinearSystem']    = np.zeros((Ns_test, maxBasis))
         computational_time['Errors']          = np.zeros((Ns_test, maxBasis))
 
         timing = Timer() 
@@ -113,7 +113,6 @@ def synt_test_error(self, test_snap: FunctionsList, mu_estimated: np.ndarray, ma
                     coeff[nn] = self.maps[nn](mu_estimated[mu].reshape(-1, n_feature))
                 computational_time['CoeffEstimation'][mu, nn] = timing.stop()
 
-            # for nn in range(maxBasis):
                 # building residual field and computing the errors
                 timing.start()
                 resid.x.array[:] = test_snap(mu) - self.PODmodes.lin_combine(coeff[:nn+1])
@@ -126,9 +125,10 @@ def synt_test_error(self, test_snap: FunctionsList, mu_estimated: np.ndarray, ma
 
         return absErr.mean(axis = 0), relErr.mean(axis = 0), computational_time
 
-    def reconstruct(self, snap: np.ndarray, mu_estimated: np.ndarray, maxBasis: int):
+    def reconstruct(self, snap: np.ndarray, mu_estimated: np.ndarray, maxBasis: int, 
+                          alpha_coeffs : np.ndarray = None):
         r"""
-        After the coefficients of the POD basis are obtained by interpolating using the maps, the `snap` is approximated using linear combination of the POD modes.
+        After the coefficients of the POD basis are obtained by interpolating using the maps or given as input, the `snap` is approximated using linear combination of the POD modes.
         
         Parameters
         ----------
@@ -138,6 +138,8 @@ def reconstruct(self, snap: np.ndarray, mu_estimated: np.ndarray, maxBasis: int)
             Arrays with the estimated parameters from the Parameter Estimation phase, it must have dimension `[1, p]` in which `p` the number of parameters.
         maxBasis : int
             Integer input indicating the maximum number of modes to use.
+        alpha_coeff : np.ndarray (optional, Default: None)
+            Array with the estimated coefficients :math:`\alpha_n`, they will be used if the input `alpha_coeffs` is not `None`.
         
         Returns
         -------
@@ -147,17 +149,23 @@ def reconstruct(self, snap: np.ndarray, mu_estimated: np.ndarray, maxBasis: int)
             Residual field using `maxBasis` POD modes.
         """
 
-        n_feature = mu_estimated.shape[1]
 
         # Variables to store the computational times
         computational_time = dict()
         timing = Timer() 
 
         # Estimate the coefficients
         timing.start()
-        coeff = np.zeros((maxBasis,))
-        for nn in range(maxBasis):
-            coeff[nn] = self.maps[nn](mu_estimated.reshape(-1, n_feature))
+
+        if mu_estimated is not None:
+            n_feature = mu_estimated.shape[1]
+            coeff = np.zeros((maxBasis,))
+            for nn in range(maxBasis):
+                coeff[nn] = self.maps[nn](mu_estimated.reshape(-1, n_feature))
+        else:
+            assert alpha_coeffs is not None, 'Both inputs mu_estimated and alpha_coeffs are None'
+            coeff = alpha_coeffs.reshape(maxBasis,)
+
         computational_time['CoeffEstimation'] = timing.stop()
 
         if isinstance(snap, Function):

pyforce/pyforce/online/pod_projection.py

@@ -76,7 +76,6 @@ def synt_test_error(self, test_snap: FunctionsList, maxBasis: int,
         # Variables to store the computational times
         computational_time = dict()
         computational_time['CoeffEstimation'] = np.zeros((Ns_test, maxBasis))
-        computational_time['LinearSystem']    = np.zeros((Ns_test, maxBasis))
         computational_time['Errors']          = np.zeros((Ns_test, maxBasis))
 
         timing = Timer()