hyperReduction.templates.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
     ██╗████████╗██╗  ██╗ █████╗  ██████╗ █████╗       ███████╗██╗   ██╗
     ██║╚══██╔══╝██║  ██║██╔══██╗██╔════╝██╔══██╗      ██╔════╝██║   ██║
     ██║   ██║   ███████║███████║██║     ███████║█████╗█████╗  ██║   ██║
     ██║   ██║   ██╔══██║██╔══██║██║     ██╔══██║╚════╝██╔══╝  ╚██╗ ██╔╝
     ██║   ██║   ██║  ██║██║  ██║╚██████╗██║  ██║      ██║      ╚████╔╝
     ╚═╝   ╚═╝   ╚═╝  ╚═╝╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝      ╚═╝       ╚═══╝
 
 * In real Time Highly Advanced Computational Applications for Finite Volumes
 * Copyright (C) 2017 by the ITHACA-FV authors
-------------------------------------------------------------------------------
 
License
    This file is part of ITHACA-FV
 
    ITHACA-FV is free software: you can redistribute it and/or modify
    it under the terms of the GN3U Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    ITHACA-FV is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public License
    along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#ifndef hyperReduction_templates_H
#define hyperReduction_templates_H
 
#include "hyperReduction.H"
 
// TODO increase modularity to avoid heavy compiling costs
 
// Template function constructor
template <typename... SnapshotsLists>
HyperReduction<SnapshotsLists...>::HyperReduction(label n_modes,
                                                  label n_nodes,
                                                  Eigen::VectorXi initialSeeds,
                                                  word problemName,
                                                  SnapshotsLists &&...snapshotsLists)
    : vectorial_dim{compute_vectorial_dim(snapshotsLists...)},
      n_modes{n_modes},
      n_nodes{n_nodes},
      initialSeeds{initialSeeds},
      problemName{problemName},
      snapshotsListTuple{std::forward_as_tuple(snapshotsLists...)}
{
  Info << "Init HyperReduction class with vectorial dim: " << vectorial_dim << endl;
  para = ITHACAparameters::getInstance();
  folderProblem = "ITHACAoutput/" + problemName;
  mkDir(folderProblem);
 
  // TODO check snapshotsLists is not empty
  // TODO check that the first snapshotsList is not empty (includes above)
 
  // TODO initialize vectorial_dim only from SnapshotsLists type
  // vectorial_dim = std::apply(compute_vectorial_dim<SnapshotsLists&&...>, snapshotsListTuple);
 
  // Get fields' names
  std::apply([this](auto& ...snapList){(..., stackNames(snapList));}, snapshotsListTuple);
 
  // Get fields' dimensions
  std::apply([this](auto& ...snapList){(..., stackDimensions(snapList));}, snapshotsListTuple);
 
  sumFieldsDim = 0;
  std::apply([this](auto& ...snapList){(..., sumDimensions(this->sumFieldsDim, snapList));}, snapshotsListTuple);
 
  Info << "Fields names: ";
  for (unsigned int ith_field = 0; ith_field < fieldNames.size(); ith_field++)
  {
    Info << fieldNames[ith_field] << " (dim=" << fieldDims[ith_field] << "); ";
  }
  Info << endl;
  
  n_snapshots = std::get<0>(snapshotsListTuple).size();
  Info << "The number of snapshots is: " << n_snapshots << endl;
  n_cells = std::get<0>(snapshotsListTuple)[0].size();
  Info << "The number of cells is: " << n_cells << endl;
  Info << "Initial seeds length: " << initialSeeds.rows() << endl;
 
  // get boundaries info
  n_boundary_patches = std::get<0>(snapshotsListTuple)[0].boundaryField().size();
 
  n_boundary_cells_list.resize(n_boundary_patches);
  n_boundary_cells = 0;
  for (unsigned int ith_boundary_patch = 0; ith_boundary_patch < n_boundary_patches; ith_boundary_patch++)
  {
    n_boundary_cells_list[ith_boundary_patch] = std::get<0>(snapshotsListTuple)[0].boundaryField()[ith_boundary_patch].size();
    n_boundary_cells = n_boundary_cells + n_boundary_cells_list[ith_boundary_patch];
  }
  Info << "Number of boundary cells: " << n_boundary_cells << endl;
}
 
 
// Template function constructor
template <typename... SnapshotsLists>
HyperReduction<SnapshotsLists...>::HyperReduction(label n_modes,
                                                  label n_nodes,
                                                  unsigned int vectorialDim,
                                                  label n_cells,
                                                  Eigen::VectorXi initialSeeds,
                                                  word problemName)
    : n_modes{n_modes},
      n_nodes{n_nodes},
      vectorial_dim{vectorialDim},
      n_cells{n_cells},
      initialSeeds{initialSeeds},
      problemName{problemName}
{
  Info << "Init HyperReduction class with vectorial dim: " << vectorial_dim << endl;
  para = ITHACAparameters::getInstance();
  folderProblem = problemName;
  mkDir(folderProblem);
 
  Info << "Initial seeds length: " << initialSeeds.rows() << endl;
 
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::getModesSVD(
    SnapshotsListTuple& snapshotsListTuple, Eigen::MatrixXd& modesSVD, Eigen::VectorXd& fieldWeights, bool saveModesFlag)
{
  // TODO move inside ITHACAPOD
  ITHACAparameters* para(ITHACAparameters::getInstance());
  word folderSVD = "ITHACAoutput/" + problemName + "/ModesSVD/";
  if (!ITHACAutilities::check_folder(folderSVD))
  {
    mkDir(folderSVD);
    // Initialize stacked snapshotsMatrix and normalizing weights
    Eigen::MatrixXd snapMatrix;
    getSnapMatrix(snapMatrix, fieldWeights);
 
    word dimensionReductionMethod = para->ITHACAdict->lookupOrDefault<word>("DimensionReductionMethod", "SVD");
 
    Eigen::MatrixXd eigenVectoreig;
 
    if (dimensionReductionMethod == "RSVD")
    {
      Info << "####### Performing RSVD for " << problemName << " #######" << endl;
 
      unsigned int r_modeSVD = para->ITHACAdict->lookupOrDefault<unsigned int>("RSVDdim", snapMatrix.cols());
      RedSVD::RedSVD svd(fieldWeights.array().cwiseInverse().matrix().asDiagonal()*snapMatrix, r_modeSVD);
 
      Info << "####### End of RSVD for " << problemName << " #######" <<
            endl;
 
      eigenValueseig = svd.singularValues().real();
      eigenVectoreig = svd.matrixU().real();
    }
    else
    {
      Info << "####### Performing SVD for " << problemName << " #######" << endl;
 
      Eigen::JacobiSVD<Eigen::MatrixXd> svd(fieldWeights.array().cwiseInverse().matrix().asDiagonal()*snapMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
 
      Info << "####### End of SVD for " << problemName << " #######" <<
            endl;
 
      eigenValueseig = svd.singularValues().real();
      eigenVectoreig = svd.matrixU().real();
    }
 
    modesSVD = eigenVectoreig;
 
    // TODO correct boundary conditions
    
    cnpy::save(eigenValueseig, "ITHACAoutput/" + problemName + "/ModesSVD/eigenvalues.npy");
 
    eigenValueseig = eigenValueseig / eigenValueseig.sum();
    Eigen::VectorXd cumEigenValues(eigenValueseig);
 
    for (label j = 1; j < cumEigenValues.size(); ++j)
    {
        cumEigenValues(j) += cumEigenValues(j - 1);
    }
 
    Eigen::saveMarketVector(eigenValueseig,
                            "ITHACAoutput/" + problemName + "/ModesSVD/Eigenvalues_" + problemName, para->precision,
                            para->outytpe);
    Eigen::saveMarketVector(cumEigenValues,
                            "ITHACAoutput/" + problemName + "/ModesSVD/CumEigenvalues_" + problemName, para->precision,
                            para->outytpe);
    cnpy::save(modesSVD, "ITHACAoutput/" + problemName + "/ModesSVD/modes.npy");
    cnpy::save(fieldWeights, "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeights.npy");
 
 
    if (saveModesFlag)
    {
      unsigned int rowIndex{0};
      for (unsigned int modeIndex = 0; modeIndex < modesSVD.cols(); modeIndex++)
      {
        std::apply([this, &modesSVD, &rowIndex, &modeIndex, &folderSVD](auto& ...snapList){(..., saveModes(snapList, modesSVD, rowIndex, modeIndex, folderSVD));}, snapshotsListTuple);
        rowIndex = 0;
      }
    }
  }
  else
  {
    Info << "Reading the existing modes" << endl;
    cnpy::load(modesSVD, folderSVD + "/modes.npy");
    cnpy::load(fieldWeights, folderSVD + "/normalizingWeights.npy");
  }
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::getModesSVD(
    SnapshotsListTuple& snapshotsListTuple, Eigen::MatrixXd& modesSVD, Eigen::VectorXd& fieldWeights, Eigen::MatrixXd& modesSVDBoundary, Eigen::VectorXd& fieldWeightsBoundary, bool saveModesFlag)
{
  // TODO move inside ITHACAPOD
  ITHACAparameters* para(ITHACAparameters::getInstance());
  word folderSVD = "ITHACAoutput/" + problemName + "/ModesSVD/";
  if (!ITHACAutilities::check_folder(folderSVD))
  {
    mkDir(folderSVD);
    // Initialize stacked snapshotsMatrix and normalizing weights
    Eigen::MatrixXd snapMatrix;
    List<Eigen::MatrixXd> snapMatrixBoundary;
    List<Eigen::VectorXd> fieldWeightsBoundaryList;
    getSnapMatrix(snapMatrix, fieldWeights, snapMatrixBoundary, fieldWeightsBoundaryList);
 
    for (unsigned int id = 0; id < n_boundary_patches; id++)
    {
      snapMatrix.conservativeResize(snapMatrix.rows()+snapMatrixBoundary[id].rows(), snapMatrix.cols());
      snapMatrix.bottomRows(snapMatrixBoundary[id].rows()) = snapMatrixBoundary[id];
 
      fieldWeightsBoundary.conservativeResize(fieldWeightsBoundary.rows()+snapMatrixBoundary[id].rows());
      fieldWeightsBoundary.tail(fieldWeightsBoundaryList[id].rows()) = fieldWeightsBoundaryList[id];
    }
 
    cnpy::save(fieldWeights, "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeights.npy");
    cnpy::save(fieldWeightsBoundary, "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeightsBoundary.npy");
 
    fieldWeights.conservativeResize(fieldWeights.rows()+fieldWeightsBoundary.rows());
    fieldWeights.tail(fieldWeightsBoundary.rows()) = fieldWeightsBoundary;
    
    word dimensionReductionMethod = para->ITHACAdict->lookupOrDefault<word>("DimensionReductionMethod", "SVD");
 
    Eigen::MatrixXd eigenVectoreig;
 
    if (dimensionReductionMethod == "RSVD")
    {
      Info << "####### Performing RSVD for " << problemName << " #######" << endl;
 
      unsigned int r_modeSVD = para->ITHACAdict->lookupOrDefault<unsigned int>("RSVDdim", snapMatrix.cols());
 
      RedSVD::RedSVD svd(fieldWeights.array().cwiseInverse().matrix().asDiagonal()*snapMatrix, r_modeSVD);
 
      Info << "####### End of RSVD for " << problemName << " #######" << endl;
 
      eigenValueseig = svd.singularValues().real();
      eigenVectoreig = svd.matrixU().real();
    }
    else
    {
      Info << "####### Performing SVD for " << problemName << " #######" << endl;
 
      Eigen::JacobiSVD<Eigen::MatrixXd> svd(fieldWeights.array().cwiseInverse().matrix().asDiagonal()*snapMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
 
      Info << "####### End of SVD for " << problemName << " #######" <<
            endl;
 
      eigenValueseig = svd.singularValues().real();
      eigenVectoreig = svd.matrixU().real();
    }
 
    modesSVD = eigenVectoreig.bottomRows(vectorial_dim*n_cells);
    modesSVDBoundary = eigenVectoreig.topRows(n_boundary_cells);
    fieldWeights.conservativeResize(fieldWeights.rows()-fieldWeightsBoundary.rows());
    
    cnpy::save(eigenValueseig, "ITHACAoutput/" + problemName + "/ModesSVD/eigenvalues.npy");
 
    eigenValueseig = eigenValueseig / eigenValueseig.sum();
    Eigen::VectorXd cumEigenValues(eigenValueseig);
 
    for (label j = 1; j < cumEigenValues.size(); ++j)
    {
        cumEigenValues(j) += cumEigenValues(j - 1);
    }
 
    Eigen::saveMarketVector(eigenValueseig,
                            "ITHACAoutput/" + problemName + "/ModesSVD/Eigenvalues_" + problemName, para->precision,
                            para->outytpe);
    Eigen::saveMarketVector(cumEigenValues,
                            "ITHACAoutput/" + problemName + "/ModesSVD/CumEigenvalues_" + problemName, para->precision,
                            para->outytpe);
    cnpy::save(modesSVD, "ITHACAoutput/" + problemName + "/ModesSVD/modes.npy");
    cnpy::save(modesSVDBoundary, "ITHACAoutput/" + problemName + "/ModesSVD/modesBoundary.npy");
 
    if (saveModesFlag)
    {
      unsigned int rowIndex{0};
      unsigned int rowIndexBoundary{0};
      for (unsigned int modeIndex = 0; modeIndex < modesSVD.cols(); modeIndex++)
      {
        std::apply([this, &modesSVD,  &modesSVDBoundary, &rowIndex, &modeIndex, &rowIndexBoundary, &folderSVD](auto& ...snapList){(..., saveModes(snapList, modesSVD, modesSVDBoundary, rowIndex, rowIndexBoundary, modeIndex, folderSVD));}, snapshotsListTuple);
        rowIndex = 0;
      }
    }
  }
  else
  {
    Info << "Reading the existing modes" << endl;
    cnpy::load(modesSVD, folderSVD + "/modes.npy");
    cnpy::load(modesSVDBoundary, "ITHACAoutput/" + problemName + "/ModesSVD/modesBoundary.npy");
    cnpy::load(fieldWeights, folderSVD + "/normalizingWeights.npy");
    cnpy::load(fieldWeightsBoundary, "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeightsBoundary.npy");
  }
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::getSnapMatrix(Eigen::MatrixXd& snapMatrix, Eigen::VectorXd& fieldWeights)
{
  std::apply([this, &fieldWeights, &snapMatrix](auto& ...snapList){(..., stackSnapshots(snapList, snapMatrix, fieldWeights));}, snapshotsListTuple);
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::getSnapMatrix(Eigen::MatrixXd& snapMatrix, Eigen::VectorXd& fieldWeights, List<Eigen::MatrixXd>& snapMatrixBoundary, List<Eigen::VectorXd>& fieldWeightsBoundary)
{
  fieldWeightsBoundary.resize(n_boundary_patches);
  snapMatrixBoundary.resize(n_boundary_patches);
 
  std::apply([this, &fieldWeights, &snapMatrix, &fieldWeightsBoundary, &snapMatrixBoundary](auto& ...snapList){
    (..., stackSnapshots(snapList, snapMatrix, fieldWeights));
    (..., stackSnapshotsBoundary(snapList, snapMatrixBoundary, fieldWeightsBoundary));},
  snapshotsListTuple);
}
 
template <typename... SnapshotsLists>
template <typename SnapshotsList>
void HyperReduction<SnapshotsLists...>::saveModes(SnapshotsList sList, Eigen::MatrixXd& snapshotsMatrix, unsigned int &rowIndex, unsigned int &modeIndex, word folder)
{
  unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
  auto fieldName = sList[0].name();
  unsigned int fieldSize = field_dim*sList[0].size();
  Eigen::VectorXd fieldBlock = snapshotsMatrix.block(rowIndex, modeIndex, fieldSize, 1);
  auto fieldOut = Foam2Eigen::Eigen2field(sList[0], fieldBlock , true);
  rowIndex +=fieldSize;
  ITHACAstream::exportSolution(fieldOut, name(modeIndex), folder);
}
 
template <typename... SnapshotsLists>
template <typename SnapshotsList>
void HyperReduction<SnapshotsLists...>::saveModes(SnapshotsList sList, Eigen::MatrixXd& snapshotsMatrix, Eigen::MatrixXd& snapshotsMatrixBoundary, unsigned int &rowIndex,  unsigned int& rowIndexBoundary, unsigned int &modeIndex, word folder)
{
  unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
  auto fieldName = sList[0].name();
  unsigned int fieldSize = field_dim*sList[0].size();
 
  Eigen::VectorXd fieldBlock = snapshotsMatrix.block(rowIndex, modeIndex, fieldSize, 1);
 
  List<Eigen::VectorXd> fieldBlockBoundary;
  fieldBlockBoundary.resize(n_boundary_patches);
  for (unsigned int id = 0; id < n_boundary_patches; id++)
  {
    unsigned int bfieldDim = field_dim*int(n_boundary_cells_list[id]);
    fieldBlockBoundary[id] = snapshotsMatrixBoundary.block(rowIndexBoundary, modeIndex, bfieldDim, 1);
    rowIndexBoundary+=bfieldDim;
  }
 
  auto fieldOut = Foam2Eigen::Eigen2field(sList[0], fieldBlock, fieldBlockBoundary);
  rowIndex +=fieldSize;
  
  ITHACAstream::exportSolution(fieldOut, name(modeIndex), folder);
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::offlineGappyDEIM(Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& weights, word folderMethodName)
{
  folderMethod = folderMethodName;
  Info << "FolderMethod : " << folderMethod << endl;
  mkDir(folderMethod);
 
  normalizingWeights = weights;
 
  nodePoints = autoPtr<IOList<label>>(new IOList<label>(
      IOobject("nodePoints", para->runTime.time().constant(), "../" + folderMethod,
                para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE))); 
 
  auto greedyMetric = para->ITHACAdict->lookupOrDefault<word>("GreedyMetric", "L2");
  
  bool offlineStage = !nodePoints().headerOk();
 
  if (offlineStage)
  {
    assert(n_modes > 0);
    assert(n_nodes >= n_modes);
 
    Eigen::VectorXd mp_not_mask = Eigen::VectorXd::Constant(n_cells * vectorial_dim, 1);
    std::set<label> nodePointsSet;
 
    // set initialSeeds
    if (initialSeeds.rows() > 0)
    {
      initSeeds(mp_not_mask, nodePointsSet);
    }
 
    int na = n_nodes - initialSeeds.rows();
    if (na > 0)
    {
      Eigen::SparseMatrix<double> reshapeMat;
      initReshapeMat(reshapeMat);
 
      if (greedyMetric=="L2")
      {
        Info << "####### Begin Greedy-L2 GappyDEIM #######" << endl;
        Info << "####### Modes=" <<  n_modes
              << ", nodePoints=" << n_nodes << " #######"
              << endl;
 
        int nb = 0;
        int nit = std::min(n_modes, na);
        int ncimin = std::floor(n_modes / nit);
        int naimin = std::floor(na / n_modes);
        Eigen::MatrixXd A;
        Eigen::VectorXd b;
        Eigen::VectorXd c;
        Eigen::VectorXd r;
        label ind_max, c1;
        double max;
 
        for (int i = 1; i <= nit; i++)
        {
          int nci = ncimin;
          // add the remaining modes in the quotient n_modes / nite
          if (i <= n_modes % nit)
          {
              nci = ncimin + 1;
          }
          int nai = naimin;
          // add the remaining nodePoints in the quotient na / n_modes
          if (i <= na % n_modes)
          {
              nai = naimin + 1;
          }
 
          Eigen::MatrixXd V;
 
          // select basis
          if (i == 1)
          {
              V = snapshotsModes.leftCols(nci);
          }
          else
          {
            for (int q = 1; q <= nci; q++)
            {
                A = P.transpose() * basisMatrix;
                b = P.transpose() * snapshotsModes.col(nb + q - 1);
                c = A.fullPivLu().solve(b);
                r = snapshotsModes.col(nb + q - 1) - basisMatrix * c;
                V.conservativeResize(snapshotsModes.rows(), q);
                V.col(q - 1) = r;
            }
          }
 
          for (int j = 1; j <= nai; j++)
          {
            if (P.cols() > 0)
            {
              max = (reshapeMat*(mp_not_mask.asDiagonal() * V)
                      .rowwise()
                      .lpNorm<2>().array().square().matrix())
                      .maxCoeff(&ind_max, &c1);
            }
            else
            {
              max = (reshapeMat*V.rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(&ind_max, &c1);
            }
            
            updateNodes(P, ind_max, mp_not_mask);
          }
 
          nb += nci;
          basisMatrix = snapshotsModes.leftCols(nb);
        }
        Info << "####### End of greedy GappyDEIM #######\n";
      }
      else if (greedyMetric=="SOPT")
      {
        Info << "####### Begin SOPT GappyDEIM #######" << endl;
        Info << "####### Modes=" <<  n_modes
              << ", nodePoints=" << n_nodes << " #######"
              << endl;
 
        label ind_max, c1;
        double max;
        Eigen::MatrixXd V = snapshotsModes.leftCols(n_modes);
 
        for (unsigned int ith_node = 0; ith_node < na; ith_node++)
        {
          if (P.cols() > 0)
          {
            // TODO make efficient
            Eigen::MatrixXd tmp = P.transpose()*V;
            Eigen::MatrixXd R = mp_not_mask.asDiagonal() * V;
            Eigen::MatrixXd inv = (tmp.transpose()*tmp).fullPivLu().inverse();
            Eigen::VectorXd num = 1 + (reshapeMat*((R * inv).array()*R.array()).rowwise().sum().matrix()).array();
            Eigen::VectorXd norma = (tmp.array().square()).colwise().sum();
            Eigen::MatrixXd adding = reshapeMat*R.array().square().matrix();
            adding.rowwise() += norma.transpose();
            Eigen::VectorXd denom = (adding.array().pow(1./V.cols()).matrix()).rowwise().prod();
 
            max = ((mp_not_mask.head(n_cells).asDiagonal()*num).array()/denom.array()).maxCoeff(&ind_max, &c1);
          }
          else
          {
            max = (reshapeMat*V.rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(&ind_max, &c1);
          }
 
          updateNodes(P, ind_max, mp_not_mask);
        }
        basisMatrix = V;
        Info << "####### End of SOPT GappyDEIM #######\n";
      }
      else if (greedyMetric=="SOPTE")
      {
        Info << "####### Begin SOPT-exact GappyDEIM #######" << endl;
        Info << "####### Modes=" <<  n_modes
              << ", nodePoints=" << n_nodes << " #######"
              << endl;
 
        label ind_max, c1;
        double max;
        Eigen::MatrixXd A;
        Eigen::VectorXd b;
        Eigen::VectorXd c;
        Eigen::VectorXd r;
        Eigen::MatrixXd V = snapshotsModes.col(0);
 
        max = (reshapeMat* (mp_not_mask.asDiagonal()*snapshotsModes.leftCols(n_modes)).rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(&ind_max, &c1);
        updateNodes(P, ind_max, mp_not_mask);
 
        for (unsigned int ith_node = 1; ith_node < na; ith_node++)
        {
          if (ith_node<n_modes)
          {
            A = P.transpose() * V;
            b = P.transpose() * snapshotsModes.col(ith_node);
            c = A.fullPivLu().solve(b);
            r = snapshotsModes.col(ith_node) - V * c;
            V.conservativeResize(snapshotsModes.rows(), ith_node + 1);
            V.col(ith_node) =  snapshotsModes.col(ith_node);
 
            max = (reshapeMat*r.rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(&ind_max, &c1);
          }
          else
          {
            // TODO optimize
            Eigen::MatrixXd tmp = P.transpose()*snapshotsModes.leftCols(n_modes);
            Info << "SOPT: " << s_optimality(tmp) << endl;
            tmp.conservativeResize(tmp.rows()+vectorial_dim, tmp.cols());
            
            Eigen::MatrixXd masked = mp_not_mask.asDiagonal()*snapshotsModes.leftCols(n_modes);
 
            Eigen::VectorXd results(n_cells);
 
            for (unsigned int ith_cell = 0; ith_cell < n_cells; ith_cell++)
            {
              for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
              {
                tmp.row(tmp.rows()-vectorial_dim+ith_field) = masked.row(ith_cell +  ith_field*n_cells);
              }
              results(ith_cell) = s_optimality(tmp);
            }
            max = results.maxCoeff(&ind_max, &c1);
          }
 
          updateNodes(P, ind_max, mp_not_mask);
        }
        basisMatrix = snapshotsModes.leftCols(n_modes);
        Info << "####### End of SOPTE GappyDEIM #######\n";
      }
    }
    else
    {
      basisMatrix = snapshotsModes.leftCols(n_modes);
      Info << "####### InitialSeeds equal to n_nodes #######\n";
    }
 
    evaluatePinv(P, basisMatrix, normalizingWeights);
    renormalizedBasisMatrix = normalizingWeights.asDiagonal()*basisMatrix;
 
    MatrixOnline = renormalizedBasisMatrix * pinvPU;
 
    P.makeCompressed();
    cnpy::save(basisMatrix, folderMethod +"/basisMatrix.npy");
    cnpy::save(P, folderMethod +"/projectionMatrix.npz");
    cnpy::save(normalizingWeights, folderMethod + "/normalizingWeights.npy");
    cnpy::save(nodes, folderMethod + "/mp.npy");
    cnpy::save(pinvPU, folderMethod + "/pinvPU.npy");
    cnpy::save(MatrixOnline, folderMethod + "/MatrixOnline.npy");
    nodePoints().write();
 
    Info << "Projection Matrix shape: " << P.rows() << " " << P.cols() << endl;
    Info << "Basis Matrix shape: " << basisMatrix.rows() << " " << basisMatrix.cols() << endl;
    Info << "Pseudo Inverse Matrix shape: " << pinvPU.rows() << " " << pinvPU.cols() << endl;
  }
  else
  {
    Info << "Read GappyDEIM\n";
    cnpy::load(basisMatrix, folderMethod +"/basisMatrix.npy");
    cnpy::load(MatrixOnline, folderMethod +"/MatrixOnline.npy");
    cnpy::load(normalizingWeights, folderMethod + "/normalizingWeights.npy");
    cnpy::load(P, folderMethod +"/projectionMatrix.npz");
    cnpy::load(nodes, folderMethod + "/mp.npy");
    cnpy::load(pinvPU, folderMethod + "/pinvPU.npy");
 
    renormalizedBasisMatrix = normalizingWeights.asDiagonal()*basisMatrix;
    
    Info << "Projection Matrix shape: " << P.rows() << " " << P.cols() << endl;
    Info << "Basis Matrix shape: " << basisMatrix.rows() << " " << basisMatrix.cols() << endl;
    Info << "Pseudo Inverse Matrix shape: " << pinvPU.rows() << " " << pinvPU.cols() << endl;
  }
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::offlineECP(Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& weights, word folderMethodName)
{
  folderMethod = folderMethodName;
  Info << "FolderMethod : " << folderMethod << endl;
  mkDir(folderMethod);
 
  n_nodes = n_modes + 1;
 
  normalizingWeights = weights;
 
  Info << "####### Begin ECP #######" << endl;
  Info << "####### Modes=" <<  n_modes
            << ", nodePoints=" << n_nodes << " #######"
            << endl;
 
  nodePoints = autoPtr<IOList<label>>(new IOList<label>(
      IOobject("nodePoints", para->runTime.time().constant(), "../" + folderMethod,
                para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE))); 
  
  bool offlineStage = !nodePoints().headerOk();
 
  if (offlineStage)
  {
    assert(n_modes > 0);
    assert(n_nodes >= n_modes);
 
    // matrices for quadratureWeights evaluation
    Eigen::MatrixXd J;//shape: (vectorial_dim * (n_modes + 1), nodePoints->size())
    Eigen::MatrixXd Jwhole(vectorial_dim * (n_modes + 1), n_cells);
    Eigen::VectorXd q(vectorial_dim * (n_modes + 1));
 
    // matrices for greedy selection of the nodes
    Eigen::MatrixXd A(vectorial_dim*n_modes, n_cells);
    Eigen::VectorXd b = Eigen::VectorXd::Constant(vectorial_dim*n_modes, 1);
 
    Eigen::VectorXd volumes = ITHACAutilities::getMassMatrixFV(std::get<0>(snapshotsListTuple)[0]);
    double volume = volumes.array().sum();
    for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
    {
      Eigen::MatrixXd block = snapshotsModes.block(ith_field*n_cells, 0, n_cells, n_modes).transpose();
 
      q.segment(ith_field*(n_modes+1), n_modes) = block.rowwise().sum();
      q(ith_field*(n_modes+1)+n_modes) = volume;
 
      Jwhole.middleRows(ith_field*(n_modes+1), n_modes) = block;
      Jwhole.row(ith_field*(n_modes+1)+n_modes) = Eigen::VectorXd::Constant(n_cells, 1);
 
      Eigen::VectorXd mean = block.rowwise().mean();
      block.colwise() -= mean;
      Eigen::VectorXd Anorm = block.colwise().lpNorm<2>();
      if (Anorm.minCoeff()> 0)
      {
        block.array().rowwise() /= Anorm.transpose().array();
      }
 
      A.middleRows(ith_field*n_modes, n_modes) = block;
    }
    Eigen::VectorXd mp_not_mask = Eigen::VectorXd::Constant(n_cells * vectorial_dim, 1);
    std::set<label> nodePointsSet;
 
    // set initialSeeds
    if (initialSeeds.rows() > 0)
    {
      initialSeeds(mp_not_mask, nodePointsSet);
      computeLS(J, Jwhole, b, q);
    }
 
    int na = n_nodes - initialSeeds.rows();
    Eigen::SparseMatrix<double> reshapeMat;
    initReshapeMat(reshapeMat);
 
    if (na > 0)
    {
      label ind_max, r1;
      double max;
 
      for (unsigned int ith_node = 0; ith_node < na; ith_node++)
      {
        max = ((b.transpose()*A*mp_not_mask.head(n_cells).asDiagonal()).colwise().sum()).maxCoeff(&r1, &ind_max);
 
        updateNodes(P, ind_max, mp_not_mask);
        computeLS(J, Jwhole, b, q);
      }
    }
 
    Info << "####### End ECP #######" << endl;
 
    basisMatrix = snapshotsModes.leftCols(n_modes);
    evaluateWPU(P, basisMatrix, normalizingWeights, quadratureWeights);
 
    cnpy::save(basisMatrix, folderMethod +"/basisMatrix.npy");
    cnpy::save(quadratureWeights, folderMethod +"/quadratureWeights.npy");
    cnpy::save(normalizingWeights, folderMethod +"/normalizingWeights.npy");
    cnpy::save(nodes, folderMethod + "/mp.npy");
    cnpy::save(wPU, folderMethod + "/wPU.npy");
    nodePoints().write();
  }
  else
  {
    Info << "Read ECP\n";
    cnpy::load(basisMatrix, folderMethod +"/basisMatrix.npy");
    cnpy::load(quadratureWeights, folderMethod +"/quadratureWeights.npy");
    cnpy::load(normalizingWeights, folderMethod +"/normalizingWeights.npy");
    cnpy::load(nodes, folderMethod + "/mp.npy");
    cnpy::load(wPU, folderMethod + "/wPU.npy");
  }
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::initSeeds(Eigen::VectorXd mp_not_mask, std::set<label> nodePointsSet)
{
  P.resize(n_cells * vectorial_dim, initialSeeds.rows() * vectorial_dim);
  P.reserve(Eigen::VectorXi::Constant(initialSeeds.rows() * vectorial_dim /*n_cols*/, 1 /*n_non_zero_elements*/));
 
  unsigned int trueSeeds{0};
  for (int i = 0; i < initialSeeds.rows(); i++)
  {
    unsigned int index = initialSeeds(i) % n_cells;
    // check that there are not repeated nodes in initialSeeds
    if (nodePointsSet.find(index) == nodePointsSet.end()){
      nodePointsSet.insert(index);
      nodePoints->append(index);
 
      for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
      {
        P.insert(index+ith_field*n_cells, vectorial_dim*trueSeeds+ith_field) = 1;
        mp_not_mask(index + ith_field * n_cells) = 0;
      }
      trueSeeds++;
    }
  }
 
  P.conservativeResize(P.rows(), nodePoints->size()*vectorial_dim);
  P.makeCompressed();
  M_Assert(nodePoints->size()<=n_nodes, "Size of 'initialSeeds' is greater than 'n_nodes'");
}
 
template <typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::updateNodes(Eigen::SparseMatrix<double>& P, label& ind, Eigen::VectorXd& mp_not_mask)
{
  nodePoints->append(ind);
 
  nodes.conservativeResize(nodes.rows()+1);
  nodes(nodes.rows()-1) = ind;
 
  unsigned int last_col{0};
 
  if (P.rows()==0)
  {
    P.resize(n_cells*vectorial_dim, vectorial_dim);
  }
  else
  {
    last_col = P.cols();
    P.resize(P.rows(), P.cols() + vectorial_dim);
  }
  int step = P.cols() / vectorial_dim - 1;
  
  for (int ith_node=0; ith_node <= step; ith_node++)
  {
    for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
    {
      P.insert(nodes[ith_node]+ith_field*n_cells,ith_node+ith_field*(step+1))=1;
    }
  }
 
  for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
  {
    mp_not_mask(ind + ith_field * n_cells) = 0;
  }
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::computeLS(Eigen::MatrixXd& J, Eigen::MatrixXd& Jwhole, Eigen::VectorXd& b, Eigen::VectorXd& q)
{
  // TODO fix with only positive weights (NNLS)
  J.resize(vectorial_dim * (n_modes + 1), nodePoints->size());
  J = Jwhole(Eigen::indexing::all, nodes);
 
  quadratureWeights.resize(nodePoints->size()*vectorial_dim);
  for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
  {
    quadratureWeights.segment(ith_field*nodePoints->size(), nodePoints->size()) = J.middleRows(ith_field*(n_modes+1), n_modes+1).colPivHouseholderQr().solve(q.segment(ith_field*(n_modes+1), n_modes+1));
    b.segment(ith_field*n_modes, n_modes) = q.segment(ith_field*(n_modes+1), n_modes)-J.middleRows(ith_field*(n_modes+1), n_modes)*quadratureWeights.segment(ith_field*nodePoints->size(), nodePoints->size());
  }
  b = b/b.lpNorm<2>();
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::initReshapeMat(Eigen::SparseMatrix<double> &reshapeMat)
{
  reshapeMat.resize(n_cells, n_cells*vectorial_dim);
  reshapeMat.reserve(Eigen::VectorXi::Constant(n_cells * vectorial_dim /*n_cols*/, 1 /*n_non_zero_elements*/));
 
  for (unsigned int i = 0; i < n_cells; i++)
  {
    for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
    {
      reshapeMat.insert(i, i+n_cells*ith_field) = 1;
    }
  }
  reshapeMat.makeCompressed();
}
 
template <typename... SnapshotsLists>
template <typename SnapshotsList>
void HyperReduction<SnapshotsLists...>::stackSnapshots(SnapshotsList sList, Eigen::MatrixXd& snapshotsMatrix, Eigen::VectorXd& fieldWeights)
{
  unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
 
  Eigen::MatrixXd tmpSnapshots = Foam2Eigen::PtrList2Eigen(sList);
  // get volumes
  Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(sList[0]);
 
 
  double maxVal = std::sqrt(tmpSnapshots.colwise().lpNorm<2>().maxCoeff());
 
  fieldWeights.conservativeResize(fieldWeights.rows()+field_dim*n_cells);
  fieldWeights.tail(n_cells * field_dim) = V.array().sqrt().cwiseInverse()*maxVal;
 
  snapshotsMatrix.conservativeResize(snapshotsMatrix.rows() + n_cells * field_dim, n_snapshots);
  snapshotsMatrix.bottomRows(n_cells * field_dim) = tmpSnapshots;
}
 
template <typename... SnapshotsLists>
template <typename SnapshotsList>
void HyperReduction<SnapshotsLists...>::stackSnapshotsBoundary(SnapshotsList sList, List<Eigen::MatrixXd>& snapshotsMatrixBoundary, List<Eigen::VectorXd>& fieldWeightsBoundary)
{
  unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
 
  List<Eigen::MatrixXd> tmpBoundarySnapshots = Foam2Eigen::PtrList2EigenBC(sList);
 
  List<double> maxVal;
  maxVal.resize(n_boundary_patches);
  
  for (label id = 0; id < n_boundary_patches; id++)
  {
    maxVal[id] = std::sqrt(tmpBoundarySnapshots[id].colwise().lpNorm<2>().maxCoeff());
 
    // get surface areas
    Eigen::VectorXd S = Foam2Eigen::field2Eigen(sList[0].mesh().magSf().boundaryField()[id]);
    S = S.replicate(field_dim, 1);
 
    unsigned int bSize = field_dim*int(n_boundary_cells_list[id]);
    fieldWeightsBoundary[id].conservativeResize(fieldWeightsBoundary[id].rows()+bSize);
    fieldWeightsBoundary[id].tail(bSize) = S.array().pow(3/4)*maxVal[id];
 
    snapshotsMatrixBoundary[id].conservativeResize(snapshotsMatrixBoundary[id].rows() + bSize, n_snapshots);
    snapshotsMatrixBoundary[id].bottomRows(bSize) = tmpBoundarySnapshots[id];
  }
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::evaluatePinv(Eigen::SparseMatrix<double> &Projector, Eigen::MatrixXd &Modes, Eigen::VectorXd &fieldWeights) {
  assert(Projector.cols() > 0);
  Eigen::MatrixXd restricted = Projector.transpose() * Modes;
  Info << "S-Optimalty: " << s_optimality(restricted) << endl;
  pinvPU = ITHACAutilities::invertMatrix(restricted, para->ITHACAdict->lookupOrDefault<word>("InversionMethod", "completeOrthogonalDecomposition"));
  pinvPU = pinvPU*(Projector.transpose() * fieldWeights.array().cwiseInverse().matrix()).asDiagonal();
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::evaluateWPU(Eigen::SparseMatrix<double> &Projector, Eigen::MatrixXd &Modes, Eigen::VectorXd &fieldWeights, Eigen::VectorXd& quadratureWeights) {
  assert(Projector.cols() > 0);
 
  Eigen::VectorXd quadratureWeightsOrderedAsProjector(quadratureWeights.rows());
  unsigned int n_weightsPerField = quadratureWeights.rows()/vectorial_dim;
  unsigned int reorderIndex{0};
 
  for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
  {  
    for (unsigned int ith_weight = 0; ith_weight < n_weightsPerField; ith_weight++)
    {
      quadratureWeightsOrderedAsProjector(reorderIndex) = quadratureWeights(ith_weight + ith_field * n_weightsPerField);
      reorderIndex++;
    }
  }
  wPU = quadratureWeightsOrderedAsProjector.transpose()*(Projector.transpose() * fieldWeights.array().cwiseInverse().matrix()).asDiagonal();
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::generateSubmesh(label layers, const fvMesh &mesh) {
 
  Info << "####### Extract submesh #######\n";
  ITHACAparameters *para(ITHACAparameters::getInstance());
 
  totalNodePoints = autoPtr<IOList<labelList>>(new IOList<labelList>(IOobject(
      "totalNodePoints", para->runTime.time().constant(),  "../" + folderMethod,
      para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)));
 
  uniqueNodePoints = autoPtr<IOList<label>>(new IOList<label>(IOobject(
      "uniqueNodePoints", para->runTime.time().constant(),  "../" + folderMethod,
      para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)));
 
  volScalarField Indici
  (
      IOobject
      (
        problemName + "_indices",
        mesh.time().timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
      ),
      mesh,
      dimensionedScalar(problemName + "_indices",
                        dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0));
 
  submesh = autoPtr<fvMeshSubset>(new fvMeshSubset(mesh));
  bool offlineStage = !totalNodePoints().headerOk();
 
  if (offlineStage) {
    List<label> indices;
    for (label i = 0; i < nodePoints().size(); i++) {
      indices = ITHACAutilities::getIndices(mesh, nodePoints()[i], layers);
      totalNodePoints().append(indices);
    }
 
    labelList a = ListListOps::combine<labelList>(totalNodePoints(),
                                                  accessOp<labelList>());
 
    inplaceUniqueSort(a);
    uniqueNodePoints() = a;
 
    scalar zerodot25 = 0.25;
    ITHACAutilities::assignIF(Indici, zerodot25,
                              uniqueNodePoints().List<label>::clone()());
    ITHACAutilities::assignONE(Indici, nodePoints());
 
    totalNodePoints().write();
    uniqueNodePoints().write();
    ITHACAstream::exportSolution(Indici, "1", folderMethod);
  }
 
  submesh->setCellSubset(uniqueNodePoints());
  submesh->subMesh().fvSchemes::readOpt() = mesh.fvSchemes::readOpt();
  submesh->subMesh().fvSolution::readOpt() = mesh.fvSolution::readOpt();
  submesh->subMesh().fvSchemes::read();
  submesh->subMesh().fvSolution::read();
  std::cout.clear();
 
  localNodePoints = global2local(nodePoints(), submesh());
 
  n_cellsSubfields = submesh().cellMap().size();
  Info << "####### End extract submesh size = " << n_cellsSubfields <<" #######\n";
 
  createMasks(offlineStage);
  Info << "####### End create masks #######\n";
}
 
template<typename... SnapshotsLists>
void HyperReduction<SnapshotsLists...>::createMasks(bool offlineStage)
{ 
  if (offlineStage)
  {
    field2submesh.resize(submesh().cellMap().size() * vectorial_dim, n_cells*vectorial_dim);
    field2submesh.reserve(Eigen::VectorXi::Constant(n_cells*vectorial_dim, 1));
 
    for (unsigned int ith_subCell{0} ; ith_subCell <submesh().cellMap().size(); ith_subCell++)
    {
        for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
        {
          field2submesh.insert(ith_subCell+ith_field*submesh().cellMap().size(), submesh().cellMap()[ith_subCell]+n_cells*ith_field) = 1;
        }
    }
    field2submesh.makeCompressed();
 
    submesh2nodes.resize(nodePoints().size() * vectorial_dim, submesh().cellMap().size()*vectorial_dim);
    submesh2nodes.reserve(Eigen::VectorXi::Constant(submesh().cellMap().size()*vectorial_dim, 1));
    submesh2nodesMask.resize(nodePoints().size() * vectorial_dim);
 
    int index_col = 0;
    for (unsigned int ith_node{0} ; ith_node <nodePoints().size(); ith_node++) {
      for (unsigned int ith_subCell{0} ; ith_subCell <submesh().cellMap().size(); ith_subCell++) {
        if (nodePoints()[ith_node] == submesh().cellMap()[ith_subCell]) {
          for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
          {
            submesh2nodes.insert(ith_node + nodePoints().size()*ith_field, index_col + ith_field*submesh().cellMap().size()) = 1;
            submesh2nodesMask(ith_node + nodePoints().size()*ith_field) = index_col + ith_field*submesh().cellMap().size();
          }
          break;
        } else {
          index_col++;
        }
      }
      index_col = 0;
    }
    submesh2nodes.makeCompressed();
 
    cnpy::save(field2submesh, folderMethod +"/field2submesh.npz");
    cnpy::save(submesh2nodes, folderMethod +"/submesh2nodes.npz");
    cnpy::save(submesh2nodesMask, folderMethod +"/submesh2nodesMask.npy");
  }
  else
  {
    cnpy::load(field2submesh, folderMethod +"/field2submesh.npz");
    cnpy::load(submesh2nodes, folderMethod +"/submesh2nodes.npz");
    cnpy::load(submesh2nodesMask, folderMethod +"/submesh2nodesMask.npy");
  }
}
 
template<typename... SnapshotsLists>
List<label> HyperReduction<SnapshotsLists...>::global2local(List<label> &points, fvMeshSubset &submesh) {
  List<label> localPoints;
 
  for (label i = 0; i < points.size(); i++) {
    for (label j = 0; j < submesh.cellMap().size(); j++) {
      if (submesh.cellMap()[j] == points[i]) {
        localPoints.append(j);
        break;
      }
    }
  }
 
  return localPoints;
}
 
#endif