text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

#include "timer.h"

void integral_recur (int nmin, int nmax, double vals[]);
void integral_gen (int nmin, int nmax, double vals[]);
double ff (double x, void *params);
int adjust_rep_count (int count, double time, double tmin, double tmax);

double ff (double x, void *params)
{
    double n = *(double *) params;
    double f = exp (-x) * pow (x, n);

    return f;
}

void integral_recur (int nmin, int nmax, double vals[])
{
    int i;

    vals[nmax] = 0.;
    for (i = nmax - 1; i >= nmin; i--)
    {
        vals[i] = vals[i + 1] / (i + 1) + 1. / ((i + 1) * M_E);
    }
}

void integral_gen (int nmin, int nmax, double vals[])
{
    int i;
    double result, error;
    double a = 0., b = 1.;
    double abserr = 1.e-9, relerr = 1.e-9;
    double n;
    size_t np = 1000;
    gsl_integration_workspace *w = gsl_integration_workspace_alloc (np);

    gsl_function F;

    F.function = &ff;

    for (i = nmax; i >= nmin; i--)
    {
        n = (double) i;
        F.params = &n;

        gsl_integration_qag (&F, a, b, abserr, relerr, np, GSL_INTEG_GAUSS15,
            w, &result, &error);

        vals[i] = result;
    }

    gsl_integration_workspace_free (w);

}

#define N 100

int main (void)
{
    int i;
    double integral1[N + 1], integral2[N + 1];
    int nmax = N;
    int nmin = 10;

    integral_recur (nmin, nmax, integral1);
    integral_gen (nmin, nmax, integral2);

    for (i = nmin; i <= nmax; i++)
    {
        printf ("%5d    % 14.6f   % 14.6f   %g\n", i, integral1[i],
            integral2[i], fabs (integral1[i] - integral2[i]));
    }

    double time, time1, time2, tmin = 1., tmax = 2.;
    int count = 1000;

    do
    {

        timer_start ();
        for (i = 0; i < count; i++)
        {
            integral_recur (nmin, nmax, integral1);
        }
        time = timer_stop ();

        time1 = time / count;
        printf ("Recur: %10.2f usec %10.6f sec %10d\n", time1 * 1.e6, time,
            count);
        /*
         * adjust count such that cpu time is between
         * tmin and tmax
         */
        count = adjust_rep_count (count, time, tmin, tmax);
    }
    while ((time > tmax) || (time < tmin));

    count = 1000;

    do
    {

        timer_start ();
        for (i = 0; i < count; i++)
        {
            integral_gen (nmin, nmax, integral2);
        }
        time = timer_stop ();

        time2 = time / count;
        printf ("Integ: %10.2f usec %10.6f sec %10d\n", time1 * 1.e6, time,
            count);
        /*
         * adjust count such that cpu time is between
         * tmin and tmax
         */
        count = adjust_rep_count (count, time, tmin, tmax);
    }
    while ((time > tmax) || (time < tmin));

    printf ("Speedup: %.2f\n", time2 / time1);

    return 0;
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_spmatrix.h>

/* parameters for functions */
struct model_params
{
  double alpha;
  gsl_spmatrix *J;
};

/* penalty function */
int
penalty_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  struct model_params *par = (struct model_params *) params;
  const double sqrt_alpha = sqrt(par->alpha);
  const size_t p = x->size;
  size_t i;
  double sum = 0.0;

  for (i = 0; i < p; ++i)
    {
      double xi = gsl_vector_get(x, i);

      gsl_vector_set(f, i, sqrt_alpha*(xi - 1.0));

      sum += xi * xi;
    }

  gsl_vector_set(f, p, sum - 0.25);

  return GSL_SUCCESS;
}

int
penalty_df (CBLAS_TRANSPOSE_t TransJ, const gsl_vector * x,
            const gsl_vector * u, void * params, gsl_vector * v,
            gsl_matrix * JTJ)
{
  struct model_params *par = (struct model_params *) params;
  const size_t p = x->size;
  size_t j;

  /* store 2*x in last row of J */
  for (j = 0; j < p; ++j)
    {
      double xj = gsl_vector_get(x, j);
      gsl_spmatrix_set(par->J, p, j, 2.0 * xj);
    }

  /* compute v = op(J) u */
  if (v)
    gsl_spblas_dgemv(TransJ, 1.0, par->J, u, 0.0, v);

  if (JTJ)
    {
      gsl_vector_view diag = gsl_matrix_diagonal(JTJ);

      /* compute J^T J = [ alpha*I_p + 4 x x^T ] */
      gsl_matrix_set_zero(JTJ);

      /* store 4 x x^T in lower half of JTJ */
      gsl_blas_dsyr(CblasLower, 4.0, x, JTJ);

      /* add alpha to diag(JTJ) */
      gsl_vector_add_constant(&diag.vector, par->alpha);
    }

  return GSL_SUCCESS;
}

int
penalty_fvv (const gsl_vector * x, const gsl_vector * v,
             void *params, gsl_vector * fvv)
{
  const size_t p = x->size;
  double normv = gsl_blas_dnrm2(v);

  gsl_vector_set_zero(fvv);
  gsl_vector_set(fvv, p, 2.0 * normv * normv);

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;
}

void
solve_system(const gsl_vector *x0, gsl_multilarge_nlinear_fdf *fdf,
             gsl_multilarge_nlinear_parameters *params)
{
  const gsl_multilarge_nlinear_type *T = gsl_multilarge_nlinear_trust;
  const size_t max_iter = 200;
  const double xtol = 1.0e-8;
  const double gtol = 1.0e-8;
  const double ftol = 1.0e-8;
  const size_t n = fdf->n;
  const size_t p = fdf->p;
  gsl_multilarge_nlinear_workspace *work =
    gsl_multilarge_nlinear_alloc(T, params, n, p);
  gsl_vector * f = gsl_multilarge_nlinear_residual(work);
  gsl_vector * x = gsl_multilarge_nlinear_position(work);
  int info;
  double chisq0, chisq, rcond, xsq;
  struct timeval tv0, tv1;

  gettimeofday(&tv0, NULL);

  /* initialize solver */
  gsl_multilarge_nlinear_init(x0, fdf, work);

  /* store initial cost */
  gsl_blas_ddot(f, f, &chisq0);

  /* iterate until convergence */
  gsl_multilarge_nlinear_driver(max_iter, xtol, gtol, ftol,
                                NULL, NULL, &info, work);

  gettimeofday(&tv1, NULL);

  /* store final cost */
  gsl_blas_ddot(f, f, &chisq);

  /* compute final ||x||^2 */
  gsl_blas_ddot(x, x, &xsq);

  /* store cond(J(x)) */
  gsl_multilarge_nlinear_rcond(&rcond, work);

  /* print summary */
  fprintf(stderr, "%-25s %-5zu %-4zu %-5zu %-6zu %-4zu %-10.4e %-10.4e %-7.2f %-11.4e %.2f\n",
          gsl_multilarge_nlinear_trs_name(work),
          gsl_multilarge_nlinear_niter(work),
          fdf->nevalf,
          fdf->nevaldfu,
          fdf->nevaldf2,
          fdf->nevalfvv,
          chisq0,
          chisq,
          1.0 / rcond,
          xsq,
          (tv1.tv_sec - tv0.tv_sec) + 1.0e-6 * (tv1.tv_usec - tv0.tv_usec));

  gsl_multilarge_nlinear_free(work);
}

int
main (void)
{
  const size_t p = 2000;
  const size_t n = p + 1;
  gsl_vector *f = gsl_vector_alloc(n);
  gsl_vector *x = gsl_vector_alloc(p);

  /* allocate sparse Jacobian matrix with 2*p non-zero elements in triplet format */
  gsl_spmatrix *J = gsl_spmatrix_alloc_nzmax(n, p, 2 * p, GSL_SPMATRIX_TRIPLET);

  gsl_multilarge_nlinear_fdf fdf;
  gsl_multilarge_nlinear_parameters fdf_params =
    gsl_multilarge_nlinear_default_parameters();
  struct model_params params;
  size_t i;

  params.alpha = 1.0e-5;
  params.J = J;

  /* define function to be minimized */
  fdf.f = penalty_f;
  fdf.df = penalty_df;
  fdf.fvv = penalty_fvv;
  fdf.n = n;
  fdf.p = p;
  fdf.params = &params;

  for (i = 0; i < p; ++i)
    {
      /* starting point */
      gsl_vector_set(x, i, i + 1.0);

      /* store sqrt(alpha)*I_p in upper p-by-p block of J */
      gsl_spmatrix_set(J, i, i, sqrt(params.alpha));
    }

  fprintf(stderr, "%-25s %-4s %-4s %-5s %-6s %-4s %-10s %-10s %-7s %-11s %-10s\n",
          "Method", "NITER", "NFEV", "NJUEV", "NJTJEV", "NAEV", "Init Cost",
          "Final cost", "cond(J)", "Final |x|^2", "Time (s)");

  fdf_params.scale = gsl_multilarge_nlinear_scale_levenberg;

  fdf_params.trs = gsl_multilarge_nlinear_trs_lm;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multilarge_nlinear_trs_lmaccel;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multilarge_nlinear_trs_dogleg;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multilarge_nlinear_trs_ddogleg;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multilarge_nlinear_trs_subspace2D;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multilarge_nlinear_trs_cgst;
  solve_system(x, &fdf, &fdf_params);

  gsl_vector_free(f);
  gsl_vector_free(x);
  gsl_spmatrix_free(J);

  return 0;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/// @file Error.h
/// @author Braxton Salyer <braxtonsalyer@gmail.com>
/// @brief Basic high-level types for communicating recoverable errors
/// @version 0.1
/// @date 2021-11-03
///
/// MIT License
/// @copyright Copyright (c) 2021 Braxton Salyer <braxtonsalyer@gmail.com>
///
/// Permission is hereby granted, free of charge, to any person obtaining a copy
/// of this software and associated documentation files (the "Software"), to
/// deal in the Software without restriction, including without limitation the
/// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
/// sell copies of the Software, and to permit persons to whom the Software is
/// furnished to do so, subject to the following conditions:
///
/// The above copyright notice and this permission notice shall be included in
/// all copies or substantial portions of the Software.
///
/// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
/// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
/// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
/// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
/// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
/// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
/// IN THE SOFTWARE.
#pragma once

#include <Hyperion/BasicTypes.h>
#include <Hyperion/HyperionDef.h>
#include <Hyperion/Ignore.h>
#include <Hyperion/Memory.h>
#include <Hyperion/error/SystemDomain.h>
#include <cstring>
#include <gsl/gsl>
#include <memory>
#include <string>
#include <system_error>
#include <type_traits>

///	@defgroup error Error
/// The Error module provides various error-handling facilities and types for
/// communicating and handling recoverable errors and for aborting on
/// irrecoverable errors in a communicable way
/// @headerfile "Hyperion/Error.h"

namespace hyperion::error {
	using namespace std::literals::string_literals;

	/// @brief Basic interface for recoverable errors. Error types may provide
	/// additional functionality beyond this, but this much is required.
	/// @ingroup error
	/// @headerfile "Hyperion/Error.h"
	IGNORE_WEAK_VTABLES_START

	class [[nodiscard]] ErrorBase {
	  public:
		constexpr ErrorBase() noexcept = default;

		constexpr ErrorBase(const ErrorBase&) noexcept = default;

		constexpr ErrorBase(ErrorBase&&) noexcept = default;

		constexpr virtual ~ErrorBase() noexcept = default;

		/// @brief Returns the value of the associated error code as an `i64`
		/// @return The value of the associated error code
		/// @ingroup error
		[[nodiscard]] virtual constexpr auto value() const noexcept -> i64 = 0;

		/// @brief Returns the message associated with the error as a `std::string`
		/// @return The message associated with the error
		/// @ingroup error
		[[nodiscard]] virtual auto message() const noexcept -> std::string = 0;

		/// @brief Converts this error into a `std::string`
		/// @return This error as a `std::string`
		/// @ingroup error
		[[nodiscard]] virtual auto to_string() const noexcept -> std::string = 0;

		constexpr auto operator=(const ErrorBase&) noexcept -> ErrorBase& = default;

		constexpr auto operator=(ErrorBase&&) noexcept -> ErrorBase& = default;
	};

	IGNORE_WEAK_VTABLES_STOP

	// class [[nodiscard]] AnyError;

	/// @brief General purpose type for communicating recoverable errors.
	///
	/// Wraps an `error::ErrorCode<Domain>` in a type-safe manner.
	/// Usually makes up the `E` variant of a `Result<T, E>`
	///
	/// @tparam Domain - The `StatusCodeDomain` associated with the this error
	/// @ingroup error
	template<StatusCodeDomain Domain = SystemDomain>
	class [[nodiscard]] Error final : public ErrorBase {
	  public:
		/// @brief The type used to represent the error (an enum, integer error code,
		/// etc)
		/// @ingroup error
		using value_type = typename ErrorCode<Domain>::value_type;

		friend class AnyError;

		/// @brief Constructs a default `Error` with an error code representing an
		/// unknown error
		/// @ingroup error
		constexpr Error() noexcept {

			if constexpr(StatusCodeEnum<value_type>) {
				m_error_code = make_error_code(static_cast<value_type>(-1));
			}
			else {
				m_error_code = make_error_code<Domain>(static_cast<value_type>(-1));
			}
		}

		/// @brief Constructs an `Error` from the given `ErrorCode<Domain>`
		///
		/// @param code - The error code
		/// @ingroup error
		constexpr Error(const ErrorCode<Domain>& code) noexcept // NOLINT
			: m_error_code(code) {
		}

		/// @brief Constructs an `Error` from the given `ErrorCode<Domain>`
		///
		/// @param code - The error code
		/// @ingroup error
		constexpr Error(ErrorCode<Domain>&& code) noexcept // NOLINT
			: m_error_code(std::move(code)) {
		}

		/// @brief Constructs an `Error` from the given error code
		///
		/// @param code - The error code
		/// @ingroup error
		constexpr Error(const value_type& code) noexcept // NOLINT
			requires(!StatusCodeEnum<value_type>)
			: m_error_code(make_error_code<Domain>(code)) {
		}

		/// @brief Constructs an `Error` from the given error code
		///
		/// @param code - The error code
		/// @ingroup error
		constexpr Error(const value_type& code) noexcept // NOLINT
			requires StatusCodeEnum<value_type> : m_error_code(make_error_code(code)) {
		}

		/// @brief Copy-constructs an `Error`
		/// @ingroup error
		constexpr Error(const Error& error) = default;

		/// @brief Move-constructs an `Error`
		/// @ingroup error
		constexpr Error(Error&& error) noexcept = default;

		/// @brief Destroys this `Error`
		/// @ingroup error
		constexpr ~Error() noexcept final = default;

		/// @brief Returns the `i64` value of the error code associated with this
		/// `Error`
		///
		/// @return The associated error code as an `i64`
		/// @ingroup error
		[[nodiscard]] auto value() const noexcept -> i64 final {

			return m_error_code.value();
		}

		/// @brief Returns the `ErrorCode` associated with this `Error`
		///
		/// @return The associated `ErrorCode`
		/// @ingroup error
		[[nodiscard]] auto code() const noexcept -> const ErrorCode<Domain>& {

			return m_error_code;
		}

		/// @brief Returns the error message associated with this `Error`
		///
		/// @return The error message
		/// @ingroup error
		[[nodiscard]] inline auto message() const noexcept -> std::string final {

			if constexpr(concepts::Same<error_code_message_type, std::string>) {
				return m_error_code.message();
			}
			else {
				return std::string(m_error_code.message());
			}
		}

		/// @brief Returns the error message associated with this `Error`, as a cstring
		///
		/// @return The error message
		/// @ingroup error
		[[nodiscard]] auto message_as_cstr() const noexcept -> const char* {

			return m_error_code.message().c_str();
		}

		/// @brief Gets the `std::string` representation of this `Error`
		///
		/// @return this `Error` formatted as a `std::string`
		/// @ingroup error
		[[nodiscard]] auto to_string() const noexcept -> std::string final {

			return fmt::format("Error: {}", m_error_code.message());
		}

		/// @brief Converts this `Error` into its underlying `ErrorCode<Domain>`
		/// @ingroup error
		constexpr operator ErrorCode<Domain>() const noexcept { // NOLINT
			return m_error_code;
		}

		/// @brief Converts this `Error` into its underlying `value_type`
		/// @ingroup error
		constexpr operator value_type() const noexcept { // NOLINT
			return static_cast<value_type>(m_error_code);
		}

		/// @brief Equality comparison operator with another `Error` from a possibly different
		/// `StatusCodeDomain`
		/// @tparam Domain2 - The `StatusCodeDomain` of `error`
		/// @ingroup error
		template<StatusCodeDomain Domain2>
		constexpr auto operator==(const Error<Domain2>& error) const noexcept -> bool {
			return m_error_code == error.m_error_code;
		}

		/// @brief Inequality comparison operator with another `Error` from a possibly different
		/// `StatusCodeDomain`
		/// @tparam Domain2 - The `StatusCodeDomain` of `error`
		/// @ingroup error
		template<StatusCodeDomain Domain2>
		constexpr auto operator!=(const Error<Domain2>& error) const noexcept -> bool {
			return m_error_code != error.m_error_code;
		}

		/// @brief Equality comparison operator with an `ErrorCode` from a possibly different
		/// `StatusCodeDomain`
		/// @tparam Domain2 - The `StatusCodeDomain` of `code`
		/// @ingroup error
		template<StatusCodeDomain Domain2>
		constexpr auto operator==(const ErrorCode<Domain2>& code) const noexcept -> bool {
			return m_error_code == code;
		}

		/// @brief Inequality comparison operator with an `ErrorCode` from a possibly different
		/// `StatusCodeDomain`
		/// @tparam Domain2 - The `StatusCodeDomain` of `code`
		/// @ingroup error
		template<StatusCodeDomain Domain2>
		constexpr auto operator!=(const ErrorCode<Domain2>& code) const noexcept -> bool {
			return m_error_code != code;
		}

		/// @brief Copy-assigns the given `error` to this `Error`
		/// @ingroup error
		constexpr auto operator=(const Error& error) noexcept -> Error& = default;

		/// @brief Move-assigns the given `error` to this `Error`
		/// @ingroup error
		constexpr auto operator=(Error&& error) noexcept -> Error& = default;

		/// @brief Copy-assigns the given `ErrorCode<Domain>` to this `Error`
		/// @ingroup error
		constexpr auto operator=(const ErrorCode<Domain>& code) noexcept -> Error& {

			m_error_code = code;

			return *this;
		}

		/// @brief Move-assigns the given `ErrorCode<Domain>` to this `Error`
		/// @ingroup error
		constexpr auto operator=(ErrorCode<Domain>&& code) noexcept -> Error& {

			m_error_code = std::move(code);

			return *this;
		}

		/// @brief Provides `std::tuple`-like structured binding support.
		///
		/// `Index == 0` returns the error code value.
		/// `Index == 1` returns the error message.
		/// Other `Index` values ar invalid
		///
		/// @return The error code value (`Index == 0`), or the error message (`Index
		/// == 1`)
		/// @ingroup error
		template<usize Index>
		auto get() const&& noexcept {

			return get_impl<Index>(*this);
		}

		/// @brief Provides `std::tuple`-like structured binding support.
		///
		/// `Index == 0` returns the error code value.
		/// `Index == 1` returns the error message.
		/// Other `Index` values ar invalid
		///
		/// @return The error code value (`Index == 0`), or the error message (`Index
		/// == 1`)
		/// @ingroup error
		template<usize Index>
		auto get() const& noexcept {

			return get_impl<Index>(*this);
		}

	  private:
		template<usize Index, typename T>
		auto get_impl(T&& t) const noexcept {

			static_assert(Index < 2, "Index out of bounds for hyperion::error::Error::get");
			if constexpr(Index == 0) {
				return std::forward<T>(t).m_error_code.value();
			}
			else {
				return std::forward<T>(t).m_error_code.message();
			}
		}

		/// error code
		ErrorCode<Domain> m_error_code;

		using error_code_message_type = decltype(m_error_code.message());
	};
} // namespace hyperion::error

/// @brief Specialize `std::tuple_size` for `error::Error<Domain>`
/// @ingroup error
template<hyperion::error::StatusCodeDomain Domain>
struct std::tuple_size<hyperion::error::Error<Domain>> {
	static constexpr hyperion::usize value = 2;
};

/// @brief Specialize `std::tuple_element<0, T>` for `error::Error<Domain>`
/// @ingroup error
template<hyperion::error::StatusCodeDomain Domain>
struct std::tuple_element<0, hyperion::error::Error<Domain>> {
	using type = typename hyperion::error::Error<Domain>::value_type;
};

/// @brief Specialize `std::tuple_element<1, T>` for `error::Error<Domain>`
/// @ingroup error
template<hyperion::error::StatusCodeDomain Domain>
struct std::tuple_element<1, hyperion::error::Error<Domain>> {
	using type = std::string;
};

namespace hyperion::error {

	/// @brief `SystemError` is an error representing the default platform/OS level
	/// errors (eg POSIX or WIN32 error codes)
	///
	/// In practice, this will usually be an alias for one of `PosixError`,
	/// `Win32Error`, or `NTError`, depending on platform and configuration
	/// @ingroup error
	using SystemError = Error<SystemDomain>;
	/// @brief `PosixError` is an error representing a platform's implementation of
	/// POSIX error codes (it includes support for additional platform specific
	/// codes not required by POSIX)
	/// @ingroup error
	using PosixError = Error<PosixDomain>;
	/// @brief `GenericError` represents the strict set of POSIX required error
	/// codes.
	///
	/// In Hyperion's error system, `GenericError` fills a role similar to a
	/// `std::error_code` of `std::generic_category`
	/// @ingroup error
	using GenericError = Error<GenericDomain>;

#if HYPERION_PLATFORM_WINDOWS
	/// @brief `Win32Error` is an error representing the Win32 error codes
	/// @ingroup error
	using Win32Error = Error<Win32Domain>;
	/// @brief `NTError` is an error representing the Windows NT error codes
	/// @ingroup error
	using NTError = Error<NTDomain>;
#endif // HYPERION_PLATFORM_WINDOWS

	/// @brief Concept specifying the requirements of a type necessary to be
	/// guaranteed compatible with Hyperion's error system.
	///
	/// Users are encouraged to fulfill this concept for their own error types,
	/// particularly when using them with Hyperion's higher-level error-handling
	/// facilities like `Result<T, E>`
	/// @ingroup error
	template<typename E>
	concept ErrorType = std::derived_from<E, ErrorBase>;

	/// @brief Concept useful for shorthanding `!ErrorType<E>`
	/// @note Removes the need for extra parentheses and enables cleaner
	/// auto-formatting (eg `SomeConcept &&(!ErrorType<E>)` vs `SomeConcept &&
	/// NotErrorType<E>`)
	/// @ingroup error
	template<typename E>
	concept NotErrorType = !ErrorType<E>;

	/// @brief `AnyError` represents a type erased error from any
	/// `error::StatusCodeDomain`
	/// @ingroup error
	class [[nodiscard]] AnyError final : public ErrorBase {
	  public:
		/// @brief Constructs an `AnyError` as an unknown error
		/// @ingroup error
		HYPERION_CONSTEXPR_STRINGS AnyError() noexcept = default;

		/// @brief Constructs an `AnyError` from the given `error::ErrorCode<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error
		///
		/// @param code - The error code this `AnyError` should represent
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS AnyError(const ErrorCode<Domain>& code) noexcept // NOLINT
			: m_error_code(code.value()), m_message(code.message()) {
		}

		/// @brief Constructs an `AnyError` from the given `error::ErrorCode<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error
		///
		/// @param code - The error code this `AnyError` should represent
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS AnyError(ErrorCode<Domain>&& code) noexcept // NOLINT
			: m_error_code(code.value()), m_message(code.message()) {
		}

		/// @brief Constructs an `AnyError` from the given `error::Error<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error
		///
		/// @param error - The error this `AnyError` should represent
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS AnyError(const Error<Domain>& error) noexcept // NOLINT
			: m_error_code(error.code().value()), m_message(error.message()) {
		}

		/// @brief Constructs an `AnyError` from the given `error::Error<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error
		///
		/// @param error - The error this `AnyError` should represent
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS AnyError(Error<Domain>&& error) noexcept // NOLINT
			: m_error_code(error.code().value()), m_message(error.message()) {
		}
		/// @brief Copy-constructs an `AnyError` from the given one
		/// @ingroup error
		HYPERION_CONSTEXPR_STRINGS AnyError(const AnyError&) noexcept = default;
		/// @brief Move-constructs an `AnyError` from the given one
		/// @ingroup error
		HYPERION_CONSTEXPR_STRINGS AnyError(AnyError&&) noexcept = default;
		/// @brief Destroys this `AnyError`
		/// @ingroup error
		HYPERION_CONSTEXPR_STRINGS ~AnyError() noexcept final = default;

		/// @brief Returns the `i64` value corresponding with the error code this
		/// represents
		/// @return The error code as an `i64`
		/// @ingroup error
		[[nodiscard]] constexpr auto value() const noexcept -> i64 final {

			return m_error_code;
		}

		/// @brief Returns the error message associated with the error code this
		/// represents
		/// @return The error code message
		/// @ingroup error
		[[nodiscard]] auto message() const noexcept -> std::string final {

			return m_message;
		}

		/// @brief Returns the string representation of this `AnyError`
		/// @return this `AnyError` as a `std::string`
		/// @ingroup error
		[[nodiscard]] HYPERION_CONSTEXPR_STRINGS auto
		to_string() const noexcept -> std::string final {

			return "Error: "s + m_message + "\n"s;
		}

		/// @brief Copy-assigns the given `AnyError` to this one
		/// @ingroup error
		HYPERION_CONSTEXPR_STRINGS auto operator=(const AnyError&) noexcept -> AnyError& = default;
		/// @brief Move-assigns the given `AnyError` to this one
		/// @ingroup error
		HYPERION_CONSTEXPR_STRINGS auto operator=(AnyError&&) noexcept -> AnyError& = default;

		/// @brief Copy-assigns this `AnyError` with the given
		/// `error::ErrorCode<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error code
		///
		/// @param code - The error code to assign to this
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS auto
		operator=(const ErrorCode<Domain>& code) noexcept -> AnyError& {

			m_error_code = code.value();
			m_message = code.message();

			return *this;
		}

		/// @brief Move-assigns this `AnyError` with the given
		/// `error::ErrorCode<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error code
		///
		/// @param code - The error code to assign to this
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS auto operator=(ErrorCode<Domain>&& code) noexcept -> AnyError& {

			m_error_code = code.value();
			m_message = code.message();

			return *this;
		}

		/// @brief Copy-assigns this `AnyError` with the given `error::Error<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error
		///
		/// @param error - The error to assign to this
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS auto
		operator=(const Error<Domain>& error) noexcept -> AnyError& {

			m_error_code = error.value();
			m_message = error.message();

			return *this;
		}

		/// @brief Move-assigns this `AnyError` with the given `error::Error<Domain>`
		///
		/// @tparam Domain - The `error::StatusCodeDomain` of the error
		///
		/// @param error - The error to assign to this
		/// @ingroup error
		template<StatusCodeDomain Domain>
		HYPERION_CONSTEXPR_STRINGS auto operator=(Error<Domain>&& error) noexcept -> AnyError& {

			m_error_code = error.value();
			m_message = error.message();

			return *this;
		}

		/// @brief Provides `std::tuple`-like structured binding support.
		///
		/// `Index == 0` returns the error code value.
		/// `Index == 1` returns the error message.
		/// Other `Index` values ar invalid
		///
		/// @return The error code value (`Index == 0`), or the error message (`Index
		/// == 1`)
		/// @ingroup error
		template<usize Index>
		auto get() const&& noexcept {

			return get_impl<Index>(*this);
		}

		/// @brief Provides `std::tuple`-like structured binding support.
		///
		/// `Index == 0` returns the error code value.
		/// `Index == 1` returns the error message.
		/// Other `Index` values ar invalid
		///
		/// @return The error code value (`Index == 0`), or the error message (`Index
		/// == 1`)
		/// @ingroup error
		template<usize Index>
		auto get() const& noexcept {

			return get_impl<Index>(*this);
		}

	  private:
		template<usize Index, typename T>
		auto get_impl(T&& t) const noexcept {

			static_assert(Index < 2, "Index out of bounds for hyperion::error::AnyError::get");
			if constexpr(Index == 0) {
				return std::forward<T>(t).m_error_code;
			}
			else {
				return std::forward<T>(t).m_message;
			}
		}

		i64 m_error_code = -1;

		std::string m_message;
	};
} // namespace hyperion::error

/// @brief Specialize `std::tuple_size` for `error::AnyError`
/// @ingroup error
template<>
struct std::tuple_size<hyperion::error::AnyError> {
	static constexpr hyperion::usize value = 2;
};

/// @brief Specialize `std::tuple_element<0, T>` for `error::AnyError`
/// @ingroup error
template<>
struct std::tuple_element<0, hyperion::error::AnyError> {
	using type = hyperion::i64;
};

/// @brief Specialize `std::tuple_element<1, T>` for `error::AnyError`
/// @ingroup error
template<>
struct std::tuple_element<1, hyperion::error::AnyError> {
	using type = std::string;
};

ext: h


text:
#pragma once

#include <cstdint>
#include <gsl/gsl>
#include "utils.h"
#include "halley/text/halleystring.h"

struct XXH64_state_s;
typedef struct XXH64_state_s XXH64_state_t;

namespace Halley {
    class Hash {
    public:
        static uint64_t hash(const Bytes& bytes);
        static uint64_t hash(gsl::span<const gsl::byte> bytes);

    	template <typename T>
    	static uint64_t hashValue(const T& v)
    	{
    		return hash(gsl::as_bytes(gsl::span<const T>(&v, 1)));
    	}

		static uint32_t compressTo32(uint64_t value);

		class Hasher
		{
		public:
			Hasher();
			~Hasher();

			void feed(const String& string)
			{
				feedBytes(gsl::as_bytes(gsl::span<const char>(string.c_str(), string.length())));
			}

			template<typename T>
			void feed(const T& data)
			{
				feedBytes(gsl::as_bytes(gsl::span<const T>(&data, 1)));
			}

			void feedBytes(gsl::span<const gsl::byte> bytes);

			[[nodiscard]] uint64_t digest();

		private:
			bool ready;
			XXH64_state_t* data;
		};
    };
}

max_stars_repo_name: code-disaster/halley

max_issues_repo_name: code-disaster/halley

max_forks_repo_name: code-disaster/halley


text:
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

const gsl_interp_type   *types[7];
const gsl_interp2d_type *types2d[2];

void init(void)
{
  types[0]   = gsl_interp_linear;
  types[1]   = gsl_interp_polynomial;
  types[2]   = gsl_interp_cspline;
  types[3]   = gsl_interp_cspline_periodic;
  types[4]   = gsl_interp_akima;
  types[5]   = gsl_interp_akima_periodic;
  types[6]   = gsl_interp_steffen;
  types2d[0] = gsl_interp2d_bilinear;
  types2d[1] = gsl_interp2d_bicubic;
}

/* 1D */
gsl_interp *mgsl_interp_alloc(int type, size_t size)
{
  if(types[0] == NULL) init();
  return gsl_interp_alloc(types[type], size);
}

unsigned int mgsl_interp_type_min_size(int type)
{
  if(types[0] == NULL) init();
  return gsl_interp_type_min_size(types[type]);
}

gsl_spline *mgsl_spline_alloc(int type, size_t size)
{
  if(types[0] == NULL) init();
  return gsl_spline_alloc(types[type], size);
}

/* 2D */
gsl_interp2d *mgsl_interp2d_alloc(int type, size_t xsize, size_t ysize)
{
  if(types[0] == NULL) init();
  return gsl_interp2d_alloc(types2d[type], xsize, ysize);
}

unsigned int mgsl_interp2d_type_min_size(int type)
{
  if(types[0] == NULL) init();
  return gsl_interp2d_type_min_size(types2d[type]);
}

gsl_spline2d *mgsl_spline2d_alloc(int type, size_t xsize, size_t ysize)
{
  if(types[0] == NULL) init();
  return gsl_spline2d_alloc(types2d[type], xsize, ysize);
}


text:
#ifdef __LOCAL_SIMULATION_ENV__

#include "model.h"

#endif

//############# START DYNAMIC_CLAMP CODE ##################

/*!
Dynamic clamp model for a current generated by multiple clusters of cooperative ion channels:
I_{coop} = -g*N_{open}*(V-E)
where the number of open channels N_{open} changes due to the stochastic opening and closing events guided by rates
which depend on the external voltage and the number of open channels among the surrounding channels in a cluster i.e.
alpha(V, o) and beta(V,o)

\par Input/Output
- V: Measured membrane potential in mV
- \f$ I_{inj} \f$ : Injected current in nA

\par Parameter
- g : conductance of a single channel in nS
- E : reversal potential of the ionic current in  mV
- cN: number of clusters of cooperative ion channels without units
- cS: number of cooperative ion channels in a cluster without units
- J : pairwise interaction strength in mV
- tau : peak reaction time scale of a single channel in ms
*/

#ifndef __KERNEL__

#include <gsl/gsl_rng.h>
#include <time.h>

#endif

#define POWER_2_EXP 28
unsigned long next = 1;
const int RAND_PERIOD = 1 << POWER_2_EXP;

// numerical recipes in c, p. 276, not recommended for usage as the period is 30k and other issues
float evil_rand(void) {
    next = next * 1103515245 + 1234;
    return (float) ((unsigned int) (next / 65536) % 32768) / 32768;
}

#if defined (__KERNEL__) || defined (DYNCLAMPMODEL)

/*! Name, by which this module is known inside Linux: */
const char *modelName;

/*! Variables used by the model */
#define MAX_CLUSTER_SIZE 8
#define DEFAULT_CLUSTER_NUMBER 100
#define DEFAULT_CLUSTER_INIT 0
int cluster_states[MAX_CLUSTER_SIZE + 1] = {DEFAULT_CLUSTER_NUMBER, 0, 0, 0, 0, 0, 0, 0, 0};
int cluster_size = MAX_CLUSTER_SIZE;
int cluster_number = DEFAULT_CLUSTER_NUMBER;
int cluster_init = DEFAULT_CLUSTER_INIT;

/*! The period length of the realtime periodic task in seconds. */
float loopInterval;
/*! One over the period length of the realtime periodic task in Hertz. */
float loopRate;

/*! Analog input that is read from the DAQ board. */

#define INPUT_N 1
/*! The \a inputNames are used to match the \a input variables with
analog input traces in Relacs. */
const char *inputNames[INPUT_N] = {"V-1"};
const char *inputUnits[INPUT_N] = {"mV"};
/*! The \a inputChannels are set automatically. */
int inputChannels[INPUT_N];
/*! \a input holds the current value that was read in from the DAQ board. */
float input[INPUT_N] = {0.0};

/*! Analog output that is written to the DAQ board. */

#define OUTPUT_N 1
const char *outputNames[OUTPUT_N] = {"Current-1"};
const char *outputUnits[OUTPUT_N] = {"nA"};
int outputChannels[OUTPUT_N];
float output[OUTPUT_N] = {0.0};

/*! Parameter that are provided by the model and can be read out. */

#define PARAMINPUT_N 2
const char *paramInputNames[PARAMINPUT_N] = {"Cooperative-Cluster-Current", "Open-Channels"};
const char *paramInputUnits[PARAMINPUT_N] = {"nA", ""};
float paramInput[PARAMINPUT_N] = {0.0, 0.0};

/*! Parameter that are read by the model and are written to the model. */

#define PARAMOUTPUT_N 8
const char *paramOutputNames[PARAMOUTPUT_N] = {"g", "E", "J", "tau", "Cluster size", "Number of clusters", "Vhalf", "Initial state"};
const char *paramOutputUnits[PARAMOUTPUT_N] = {"pS", "mV", "mV", "ms", "", "", "mV", ""};
float paramOutput[PARAMOUTPUT_N] = {0.0, 0.0, 0.0, 15, MAX_CLUSTER_SIZE, DEFAULT_CLUSTER_NUMBER, -1.0, 0};

/*! Functions used by the model */
double opening_rate(double voltage) {
    double oR;
    float vhalf = paramOutput[6];
    if (paramOutput[3] != 0.0) {
        oR = cosh((voltage - vhalf) / 30.0) * (1.0 + tanh((voltage - vhalf) / 15.0)) / paramOutput[3] / 2.0;
    } else {
        oR = 0;
    }
    return oR;
}

double closing_rate(double voltage) {
    double cR;
    float vhalf = paramOutput[6];
    if (paramOutput[3] != 0.0) {
        cR = cosh((voltage - vhalf) / 30.0) * (1.0 - tanh((voltage - vhalf) / 15.0)) / paramOutput[3] / 2.0;
    } else {
        cR = 0;
    }
    return cR;
}

double get_rand(void) {
    double u;

#ifdef ENABLE_LOOKUPTABLES
    // random number from lookuptable:
    next = (next >= (RAND_PERIOD - 1)) ? 0 : (next + 1);
    u = (double) lookupx[0][next];
#else
    u = (double) evil_rand();
#endif
    return u;
}

int get_number_of_open_channels(void) {
    int sum = 0;
    int size = cluster_size;
    int i;
    for (i = 0; i <= size; i++) {
        sum += cluster_states[i] * i;
    }
    return sum;
}

/*
brute_force on cluster state basis: rates based, does not need any internal state
- needs 2*cluster_size random numbers per time step
- problem: as multiple reactions might take place in a single time step an opening and a closing reaction might take
place simultaneously leaving behind an inconsistent state e.g. [0 1 0] -> [1 -1 1]
*/
void evolve_cluster_states(void) {
    double rand;
    double dt = loopInterval * 1000;
    double voltage = input[0];
    int cS = cluster_size;
    float cluster_interaction = paramOutput[2];

    //helpers
    int i;
    double opening_rates[cS];
    double closing_rates[cS];

    // calculate opening and closing rates
    for (i = 0; i < cS; i++) {
        // opening_rate[i]: rate from (i channels open) to (i+1 channels open)
        opening_rates[i] = cluster_states[i] * (cS - i) * opening_rate(voltage + i * cluster_interaction);
        // closing_rate[i]: rate from (i+1 channels open) to (i channels open)
        closing_rates[i] = cluster_states[i + 1] * (i + 1) * closing_rate(voltage + i * cluster_interaction);
    }

    // update cluster states
    for (i = 0; i < cS; i++) {
        //opening reaction from (i channels open) to (i+1 channels open)
        rand = get_rand();
        if (rand <= opening_rates[i] * dt) {
            //DEBUG_MSG("C -> Likelihood %.16f, Threshold:%.16f\n", opening_rates[i]*dt, rand);
            //DEBUG_MSG("Opening %d:%d\n", i, cluster_states[i]);
            if (cluster_states[i] != 0) {
                cluster_states[i] -= 1;
                cluster_states[i + 1] += 1;
            }
        }

        // closing reaction from (i+1 channels open) to (i channels open)
        rand = get_rand();
        if (rand <= closing_rates[i] * dt) {
            //DEBUG_MSG("C -> Likelihood %.16f, Threshold:%.16f\n", closing_rates[i]*dt, rand);
            //DEBUG_MSG("Closing %d:%d\n", i+1,cluster_states[i+1]);
            if (cluster_states[i + 1] != 0) {
                cluster_states[i + 1] -= 1;
                cluster_states[i] += 1;
            }
        }
    }
}

void populate_state(int state_idx) {
    int idx;
    cluster_states[state_idx] = cluster_number;
    for (idx = 0; idx <= cluster_size; idx++) {
            if(idx != state_idx) {
                cluster_states[idx] = 0;
            }
        }
}

void process_cluster_property_changes(void) {
    int cluster_size_idx = 4;
    int cluster_number_idx = 5;
    int cluster_init_idx = 7;

    if (cluster_size != (int) paramOutput[cluster_size_idx] ||
        cluster_number != (int) paramOutput[cluster_number_idx] || cluster_init != (int) paramOutput[cluster_init_idx]) {
        if (paramOutput[cluster_size_idx] > MAX_CLUSTER_SIZE) {
            cluster_size = MAX_CLUSTER_SIZE;
        } else {
            cluster_size = (paramOutput[cluster_size_idx] >= 1) ? (int) paramOutput[cluster_size_idx] : 1;
        };

        if (paramOutput[cluster_init_idx] > MAX_CLUSTER_SIZE) {
            cluster_init = MAX_CLUSTER_SIZE;
        } else {
            cluster_init = (paramOutput[cluster_init_idx] >= 0) ? (int) paramOutput[cluster_init_idx] : 0;
        };

        cluster_number = (paramOutput[cluster_number_idx] > 0) ? (int) paramOutput[cluster_number_idx] : 0;

        populate_state(cluster_init);
    }
}

/*! Functions provided by the model for the dynamic clamp module */

void initModel(void) {
    modelName = "cooperative_clusters";
    process_cluster_property_changes();
}

#define CONDUCTANCE_IDX 0
#define MEMBRANE_POTENTIAL_IDX 0
#define REVERSAL_POTENTIAL_IDX 1

void computeModel(void) {

    int number_of_open_channels = 0;
    process_cluster_property_changes();
    evolve_cluster_states();
    // number of open channels
    number_of_open_channels = get_number_of_open_channels();
    paramInput[0] = (float) (-0.000001 * paramOutput[CONDUCTANCE_IDX] * number_of_open_channels *
            (input[MEMBRANE_POTENTIAL_IDX] - paramOutput[REVERSAL_POTENTIAL_IDX]));
    // current flowing through cooperative clusters :
    paramInput[1] = number_of_open_channels;
    // total injected current:
    output[0] = paramInput[0];
}

#endif

#ifndef __KERNEL__
#ifdef ENABLE_LOOKUPTABLES

/*! This function is called from DynClampAnalogOutput in user
    space/C++ context and can be used to create some lookuptables for
    nonlinear functions to be used by computeModel(). The implementation of this
    functions has to allocate an \a x and \a y array of floats of a sensible size \a n.
    \param[in] \a k : the index for the lookup table to be generated.
    \param[out] \a n : the size of the lookup table (the number of elements in \a x and \a y).
    \param[out] \a x : the x-values.
    \param[out] \a y : the corresponding y-values.
    \return: 0 if a lookuptable was generated, -1 otherwise.
*/

gsl_rng *rng_state;

unsigned long get_seed_from_clock();

void initialize_random_number_generator() {
    const gsl_rng_type *T;
    unsigned long int seed = get_seed_from_clock();

    gsl_rng_env_setup();
    T = gsl_rng_default;

    rng_state = gsl_rng_alloc(T);
    gsl_rng_set(rng_state, seed);
}

unsigned long get_seed_from_clock() {
    clock_t current;
    current = clock();
    return (unsigned long int) current;
}

int generateLookupTable(int k, float **x, float **y, int *n) {
    initialize_random_number_generator();
    if (k == 0) {
        /* Lookup-table for the Random number generator: */
        const int nn = RAND_PERIOD;
        *n = nn;
#ifdef __LOCAL_SIMULATION_ENV__
        *x = (float *) malloc(nn * sizeof(float));
#else
        *x = new float[nn];
#endif
        for (int j = 0; j < nn; j++) {
            float xx = gsl_rng_uniform(rng_state);
            (*x)[j] = xx;
        }
        return 0;
    }

    return -1;
}

#endif
#endif

//############# END DYNAMIC_CLAMP CODE ##################

#ifdef __LOCAL_SIMULATION_ENV__

void setVoltage(float voltage) {
    input[0] = voltage;
}

void setConductance(float conductance) {
    paramOutput[0] = conductance;
}

void setReversalPotential(float reversalPotential) {
    paramOutput[1] = reversalPotential;
}

void setClusterInteraction(float clusterInteraction) {
    paramOutput[2] = clusterInteraction;
}

void setTimeConstant(float timeConstant) {
    paramOutput[3] = timeConstant;
}

void setClusterSize(float clusterSize) {
    paramOutput[4] = clusterSize;
}

void setClusterNumber(float clusterNumber) {
    paramOutput[5] = clusterNumber;
}

void setVhalf(float vhalf) {
    paramOutput[6] = vhalf;
}

void setClusterInit(float clusterInit) {
    paramOutput[7] = clusterInit;
}

void setLoopInterval(float interval) {
    loopInterval = interval;
}

float getCurrentOutput(void) {
    return output[0];
}

int getClusterSize(void) {
    return cluster_size;
}

int getClusterState(int i) {
    return cluster_states[i];
}

#endif


text:
#pragma once
#include "Buffer.h"
#include "Enum.h"
#include "Error.h"
#include "detail/Defaulted.h"
#include <boost/hana/define_struct.hpp>
#include <gsl/span>
#include <json.hpp>

namespace gltfpp {
	inline namespace v1 {
		BETTER_ENUM(BufferViewTarget, int32_t, ARRAY_BUFFER = 34962, ELEMENT_ARRAY_BUFFER = 34963)

		struct BufferView {
			BufferView()
			    : byteStride{0} {
			}

			BOOST_HANA_DEFINE_STRUCT(BufferView,
			                         (option<BufferViewTarget>, target),
			                         (option<std::string>, name),
			                         (defaulted<uint8_t>, byteStride),
			                         (option<nlohmann::json>, extensions),
			                         (option<nlohmann::json>, extras));
			gsl::span<byte const> span;
			Buffer const *buffer;
		};

		auto parse(BufferView &) noexcept;
	}    // namespace v1
}    // namespace gltfpp

max_stars_repo_name: mmha/gltfpp

max_issues_repo_name: mmha/gltfpp

max_forks_repo_name: mmha/gltfpp


text:
/**
 * @file batchv_zdotu_sub.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date
 *
 * @precisions normal z -> c d s
 *
 **/

#include <cblas.h>
#include "bblas.h"

#define COMPLEX

void batchv_zdotu_sub(
		  const int *n,
		  BBLAS_Complex64_t const * const * x,
		  const int *incx,
		  BBLAS_Complex64_t const * const * y,
		  const int *incy,
		  BBLAS_Complex64_t *dotu,
		  const int batch_count, int* info)
{
	/* Local variables */
	int batch_iter = 0;
	char func_name[15] = "batchv_zdotu";

	/*initialize the result */
	for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
	  {
	    dotu[batch_iter] = (BBLAS_Complex64_t)0.0;
	  }

	/* Check input arguments */
	if (batch_count < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_BATCH_COUNT, -1);
	}

	for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
	  {
	    if (n[batch_iter] < 0)
	      {
		xerbla_batch(func_name, BBLAS_ERR_N, batch_iter);
		info[batch_iter] = BBLAS_ERR_N;
	      }
	    if (incx[batch_iter] < 1)
	      {
		xerbla_batch(func_name, BBLAS_ERR_INCX, batch_iter);
		info[batch_iter] = BBLAS_ERR_INCX;
	      }
	    if (incy[batch_iter] < 1)
	      {
		xerbla_batch(func_name, BBLAS_ERR_INCY, batch_iter);
		info[batch_iter] = BBLAS_ERR_INCY;
	      }

	    cblas_zdotu_sub(
		      n[batch_iter],
		      (void *)x[batch_iter],
		      incx[batch_iter],
		      (void *)y[batch_iter],
		      incy[batch_iter],
		      &dotu[batch_iter]);
	    /* Successful */
	    info[batch_iter] = BBLAS_SUCCESS;
	  } /* End fixed size for loop */
}
#undef COMPLEX

max_stars_repo_name: mawussi/BBLAS-group


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code for EOBNRv2HM reduced order model (non-spinning version).
 * See CQG 31 195010, 2014, arXiv:1402.4146 for details on the reduced order method.
 * See arXiv:1106.1021 for the EOBNRv2HM model.
 *
 * Borrows from the SEOBNR ROM LAL code written by Michael Puerrer and John Veitch.
 *
 * Put the untared data in the directory designated by the environment variable ROM_DATA_PATH.
 *
 * Parameter ranges:
 *   q = 1-12 (almost)
 *   No spin
 *   Mtot >= 10Msun for fstart=8Hz
 *
 */

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "EOBNRv2HMROMstruct.h"
#include "EOBNRv2HMROM.h"

/*************************************************/
/********* Some general definitions **************/

/* Number and list of modes to generate - to be modified ultimately to allow for a selection of the desired mode(s) */
/* By convention the first mode of the list is used to set phiRef */
/* nbmodemax = 5 has been defined in the header */
const int listmode[nbmodemax][2] = { {2,2}, {2,1}, {3,3}, {4,4}, {5,5} };

/* Maximal mass ratio covered by the model - q=12 (almost) for now */
static const double q_max = 11.9894197212;
/* Minimal geometric frequency covered by the model - f=8Hz for M=10Msol for now */
static const double Mf_ROM_min = 0.0003940393857519091;
/* Maximal geometric frequency covered by the model - Mf=0.285 for the 55 mode - used as default */
static const double Mf_ROM_max = 0.285;
/* Total mass (in units of solar mass) used to generate the ROM - used to restore the correct amplitude (global factor M) when coming back to physical units */
static const double M_ROM = 10.;

/* Define the number of points in frequency used by the SVD, identical for all modes */
static const int nbfreq = 300;
/* Define the number of training waveforms used by the SVD, identical for all modes */
static const int nbwf = 301;

/******************************************************************/
/********* Definitions for the persistent structures **************/

/* SEOBNR-ROM structures are generalized to lists */
ListmodesEOBNRHMROMdata* __EOBNRv2HMROM_data_init = NULL; /* for initialization only */
ListmodesEOBNRHMROMdata** const __EOBNRv2HMROM_data = &__EOBNRv2HMROM_data_init;
ListmodesEOBNRHMROMdata_interp* __EOBNRv2HMROM_interp_init = NULL; /* for initialization only */
ListmodesEOBNRHMROMdata_interp** const __EOBNRv2HMROM_interp = &__EOBNRv2HMROM_interp_init;
/* Tag indicating whether the data has been loaded and interpolated */
int __EOBNRv2HMROM_setup = FAILURE; /* To be set to SUCCESS after initialization*/

/********************* Miscellaneous ********************/

/* Return the closest higher power of 2 */
static size_t NextPow2(const size_t n) {
  return 1 << (size_t) ceil(log2(n));
}

/* Arbitrary tuned q-dependent functions by which the frequencies for the 44 and 55 modes have been multiplied (to put the ringdowns at the same Mf). The frequencies stored in the data for the 44 and 55 modes are the rescaled ones, not the original ones. */
static double Scaling44(const double q) {
  return 1.-1./4.*exp(-(q-1.)/5.);
}
static double Scaling55(const double q) {
  return 1.-1./3.*exp(-(q-1.)/5.);
}

/* Function evaluating eta as a function of q */
static double EtaOfq(const double q) {
  return q/(1.+q)/(1.+q);
}
/* Function evaluating delta m/m = (m1-m2)/(m1+m2) as a function of q */
static double DeltaOfq(const double q) {
  return( (q-1.)/(q+1.) );
}

/* Fit of the frequency of the 22 mode at the peak amplitude - from table III in the EOBNRv2HM paper, Pan&al 1106 */
static double omega22peakOfq(const double q) {
  double eta = EtaOfq(q);
  return 0.2733 + 0.2316*eta + 0.4463*eta*eta;
}

/* Amplitude factors scaled out for each mode; this does not include the global amplitude factor scaled out from all modes. */
static double ModeAmpFactor(const int l, const int m, const double q) {
  double eta = EtaOfq(q);
  if( l==2 && m==2 ) return(sqrt(eta));
  else if( l==2 && m==1 ) return( sqrt(eta)*DeltaOfq(q) );
  else if( l==3 && m==3 ) return( sqrt(eta)*DeltaOfq(q) );
  else if( l==4 && m==4 ) return( sqrt(Scaling44(q))*sqrt(eta)*(1.-3.*eta) );
  else if( l==5 && m==5 ) return( pow(Scaling55(q), 1./6)*sqrt(eta)*DeltaOfq(q)*(1.-2.*eta) );
  else {
    fprintf(stderr, "Unknown mode (%d,%d) for the amplitude factor.\n", l, m);
    return(FAILURE);
  }
}

/********************* Functions to initialize and cleanup contents of data structures ********************/
void EOBNRHMROMdata_Init(EOBNRHMROMdata **data) {
  if(!data) exit(1);
  /* Create storage for structures */
  if(!*data) *data=malloc(sizeof(EOBNRHMROMdata));
  else
  {
    EOBNRHMROMdata_Cleanup(*data);
  }
  gsl_set_error_handler(&Err_Handler);
  (*data)->q = gsl_vector_alloc(nbwf);
  (*data)->freq = gsl_vector_alloc(nbfreq);
  (*data)->Camp = gsl_matrix_alloc(nk_amp,nbwf);
  (*data)->Cphi = gsl_matrix_alloc(nk_phi,nbwf);
  (*data)->Bamp = gsl_matrix_alloc(nbfreq,nk_amp);
  (*data)->Bphi = gsl_matrix_alloc(nbfreq,nk_phi);
  (*data)->shifttime = gsl_vector_alloc(nbwf);
  (*data)->shiftphase = gsl_vector_alloc(nbwf);
}
void EOBNRHMROMdata_interp_Init(EOBNRHMROMdata_interp **data_interp) {
  if(!data_interp) exit(1);
  /* Create storage for structures */
  if(!*data_interp) *data_interp=malloc(sizeof(EOBNRHMROMdata_interp));
  else
  {
    EOBNRHMROMdata_interp_Cleanup(*data_interp);
  }
  (*data_interp)->Camp_interp = NULL;
  (*data_interp)->Cphi_interp = NULL;
  (*data_interp)->shifttime_interp = NULL;
  (*data_interp)->shiftphase_interp = NULL;
}
void EOBNRHMROMdata_coeff_Init(EOBNRHMROMdata_coeff **data_coeff) {
  if(!data_coeff) exit(1);
  /* Create storage for structures */
  if(!*data_coeff) *data_coeff=malloc(sizeof(EOBNRHMROMdata_coeff));
  else
  {
    EOBNRHMROMdata_coeff_Cleanup(*data_coeff);
  }
  gsl_set_error_handler(&Err_Handler);
  (*data_coeff)->Camp_coeff = gsl_vector_alloc(nk_amp);
  (*data_coeff)->Cphi_coeff = gsl_vector_alloc(nk_phi);
  (*data_coeff)->shifttime_coeff = malloc(sizeof(double));
  (*data_coeff)->shiftphase_coeff = malloc(sizeof(double));
}
void EOBNRHMROMdata_Cleanup(EOBNRHMROMdata *data /* data to destroy */) {
  if(data->q) gsl_vector_free(data->q);
  if(data->freq) gsl_vector_free(data->freq);
  if(data->Camp) gsl_matrix_free(data->Camp);
  if(data->Cphi) gsl_matrix_free(data->Cphi);
  if(data->Bamp) gsl_matrix_free(data->Bamp);
  if(data->Bphi) gsl_matrix_free(data->Bphi);
  if(data->shifttime) gsl_vector_free(data->shifttime);
  if(data->shiftphase) gsl_vector_free(data->shiftphase);
  free(data);
}
void EOBNRHMROMdata_coeff_Cleanup(EOBNRHMROMdata_coeff *data_coeff) {
  if(data_coeff->Camp_coeff) gsl_vector_free(data_coeff->Camp_coeff);
  if(data_coeff->Cphi_coeff) gsl_vector_free(data_coeff->Cphi_coeff);
  if(data_coeff->shifttime_coeff) free(data_coeff->shifttime_coeff);
  if(data_coeff->shiftphase_coeff) free(data_coeff->shiftphase_coeff);
  free(data_coeff);
}
void EOBNRHMROMdata_interp_Cleanup(EOBNRHMROMdata_interp *data_interp) {
  if(data_interp->Camp_interp) SplineList_Destroy(data_interp->Camp_interp);
  if(data_interp->Cphi_interp) SplineList_Destroy(data_interp->Cphi_interp);
  if(data_interp->shifttime_interp) SplineList_Destroy(data_interp->shifttime_interp);
  if(data_interp->shiftphase_interp) SplineList_Destroy(data_interp->shiftphase_interp);
  free(data_interp);
}

/* Read binary ROM data for frequency vectors, coefficients matrices, basis functions matrices, and shiftvectors in time and phase */
int Read_Data_Mode(const char dir[], const int mode[2], EOBNRHMROMdata *data) {
  /* Load binary data for amplitude and phase spline coefficients as computed in Mathematica */
  int ret = SUCCESS;
  char* file_q = malloc(strlen(dir)+64);
  char* file_freq = malloc(strlen(dir)+64);
  char* file_Camp = malloc(strlen(dir)+64);
  char* file_Cphi = malloc(strlen(dir)+64);
  char* file_Bamp = malloc(strlen(dir)+64);
  char* file_Bphi = malloc(strlen(dir)+64);
  char* file_shifttime = malloc(strlen(dir)+64);
  char* file_shiftphase = malloc(strlen(dir)+64);
  sprintf(file_q, "%s", "EOBNRv2HMROM_q.dat"); /* The q vector is the same for all modes */
  sprintf(file_freq, "%s%d%d%s", "EOBNRv2HMROM_freq_", mode[0], mode[1], ".dat");
  sprintf(file_Camp, "%s%d%d%s", "EOBNRv2HMROM_Camp_", mode[0], mode[1], ".dat");
  sprintf(file_Cphi, "%s%d%d%s", "EOBNRv2HMROM_Cphi_", mode[0], mode[1], ".dat");
  sprintf(file_Bamp, "%s%d%d%s", "EOBNRv2HMROM_Bamp_", mode[0], mode[1], ".dat");
  sprintf(file_Bphi, "%s%d%d%s", "EOBNRv2HMROM_Bphi_", mode[0], mode[1], ".dat");
  sprintf(file_shifttime, "%s%d%d%s", "EOBNRv2HMROM_shifttime_", mode[0], mode[1], ".dat");
  sprintf(file_shiftphase, "%s%d%d%s", "EOBNRv2HMROM_shiftphase_", mode[0], mode[1], ".dat");
  ret |= Read_Vector(dir, file_q, data->q);
  ret |= Read_Vector(dir, file_freq, data->freq);
  ret |= Read_Matrix(dir, file_Camp, data->Camp);
  ret |= Read_Matrix(dir, file_Cphi, data->Cphi);
  ret |= Read_Matrix(dir, file_Bamp, data->Bamp);
  ret |= Read_Matrix(dir, file_Bphi, data->Bphi);
  ret |= Read_Vector(dir, file_shifttime, data->shifttime);
  ret |= Read_Vector(dir, file_shiftphase, data->shiftphase);
  free(file_q);
  free(file_freq);
  free(file_Camp);
  free(file_Cphi);
  free(file_Bamp);
  free(file_Bphi);
  free(file_shifttime);
  free(file_shiftphase);
  return(ret);
}

/* Function interpolating the data in matrix/vector form produces an interpolated data in the form of SplineLists */
int Interpolate_Spline_Data(const EOBNRHMROMdata *data, EOBNRHMROMdata_interp *data_interp) {

  gsl_set_error_handler(&Err_Handler);
  SplineList* splinelist;
  gsl_spline* spline;
  gsl_interp_accel* accel;
  gsl_vector* matrixline;
  gsl_vector* vector;
  int j;

  /* Interpolating Camp */
  splinelist = data_interp->Camp_interp;
  for (j = 0; j<nk_amp; j++) {
    matrixline = gsl_vector_alloc(nbwf);
    gsl_matrix_get_row(matrixline, data->Camp, j);

    accel = gsl_interp_accel_alloc();
    spline = gsl_spline_alloc(gsl_interp_cspline, nbwf);
    gsl_spline_init(spline, gsl_vector_const_ptr(data->q, 0), gsl_vector_const_ptr(matrixline, 0), nbwf);

    splinelist = SplineList_AddElementNoCopy(splinelist, spline,  accel, j);
    gsl_vector_free(matrixline);
  }
  data_interp->Camp_interp = splinelist;

  /* Interpolating Cphi */
  splinelist = data_interp->Cphi_interp;
  for (j = 0; j<nk_phi; j++) {
    matrixline = gsl_vector_alloc(nbwf);
    gsl_matrix_get_row(matrixline, data->Cphi, j);

    accel = gsl_interp_accel_alloc();
    spline = gsl_spline_alloc(gsl_interp_cspline, nbwf);
    gsl_spline_init(spline, gsl_vector_const_ptr(data->q, 0), gsl_vector_const_ptr(matrixline, 0), nbwf);

    splinelist = SplineList_AddElementNoCopy(splinelist, spline,  accel, j);
    gsl_vector_free(matrixline);
  }
  data_interp->Cphi_interp = splinelist;

  /* Interpolating shifttime */
  splinelist = data_interp->shifttime_interp;
  vector = data->shifttime;

  accel = gsl_interp_accel_alloc();
  spline = gsl_spline_alloc(gsl_interp_cspline, nbwf);
  gsl_spline_init(spline, gsl_vector_const_ptr(data->q, 0), gsl_vector_const_ptr(vector, 0), nbwf);

  splinelist = SplineList_AddElementNoCopy(NULL, spline,  accel, 0); /* This SplineList has only 1 element */
  data_interp->shifttime_interp = splinelist;

  /* Interpolating shiftphase */
  splinelist = data_interp->shiftphase_interp;
  vector = data->shiftphase;

  accel = gsl_interp_accel_alloc();
  spline = gsl_spline_alloc(gsl_interp_cspline, nbwf);
  gsl_spline_init(spline, gsl_vector_const_ptr(data->q, 0), gsl_vector_const_ptr(vector, 0), nbwf);

  splinelist = SplineList_AddElementNoCopy(NULL, spline,  accel, 0); /* This SplineList has only 1 element */
  data_interp->shiftphase_interp = splinelist;

  return SUCCESS;
}

/* Function taking as input interpolated data in the form of SplineLists
 * evaluates for a given q the projection coefficients and shifts in time and phase
*/
int Evaluate_Spline_Data(const double q, const EOBNRHMROMdata_interp* data_interp, EOBNRHMROMdata_coeff* data_coeff){

  SplineList* splinelist;
  /* Evaluating the vector of projection coefficients for the amplitude */
  for (int j=0; j<nk_amp; j++) {
    splinelist = SplineList_GetElement(data_interp->Camp_interp, j);
    gsl_vector_set(data_coeff->Camp_coeff, j, gsl_spline_eval(splinelist->spline, q, splinelist->accel));
  }
  /* Evaluating the vector of projection coefficients for the phase */
  for (int j=0; j<nk_phi; j++) {
    splinelist = SplineList_GetElement(data_interp->Cphi_interp, j);
    gsl_vector_set(data_coeff->Cphi_coeff, j, gsl_spline_eval(splinelist->spline, q, splinelist->accel));
  }
  /* Evaluating the shift in time */
  splinelist = SplineList_GetElement(data_interp->shifttime_interp, 0); /* This SplineList has only one element */
  *(data_coeff->shifttime_coeff) = gsl_spline_eval(splinelist->spline, q, splinelist->accel);
  /* Evaluating the shift in phase */
  splinelist = SplineList_GetElement(data_interp->shiftphase_interp, 0); /* This SplineList has only one element */
  *(data_coeff->shiftphase_coeff) = gsl_spline_eval(splinelist->spline, q, splinelist->accel);

  return SUCCESS;
}

/*
 * Core function for computing the ROM waveform.
 * Evaluates projection coefficients and shifts in time and phase at desired q.
 * Construct 1D splines for amplitude and phase.
 * Compute strain waveform from amplitude and phase.
*/
int EOBNRv2HMROMCore(
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,
  int nbmode,
  double deltatRef,
  double phiRef,
  double fRef,
  double Mtot_sec,
  double q,
  double distance,
  int setphiRefatfRef)
{
  int ret = SUCCESS;
  int j;
  double tpeak22estimate = 0;
  /* Check output arrays */
  if(!listhlm) exit(1);
  if(*listhlm)
  {
    printf("Error: (*listhlm) is supposed to be NULL, but got %p\n",(*listhlm));
    exit(1);
  }
  /* Check number of modes */
  if(nbmode<1 || nbmode>nbmodemax) {
    printf("Error: incorrect number of modes: %d", nbmode);
    exit(1);
  }

  /* Check if the data has been set up */
  if(__EOBNRv2HMROM_setup) {
    printf("Error: the ROM data has not been set up\n");
    exit(1);
  }
  /* Set the global pointers to data */
  ListmodesEOBNRHMROMdata* listdata = *__EOBNRv2HMROM_data;
  ListmodesEOBNRHMROMdata_interp* listdata_interp = *__EOBNRv2HMROM_interp;

  /* Global amplitude prefactor - includes total mass scaling, Fourier scaling, distance scaling, and undoing an additional arbitrary scaling */
  double Mtot_msol = Mtot_sec / MTSUN_SI; /* Mtot_msol and M_ROM in units of solar mass */
  double amp0 = (Mtot_msol/M_ROM) * Mtot_sec * 1.E-16 * 1.E6 * PC_SI / distance;

  /* Highest allowed geometric frequency for the first mode of listmode in the ROM - used for fRef
   * by convention, we use the first mode of listmode (presumably the 22) for phiref */
  ListmodesEOBNRHMROMdata* listdata_ref = ListmodesEOBNRHMROMdata_GetMode(listdata, listmode[0][0], listmode[0][1]);
  EOBNRHMROMdata* data_ref = listdata_ref->data;
  double Mf_ROM_max_ref = gsl_vector_get(data_ref->freq, nbfreq-1);
  /* Convert to geometric units the reference time and frequency */
  double deltatRef_geom = deltatRef / Mtot_sec;
  double fRef_geom = fRef * Mtot_sec;

  /* Enforce allowed geometric frequency range for fRef */
  /* In case the user asks for a reference frequency higher than covered by the ROM, we keep it that way as we will just 0-pad the waveform (and do it anyway for some modes) */
  if (fRef_geom > Mf_ROM_max_ref || fRef_geom == 0)
    fRef_geom = Mf_ROM_max_ref; /* If fRef > fhigh or 0 we reset fRef to default value of cutoff frequency for the first mode of the list (presumably the 22 mode) */
  if (0 < fRef_geom && fRef_geom < Mf_ROM_min) {
    //printf("WARNING: Reference frequency Mf_ref=%g is smaller than lowest frequency in ROM Mf=%g. Setting it to the lowest frequency in ROM.\n", fRef_geom, Mf_ROM_min);
    fRef_geom = Mf_ROM_min;
  }

  /* Internal storage for the projection coefficients and shifts in time and phase */
  /* Initialized only once, and reused for the different modes */
  EOBNRHMROMdata_coeff *data_coeff = NULL;
  EOBNRHMROMdata_coeff_Init(&data_coeff);

  /* The phase change imposed by phiref, from the phase of the first mode in the list - to be set in the first step of the loop on the modes */
  double phase_change_ref = 0;

  /* Main loop over the modes */
  for(int i=0; i<nbmode; i++ ){
    int l = listmode[i][0];
    int m = listmode[i][1];

    /* Getting the relevant modes in the lists of data */
    ListmodesEOBNRHMROMdata* listdata_mode = ListmodesEOBNRHMROMdata_GetMode(listdata, l, m);
    ListmodesEOBNRHMROMdata_interp* listdata_interp_mode = ListmodesEOBNRHMROMdata_interp_GetMode(listdata_interp, l, m);

    /* Evaluating the projection coefficients and shift in time and phase */
    ret |= Evaluate_Spline_Data(q, listdata_interp_mode->data_interp, data_coeff);

    /* Evaluating the unnormalized amplitude and unshifted phase vectors for the mode */
    /* Notice a change in convention: B matrices are transposed with respect to the B matrices in SEOBNRROM */
    /* amp_pts = Bamp . Camp_coeff */
    /* phi_pts = Bphi . Cphi_coeff */
    gsl_vector* amp_f = gsl_vector_alloc(nbfreq);
    gsl_vector* phi_f = gsl_vector_alloc(nbfreq);
    //clock_t begblas = clock();
    gsl_blas_dgemv(CblasNoTrans, 1.0, listdata_mode->data->Bamp, data_coeff->Camp_coeff, 0.0, amp_f);
    gsl_blas_dgemv(CblasNoTrans, 1.0, listdata_mode->data->Bphi, data_coeff->Cphi_coeff, 0.0, phi_f);
    //clock_t endblas = clock();
    //printf("Mode (%d,%d) Blas time: %g s\n", l, m, (double)(endblas - begblas) / CLOCKS_PER_SEC);

     /* The downsampled frequencies for the mode - we undo the rescaling of the frequency for the 44 and 55 modes */
    gsl_vector* freq_ds = gsl_vector_alloc(nbfreq);
    gsl_vector_memcpy(freq_ds, listdata_mode->data->freq);
    if ( l==4 && m==4) gsl_vector_scale( freq_ds, 1./Scaling44(q));
    if ( l==5 && m==5) gsl_vector_scale( freq_ds, 1./Scaling55(q));

    /* Evaluating the shifts in time and phase - conditional scaling for the 44 and 55 modes */
    /* Note: the stored values of 'shifttime' correspond actually to 2pi*Deltat */
    SplineList* shifttime_splinelist = listdata_interp_mode->data_interp->shifttime_interp;
    SplineList* shiftphase_splinelist = listdata_interp_mode->data_interp->shiftphase_interp;
    double twopishifttime;
    if( l==4 && m==4) {
      twopishifttime = gsl_spline_eval(shifttime_splinelist->spline, q, shifttime_splinelist->accel) * Scaling44(q);
    }
    else if( l==5 && m==5) {
      twopishifttime = gsl_spline_eval(shifttime_splinelist->spline, q, shifttime_splinelist->accel) * Scaling55(q);
    }
    else {
      twopishifttime = gsl_spline_eval(shifttime_splinelist->spline, q, shifttime_splinelist->accel);
    }
    double shiftphase = gsl_spline_eval(shiftphase_splinelist->spline, q, shiftphase_splinelist->accel);

    /* If first mode in the list, assumed to be the 22 mode, set totalshifttime and phase_change_ref */
    if( i==0 ) {
      if(l==2 && m==2) {
      /* Setup 1d cubic spline for the phase of the 22 mode */
      gsl_interp_accel* accel_phi22 = gsl_interp_accel_alloc();
      gsl_spline* spline_phi22 = gsl_spline_alloc(gsl_interp_cspline, nbfreq);
      gsl_spline_init(spline_phi22, gsl_vector_const_ptr(freq_ds,0), gsl_vector_const_ptr(phi_f,0), nbfreq);
      /* Compute the shift in time needed to set the peak of the 22 mode roughly at deltatRef */
      /* We use the SPA formula tf = -(1/2pi)*dPsi/df to estimate the correspondence between frequency and time */
      /* The frequency corresponding to the 22 peak is omega22peak/2pi, with omega22peak taken from the fit to NR in Pan&al 1106 EOBNRv2HM paper */
      double f22peak = fmin(omega22peakOfq(q)/(2*PI), Mf_ROM_max_ref); /* We ensure we evaluate the spline within its range */
      /* Note : twopishifttime is almost 0 (order 1e-8) by construction for the 22 mode, so it does not intervene here */
      tpeak22estimate = -1./(2*PI) * gsl_spline_eval_deriv(spline_phi22, f22peak, accel_phi22);
      /* Determine the change in phase (to be propagated to all modes) required to have phi22(fRef) = 2*phiRef */
      if(setphiRefatfRef) {
        phase_change_ref = 2*phiRef + (gsl_spline_eval(spline_phi22, fRef_geom, accel_phi22) - (twopishifttime - 2*PI*tpeak22estimate + 2*PI*deltatRef_geom) * fRef_geom - shiftphase);
      }
      else {
        phase_change_ref = 2*phiRef;
      }
      gsl_spline_free(spline_phi22);
      gsl_interp_accel_free(accel_phi22);
      }
      else {
      	printf("Error: the first mode in listmode must be the 22 mode to set the changes in phase and time \n");
      	return FAILURE;
      }
    }
    /* Total shift in time, and total change in phase for this mode */
    double totaltwopishifttime = twopishifttime - 2*PI*tpeak22estimate + 2*PI*deltatRef_geom;
    double constphaseshift = m/2. * phase_change_ref + shiftphase;

    //
    //printf("deltatRef_geom: %g\n", deltatRef_geom);
    //printf("freq_ds 0: %g\n", gsl_vector_get(freq_ds, 0));

    //printf("2*PI*deltatRef_geom * gsl_vector_get(freq_ds, 0), phase_change_ref: %g, %g\n", 2*PI*deltatRef_geom * gsl_vector_get(freq_ds, 0), phase_change_ref);

    /* Initialize the complex series for the mode */
    CAmpPhaseFrequencySeries *modefreqseries = NULL;
    int len = (int) freq_ds->size;
    CAmpPhaseFrequencySeries_Init(&modefreqseries, len);

    /* Mode-dependent complete amplitude prefactor */
    double amp_pre = amp0 * ModeAmpFactor( l, m, q);

    /* Final result for the mode */
    /* Scale and set the amplitudes (amplitudes are real at this stage)*/
    gsl_vector_scale(amp_f, amp_pre);
    gsl_vector_memcpy(modefreqseries->amp_real, amp_f);
    gsl_vector_set_zero(modefreqseries->amp_imag); /* Amplitudes are real at this stage */
    /* Add the linear term and the constant (including the shift to phiRef), and set the phases */
    gsl_vector_scale(phi_f, -1.); /* Change the sign of the phases: ROM convention Psi=-phase */
    gsl_blas_daxpy(totaltwopishifttime, freq_ds, phi_f); /*Beware: here freq_ds must still be in geometric units*/
    gsl_vector_add_constant(phi_f, constphaseshift);
    gsl_vector_memcpy(modefreqseries->phase, phi_f);

//
//printf("first phi_f: %g\n", gsl_vector_get(phi_f, 0 ));
//printf("last phi_f: %g\n", gsl_vector_get(phi_f, phi_f->size -1 ));

    /* Scale (to physical units) and set the frequencies */
    gsl_vector_scale(freq_ds, 1./Mtot_sec);
    gsl_vector_memcpy(modefreqseries->freq, freq_ds);

    /* Append the computed mode to the ListmodesCAmpPhaseFrequencySeries structure */
    *listhlm = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*listhlm, modefreqseries, l, m);

    //
    //printf("Mode (%d,%d) generated:\n", l, m);
    //for( int j=0; j<len; j++ ){
    //  printf("%g %g %g %g\n", gsl_vector_get(modefreqseries->freq, j), gsl_vector_get(modefreqseries->amp_real, j), gsl_vector_get(modefreqseries->amp_imag, j), gsl_vector_get(modefreqseries->phase, j));
    //}

    /* Cleanup for the mode */
    gsl_vector_free(freq_ds);
    gsl_vector_free(amp_f);
    gsl_vector_free(phi_f);
  }

  /* Cleanup of the coefficients data structure */
  EOBNRHMROMdata_coeff_Cleanup(data_coeff);

  return(SUCCESS);
}

/* Compute waveform in downsampled frequency-amplitude-phase format */
int SimEOBNRv2HMROM(
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,  /* Output: list of modes in Frequency-domain amplitude and phase form */
  int nbmode,                                    /* Number of modes to generate (starting with the 22) */
  double deltatRef,                              /* Time shift so that the peak of the 22 mode occurs at deltatRef */
  double phiRef,                                 /* Phase at reference frequency */
  double fRef,                                   /* Reference frequency (Hz); 0 defaults to fLow */
  double m1SI,                                   /* Mass of companion 1 (kg) */
  double m2SI,                                   /* Mass of companion 2 (kg) */
  double distance,                               /* Distance of source (m) */
  int setphiRefatfRef)                           /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) */
{
  /* Get masses in terms of solar mass */
  double mass1 = m1SI / MSUN_SI;
  double mass2 = m2SI / MSUN_SI;
  double Mtot = mass1 + mass2;
  double q = fmax(mass1/mass2, mass2/mass1);    /* Mass-ratio >1 by convention*/
  double Mtot_sec = Mtot * MTSUN_SI; /* Total mass in seconds */

  if ( q > q_max ) {
    //printf( "Error - %s: q out of range!\nEOBNRv2HMROM is only available for a mass ratio in the range q <= %g.\n", __func__, q_max);
    return FAILURE;
  }

  /* Set up (load and build interpolation) ROM data if not setup already */
  //clock_t beg = clock();
  EOBNRv2HMROM_Init_DATA();
  //clock_t end = clock();
  //printf("Initialization time: %g s\n", (double)(end - beg) / CLOCKS_PER_SEC);

  //beg = clock();
  int retcode = EOBNRv2HMROMCore(listhlm, nbmode, deltatRef, phiRef, fRef, Mtot_sec, q, distance, setphiRefatfRef);
  //end = clock();
  //printf("ROM evaluation time: %g s\n", (double)(end - beg) / CLOCKS_PER_SEC);

  return(retcode);
}

/* Setup EOBNRv2HMROM model using data files installed in $ROM_DATA_PATH */
int EOBNRv2HMROM_Init_DATA(void) {
  if (!__EOBNRv2HMROM_setup) return SUCCESS;

  int ret=FAILURE;
  char *envpath=NULL;
  char path[32768];
  char *brkt,*word;
  envpath=getenv("ROM_DATA_PATH");
  if(!envpath) {
    printf("Error: the environment variable ROM_DATA_PATH, giving the path to the ROM data, seems undefined\n");
    return(FAILURE);
  }
  strncpy(path,envpath,sizeof(path));

#pragma omp critical
  {
    for(word=strtok_r(path,":",&brkt); word; word=strtok_r(NULL,":",&brkt))
      {
	ret = EOBNRv2HMROM_Init(word);
	if(ret == SUCCESS) break;
      }
    if(ret!=SUCCESS) {
      printf("Error: unable to find EOBNRv2HMROM data files in $ROM_DATA_PATH\n");
      exit(FAILURE);
    }
    __EOBNRv2HMROM_setup = ret;
  }
  return(ret);
}

/* Setup EOBNRv2HMROM model using data files installed in dir */
int EOBNRv2HMROM_Init(const char dir[]) {
  if(!__EOBNRv2HMROM_setup) {
    printf("Error: EOBNRHMROMdata was already set up!");
    exit(1);
  }

  int ret = SUCCESS;
  ListmodesEOBNRHMROMdata* listdata = *__EOBNRv2HMROM_data;
  ListmodesEOBNRHMROMdata_interp* listdata_interp = *__EOBNRv2HMROM_interp;

  for (int j=0; j<nbmodemax; j++) { /* At setup, we initialize all available modes anyway */

    EOBNRHMROMdata* data = NULL;
    EOBNRHMROMdata_Init(&data);
    ret |= Read_Data_Mode( dir, listmode[j], data);
    if (!ret) {
      listdata = ListmodesEOBNRHMROMdata_AddModeNoCopy( listdata, data, listmode[j][0], listmode[j][1]);

      EOBNRHMROMdata_interp* data_interp = NULL;
      EOBNRHMROMdata_interp_Init(&data_interp);
      ret |= Interpolate_Spline_Data(data, data_interp);
      if (!ret) listdata_interp = ListmodesEOBNRHMROMdata_interp_AddModeNoCopy( listdata_interp, data_interp, listmode[j][0], listmode[j][1]);
    }
  }

  __EOBNRv2HMROM_setup = ret;
  if (!ret) {
    *__EOBNRv2HMROM_data = listdata;
    *__EOBNRv2HMROM_interp = listdata_interp;
  }
  return(ret);
}

/* Non-spinning merger TaylorF2 waveform, copied and condensed from LAL */
/* Used by SimEOBNRv2HMROMExtTF2 to extend the signal to arbitrarily low frequencies */
static void TaylorF2nonspin(
		double *amp,                            /**< FD waveform amplitude (modulus)*/
		double *phase,                          /**< FD waveform phase */
		const double *freqs,                    /**< frequency points at which to evaluate the waveform (Hz) */
		const int size,                         /** number of freq samples */
		const double m1_SI,                     /**< mass of companion 1 (kg) */
		const double m2_SI,                     /**< mass of companion 2 (kg) */
		const double distance,                  /** distance (m) */
		const double imatch                     /**< index at which to match phase;
							   assumes arrays are preloaded at imatch and imatch+1
							   with the required result */
		     )
{
  //The meat of this computation is copied from LAL: XLALSimInspiralPNPhasing_F2
  //We dont need the spin terms
  double m1 = m1_SI / MSUN_SI;
  double m2 = m2_SI / MSUN_SI;
  double mtot = m1 + m2;
  double d = (m1 - m2) / (m1 + m2);
  double eta = m1*m2/mtot/mtot;
  double m1M = m1/mtot;
  double m2M = m2/mtot;
  double m_sec = mtot * MTSUN_SI;
  double piM = PI * m_sec;

  double pfaN = 3.L/(128.L * eta);

  /* Non-spin phasing terms - see arXiv:0907.0700, Eq. 3.18 */
  double pfav0 = 1.L;
  double pfav2 = 5.L*(743.L/84.L + 11.L * eta)/9.L;
  double pfav3 = -16.L*PI;
  double pfav4 = 5.L*(3058.673L/7.056L + 5429.L/7.L * eta
		+ 617.L * eta*eta)/72.L;
  double pfav5 = 5.L/9.L * (7729.L/84.L - 13.L * eta) * PI;
  double pfalogv5 = 5.L/3.L * (7729.L/84.L - 13.L * eta) * PI;
  double pfav6 = (11583.231236531L/4.694215680L
	    - 640.L/3.L * PI * PI - 6848.L/21.L*GAMMA)
    + eta * (-15737.765635L/3.048192L
	     + 2255./12. * PI * PI)
    + eta*eta * 76055.L/1728.L
    - eta*eta*eta * 127825.L/1296.L;
  pfav6 += (-6848.L/21.L)*log(4.);
  double pfalogv6 = -6848.L/21.L;
  double pfav7 = PI * ( 77096675.L/254016.L
		      + 378515.L/1512.L * eta - 74045.L/756.L * eta*eta);

  /* Non-spin 2-2 amplitude terms (Blanchet LRR)*/
  double a2 = ( -107 + 55*eta ) / 42.;
  double a3 = 2*PI;
  double a4 = ( ( 2047.*eta - 7483. ) * eta - 2173. ) / 1512.;
  /* Blanchet has more terms, but there should diminishing returns:
     expect v^5 ~ 1e-5 and the higher terms are more complicated and, indeed, complex */

  //Lead coefficients
  //double amp0 = -4. * m1 * m2 / distance * C_SI * MTSUN_SI * MTSUN_SI * sqrt(PI/12.L); //(from LAL)
  double amp0B = 2. * m1 * m2 / distance * C_SI * MTSUN_SI * MTSUN_SI * sqrt(16*PI/5.L); //Based on Blanchet-LRR (327)
  //Note: amp0B = -4 * sqrt( 3/5) * amp0;
  double FTaN =  32.0 * eta*eta / 5.0;
  //printf("eta=%g\n",eta);
  //Compute raw TaylorF2
  int i;
  for (i = 0; i < size; i++) {
    double f = freqs[i];
    double v = cbrt(piM*f);
    double logv = log(v);
    double v2 = v*v;
    double v5 = v2*v2*v;
    double v10 = v5*v5;

    //printf("taylorf2: f=%g  v=%g\n",f,v);
    double phasing=0;
    phasing = pfav7 * v;
    phasing = (phasing + pfav6 + pfalogv6 * logv) * v;
    phasing = (phasing + pfav5 + pfalogv5 * logv) * v;
    phasing = (phasing + pfav4) * v;
    phasing = (phasing + pfav3) * v;
    phasing = (phasing + pfav2) * v2;
    phasing += 1;
    phasing *=pfaN;

    double amp22fac;
    amp22fac = a4*v;
    amp22fac = ( amp22fac + a3 ) * v;
    amp22fac = ( amp22fac + a2 ) * v2;
    amp22fac += 1.0;

    phasing /= v5;
    double flux = FTaN * v10;
    double dEnergy = -eta * v;
    phase[i]=phasing;
    //Notes for amplitude: Blanchet at leading order:
    /* mf=x^(3/2); fdot=3/2/m x^(1/2) xdot ~ 3/2/m x^(1/2) * (-1/16)*(4x)^5(-eta/5/m) = 96/5*eta/m^2 * x^(11/2) */
    /*-flux/dEnergy = 32.0 * eta*eta / 5.0 / eta *v^9 */
    /*--> -flux/dEnergy =  fdot / (3*v^2) [using x=v^2]*/
    //amp[i] = amp0 * sqrt(-dEnergy/flux) * v;  (Based on LAL)
    amp[i] = amp0B * amp22fac * v2 / sqrt(-flux/dEnergy * 3 * v2 );
    //printf("v=%g: a=%g ph=%g;  amp22fac=%g sqrt(v^-9)=%g\n",v,amp[i],phase[i],amp22fac,sqrt(v/v10));
    //Note ampB = - amp * 4*sqrt(3/5) * / sqrt(3) + higher-order = -4/sqrt(5)*( 1 + higher-order )
    //  ...possibly related to sph-harm normalization
    // HACK: Strangely it seems that an additional factor of 4/sqrt(5) is just right to nearly match the EOB wf FT
    amp[i] *= 4/sqrt(5);   //HACK

    //Here we depart from LAL, referencing phase and time-shift to two neighboring freq points
    //First we match the freq derivative
  }
}

/*Wrapper for waveform generation with possibly a combination of EOBNRv2HMROM and TaylorF2*/
/* Note: GenerateWaveform accepts masses and distances in SI units, whereas LISA params is in solar masses and Mpc */
/* Note: the extended waveform will now a different number of frequency points for each mode */
/* Note: phiRef is readjusted after the extension -- for the case where fRef is below the band covered by ROM, in which case the core function defaults fRef to the max geometric freq of the ROM */
int SimEOBNRv2HMROMExtTF2(
  ListmodesCAmpPhaseFrequencySeries **listhlm,   /* Output: list of modes in Frequency-domain amplitude and phase form */
  int nbmode,                                    /* Number of modes to generate (starting with the 22) */
  double Mf_match,                               /* Minimum frequency using EOBNRv2HMROM in inverse total mass units*/
  double minf,                                   /* Minimum frequency required */
  int tagexthm,                                  /* Tag to decide whether or not to extend the higher modes as well */
  double deltatRef,                              /* Time shift so that the peak of the 22 mode occurs at deltatRef */
  double phiRef,                                 /* Phase at reference frequency */
  double fRef,                                   /* Reference frequency (Hz); 0 defaults to fLow */
  double m1SI,                                   /* Mass of companion 1 (kg) */
  double m2SI,                                   /* Mass of companion 2 (kg) */
  double distance,                               /* Distance of source (m) */
  int setphiRefatfRef)                           /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) */
{//
  //printf("calling SimEOBNRv2HMROMExtTF2 with minf=%g\n", minf);
  //printf("params: %d %g %g %g %g %g %g %g %g\n", nbmode, Mf_match, minf, deltatRef, phiRef, fRef, m1SI, m2SI, distance);

  int ret,i;
  ListmodesCAmpPhaseFrequencySeries* listROM = NULL;
  //int lout=-1,mout=-1;
  int lout=-1,mout=-1;

  /* Generate the waveform with the ROM */
  ret = SimEOBNRv2HMROM(&listROM, nbmode, deltatRef, phiRef, fRef, m1SI, m2SI, distance, setphiRefatfRef);

  /* If the ROM waveform generation failed (e.g. parameters were out of bounds) return FAILURE */
  //if(ret==FAILURE)printf("SimEOBNRv2HMROMExtTF2: Generation of ROM for injection failed!\n");
  if(ret==FAILURE) return FAILURE;

  /* Main loop over the modes (as linked list) to perform the extension */
  /* The 2-2 mode will be extended by TaylorF2 model with the phase and time offset
  determined by matching conditions. All other modes will be extended as some
  sort of power-law fall-off in amplitude and power-law growth in phase.          */
  ListmodesCAmpPhaseFrequencySeries* listelement = listROM;
  while(listelement) {    // For each l-m (ie each listelement)

    /* Definitions: l,m, frequency series and length */
    int l = listelement->l;
    int m = listelement->m;

    //NOTE: temporary hack to allow avoiding power-law extension of higher modes which seems problematic
    //    if((l==2&&m==2) || tagexthm) {

    /* First we must compute a new frequency grid including a possible extension to lower frequencies*/
    gsl_vector *freq_new;
    gsl_vector* freq = listelement->freqseries->freq;
    int len = (int) freq->size;
    // Construct frequency grid extension on the geometric mean of the lowest few ROM frequencies after the matching point
    const int Navg=3;
    int imatch=-1;
    double f_match;
    if(Mf_match<=0) {
      f_match = Mf_ROM_min/(m1SI+m2SI)*MSUN_SI/MTSUN_SI;
    }
    else if(Mf_match>Mf_ROM_min) {
      f_match = Mf_match/(m1SI+m2SI)*MSUN_SI/MTSUN_SI;
    }
    else {
      printf("WARNING: Mf_match lower than the lowest geometric frequency of the ROM model - used instead\n");
      f_match = Mf_ROM_min/(m1SI+m2SI)*MSUN_SI/MTSUN_SI;
    }
    f_match = f_match*m/2; //Shift the matching frequency for non-22 modes to correspond to a comparable orbital freq.
    //printf("f_match=%g\n",f_match);
    for(i=0;i<len-Navg;i++){
      if(gsl_vector_get(freq,i)>f_match){
        imatch=i;
        break;
      }
    }
    if(imatch<0){
      printf("WARNING: f_match exceeds high-freq range of the ROM model\n");
      imatch=len-Navg-1;
    }
    double lnfmatch=log(gsl_vector_get(freq,imatch));
    double lnfmin=log(minf);
    double dlnf=(log(gsl_vector_get(freq,imatch+Navg))-lnfmatch)/Navg;
    double dffac=exp(dlnf);
    //The new grid will include the old grid, from imatch and after, plus adequate lower-freq
    int len_add = (lnfmatch-lnfmin)/dlnf;
    if(len_add<0) len_add=0;
    int len_new = len - imatch + len_add;
    //printf("extending waveform: len_add + len - imatch = len_new: %i + %i - %i = %i\n",len_add,len,imatch,len_new);
    /* construct the extended freq grid */
    freq_new=gsl_vector_alloc(len_new);
    for(i=len_add;i<len_new;i++)gsl_vector_set(freq_new,i,gsl_vector_get(freq,len-len_new+i));
    for(i=len_add;i>0;i--)gsl_vector_set(freq_new,i-1,gsl_vector_get(freq_new,i)/dffac);
    //for(i=0;i<len_new;i++)printf("%i: %g %g\n",i,(i>=len_new-len?freq->data[i-len_new+len]:0),freq_new->data[i]);

    //copy the old freqseries data to a new one and extend with power-law
    CAmpPhaseFrequencySeries* freqseries = listelement->freqseries;
    if(l==lout&&m==mout){  //write to file for debugging
      printf("Writing waveform-before for (%i,%i)\n",l,m);
      FILE *f = fopen("waveform-before.dat", "w");
      for(i=0;i<freqseries->freq->size;i++){
        fprintf(f,"%i %i %g  %g  %g  %g  %i\n",l,m,
        gsl_vector_get(freqseries->freq,i),
        gsl_vector_get(freqseries->amp_real,i),
        gsl_vector_get(freqseries->amp_imag,i),
        gsl_vector_get(freqseries->phase,i),
        i);
      }
      fclose(f);
    }
    CAmpPhaseFrequencySeries* freqseries_new=0;     //New result will be assembled here
    CAmpPhaseFrequencySeries_Init(&freqseries_new,len_new);
    //set the new freqs
    for(i=0;i<len_new;i++){
      gsl_vector_set(freqseries_new->freq,i,gsl_vector_get(freq_new,i));
    }
    //copy in the high-freq ROM-model data
    //printf("l,m = %i,%i;  lenghts=%i,%i\n",l,m,freqseries->freq->size, freqseries_new->freq->size);
    for(i=len_add;i<len_new;i++){
      //printf("i, len-len_new+i: %i, %i\n",i,len-len_new+i);
      gsl_vector_set(freqseries_new->amp_real,i,gsl_vector_get(freqseries->amp_real,len-len_new+i));
      gsl_vector_set(freqseries_new->amp_imag,i,gsl_vector_get(freqseries->amp_imag,len-len_new+i));
      gsl_vector_set(freqseries_new->phase,i,gsl_vector_get(freqseries->phase,len-len_new+i));
      //printf("%i: copying %g  %g  %g  %g\n",i,freqseries_new->freq->data[i],freqseries_new->amp_real->data[i],freqseries_new->amp_imag->data[i],freqseries_new->phase->data[i]);
    }
    //extend
    if(l==2&&m==2&&len_add>0) {//extend 2-2 with TaylorF2
      //Assemble data for matching
      double f0=freq_new->data[len_add],f1=freq_new->data[len_add+1];
      double ph0=freqseries_new->phase->data[len_add],ph1=freqseries_new->phase->data[len_add+1];
      double amp=sqrt(freqseries_new->amp_real->data[len_add]*freqseries_new->amp_real->data[len_add]
        +freqseries_new->amp_imag->data[len_add]*freqseries_new->amp_imag->data[len_add]);
      double amp1=sqrt(freqseries_new->amp_real->data[len_add+1]*freqseries_new->amp_real->data[len_add+1]
        +freqseries_new->amp_imag->data[len_add+1]*freqseries_new->amp_imag->data[len_add+1]);
      double amprfac=freqseries_new->amp_real->data[len_add]/amp,ampifac=freqseries_new->amp_imag->data[len_add]/amp;
      //Compute raw TaylorF2
      TaylorF2nonspin(freqseries_new->amp_real->data,freqseries_new->phase->data,freq_new->data,len_add+2,m1SI,m2SI,distance,imatch);
      //Compute offsets in phase, first phase derivative, and amplitude argument
      double dphase0tf2=(freqseries_new->phase->data[len_add+1]-freqseries_new->phase->data[len_add])/(f1-f0);
      double phase0tf2=freqseries_new->phase->data[len_add];
      double dphase0eob=(ph1-ph0)/(f1-f0);
      double dphase0=dphase0eob-dphase0tf2;
      double phase0=ph0 - phase0tf2 - f0*dphase0;
      double amp0bcoeff = amp / freqseries_new->amp_real->data[len_add] - 1.0; //Compute correction for continuity matching.
      double amp0ccoeff = ((amp1/freqseries_new->amp_real->data[len_add+1]-1.0)/amp0bcoeff/(f1*f1/f0/f0)-1.0)/(f1/f0-1.0);
      //TESTING
      //printf("%g, %g\n", amp0bcoeff, amp0ccoeff);
      //Compute correction for continuity matching.
      /*
      printf("ph0eob,dph0eob= %g,  %g\n",ph0,dphase0eob);
      printf("ph0tf2,dph0tf2= %g,  %g\n",phase0tf2,dphase0tf2);
      printf("ph0,dph0= %g,  %g\n",phase0,dphase0);
      printf("f0,f0*dph0= %g,  %g\n",f0,dphase0*f0);
      printf("imatch=%i\n",imatch);
      */
      //Apply offsets
      for (i = 0; i < len_add+2; i++){
        double f=freqseries_new->freq->data[i];
        //printf("%i<%i,%g\n",i,len_add,len_add-i);
        //printf("f,ph0+f*dph0= %g, %g\n",freqseries_new->freq->data[i],phase0 + dphase0*f);
        freqseries_new->phase->data[i] += phase0 + dphase0*f;
        //First apply continuity matching
        // amp -> amp * ( 1 + b*f^2/f0^2 * ( 1 + c*(f/f0 -1) )
        //(starts at order f^2 since we only keep 2PN order ampl corrections in TaylorF2 code below; could change to f^3 if higher order terms are used)
        freqseries_new->amp_real->data[i] *= 1.0 + amp0bcoeff*f*f/f0/f0*(1+amp0ccoeff*(f/f0-1));
        freqseries_new->amp_imag->data[i] = freqseries_new->amp_real->data[i]*ampifac;
        freqseries_new->amp_real->data[i] *= amprfac;
        //printf("%i: extending TF2 %g  %g  %g  %g\n",i,freqseries_new->freq->data[i],freqseries_new->amp_real->data[i],freqseries_new->amp_imag->data[i],freqseries_new->phase->data[i]);
      }
    } else { //extend other modes with power-law
      //The results are many cycles out of phase almost immediately, so this definitely is not an accurate
      //waveform, but the results are reasonably smooth and of plausible structure.
      //Alternatively, we could also extend these with TaylorF2, btu we are mostly assuming this part of the WF is small
      double phref=freqseries->phase->data[len-1];//For phase we extend by a power-law referenced to zero phase at end of ROM
      for(i=0;i<len;i++)if(phref>freqseries->phase->data[i])phref=freqseries->phase->data[i];//get the smallest value of phi to use as ref.
      phref=phref+1.0;//add one more
      double ldArfac = 0;
      if(gsl_vector_get(freqseries->amp_real,imatch)>0) //avoid div0 in trivial cases
      //dArfac=exp(-log( gsl_vector_get(freqseries->amp_real,imatch+Navg)
      //		 /gsl_vector_get(freqseries->amp_real,imatch) ) / Navg);
      ldArfac=(log( gsl_vector_get(freqseries->amp_real,imatch+Navg)
      /gsl_vector_get(freqseries->amp_real,imatch) ) /
      log( gsl_vector_get(freqseries->freq,imatch+Navg)
      /gsl_vector_get(freqseries->freq,imatch) )  );
      //double dphfac = exp(-log( (gsl_vector_get(freqseries->phase,imatch+Navg)-phref)
      //			/(gsl_vector_get(freqseries->phase,imatch) -phref)) / Navg);
      double ldphfac = (log( (gsl_vector_get(freqseries->phase,imatch+Navg)-phref)
      /(gsl_vector_get(freqseries->phase,imatch) -phref)) /
      log( gsl_vector_get(freqseries->freq,imatch+Navg)
      /gsl_vector_get(freqseries->freq,imatch) )  );
      if(1&&l==lout&&m==mout)printf("ldphfac(%i,%i)=%g\n",l,m,ldphfac);
      //double f0=gsl_vector_get(freqseries->freq,imatch);
      for(i=len_add;i>0;i--){
        double fratio=gsl_vector_get(freqseries_new->freq,i-1)/gsl_vector_get(freqseries_new->freq,i);
        double dArfac=pow(fratio,ldArfac);
        gsl_vector_set(freqseries_new->amp_real,i-1,gsl_vector_get(freqseries_new->amp_real,i)*dArfac);
        gsl_vector_set(freqseries_new->amp_imag,i-1,gsl_vector_get(freqseries_new->amp_imag,i)*dArfac);
        double dphfac=pow(fratio,ldphfac);
        gsl_vector_set(freqseries_new->phase,i-1,(gsl_vector_get(freqseries_new->phase,i)-phref)*dphfac+phref);
        //printf("%i: extending %g  %g  %g  %g\n",i,freqseries_new->freq->data[i-1],freqseries_new->amp_real->data[i-1],freqseries_new->amp_imag->data[i-1],freqseries_new->phase->data[i-1]);
      }
      //printf("Extended (%i,%i) down to f=%g, ampR=%g, ampI=%g, phase=%g\n",l,m,freqseries_new->freq->data[0],freqseries_new->amp_real->data[0],freqseries_new->amp_imag->data[0],freqseries_new->phase->data[0]);
    }
    //delete the old content data and replace with the new
    CAmpPhaseFrequencySeries_Cleanup(freqseries);
    listelement->freqseries=freqseries_new;

    //TESTING
    //printf("In ext: len freqseries: %d\n", freqseries_new->freq->size);
    //for(int i=0; i<freqseries_new->freq->size; i++) {
    //printf("%g %g %g %g\n", gsl_vector_get(freqseries_new->freq, i), gsl_vector_get(freqseries_new->amp_real, i), gsl_vector_get(freqseries_new->amp_imag, i), gsl_vector_get(freqseries_new->phase, i));
    //}

    //debugging
    //   freqseries=listelement->freqseries;
    //   if(1&&l==lout&&m==mout){  //write to file for debugging
    //     FILE *f = fopen("waveform.dat", "w");
    //     for(i=0;i<freqseries->freq->size;i++){
    // fprintf(f,"%i %i %g  %g %g  %g %i\n",l,m,
    // 	freqseries->freq->data[i],
    // 	freqseries->amp_real->data[i],
    // 	freqseries->amp_imag->data[i],
    // 	freqseries->phase->data[i],
    // 	i);
    //     }
    //     fclose(f);
    //   }

    gsl_vector_free(freq_new);
    listelement=listelement->next;
  }

  /* Compute phase shift to set phiRef at fRef */
  /* Covers the case where input fRef is outside the range of the ROM, in which case the Core function defaulted to Mfmax_ROM */
  /* Not very clean and a bit redundant */
  if (setphiRefatfRef) {
    CAmpPhaseFrequencySeries* h22 = ListmodesCAmpPhaseFrequencySeries_GetMode(listROM, 2, 2)->freqseries;
    gsl_vector* freq22 = h22->freq;
    gsl_vector* phase22 = h22->phase;
    int nbfreq = freq22->size;
    gsl_interp_accel* accel_phi22 = gsl_interp_accel_alloc();
    gsl_spline* spline_phi22 = gsl_spline_alloc(gsl_interp_cspline, nbfreq);
    gsl_spline_init(spline_phi22, gsl_vector_const_ptr(freq22,0), gsl_vector_const_ptr(phase22,0), nbfreq);
    /* If fRef was not set (fRef=0), use the last frequency generated for 22 mode -- as is done internally in Core */
    if (fRef==0.) fRef = freq22->data[freq22->size - 1];
    /* Compute 22 phase at fRef before adjustment, check the extended range */
    if ( (fRef<freq22->data[0]) || (fRef>freq22->data[freq22->size - 1]) ) {
      printf("Error: fRef is not covered by the frequency range of the waveform after extension.\n");
      return FAILURE;
    }
    double phi22atfRef = gsl_spline_eval(spline_phi22, fRef, accel_phi22);
    /* Phase shift, as an orbital or observer phase shift (change in 22 is 2*phaseshift) */
    double phaseshift = (2*phiRef - phi22atfRef)/2.;

    /* Propagate phase shift to full list of modes */
    listelement = listROM;
    while(listelement) {
      int l = listelement->l;
      int m = listelement->m;

      double phaseshiftlm = m/2. * phaseshift;
      gsl_vector* philm = listelement->freqseries->phase;
      gsl_vector_add_constant(philm, phaseshiftlm);

      listelement = listelement->next;
    }

    /* Cleanup */
    gsl_spline_free(spline_phi22);
    gsl_interp_accel_free(accel_phi22);
  }

  /* Output */
  *listhlm = listROM;

      /*
      printf("generated listROM: n=%i l=%i m=%i\n",(*listhlm)->freqseries->amp_real->size,listROM->l,listROM->m);
      for(i=0;i<(*listhlm)->freqseries->freq->size;i++){
      printf("%i:  %g  %g  %g  %g\n",i,(*listhlm)->freqseries->freq->data[i],(*listhlm)->freqseries->amp_real->data[i],(*listhlm)->freqseries->amp_imag->data[i],(*listhlm)->freqseries->phase->data[i]);
    }
    printf("listlhm=%x\n",listhlm);
    printf("*listlhm=%x\n",*listhlm);
    */
  return SUCCESS;
}

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare