gm_to_dirac_even_optimization_params.h

#ifndef GM_TO_DIRAC_EVEN_OPTIMIZATION_PARAMS_H
#define GM_TO_DIRAC_EVEN_OPTIMIZATION_PARAMS_H
 
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
 
#include <cassert>
#include <vector>
 
#include "gsl_minimizer.h"
#include "lcd_even_closed_form.h"
#include "squared_euclidean_distance_utils.h"
 
struct GMToDiracEvenOptimizationParams : public GslMinimizerOptimizationParams {
 public:
  GMToDiracEvenOptimizationParams(const gsl_vector* wX, size_t N, size_t L,
                                  double bMax, bool includeD1)
      : GslMinimizerOptimizationParams(L, N),
        wX(wX),
        bMax(bMax),
        k(N / 2),
        reportedOffset(
            lcd_reported_offset(N, bMax, weightSumOf(wX, L))),
        objectiveOffset(includeD1 ? reportedOffset : 0.00),
        cSqrdNorm(L, 0.00) {
    assert(wX != nullptr);
    assert(wX->size == L);
    assert(N % 2 == 0);
    assert(N >= 2);
    assert(bMax > 0.00);
    assert(L >= 1);
 
    if (L > 1)
      squaredEuclideanDistanceUtilLL =
          new SquaredEuclideanDistance_LL_vectorized(L, N);
  }
 
  ~GMToDiracEvenOptimizationParams() {
    if (squaredEuclideanDistanceUtilLL) delete squaredEuclideanDistanceUtilLL;
  }
 
  inline void update(const gsl_vector* x) {
    assert(x->size == L * N);
 
    for (size_t i = 0; i < L; ++i) {
      double sum = 0.00;
      for (size_t d = 0; d < N; ++d) {
        const double xid = x->data[i * N + d];
        sum += xid * xid;
      }
      cSqrdNorm[i] = sum;
    }
 
    if (!squaredEuclideanDistanceUtilLL) return;
 
    const gsl_matrix xMatrix = gsl_matrix_view_array(x->data, L, N).matrix;
    squaredEuclideanDistanceUtilLL->calculateDistance(&xMatrix, nullptr);
  }
 
  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 double bMax;
  const size_t k;  
 
  const double reportedOffset;
 
  const double objectiveOffset;
 
  std::vector<double> cSqrdNorm;   
  SquaredEuclideanDistanceUtilsLL* squaredEuclideanDistanceUtilLL = nullptr;
 
 private:
  static double weightSumOf(const gsl_vector* wX, size_t L) {
    double sum = 0.00;
    for (size_t i = 0; i < L; ++i) sum += wX->data[i];
    return sum;
  }
};
 
#endif  // GM_TO_DIRAC_EVEN_OPTIMIZATION_PARAMS_H