gm_to_dirac_even_binf_optimization_params.h

#ifndef GM_TO_DIRAC_EVEN_BINF_OPTIMIZATION_PARAMS_H
#define GM_TO_DIRAC_EVEN_BINF_OPTIMIZATION_PARAMS_H
 
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
 
#include <cassert>
#include <vector>
 
#include "gsl_minimizer.h"
#include "lcd_even_binf_closed_form.h"
#include "squared_euclidean_distance_utils.h"
 
struct GMToDiracEvenBinfOptimizationParams
    : public GslMinimizerOptimizationParams {
 public:
  GMToDiracEvenBinfOptimizationParams(const gsl_vector* wX, size_t N, size_t L)
      : GslMinimizerOptimizationParams(L, N),
        wX(wX),
        k(N / 2),
        constant(lcd_binf_constant(N / 2)),
        sProjected(L * N, 0.00),
        cSqrdNorm(L, 0.00) {
    assert(wX != nullptr);
    assert(wX->size == L);
    assert(N % 2 == 0);
    assert(N >= 2);
    assert(L >= 1);
 
    if (L > 1)
      squaredEuclideanDistanceUtilLL =
          new SquaredEuclideanDistance_LL_vectorized(L, N);
  }
 
  ~GMToDiracEvenBinfOptimizationParams() {
    if (squaredEuclideanDistanceUtilLL) delete squaredEuclideanDistanceUtilLL;
  }
 
  inline void update(const gsl_vector* t) {
    assert(t->size == L * N);
 
    // s = t - 1 * (w^T t)
    std::vector<double> mean(N, 0.00);
    for (size_t i = 0; i < L; ++i) {
      const double wXi = wX->data[i];
      for (size_t d = 0; d < N; ++d) mean[d] += wXi * t->data[i * N + d];
    }
 
    for (size_t i = 0; i < L; ++i) {
      double sum = 0.00;
      for (size_t d = 0; d < N; ++d) {
        const double sid = t->data[i * N + d] - mean[d];
        sProjected[i * N + d] = sid;
        sum += sid * sid;
      }
      cSqrdNorm[i] = sum;
    }
 
    if (!squaredEuclideanDistanceUtilLL) return;
 
    const gsl_matrix sMatrix =
        gsl_matrix_view_array(sProjected.data(), L, N).matrix;
    squaredEuclideanDistanceUtilLL->calculateDistance(&sMatrix, nullptr);
  }
 
  inline bool hasDegenerateSample() const {
    for (size_t i = 0; i < L; ++i)
      if (!(cSqrdNorm[i] > 0.00)) return true;
    return false;
  }
 
  inline double sample(size_t i, size_t d) const {
    return sProjected[i * N + d];
  }
 
  inline double distanceSq(size_t i, size_t j) const {
    if (!squaredEuclideanDistanceUtilLL) return 0.00;  // L == 1
    return squaredEuclideanDistanceUtilLL->getDistance(i, j);
  }
 
  const gsl_vector* wX;
  const size_t k;         
  const double constant;  
 
  std::vector<double> sProjected;  
  std::vector<double> cSqrdNorm;   
  SquaredEuclideanDistanceUtilsLL* squaredEuclideanDistanceUtilLL = nullptr;
};
 
#endif  // GM_TO_DIRAC_EVEN_BINF_OPTIMIZATION_PARAMS_H