<commit_before>/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkSymmetricEigenAnalysis_h
#define itkSymmetricEigenAnalysis_h

#include "itkMacro.h"
#include "itk_eigen.h"
#include ITK_EIGEN(Eigenvalues)
#include <numeric>
#include <vector>
// For GetPointerToMatrixData
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_matrix_fixed.h"
#include "itkMatrix.h"

namespace itk
{
namespace detail
{
  /* Helper functions returning pointer to matrix data for different types.  */
  template<typename TValueType, unsigned int NRows, unsigned int NCols>
    const TValueType* GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & inputMatrix)
    {
      return inputMatrix.data_block();
    };
  template<typename TValueType>
    const TValueType* GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)
    {
      return inputMatrix.data_block();
    };

  template<typename TValueType, unsigned int NRows, unsigned int NCols>
    const TValueType* GetPointerToMatrixData(const itk::Matrix<TValueType, NRows, NCols> & inputMatrix)
    {
      return inputMatrix.GetVnlMatrix().data_block();
    };

  /** Sort input to be ordered by magnitude, and returns container with the
   * permutations required for the sorting.
   *
   * For example, if input eigenValues = {10, 0, 40}, the output would be: {2,0,1}
   * and the eigenValues would be modified in-place: {40, 10, 0}.
   *
   * The permutations indices is used to order the matrix of eigenVectors.
   * \sa permuteEigenVectorsWithSortPermutations
   *
   * @tparam TArray  array type with operator []
   * @param eigenValues input array, requires operator []
   * @param numberOfElements size of array
   *
   * @return the permutations needed to sort the input array
   */
  template<typename TArray>
    std::vector<int> sortEigenValuesByMagnitude(TArray & eigenValues, const unsigned int numberOfElements)
    {
      std::vector<int> indicesSortPermutations(numberOfElements, 0);
      std::iota(std::begin(indicesSortPermutations), std::end(indicesSortPermutations), 0);

      std::sort(
          std::begin(indicesSortPermutations),
          std::end(indicesSortPermutations),
          [&eigenValues](unsigned int a, unsigned int b)
          { return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); }
          );
      auto tmpCopy = eigenValues;
      for (unsigned int i = 0; i < numberOfElements; ++i)
      {
        eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
      }
      return indicesSortPermutations;
    }

  /** Permute a eigenVectors matrix according to the permutation indices
   * computed from the output of a sort function like \sa detail::sortEigenValuesByMagnitude
   *
   * @tparam QMatrix a Eigen3 matrix
   * @param eigenVectors stored in columns
   * @param indicesSortPermutations container with the permutations from the output of
   * a sort function.
   */
  template<typename QMatrix>
    void permuteColumnsWithSortIndices( QMatrix & eigenVectors,
        const std::vector<int> & indicesSortPermutations)
    {
      using EigenLibPermutationMatrix =
        Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
      auto numberOfElements = indicesSortPermutations.size();
      // Creates a NxN permutation matrix copying our permutation to the matrix indices.
      // Which holds the 1D array representation of a permutation.
      EigenLibPermutationMatrix perm(numberOfElements);
      perm.setIdentity();
      std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(),
          perm.indices().data());
      // Apply it
      eigenVectors = eigenVectors * perm;
    }
} // end namespace detail

/** \class SymmetricEigenAnalysis
 * \brief Find Eigen values of a real 2D symmetric matrix. It
 * serves as a thread-safe alternative to the class:
 * vnl_symmetric_eigensystem, which uses netlib routines.
 *
 * The class is templated over the input matrix (which is expected to provide
 * access to its elements with the [][] operator), matrix to store eigen
 * values (must provide write operations on its elements with the [] operator), and
 * EigenMatrix to store eigen vectors (must provide write access to its elements
 * with the [][] operator).
 *
 * The SetOrderEigenValues() method can be used to order eigen values (and their
 * corresponding eigen vectors if computed) in ascending order. This is the
 * default ordering scheme. Eigen vectors and values can be obtained without
 * ordering by calling SetOrderEigenValues(false).
 *
 * The SetOrderEigenMagnitudes() method can be used to order eigen values (and
 * their corresponding eigen vectors if computed) by magnitude in ascending order.
 *
 * The user of this class is explicitly supposed to set the dimension of the
 * 2D matrix using the SetDimension() method.
 *
 * The class contains routines taken from netlib sources (www.netlib.org).
 * netlib/tql1.c
 * netlib/tql2.c
 * netlib/tred1.c
 * netlib/tred2.c
 *
 * Reference:
 *     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
 *     wilkinson.
 *     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
 * \ingroup ITKCommon
 */

template< typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix >
class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
{
public:
  typedef enum
  {
    OrderByValue = 1,
    OrderByMagnitude,
    DoNotOrder
  }
  EigenValueOrderType;

  SymmetricEigenAnalysis() :
  m_OrderEigenValues(OrderByValue) {}

  SymmetricEigenAnalysis(const unsigned int dimension):
  m_UseEigenLibrary(false),
  m_Dimension(dimension),
  m_Order(dimension),
  m_OrderEigenValues(OrderByValue) {}

  ~SymmetricEigenAnalysis() = default;

  using MatrixType = TMatrix;
  using EigenMatrixType = TEigenMatrix;
  using VectorType = TVector;

  /** Compute Eigen values of A
   * A is any type that overloads the [][] operator and contains the
   * symmetric matrix. In practice only the upper triangle of the
   * matrix will be accessed. (Both itk::Matrix and vnl_matrix
   * overload [][] operator.)
   *
   * 'EigenValues' is any type that overloads the [][] operator and will contain<commit_msg>STYLE: Prevent double initialization<commit_after>/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkSymmetricEigenAnalysis_h
#define itkSymmetricEigenAnalysis_h

#include "itkMacro.h"
#include "itk_eigen.h"
#include ITK_EIGEN(Eigenvalues)
#include <numeric>
#include <vector>
// For GetPointerToMatrixData
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_matrix_fixed.h"
#include "itkMatrix.h"

namespace itk
{
namespace detail
{
  /* Helper functions returning pointer to matrix data for different types.  */
  template<typename TValueType, unsigned int NRows, unsigned int NCols>
    const TValueType* GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & inputMatrix)
    {
      return inputMatrix.data_block();
    };
  template<typename TValueType>
    const TValueType* GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)
    {
      return inputMatrix.data_block();
    };

  template<typename TValueType, unsigned int NRows, unsigned int NCols>
    const TValueType* GetPointerToMatrixData(const itk::Matrix<TValueType, NRows, NCols> & inputMatrix)
    {
      return inputMatrix.GetVnlMatrix().data_block();
    };

  /** Sort input to be ordered by magnitude, and returns container with the
   * permutations required for the sorting.
   *
   * For example, if input eigenValues = {10, 0, 40}, the output would be: {2,0,1}
   * and the eigenValues would be modified in-place: {40, 10, 0}.
   *
   * The permutations indices is used to order the matrix of eigenVectors.
   * \sa permuteEigenVectorsWithSortPermutations
   *
   * @tparam TArray  array type with operator []
   * @param eigenValues input array, requires operator []
   * @param numberOfElements size of array
   *
   * @return the permutations needed to sort the input array
   */
  template<typename TArray>
    std::vector<int> sortEigenValuesByMagnitude(TArray & eigenValues, const unsigned int numberOfElements)
    {
      std::vector<int> indicesSortPermutations(numberOfElements, 0);
      std::iota(std::begin(indicesSortPermutations), std::end(indicesSortPermutations), 0);

      std::sort(
          std::begin(indicesSortPermutations),
          std::end(indicesSortPermutations),
          [&eigenValues](unsigned int a, unsigned int b)
          { return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); }
          );
      auto tmpCopy = eigenValues;
      for (unsigned int i = 0; i < numberOfElements; ++i)
      {
        eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
      }
      return indicesSortPermutations;
    }

  /** Permute a eigenVectors matrix according to the permutation indices
   * computed from the output of a sort function like \sa detail::sortEigenValuesByMagnitude
   *
   * @tparam QMatrix a Eigen3 matrix
   * @param eigenVectors stored in columns
   * @param indicesSortPermutations container with the permutations from the output of
   * a sort function.
   */
  template<typename QMatrix>
    void permuteColumnsWithSortIndices( QMatrix & eigenVectors,
        const std::vector<int> & indicesSortPermutations)
    {
      using EigenLibPermutationMatrix =
        Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
      auto numberOfElements = indicesSortPermutations.size();
      // Creates a NxN permutation matrix copying our permutation to the matrix indices.
      // Which holds the 1D array representation of a permutation.
      EigenLibPermutationMatrix perm(numberOfElements);
      perm.setIdentity();
      std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(),
          perm.indices().data());
      // Apply it
      eigenVectors = eigenVectors * perm;
    }
} // end namespace detail

/** \class SymmetricEigenAnalysis
 * \brief Find Eigen values of a real 2D symmetric matrix. It
 * serves as a thread-safe alternative to the class:
 * vnl_symmetric_eigensystem, which uses netlib routines.
 *
 * The class is templated over the input matrix (which is expected to provide
 * access to its elements with the [][] operator), matrix to store eigen
 * values (must provide write operations on its elements with the [] operator), and
 * EigenMatrix to store eigen vectors (must provide write access to its elements
 * with the [][] operator).
 *
 * The SetOrderEigenValues() method can be used to order eigen values (and their
 * corresponding eigen vectors if computed) in ascending order. This is the
 * default ordering scheme. Eigen vectors and values can be obtained without
 * ordering by calling SetOrderEigenValues(false).
 *
 * The SetOrderEigenMagnitudes() method can be used to order eigen values (and
 * their corresponding eigen vectors if computed) by magnitude in ascending order.
 *
 * The user of this class is explicitly supposed to set the dimension of the
 * 2D matrix using the SetDimension() method.
 *
 * The class contains routines taken from netlib sources (www.netlib.org).
 * netlib/tql1.c
 * netlib/tql2.c
 * netlib/tred1.c
 * netlib/tred2.c
 *
 * Reference:
 *     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
 *     wilkinson.
 *     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
 * \ingroup ITKCommon
 */

template< typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix >
class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
{
public:
  typedef enum
  {
    OrderByValue = 1,
    OrderByMagnitude,
    DoNotOrder
  }
  EigenValueOrderType;

  SymmetricEigenAnalysis() :
  m_OrderEigenValues(OrderByValue) {}

  SymmetricEigenAnalysis(const unsigned int dimension):

  m_Dimension(dimension),
  m_Order(dimension),
  m_OrderEigenValues(OrderByValue) {}

  ~SymmetricEigenAnalysis() = default;

  using MatrixType = TMatrix;
  using EigenMatrixType = TEigenMatrix;
  using VectorType = TVector;

  /** Compute Eigen values of A
   * A is any type that overloads the [][] operator and contains the
   * symmetric matrix. In practice only the upper triangle of the
   * matrix will be accessed. (Both itk::Matrix and vnl_matrix
   * overload [][] operator.)
   *
   * 'EigenValues' is any type that overloads the [][] operator and will contain
{'fraction_non_alphanumeric': 0.10422740524781342, 'fraction_numerical': 0.0074344023323615165, 'mean_word_length': 3.9320632638389648, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3135352', 'n_tokens_mistral': 4004, 'n_tokens_neox': 3780, 'n_words': 1309}