Spacing.hpp

Go to the documentation of this file.

#pragma once
 
#include "janus/core/JanusConcepts.hpp"
#include "janus/core/JanusError.hpp"
#include "janus/core/JanusTypes.hpp"
#include "janus/math/Arithmetic.hpp"
#include "janus/math/Trig.hpp"
#include <Eigen/Dense>
#include <numbers>
 
namespace janus {
 
// --- linspace ---
template <typename T> JanusVector<T> linspace(const T &start, const T &end, int n) {
    if (n < 1) {
        throw InvalidArgument("linspace: n must be >= 1");
    }
    if (n == 1) {
        JanusVector<T> ret(1);
        ret(0) = start;
        return ret;
    }
 
    JanusVector<T> ret(n);
    // Explicit loop to support symbolic types (cannot use .setLinSpaced)
    T step = (end - start) / static_cast<double>(n - 1);
 
    for (int i = 0; i < n; ++i) {
        ret(i) = start + static_cast<double>(i) * step;
    }
    // Ensure exact end point
    ret(n - 1) = end;
    return ret;
}
 
// --- cosine_spacing ---
template <typename T> JanusVector<T> cosine_spacing(const T &start, const T &end, int n) {
    if (n < 1) {
        throw InvalidArgument("cosine_spacing: n must be >= 1");
    }
    if (n == 1) {
        JanusVector<T> ret(1);
        ret(0) = start;
        return ret;
    }
 
    JanusVector<T> ret(n);
    T center = 0.5 * (start + end);
    T radius = 0.5 * (end - start);
    double pi = std::numbers::pi_v<double>;
 
    for (int i = 0; i < n; ++i) {
        double angle = pi * static_cast<double>(i) / static_cast<double>(n - 1);
        ret(i) = center - radius * std::cos(angle);
    }
    return ret;
}
 
// --- sinspace ---
template <typename T>
JanusVector<T> sinspace(const T &start, const T &end, int n, bool reverse_spacing = false) {
    if (n < 1) {
        throw InvalidArgument("sinspace: n must be >= 1");
    }
    if (n == 1) {
        JanusVector<T> ret(1);
        ret(0) = start;
        return ret;
    }
 
    if (reverse_spacing) {
        auto ret = sinspace(end, start, n, false);
        return ret.reverse();
    }
 
    JanusVector<T> ret(n);
    double pi_half = std::numbers::pi_v<double> / 2.0;
 
    for (int i = 0; i < n; ++i) {
        // angle goes from 0 to pi/2
        double angle = pi_half * static_cast<double>(i) / static_cast<double>(n - 1);
        double factor = 1.0 - std::cos(angle); // 0 to 1
        ret(i) = start + (end - start) * factor;
    }
    // Ensure exact endpoints
    ret(0) = start;   // cos(0) = 1, factor=0
    ret(n - 1) = end; // cos(pi/2) = 0, factor=1
    return ret;
}
 
// --- logspace ---
template <typename T> JanusVector<T> logspace(const T &start, const T &end, int n) {
    if (n < 1) {
        throw InvalidArgument("logspace: n must be >= 1");
    }
    if (n == 1) {
        JanusVector<T> ret(1);
        ret(0) = janus::pow(10.0, start);
        return ret;
    }
 
    JanusVector<T> ret(n);
    // Linear spacing in exponent
    for (int i = 0; i < n; ++i) {
        double fraction = static_cast<double>(i) / static_cast<double>(n - 1);
        T exponents = start + (end - start) * fraction;
        ret(i) = janus::pow(10.0, exponents);
    }
    // ret(0) is 10^start, ret(n-1) is 10^end
    return ret;
}
 
// --- geomspace ---
template <typename T> JanusVector<T> geomspace(const T &start, const T &end, int n) {
    if (n < 1) {
        throw InvalidArgument("geomspace: n must be >= 1");
    }
    T log_start = janus::log10(start);
    T log_end = janus::log10(end);
    return logspace(log_start, log_end, n);
}
 
} // namespace janus