least_squares_fit.hpp

Go to the documentation of this file.

#pragma once
 
#include <array>
#include <ranges>
#include <span>
#include <stdexcept>
#include <vector>
 
#include <Eigen/Dense>
 
#include "interpolation.hpp"
 
namespace kynema::math {
 
inline std::vector<double> MapGeometricLocations(std::span<const double> geom_locations) {
    // Get first and last points of the input domain (assumed to be sorted)
    const double domain_start = geom_locations.front();
    const double domain_end = geom_locations.back();
    if (domain_end == domain_start) {
        throw std::invalid_argument(
            "Invalid geometric locations: domain start and end points are equal."
        );
    }
 
    // Map each point from domain -> [-1, 1]
    std::vector<double> mapped_locations(geom_locations.size());
    const auto domain_span = domain_end - domain_start;
    std::ranges::transform(
        geom_locations, std::begin(mapped_locations),
        [domain_start, domain_span](auto geom_location) {
            return (2. * (geom_location - domain_start) / domain_span) - 1.;
        }
    );
    return mapped_locations;
}
 
inline std::vector<std::vector<double>> ComputeShapeFunctionValues(
    std::span<const double> input_points, std::span<const double> output_points
) {
    // Number of points in input and output arrays
    const auto num_input_points = input_points.size();
    const auto num_output_points = output_points.size();
 
    // Compute weights for the shape functions and its derivatives at the evaluation points
    std::vector<double> weights(num_output_points, 0.);
 
    // Create shape function interpolation matrix
    auto shape_functions = std::vector<std::vector<double>>(
        num_output_points, std::vector<double>(num_input_points, 0.)
    );
    for (auto input_point : std::views::iota(0U, num_input_points)) {
        math::LagrangePolynomialInterpWeights(input_points[input_point], output_points, weights);
        for (auto output_point : std::views::iota(0U, num_output_points)) {
            shape_functions[output_point][input_point] = weights[output_point];
        }
    }
 
    return shape_functions;
}
 
inline std::vector<std::vector<double>> ComputeShapeFunctionDerivatives(
    std::span<const double> input_points, std::span<const double> output_points
) {
    // Number of points in input and output arrays
    const auto num_input_points = input_points.size();
    const auto num_output_points = output_points.size();
 
    // Compute weights for the shape functions and its derivatives at the evaluation points
    std::vector<double> weights(num_output_points, 0.);
 
    // Create shape function derivative matrix
    auto derivative_functions = std::vector<std::vector<double>>(
        num_output_points, std::vector<double>(num_input_points, 0.)
    );
    for (auto input_point : std::views::iota(0U, num_input_points)) {
        math::LagrangePolynomialDerivWeights(input_points[input_point], output_points, weights);
        for (auto output_point : std::views::iota(0U, num_output_points)) {
            derivative_functions[output_point][input_point] = weights[output_point];
        }
    }
 
    return derivative_functions;
}
 
std::vector<std::array<double, 3>> PerformLeastSquaresFitting(
    std::span<const std::vector<double>> shape_functions,
    std::span<const std::array<double, 3>> points_to_fit
);
}  // namespace kynema::math