text:
#pragma once

#include <cstddef>
#include <gsl/gsl>
#include <iterator>
#include <limits>

#include "../math/General.h"
#include "Concepts.h"

namespace utils {
	using concepts::DefaultConstructible, concepts::Integral, concepts::Copyable, concepts::Movable,
		concepts::NotMovable, concepts::ConstructibleFrom;
	/// @brief A simple Ring Buffer implementation.
	/// Supports resizing, writing, reading, erasing, and provides mutable and immutable
	/// random access iterators.
	///
	/// # Iterator Invalidation
	/// * Iterators are lazily evaluated, so will only ever be invalidated at their current state.
	/// Performing any mutating operation (mutating the iterator, not the underlying data) on them
	/// will re-sync them with their associated `RingBuffer`.
	/// The following operations will invalidate an iterator's current state:
	/// - Read-only operations: never
	/// - clear: always
	/// - reserve: only if the `RingBuffer` changed capacity
	/// - erase: Erased elements and all following elements
	/// - push_back, emplace_back: only `end()`
	/// - insert, emplace: only the element at the position inserted/emplaced
	/// - pop_back: the element removed and `end()`
	///
	/// @tparam T - The type to store in the `RingBuffer`; Must be Default Constructible
	template<DefaultConstructible T>
	class RingBuffer {
	  public:
		/// Default capacity of `RingBuffer`
		static const constexpr size_t DEFAULT_CAPACITY = 16;

		/// @brief Random-Access Bidirectional iterator for `RingBuffer`
		/// @note All navigation operators are checked such that any movement past `begin()` or
		/// `end()` is ignored.
		class Iterator {
		  public:
			using iterator_category = std::random_access_iterator_tag;
			using difference_type = std::ptrdiff_t;
			using value_type = T;
			using pointer = value_type*;
			using reference = value_type&;

			constexpr explicit Iterator(pointer ptr,
										RingBuffer* containerPtr,
										size_t currentIndex) noexcept
				: mPtr(ptr), mContainerPtr(containerPtr), mCurrentIndex(currentIndex) {
			}
			constexpr Iterator(const Iterator& iter) noexcept = default;
			constexpr Iterator(Iterator&& iter) noexcept = default;
			~Iterator() noexcept = default;

			/// @brief Returns the index in the `RingBuffer` that corresponds
			/// to the element this iterator points to
			///
			/// @return The index corresponding with the element this points to
			[[nodiscard]] constexpr inline auto getIndex() const noexcept -> size_t {
				return mCurrentIndex;
			}

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

			constexpr inline auto operator==(const Iterator& rhs) const noexcept -> bool {
				return mPtr == rhs.mPtr;
			}

			constexpr inline auto operator!=(const Iterator& rhs) const noexcept -> bool {
				return mPtr != rhs.mPtr;
			}

			constexpr inline auto operator*() const noexcept -> reference {
				return *mPtr;
			}

			constexpr inline auto operator->() noexcept -> pointer {
				return mPtr;
			}

			constexpr inline auto operator++() noexcept -> Iterator& {
				mCurrentIndex++;
				if(mCurrentIndex >= mContainerPtr->capacity()) {
					mCurrentIndex = mContainerPtr->capacity();
					mPtr = mContainerPtr->end().mPtr;
				}
				else {
					mPtr = &(*mContainerPtr)[mCurrentIndex];
				}
				return *this;
			}

			constexpr inline auto operator++(int) noexcept -> Iterator {
				Iterator temp = *this;
				++(*this);
				return temp;
			}

			constexpr inline auto operator--() noexcept -> Iterator& {
				if(mCurrentIndex == 0) {
					return *this;
				}
				else {
					mCurrentIndex--;
					mPtr = &(*mContainerPtr)[mCurrentIndex];
				}
				return *this;
			}

			constexpr inline auto operator--(int) noexcept -> Iterator {
				Iterator temp = *this;
				--(*this);
				return temp;
			}

			constexpr inline auto operator+(Integral auto rhs) const noexcept -> Iterator {
				const auto diff = static_cast<size_t>(rhs);
				if(rhs < 0) {
					return std::move(*this - -rhs);
				}

				auto temp = *this;
				temp.mCurrentIndex += diff;
				if(temp.mCurrentIndex > temp.mContainerPtr->capacity()) {
					temp.mCurrentIndex = temp.mContainerPtr->capacity();
					temp.mPtr = temp.mContainerPtr->end().mPtr;
				}
				else {
					temp.mPtr = &(*temp.mContainerPtr)[temp.mCurrentIndex];
				}
				return temp;
			}

			constexpr inline auto operator+=(Integral auto rhs) noexcept -> Iterator& {
				*this = std::move(*this + rhs);
				return *this;
			}

			constexpr inline auto operator-(Integral auto rhs) const noexcept -> Iterator {
				const auto diff = static_cast<size_t>(rhs);
				if(rhs < 0) {
					return std::move(*this + -rhs);
				}

				auto temp = *this;
				if(diff > temp.mCurrentIndex) {
					temp.mPtr = temp.mContainerPtr->begin().mPtr;
					temp.mCurrentIndex = 0;
				}
				else {
					temp.mCurrentIndex -= diff;
					temp.mPtr = &(*temp.mContainerPtr)[temp.mCurrentIndex];
				}
				return temp;
			}

			constexpr inline auto operator-=(Integral auto rhs) noexcept -> Iterator& {
				*this = std::move(*this - rhs);
				return *this;
			}

			constexpr inline auto operator-(const Iterator& rhs) const noexcept -> difference_type {
				return mPtr - rhs.mPtr;
			}

			constexpr inline auto operator[](Integral auto index) noexcept -> Iterator {
				return std::move(*this + index);
			}

			constexpr inline auto operator>(const Iterator& rhs) const noexcept -> bool {
				return mCurrentIndex > rhs.mCurrentIndex;
			}

			constexpr inline auto operator<(const Iterator& rhs) const noexcept -> bool {
				return mCurrentIndex < rhs.mCurrentIndex;
			}

			constexpr inline auto operator>=(const Iterator& rhs) const noexcept -> bool {
				return mCurrentIndex >= rhs.mCurrentIndex;
			}

			constexpr inline auto operator<=(const Iterator& rhs) const noexcept -> bool {
				return mCurrentIndex <= rhs.mCurrentIndex;
			}

		  private:
			pointer mPtr;
			RingBuffer* mContainerPtr = nullptr;
			size_t mCurrentIndex = 0;
		};

		/// @brief Read-only Random-Access Bidirectional iterator for `RingBuffer`
		/// @note All navigation operators are checked such that any movement past `begin()` or
		/// `end()` is ignored.
		class ConstIterator {
		  public:
			using iterator_category = std::random_access_iterator_tag;
			using difference_type = std::ptrdiff_t;
			using value_type = T;
			using pointer = const value_type*;
			using reference = const value_type&;

			constexpr explicit ConstIterator(pointer ptr,
											 RingBuffer* containerPtr,
											 size_t currentIndex) noexcept
				: mPtr(ptr), mContainerPtr(containerPtr), mCurrentIndex(currentIndex) {
			}
			constexpr ConstIterator(const ConstIterator& iter) noexcept = default;
			constexpr ConstIterator(ConstIterator&& iter) noexcept = default;
			~ConstIterator() noexcept = default;

			/// @brief Returns the index in the `RingBuffer` that corresponds
			/// to the element this iterator points to
			///
			/// @return The index corresponding with the element this points to
			[[nodiscard]] constexpr inline auto getIndex() const noexcept -> size_t {
				return mCurrentIndex;
			}

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

			constexpr inline auto operator==(const ConstIterator& rhs) const noexcept -> bool {
				return mPtr == rhs.mPtr;
			}

			constexpr inline auto operator!=(const ConstIterator& rhs) const noexcept -> bool {
				return mPtr != rhs.mPtr;
			}

			constexpr inline auto operator*() const noexcept -> reference {
				return *mPtr;
			}

			constexpr inline auto operator->() const noexcept -> pointer {
				return mPtr;
			}

			constexpr inline auto operator++() noexcept -> ConstIterator& {
				mCurrentIndex++;
				if(mCurrentIndex >= mContainerPtr->capacity()) {
					mCurrentIndex = mContainerPtr->capacity();
					mPtr = mContainerPtr->end().mPtr;
				}
				else {
					mPtr = &(*mContainerPtr)[mCurrentIndex];
				}
				return *this;
			}

			constexpr inline auto operator++(int) noexcept -> ConstIterator {
				ConstIterator temp = *this;
				++(*this);
				return temp;
			}

			constexpr inline auto operator--() noexcept -> ConstIterator& {
				if(mCurrentIndex == 0) {
					return *this;
				}
				else {
					mCurrentIndex--;
					mPtr = &(*mContainerPtr)[mCurrentIndex];
				}
				return *this;
			}

			constexpr inline auto operator--(int) noexcept -> ConstIterator {
				ConstIterator temp = *this;
				--(*this);
				return temp;
			}

			constexpr inline auto operator+(Integral auto rhs) const noexcept -> ConstIterator {
				const auto diff = static_cast<size_t>(rhs);
				if(rhs < 0) {
					return std::move(*this - -rhs);
				}

				auto temp = *this;
				temp.mCurrentIndex += diff;
				if(temp.mCurrentIndex > temp.mContainerPtr->capacity()) {
					temp.mCurrentIndex = temp.mContainerPtr->capacity();
					temp.mPtr = temp.mContainerPtr->end().mPtr;
				}
				else {
					temp.mPtr = &(*temp.mContainerPtr)[temp.mCurrentIndex];
				}
				return temp;
			}

			constexpr inline auto operator+=(Integral auto rhs) noexcept -> ConstIterator& {
				*this = std::move(*this + rhs);
				return *this;
			}

			constexpr inline auto operator-(Integral auto rhs) const noexcept -> ConstIterator {
				const auto diff = static_cast<size_t>(rhs);
				if(rhs < 0) {
					return std::move(*this + -rhs);
				}

				auto temp = *this;
				if(diff > temp.mCurrentIndex) {
					temp.mPtr = temp.mContainerPtr->begin().mPtr;
					temp.mCurrentIndex = 0;
				}
				else {
					temp.mCurrentIndex -= diff;
					temp.mPtr = &(*temp.mContainerPtr)[temp.mCurrentIndex];
				}
				return temp;
			}

			constexpr inline auto operator-=(Integral auto rhs) noexcept -> ConstIterator& {
				*this = std::move(*this - rhs);
				return *this;
			}

			constexpr inline auto
			operator-(const ConstIterator& rhs) const noexcept -> difference_type {
				return mPtr - rhs.mPtr;
			}

			constexpr inline auto operator[](Integral auto index) const noexcept -> ConstIterator {
				return std::move(*this + index);
			}

			constexpr inline auto operator>(const ConstIterator& rhs) const noexcept -> bool {
				return mCurrentIndex > rhs.mCurrentIndex;
			}

			constexpr inline auto operator<(const ConstIterator& rhs) const noexcept -> bool {
				return mCurrentIndex < rhs.mCurrentIndex;
			}

			constexpr inline auto operator>=(const ConstIterator& rhs) const noexcept -> bool {
				return mCurrentIndex >= rhs.mCurrentIndex;
			}

			constexpr inline auto operator<=(const ConstIterator& rhs) const noexcept -> bool {
				return mCurrentIndex <= rhs.mCurrentIndex;
			}

		  private:
			pointer mPtr;
			RingBuffer* mContainerPtr = nullptr;
			size_t mCurrentIndex = 0;
		};

		/// @brief Creates a `RingBuffer` with default capacity
		RingBuffer() noexcept : mBuffer(new T[DEFAULT_CAPACITY_INTERNAL]) {
		}

		/// @brief Creates a `RingBuffer` with (at least) the given initial capacity
		///
		/// @param initialCapacity - The initial capacity of the `RingBuffer`
		constexpr explicit RingBuffer(size_t initialCapacity) noexcept
			: mBuffer(new T[initialCapacity + 1]), mLoopIndex(initialCapacity),
			  mCapacity(initialCapacity) {
		}

		/// @brief Constructs a new `RingBuffer` with the given initial capacity and
		/// fills it with `defaultValue`
		///
		/// @param initialCapacity - The initial capacity of the `RingBuffer`
		/// @param defaultValue - The value to fill the `RingBuffer` with
		constexpr RingBuffer(size_t initialCapacity,
							 const T& defaultValue) noexcept requires Copyable<T>
			: mBuffer(new T[initialCapacity + 1]),
			  mLoopIndex(initialCapacity),
			  mCapacity(initialCapacity) {
			for(auto i = 0ULL; i < mCapacity; ++i) {
				push_back(defaultValue);
			}
		}

		constexpr RingBuffer(const RingBuffer& buffer) noexcept requires Copyable<T>
			: mBuffer(new T[buffer.mCapacity + 1]),
			  mWriteIndex(0ULL), // NOLINT
			  mStartIndex(0ULL), // NOLINT
			  mLoopIndex(buffer.mLoopIndex),
			  mCapacity(buffer.mCapacity) {
			for(auto i = 0ULL; i < mCapacity; ++i) {
				push_back(buffer.mBuffer[i]);
			}
			mStartIndex = buffer.mStartIndex; // NOLINT(cppcoreguidelines-prefer-member-initializer)
			mWriteIndex = buffer.mWriteIndex; // NOLINT(cppcoreguidelines-prefer-member-initializer)
			mSize = buffer.mSize;			  // NOLINT(cppcoreguidelines-prefer-member-initializer)
		}

		constexpr RingBuffer(RingBuffer&& buffer) noexcept
			: mBuffer(std::move(buffer.mBuffer)), mStartIndex(buffer.mStartIndex),
			  mLoopIndex(buffer.mLoopIndex), mCapacity(buffer.mCapacity), mSize(buffer.mSize) {
			buffer.mCapacity = 0ULL;
			buffer.mLoopIndex = 0ULL;
			buffer.mSize = 0ULL;
			buffer.mWriteIndex = 0ULL;
			buffer.mStartIndex = 0ULL;
			buffer.mBuffer = nullptr;
		}

		~RingBuffer() noexcept = default;

		/// @brief Returns the element at the given index.
		/// @note This is not checked in the same manner as STL containers:
		/// if index >= capacity, the element at capacity - 1 is returned.
		///
		/// @param index - The index of the desired element
		///
		/// @return The element at the given index, or at capacity - 1 if index >= capacity
		[[nodiscard]] constexpr inline auto at(Integral auto index) noexcept -> T& {
			auto i = getAdjustedInternalIndex(index);

			if(i > mLoopIndex) {
				i = mLoopIndex;
			}

			return mBuffer[i]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
		}

		/// @brief Returns the first element in the `RingBuffer`
		///
		/// @return The first element
		[[nodiscard]] constexpr inline auto front() noexcept -> T& {
			return mBuffer[mStartIndex]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
		}

		/// @brief Returns the last element in the `RingBuffer`
		/// @note If <= 1 elements are in the `RingBuffer`, this will be the same as `front`
		///
		/// @return The last element
		[[nodiscard]] constexpr inline auto back() noexcept -> T& {
			auto index = mWriteIndex - 1;
			if(mWriteIndex == 0) {
				if(mStartIndex == 0) {
					index = 0;
				}
				else {
					index = mLoopIndex;
				}
			}
			return mBuffer[index]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
		}

		/// @brief Returns a pointer to the underlying data in the `RingBuffer`.
		/// @note This is not sorted in any way to match the representation used by the `RingBuffer`
		///
		/// @return A pointer to the underlying data
		[[nodiscard]] constexpr inline auto data() noexcept -> T* {
			return mBuffer;
		}

		/// @brief Returns whether the `RingBuffer` is empty
		///
		/// @return `true` if the `RingBuffer` is empty, `false` otherwise
		[[nodiscard]] constexpr inline auto empty() noexcept -> bool {
			return mSize == 0;
		}

		/// @brief Returns the current number of elements in the `RingBuffer`
		///
		/// @return The current number of elements
		[[nodiscard]] constexpr inline auto size() noexcept -> size_t {
			return mSize;
		}

		/// @brief Returns the maximum possible number of elements this `RingBuffer` could store
		/// if grown to maximum possible capacity
		///
		/// @return The maximum possible number of storable elements
		[[nodiscard]] constexpr inline auto max_size() noexcept -> size_t {
			return std::numeric_limits<size_t>::max();
		}

		/// @brief Returns the current capacity of the `RingBuffer`;
		/// the number of elements it can currently store
		///
		/// @return The current capacity
		[[nodiscard]] constexpr inline auto capacity() noexcept -> size_t {
			return mCapacity;
		}

		/// @brief Reserves more storage for the `RingBuffer`. If `newCapacity` is > capacity,
		/// then the capacity of the `RingBuffer` will be extended until at least `newCapacity`
		/// elements can be stored.
		/// @note Memory contiguity is maintained, so no **elements** will be lost or invalidated.
		/// However, all iterators and references to elements will be invalidated.
		///
		/// @param newCapacity - The new capacity of the `RingBuffer`
		constexpr inline auto reserve(size_t newCapacity) noexcept -> void {
			// we only need to do anything if `newCapacity` is actually larger than `mCapacity`
			if(newCapacity > mCapacity) {
				gsl::owner<T*> temp = new T[newCapacity]; // NOLINT
				auto span = gsl::make_span(temp, newCapacity);
				std::copy(begin(), end(), span.begin());
				mBuffer.reset(temp);
				mStartIndex = 0;
				mWriteIndex = mLoopIndex + 1;
				mLoopIndex = newCapacity;
				mCapacity = newCapacity;
			}
		}

		/// @brief Erases all elements from the `RingBuffer`
		constexpr inline auto clear() noexcept -> void {
			mStartIndex = 0;
			mWriteIndex = 0;
			mSize = 0;
		}

		/// @brief Inserts the given element at the end of the `RingBuffer`
		/// @note if `size() == capacity()` then this loops and overwrites `front()`
		///
		/// @param value - the element to insert
		constexpr inline auto push_back(const T& value) noexcept -> void requires Copyable<T> {
			// clang-format off
		mBuffer[mWriteIndex] = value; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
			// clang-format on

			incrementIndices();
		}

		/// @brief Inserts the given element at the end of the `RingBuffer`
		/// @note if `size() == capacity()` then this loops and overwrites `front()`
		///
		/// @param value - the element to insert
		constexpr inline auto push_back(T&& value) noexcept -> void {
			// clang-format off
			mBuffer[mWriteIndex] = std::forward<T>(value); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
			// clang-format on

			incrementIndices();
		}

		/// @brief Constructs the given element in place at the end of the `RingBuffer`
		/// @note if `size() == capacity()` then this loops and overwrites `front()`
		///
		/// @tparam Args - The types of the element's constructor arguments
		/// @param args - The constructor arguments for the element
		///
		/// @return A reference to the element constructed at the end of the `RingBuffer`
		template<typename... Args>
		requires ConstructibleFrom<T, Args...>
		constexpr inline auto emplace_back(Args&&... args) noexcept -> T& requires Movable<T> {
			// clang-format off
			new (&mBuffer[mWriteIndex]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
			// clang-format on

			incrementIndices();

			auto index = mWriteIndex == 0 ? mLoopIndex : mWriteIndex - 1;
			return mBuffer[index]; // NOLINT
		}

		/// @brief Constructs the given element in place at the location
		/// indicated by the `Iterator` `position`
		///
		/// @tparam Args - The types of the element's constructor arguments
		/// @param position - `Iterator` indicating where in the `RingBuffer` to construct the
		/// element
		/// @param args - The constructor arguments for the element
		///
		/// @return A reference to the element constructed at the location indicated by `position`
		template<typename... Args>
		requires ConstructibleFrom<T, Args...>
		constexpr inline auto emplace(const Iterator& position, Args&&... args) noexcept -> T& {
			auto index = getAdjustedInternalIndex(position.getIndex());

			// clang-format off
			new (&mBuffer[index]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
			// clang-format on

			return mBuffer[index]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
		}

		/// @brief Constructs the given element in place at the location
		/// indicated by the `ConstIterator` `position`
		///
		/// @tparam Args - The types of the element's constructor arguments
		/// @param position - `ConstIterator` indicating where in the `RingBuffer` to construct the
		/// element
		/// @param args - The constructor arguments for the element
		///
		/// @return A reference to the element constructed at the location indicated by `position`
		template<typename... Args>
		constexpr inline auto
		emplace(const ConstIterator& position, Args&&... args) noexcept -> T& {
			auto index = getAdjustedInternalIndex(position.getIndex());

			// clang-format off
			new (&mBuffer[index]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
			// clang-format on

			return mBuffer[index]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
		}

		/// @brief Assigns the given element to the position indicated
		/// by the `Iterator` `position`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param position - `Iterator` indicating where in the `RingBuffer` to place the element
		/// @param element - The element to store in the `RingBuffer`
		constexpr inline auto
		insert(const Iterator& position, const T& element) noexcept -> void requires Copyable<T> {
			insert_internal(position.getIndex(), element);
		}

		/// @brief Assigns the given element to the position indicated
		/// by the `Iterator` `position`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param position - `Iterator` indicating where in the `RingBuffer` to place the element
		/// @param element - The element to store in the `RingBuffer`
		constexpr inline auto insert(const Iterator& position, T&& element) noexcept -> void {
			insert_internal(position.getIndex(), std::forward<T>(element));
		}

		/// @brief Assigns the given element to the position indicated
		/// by the `ConstIterator` `position`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param position - `ConstIterator` indicating where in the `RingBuffer` to place the
		/// element
		/// @param element - The element to store in the `RingBuffer`
		constexpr inline auto insert(const ConstIterator& position, const T& element) noexcept
			-> void requires Copyable<T> {
			insert_internal(position.getIndex(), element);
		}

		/// @brief Assigns the given element to the position indicated
		/// by the `ConstIterator` `position`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param position - `ConstIterator` indicating where in the `RingBuffer` to place the
		/// element
		/// @param element - The element to store in the `RingBuffer`
		constexpr inline auto insert(const ConstIterator& position, T&& element) noexcept -> void {
			insert_internal(position.getIndex(), std::forward<T>(element));
		}

		/// @brief Constructs the given element at the insertion position indicated
		/// by the `ConstIterator` `position`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param position - `ConstIterator` indicating where in the `RingBuffer` to place the
		/// element
		/// @param args - The arguments to the constructor for
		/// the element to store in the `RingBuffer`
		template<typename... Args>
		requires ConstructibleFrom<T, Args...>
		constexpr inline auto
		insert_emplace(const Iterator& position, Args&&... args) noexcept -> T& {
			return insert_emplace_internal(position.getIndex(), args...);
		}

		/// @brief Constructs the given element at the insertion position indicated
		/// by the `ConstIterator` `position`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param position - `ConstIterator` indicating where in the `RingBuffer` to place the
		/// element
		/// @param args - The arguments to the constructor for
		/// the element to store in the `RingBuffer`
		template<typename... Args>
		requires ConstructibleFrom<T, Args...>
		constexpr inline auto
		insert_emplace(const ConstIterator& position, Args&&... args) noexcept -> T& {
			return insert_emplace_internal(position.getIndex(), args...);
		}

		/// @brief Erases the element at the given `position`, moving other elements backward
		/// in the buffer to maintain contiguity
		///
		/// @param position - The element to erase
		///
		/// @return `Iterator` pointing to the element after the one erased
		constexpr inline auto erase(const Iterator& position) noexcept -> Iterator {
			return erase_internal(position.getIndex());
		}

		/// @brief Erases the element at the given `position`, moving other elements backward
		/// in the buffer to maintain contiguity
		///
		/// @param position - The element to erase
		///
		/// @return `Iterator` pointing to the element after the one erased
		constexpr inline auto erase(const ConstIterator& position) noexcept -> Iterator {
			return erase_internal(position.getIndex());
		}

		/// @brief Erases the range of elements in [`first`, `last`)
		/// Returns an `Iterator` to the element after the last one erased
		/// @note In the case `first` >= `last`, no elements are erased and `last` is returned;
		///
		/// @param first - The first element in the range to erase
		/// @param last - The last element in the range
		///
		/// @return `Iterator` pointing to the element after the last one erased
		constexpr inline auto
		erase(const Iterator& first, const Iterator& last) noexcept -> Iterator {
			if(first >= last) {
				return last;
			}

			return erase_internal(first.getIndex(), last.getIndex());
		}

		/// @brief Erases the range of elements in [`first`, `last`)
		/// Returns an `Iterator` to the element after the last one erased
		/// @note In the case `first` >= `last`, no elements are erased and `last` is returned;
		///
		/// @param first - The first element in the range to erase
		/// @param last - The last element in the range
		///
		/// @return `Iterator` pointing to the element after the last one erased
		constexpr inline auto
		erase(const ConstIterator& first, const ConstIterator& last) noexcept -> Iterator {
			if(first >= last) {
				return last;
			}

			return erase_internal(first.getIndex(), last.getIndex());
		}

		/// @brief Removes the last element in the `RingBuffer` and returns it
		///
		/// @return The last element in the `RingBuffer`
		[[nodiscard]] constexpr inline auto pop_back() noexcept -> T requires Copyable<T> {
			T _back = back();
			erase(--end());
			return _back;
		}

		/// @brief Returns a Random Access Bidirectional iterator over the `RingBuffer`,
		/// at the beginning
		///
		/// @return The iterator, at the beginning
		[[nodiscard]] constexpr inline auto begin() -> Iterator {
			T* p = &mBuffer[mStartIndex]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)

			return Iterator(p, this, 0ULL);
		}

		/// @brief Returns a Random Access Bidirectional iterator over the `RingBuffer`,
		/// at the end
		///
		/// @return The iterator, at the end
		[[nodiscard]] constexpr inline auto end() -> Iterator {
			T* p = &mBuffer[mWriteIndex]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)

			return Iterator(p, this, mSize);
		}

		/// @brief Returns a Random Access Bidirectional read-only iterator over the `RingBuffer`,
		/// at the beginning
		///
		/// @return The iterator, at the beginning
		[[nodiscard]] constexpr inline auto cbegin() -> ConstIterator {
			T* p = &mBuffer[mStartIndex]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)

			return ConstIterator(p, this, 0ULL);
		}

		/// @brief Returns a Random Access Bidirectional read-only iterator over the `RingBuffer`,
		/// at the end
		///
		/// @return The iterator, at the end
		[[nodiscard]] constexpr inline auto cend() -> ConstIterator {
			T* p = &mBuffer[mWriteIndex]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
			return ConstIterator(p, this, mSize);
		}

		/// @brief Unchecked access-by-index operator
		///
		/// @param index - The index to get the corresponding element for
		///
		/// @return - The element at index
		[[nodiscard]] constexpr inline auto operator[](Integral auto index) noexcept -> T& {
			auto i = getAdjustedInternalIndex(index);

			return mBuffer[i]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
		}

		constexpr auto
		operator=(const RingBuffer& buffer) noexcept -> RingBuffer& requires Copyable<T> {
			if(this == &buffer) {
				return *this;
			}
			mBuffer.reset(new T[buffer.mCapacity + 1]); // NOLINT
			mStartIndex = 0ULL;
			mWriteIndex = 0ULL;
			mCapacity = buffer.mCapacity;
			for(auto i = 0ULL; i < mCapacity; ++i) {
				push_back(buffer.mBuffer[i]);
			}
			mStartIndex = buffer.mStartIndex;
			mWriteIndex = buffer.mWriteIndex;
			mLoopIndex = buffer.mLoopIndex;
			mSize = buffer.mSize;
			return *this;
		}
		constexpr auto operator=(RingBuffer&& buffer) noexcept -> RingBuffer& {
			mBuffer = std::move(buffer.mBuffer);
			mWriteIndex = buffer.mWriteIndex;
			mStartIndex = buffer.mStartIndex;
			mLoopIndex = buffer.mLoopIndex;
			mCapacity = buffer.mCapacity;
			mSize = buffer.mSize;
			buffer.mBuffer = nullptr;
			buffer.mWriteIndex = 0ULL;
			buffer.mStartIndex = 0ULL;
			buffer.mCapacity = 0ULL;
			buffer.mSize = 0ULL;
			return *this;
		}

	  private:
		static const constexpr size_t DEFAULT_CAPACITY_INTERNAL = DEFAULT_CAPACITY + 1;
		std::unique_ptr<T[]> mBuffer = new T[DEFAULT_CAPACITY_INTERNAL]; // NOLINT
		size_t mWriteIndex = 0;
		size_t mStartIndex = 0;
		size_t mLoopIndex = DEFAULT_CAPACITY;
		size_t mCapacity = DEFAULT_CAPACITY;
		size_t mSize = 0;

		/// @brief Converts the given `RingBuffer` index into the corresponding index into then
		/// underlying `T` array
		///
		/// @param index - The `RingBuffer` index to convert
		///
		/// @return The corresponding index into the underlying `T` array
		[[nodiscard]] constexpr inline auto
		getAdjustedInternalIndex(Integral auto index) const noexcept -> size_t {
			auto i = static_cast<size_t>(index);
			if(mStartIndex + i > mLoopIndex) {
				i = (mStartIndex + i) - (mLoopIndex + 1);
			}
			else {
				i += mStartIndex;
			}

			return i;
		}

		/// @brief Converts the given index into the underlying `T` array into
		/// a using facing index into the `RingBuffer`
		///
		/// @param index - The internal index
		///
		/// @return The corresponding user-facing index
		[[nodiscard]] constexpr inline auto
		getExternalIndexFromInternal(Integral auto index) const noexcept -> size_t {
			auto i = static_cast<size_t>(index);
			if(i > mStartIndex && i <= mLoopIndex) {
				return i - mStartIndex;
			}
			else if(i < mStartIndex) {
				return (mLoopIndex - mStartIndex) + i + 1;
			}
			else if(i == mStartIndex) {
				return 0;
			}
			else {
				return mSize;
			}
		}

		/// @brief Increments the start and write indices into the underlying `T` array,
		/// and the size property, maintaining the logical `RingBuffer` structure
		constexpr inline auto incrementIndices() noexcept -> void {
			mWriteIndex++;
			mSize = math::General::min(mSize + 1, mCapacity);

			// if write index is at start - 1, we need to push start forward to maintain
			// the "invalid" spacer element for this.end()
			if(mWriteIndex > mLoopIndex && mStartIndex == 0) {
				mWriteIndex = 0;
				mStartIndex = 1;
			}
			else if(mWriteIndex == mStartIndex) {
				mStartIndex++;
				if(mStartIndex > mLoopIndex) {
					mStartIndex = 0;
				}
			}
		}

		/// @brief Inserts the given element at the position indicated
		/// by the `externalIndex`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param externalIndex - The user-facing index into the `RingBuffer` to insert the element
		/// at
		/// @param elem - The element to store in the `RingBuffer`
		constexpr inline auto
		insert_internal(size_t externalIndex, const T& elem) noexcept -> void requires Movable<T> {
			const auto index = getAdjustedInternalIndex(externalIndex);

			if(index == mWriteIndex) {
				emplace_back(elem);
			}
			else {
				auto numToMove = mSize - externalIndex;
				auto iter = begin() + externalIndex;
				if(mSize == mCapacity) [[likely]] { // NOLINT
					for(size_t i = 1ULL, j = numToMove - 2; i < numToMove; ++i, --j) {
						*((end() - i) + 1) = std::move(*(iter + j));
					}
				}
				else {
					for(size_t i = 0ULL, j = numToMove - 1; i < numToMove; ++i, --j) {
						*(end() - i) = std::move(*(iter + j));
					}
				}
				incrementIndices();
				if(externalIndex == 0) {
					// clang-format off
					mBuffer[mStartIndex] = elem; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
				else {
					// clang-format off
					mBuffer[index] = elem; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
			}
		}

		/// @brief Inserts the given element at the position indicated
		/// by the `externalIndex`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param externalIndex - The user-facing index into the `RingBuffer` to insert the element
		/// at
		/// @param elem - The element to store in the `RingBuffer`
		constexpr inline auto insert_internal(size_t externalIndex, const T& elem) noexcept
			-> void requires NotMovable<T> {
			const auto index = getAdjustedInternalIndex(externalIndex);

			if(index == mWriteIndex) {
				emplace_back(elem);
			}
			else {
				auto numToMove = mSize - externalIndex;
				auto iter = begin() + externalIndex;
				if(mSize == mCapacity) [[likely]] { // NOLINT
					for(size_t i = 1ULL, j = numToMove - 2; i < numToMove; ++i, --j) {
						*((end() - i) + 1) = *(iter + j);
					}
				}
				else {
					for(size_t i = 0ULL, j = numToMove - 1; i < numToMove; ++i, --j) {
						*(end() - i) = *(iter + j);
					}
				}
				incrementIndices();
				if(externalIndex == 0) {
					// clang-format off
					mBuffer[mStartIndex] = elem; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
				else {
					// clang-format off
					mBuffer[index] = elem; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
			}
		}

		/// @brief Inserts the given element at the position indicated
		/// by the `externalIndex`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param externalIndex - The user-facing index into the `RingBuffer` to insert the element
		/// at
		/// @param elem - The element to store in the `RingBuffer`
		constexpr inline auto
		insert_internal(size_t externalIndex, T&& elem) noexcept -> void requires Movable<T> {
			const auto index = getAdjustedInternalIndex(externalIndex);

			if(index == mWriteIndex) {
				emplace_back(elem);
			}
			else {
				auto numToMove = mSize - externalIndex;
				auto iter = begin() + externalIndex;
				if(mSize == mCapacity) [[likely]] { // NOLINT
					for(size_t i = 1ULL, j = numToMove - 2; i < numToMove; ++i, --j) {
						*((end() - i) + 1) = std::move(*(iter + j));
					}
				}
				else {
					for(size_t i = 0ULL, j = numToMove - 1; i < numToMove; ++i, --j) {
						*(end() - i) = std::move(*(iter + j));
					}
				}
				incrementIndices();
				if(externalIndex == 0) {
					// clang-format off
					mBuffer[mStartIndex] = std::move(elem); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
				else {
					// clang-format off
					mBuffer[index] = std::move(elem); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
			}
		}

		/// @brief Inserts the given element at the position indicated
		/// by the `externalIndex`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param externalIndex - The user-facing index into the `RingBuffer` to insert the element
		/// at
		/// @param elem - The element to store in the `RingBuffer`
		constexpr inline auto
		insert_internal(size_t externalIndex, T&& elem) noexcept -> void requires NotMovable<T> {
			const auto index = getAdjustedInternalIndex(externalIndex);

			if(index == mWriteIndex) {
				emplace_back(elem);
			}
			else {
				auto numToMove = mSize - externalIndex;
				auto iter = begin() + externalIndex;
				if(mSize == mCapacity) [[likely]] { // NOLINT
					for(size_t i = 1ULL, j = numToMove - 2; i < numToMove; ++i, --j) {
						*((end() - i) + 1) = *(iter + j);
					}
				}
				else {
					for(size_t i = 0ULL, j = numToMove - 1; i < numToMove; ++i, --j) {
						*(end() - i) = *(iter + j);
					}
				}
				incrementIndices();
				if(externalIndex == 0) {
					// clang-format off
					mBuffer[mStartIndex] = elem; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
				else {
					// clang-format off
					mBuffer[index] = elem; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
			}
		}

		/// @brief Constructs the given element at the insertion position indicated
		/// by the `externalIndex`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param externalIndex - The user-facing index into the `RingBuffer` to insert the element
		/// at
		/// @param args - The arguments to the constructor for
		/// the element to store in the `RingBuffer`
		template<typename... Args>
		requires ConstructibleFrom<T, Args...>
		constexpr inline auto insert_emplace_internal(size_t externalIndex, Args&&... args) noexcept
			-> T& requires Movable<T> {
			const auto index = getAdjustedInternalIndex(externalIndex);

			if(index == mWriteIndex) {
				return emplace_back(args...);
			}
			else {
				auto numToMove = mSize - externalIndex;
				auto iter = begin() + externalIndex;
				if(mSize == mCapacity) [[likely]] { // NOLINT
					for(size_t i = 1ULL, j = numToMove - 2; i < numToMove; ++i, --j) {
						*((end() - i) + 1) = std::move(*(iter + j));
					}
				}
				else {
					for(size_t i = 0ULL, j = numToMove - 1; i < numToMove; ++i, --j) {
						*(end() - i) = std::move(*(iter + j));
					}
				}
				incrementIndices();
				if(externalIndex == 0) {
					// clang-format off
					new (&mBuffer[mStartIndex]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					return mBuffer[mStartIndex];// NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
				else {
					// clang-format off
					new(&mBuffer[index]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					return mBuffer[index];// NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
			}
		}

		/// @brief Constructs the given element at the insertion position indicated
		/// by the `externalIndex`
		/// @note if `size() == capacity()` this drops the last element out of the `RingBuffer`
		///
		/// @param externalIndex - The user-facing index into the `RingBuffer` to insert the element
		/// at
		/// @param args - The arguments to the constructor for
		/// the element to store in the `RingBuffer`
		template<typename... Args>
		requires ConstructibleFrom<T, Args...>
		constexpr inline auto insert_emplace_internal(size_t externalIndex, Args&&... args) noexcept
			-> T& requires NotMovable<T> {
			const auto index = getAdjustedInternalIndex(externalIndex);

			if(index == mWriteIndex) {
				return emplace_back(args...);
			}
			else {
				auto numToMove = mSize - externalIndex;
				auto iter = begin() + externalIndex;
				if(mSize == mCapacity) [[likely]] { // NOLINT
					for(size_t i = 1ULL, j = numToMove - 2; i < numToMove; ++i, --j) {
						*((end() - i) + 1) = *(iter + j);
					}
				}
				else {
					for(size_t i = 0ULL, j = numToMove - 1; i < numToMove; ++i, --j) {
						*(end() - i) = *(iter + j);
					}
				}
				incrementIndices();
				if(externalIndex == 0) {
					// clang-format off
					new (&mBuffer[mStartIndex]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					return mBuffer[mStartIndex];// NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
				else {
					// clang-format off
					new(&mBuffer[index]) T(args...); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					return mBuffer[index];// NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}
			}
		}

		/// @brief Erases the element at the given index, returning an `Iterator` to the element
		/// after the removed one
		///
		/// @param index - The index to the element to remove. This should be a `RingBuffer` index:
		/// IE, not an interal one into the `T` array
		///
		/// @return `Iterator` pointing to the element after the one removed
		[[nodiscard]] constexpr inline auto erase_internal(size_t index) noexcept -> Iterator {
			auto indexInternal = getAdjustedInternalIndex(index);
			if(indexInternal == mWriteIndex) [[unlikely]] { // NOLINT
				return end();
			}
			else {
				auto numToMove = (mSize - 1) - index;
				mWriteIndex = indexInternal;
				mSize -= numToMove + 1;
				auto posToMove = index + 1;
				for(auto i = 0ULL; i < numToMove; ++i) {
					// clang-format off
					emplace_back(mBuffer[getAdjustedInternalIndex(posToMove + i)]); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}

				auto iter = begin() + getExternalIndexFromInternal(mWriteIndex);
				return iter;
			}
		}

		/// @brief Erases the range of elements in [`first`, `last`)
		/// Returns an `Iterator` to the element after the last one erased
		///
		/// @param first - The first index in the range to erase. This should be a `RingBuffer`
		/// index: IE, not an internal one into the `T` array
		/// @param last - The last index` in the range to erase. This should be a `RingBuffer`
		/// index: IE, not an internal one into the `T` array
		///
		/// @return `Iterator` pointing to the element after the last one erased
		[[nodiscard]] constexpr inline auto
		erase_internal(size_t first, size_t last) noexcept -> Iterator {
			auto firstInternal = getAdjustedInternalIndex(first);
			auto lastInternal = getAdjustedInternalIndex(last);
			if(lastInternal == mWriteIndex) {
				if(mWriteIndex == mLoopIndex) {
					mStartIndex--;
					mWriteIndex -= last - first;
				}
				else if(mWriteIndex < mStartIndex) {
					auto numToRemove = (last - first);
					auto numBeforeZero = numToRemove - mWriteIndex;
					auto numAfterZero = numToRemove - numBeforeZero;
					mStartIndex -= numBeforeZero;
					if(numAfterZero > 0) {
						mWriteIndex = mLoopIndex;
						numAfterZero--;
						mWriteIndex -= numAfterZero;
					}
					else {
						mWriteIndex -= numBeforeZero;
					}
				}
				else if(mWriteIndex > mStartIndex) {
					auto numToRemove = (last - first);
					mWriteIndex -= numToRemove;
				}
				mSize -= last - first;
				return end();
			}
			else {
				auto numToMove = mSize - last;
				mWriteIndex = firstInternal;
				mSize -= numToMove + (last - first);
				auto posToMove = last;
				for(auto i = 0ULL; i < numToMove; ++i) {
					// clang-format off
					emplace_back(mBuffer[getAdjustedInternalIndex(posToMove + i)]); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
					// clang-format on
				}

				auto iter = begin() + first;
				return iter;
			}
		}
	};

} // namespace utils


text:
/* ******************************************************************************
 *
 *
 * This program and the accompanying materials are made available under the
 * terms of the Apache License, Version 2.0 which is available at
 * https://www.apache.org/licenses/LICENSE-2.0.
 *
 *  See the NOTICE file distributed with this work for additional
 *  information regarding copyright ownership.
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
 * License for the specific language governing permissions and limitations
 * under the License.
 *
 * SPDX-License-Identifier: Apache-2.0
 ******************************************************************************/

//
// @author raver119@gmail.com
//

#ifndef LIBND4J_GEMM_H
#define LIBND4J_GEMM_H

// work around conflict with OpenBLAS
struct bfloat16;
#define BFLOAT16 BFLOAT16
#include <cblas.h>
#include <math/templatemath.h>
#include <system/op_boilerplate.h>

namespace sd {
namespace blas {
template <typename T>
static void *transpose(int orderSource, int orderTarget, int rows, int cols, void *source);

static inline int linearIndexC(int rows, int cols, int r, int c);
static inline int linearIndexF(int rows, int cols, int r, int c);

template <typename X, typename Y, typename Z>
class GEMM {
 protected:
 public:
  static void op(int Order, int TransA, int TransB, int M, int N, int K, double alpha, void *A, int lda, void *B,
                 int ldb, double beta, void *C, int ldc);
};

template <typename X, typename Y, typename Z>
class GEMV : public sd::blas::GEMM<X, Y, Z> {
 public:
  static void op(int TRANS, int M, int N, double alpha, void *vA, int lda, void *vX, int incx, double beta, void *vY,
                 int incy);
};

int SD_INLINE linearIndexC(int rows, int cols, int r, int c) { return (r * cols + c); }

int SD_INLINE linearIndexF(int rows, int cols, int r, int c) { return (c * rows + r); }

}  // namespace blas
}  // namespace sd

#endif  // LIBND4J_GEMM_H

ext: h


text:
#pragma once
//=============================================================================
// EXTERNAL DECLARATIONS
//=============================================================================
#include "Core/Game/IManagerComponent.h"
#include "Core/Game/IManager.h"
#include "Core/Game/IGameEngine.h"
#include "Core/Components/ComponentBase.h"
#include <gsl/span>

namespace engine {

//=============================================================================
// CLASS GameManager
//=============================================================================
template<class TDerived>
class GameManager
	: virtual public IManagerComponent
	, virtual public IManager
{
public:
	virtual void onAttached(const GameEngineRef &iGameEngine) override {}
	virtual void onDetached(const GameEngineRef &iGameEngine) override {}

	virtual gsl::span<IdType> interfaces() override
	{
		return TDerived::Interfaces;
	}
};
} // namespace engine


text:
/*
  NAME:
      read_data
  PURPOSE:
      read the data from the data input file
  CALLING SEQUENCE:
      read_data(inputfilename)
  INPUT:
      inputfilename - name of the data file
  OUPUT:
      reads everything into the right structures
  REVISION HISTORY:
     2008-09-21 - Written Bovy
*/
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <proj_gauss_main.h>
#include <proj_gauss_mixtures.h>

bool read_data(char inputfilename[]){

  FILE *inputfile;
  if ( (inputfile= fopen(inputfilename,"r")) == NULL){
    printf ("Opening the data file failed...\n");
    return false;
  }

  //First count the number of datapoints
  N=0;
  while (feof(inputfile) ==0) if (getc(inputfile) == '\n') N++;

  printf("%i datapoints found in ",N);
  printf("%s", inputfilename);
  printf("\n");

  //Then read the actual data
  fseek(inputfile,0,SEEK_SET);

  //allocate data
  data = (struct datapoint *) malloc(N*sizeof (struct datapoint));
  if (data == NULL){
    printf("Allocation of arrays failed, not enough free memory?\n");
    exit(-1);
  }
  startdata = data;

  //current line holds a maximum of d+dV+d*d elements
  //First read the whole line
  double *curr_line = (double *) calloc(d+dV+d*d, sizeof (double) );
  bool end_line;
  int ii=0,dd=0,di;//di is the dimension of the individual datapoints
  while (feof(inputfile) == 0){
    char curr_value[20]="";
    end_line = read_till_sep(curr_value,inputfile,'|');
    *(curr_line++) = atof(curr_value);
    ++dd;
    if (end_line){
      ++ii;
      curr_line -= dd;
      //Determine the dimension of the datapoint, di, from the equation di+di*(di+1)/2+d*di=dd, or,
      di = (int) ((-(3 + 2 * d) + sqrt((3 + 2 * d)*(3 + 2 * d)+8 * dd))/2);
      //then write data values to memory
      //first allocate the space for the data
      data->ww = gsl_vector_alloc(di);
      data->SS = gsl_matrix_alloc(di,di);
      data->RR = gsl_matrix_alloc(di,d);
      int dd1,dd2;
      for (dd1=0; dd1 != di; ++dd1)
	gsl_vector_set(data->ww,dd1,*(curr_line++));
      for (dd1=0; dd1 != di; ++dd1)
	gsl_matrix_set(data->SS,dd1,dd1,*(curr_line++));
      for (dd1=0; dd1 != di-1; ++dd1)
	for (dd2=dd1+1; dd2 != di; ++dd2){
	  gsl_matrix_set(data->SS,dd1,dd2,*curr_line);
	  gsl_matrix_set(data->SS,dd2,dd1,*curr_line);
	  curr_line++;
	}
      for (dd1=0; dd1 != di; ++dd1)
	for (dd2=0; dd2 != d; ++dd2)
	  gsl_matrix_set(data->RR,dd1,dd2,*(curr_line++));
      curr_line -= dd;
      dd = 0;
      data++;
      if (ii == N) break;
    }
  }

  data = startdata;

  fclose(inputfile);
  free(curr_line);

  return true;
}

ext: c


text:
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>

#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_linalg.h>

#define maxnd 10000
#define H 0.5

/*
  PROGRAM: simSFS
  AUTHOR: Athanasios Kousathanas
  VERSION: 1.0

  Generate allele frequency vectors (AFVs) using Wright-Fisher transition matrix methods and sample
  to obtain site frequency spectra (SFS). Code implements methods described by Keightley and Eyre-Walker(2007).

  Features:
  Demography: 1-epoch or 2-epoch population changes
  Selection: only single s

  Compilation:
  use makefile or compile as
  gcc -O2 -o simSFS simSFS.v1.0.c -lm -lgsl -lgslcblas -w
*/

////////////////////////////////////////
//calculate allele frequency change, delta(q)
////////////////////////////////////////
double selfn(double q, double s, double h)
{
  double num, denom;
  num = s*q*(1-q)*(q + h*(1-2*q));
  denom = 1 + 2.0*h*s*q*(1-q) + s*q*q;
  if (denom== 0) denom = .00000000000001;
  return( q + (num/denom));
}
////////////////////////////////////////
//initialise rows of matrix (used by setupmatrix)
////////////////////////////////////////
setuprow(double **a, double p, int row, int n){
  double z;         /*probability of zero failures*/
  double x, y, temp1, temp2;
  int j;

  temp1 = 1.0 - p;
  if (temp1 == 0.0)  temp1 = .000000000001;            /*prevent zero divide*/
  temp1 = p/temp1;       /*First prevent negative log*/
  if (temp1 <=0)   temp1 = .000000000001;
  x = log(temp1);    /*used for later multiplication of series*/
  temp2 = 1.0-p;
  if (temp2<=0)   temp2 = .000000000001;
  z = (double)n*log(temp2);
  if (z > 150)   z = 150;
  a[row][0] = exp(z);
  y = a[row][0];        /*add up for later*/
  for (j = 1;j<=n; j++)
    {

      z = z + log(((double)n + 1.0 - (double)j)/(double)j) + x;
      /*             cumulative binomial coeff.*/
      a[row][j] = exp(z);
      y = y + a[row][j];          /*add up*/
    }
  if (y == 0.0)   y = .0000000000001;
  for (j=0; j<=n; j++) a[row][j] = a[row][j]/y;
}
////////////////////////////////////////
//initialise transition matrix
////////////////////////////////////////
setupmatrix(double **a, int n1, int n2, double s, double h){
  int i;
  double p;

  for (i = 1; i<=n1 - 1; i++)             /*do variable rows*/
    {
      p = (double)i/(double)n1;
      p = selfn(p, s,h);
      setuprow(a, p, i, n2);
    }
  for (i = 0; i<=n2; i++)
    {              /*do constant edges*/
      a[0][i] = 0;
      a[n1][i] = 0;
    }
  a[0][0] = 1.0;                          /*do absorbing corners*/
  a[n1][n2] = 1.0;

}
////////////////////////////////////////
//Perform matrix inversion
////////////////////////////////////////
matrixinvert(int N, double **a,double *inverter){

  int i, j, lotkin_signum;

  gsl_matrix  *lotkin_a,*lotkin_inv;
  gsl_vector *x, *lotkin_b, *lotkin_x;

  gsl_permutation *lotkin_perm;

  /* allocate a, x, b */
  lotkin_a = gsl_matrix_alloc(N-1, N-1);
  lotkin_inv = gsl_matrix_alloc(N-1, N-1);
  x = gsl_vector_alloc(N-1);
  lotkin_b = gsl_vector_alloc(N-1);

  lotkin_x = gsl_vector_alloc(N-1);

  gsl_matrix *identity=gsl_matrix_alloc(N-1, N-1);

  for (i = 0; i<=N-2; i++)
    {

      for (j = 0; j<=N-2; j++)
	{
          gsl_matrix_set( identity, i, i,1);
	}

    }

  for (i = 0; i<=N-2; i++)
    {
      for (j = 0; j<=N-2; j++)
	{
	  gsl_matrix_set(lotkin_a,i,j,gsl_matrix_get(identity,i,j)-a[i+1][j+1]);
	}

    }

  /* LU decomposition and forward&backward substition */

  //gsl_blas_dgemv(CblasNoTrans, 1.0, lotkin_a, x, 0.0, lotkin_b);

  lotkin_perm = gsl_permutation_alloc(N-1);
  gsl_linalg_LU_decomp(lotkin_a, lotkin_perm, &lotkin_signum);
  gsl_linalg_LU_solve(lotkin_a, lotkin_perm, lotkin_b, lotkin_x);
  gsl_linalg_LU_invert (lotkin_a,lotkin_perm, lotkin_inv);

  /*apparently the inversion adds +1 to the first element of the vector
    so I substract it*/

  gsl_matrix_set(lotkin_inv,0,0,gsl_matrix_get(lotkin_inv,0,0)-1);

  for (i=0;i<=N-2;i++){

    inverter[i+1]=gsl_matrix_get(lotkin_inv,0,i);

  }

  gsl_matrix_free(lotkin_a);
  gsl_matrix_free(lotkin_inv);
  gsl_vector_free(x);
  gsl_vector_free(lotkin_b);
  gsl_vector_free(lotkin_x);
  gsl_matrix_free(identity);
  gsl_permutation_free(lotkin_perm);

}
////////////////////////////////////////
//perform TM iteration for t generations
////////////////////////////////////////
tmiterate(double **a,double *mut1,
	  int t, int n1, int n2,int decay,double *sumf){

  int k, i, j = 0;
  double z;

  double * steadystatefreq = (double*) calloc (maxnd+1, sizeof(double));
  double * mut2 = (double*) calloc (maxnd+1, sizeof(double));

  for (i=0; i<=n1; i++) steadystatefreq[i] = 0;
  for (k = 1; k<=t; k++)
    {                 /*t iterations*/

      /* Increment steadystatefreq. distribution*/
      for (i = 0; i<=n1; i++)
	{
	  steadystatefreq[i] = steadystatefreq[i] + mut1[i];
	}

      /* Perform matrix multiplication*/
      for (i = 0; i<=n2; i++)
	{
	  z = 0;
	  for (j = 0; j<=n1; j++)
	    {
	      z = z + a[j][i]*mut1[j];
	    }
	  mut2[i] = z;
	}

      /* Copy result of multiplication*/
      for (i=0; i<=n2; i++) {mut1[i] = mut2[i];

	if (decay==0){
	  sumf[i]+=mut1[i];}
	if(decay==1){sumf[i]=mut1[i];}
      }

    }

  free(steadystatefreq);
  free(mut2);
}

eqf_using_matrix_inversion(int n1,double s,double **a,double *egf_out){

  setupmatrix(a, n1, n1, s, H);
  matrixinvert(n1,a,egf_out);

}
////////////////////////////////////////
//scale AFV for mutations eliminated by selection/not mutated
////////////////////////////////////////

egf_scaling_s(int n1,double *v_s_in, double *v_s_out,double f0){
  int n1d=2*n1;
  double sumx,sumy;
  int i;
  sumx=0;
  sumy=0;

  for (i=1;i<n1d;i++){
    sumx+=v_s_in[i];
  }

  for (i=1;i<n1d;i++){
    v_s_out[i]/=sumx;
    sumy+=v_s_out[i];
  }

  v_s_out[0]=(1-sumy)*(1-f0);
}
////////////////////////////////////////
//scale AFV for sites not mutated
////////////////////////////////////////
egf_scaling_f0(int n1,double *fv,double f0){

  int i;
  int n1d=2*n1;
  for (i=1;i<n1d;i++){
    fv[i]*=(1-f0);
  }
  fv[0]+=f0;
}
////////////////////////////////////////
//binomial sampling from AFV to generate SFS
////////////////////////////////////////
binomial_sampling(int n1,int n1b,int sample, double *invy,int *discrete,gsl_rng *rgen){

  int s1,success;
  int i;

  double prob;

  /*
    based on equilibrium frequency vector (use it as a probability vector),
    I randomly generate numbers from 0 to n1d-1 for sample sites.
  */

  gsl_ran_discrete_t *r= gsl_ran_discrete_preproc (n1,invy);
  for (i=0;i<sample;i++){
    s1=gsl_ran_discrete (rgen, r);
    prob=(double)(s1)/n1;

    success=gsl_ran_binomial(rgen,prob,n1b+1);
    discrete[success]++;
  }
  gsl_ran_discrete_free(r);
}

output_sfs_to_file(int n1,int *sfs1,int *sfs2,char *filename,long seed){
  int i;
  FILE *file;

  file = fopen(filename,"w");
 fprintf(file,"#seed used: %d\n",seed);
 fprintf(file,"#neutral\n");
  for (i=0;i<=n1;i++){
    fprintf(file,"%d ",sfs1[i]);
  }
 fprintf(file,"\n#selected\n");

  for (i=0;i<=n1;i++){
    fprintf(file,"%d ",sfs2[i]);
  }

  fclose(file);
}
////////////////////////////////////////
//Calculate weighted average of two vectors
// Normally vectors w(s) x(s) by N1 and N2 respectively
////////////////////////////////////////

vector_average(int n1,int n2,double *fv1,double *fv2, double *fv_out){
  int i;
  int n1d=n1*2;
  int n2d=n2*2;
  for (i=0;i<=n2d;i++){

    fv_out[i]=(n1*fv1[i]+n2*fv2[i])/(n1+n2);

  }
}
//////////////////////////////////////
//Obtain AFV for 1-epoch model
////////////////////////////////////////
calculate_FV_one_epoch(int n1,double s,double *FV){
  int i=0;
  int n1d=2*n1;
  double **FVW = calloc(maxnd+1, sizeof(double *));

  for(i = 0; i < maxnd+1; i++){
    FVW[i] = calloc(maxnd+1,  sizeof(double));}
  eqf_using_matrix_inversion(n1d,s,FVW,FV);
  FV[1]+=1;
  for(i = 0; i < maxnd+1; i++)
    free(FVW[i]);
  free(FVW);
}
////////////////////////////////////////
//Obtain AFV for 2-epoch model
////////////////////////////////////////
calculate_FV_two_epoch(int n1,int n2,int t,double s,double *FV){
  int i=0;
  double **egf_0 = calloc(maxnd+1, sizeof(double *));
  for(i = 0; i < maxnd+1; i++){
    egf_0[i] = calloc(maxnd+1,  sizeof(double));
  }

  double * FVW = (double*) calloc (maxnd+1, sizeof(double));
  double * FVX = (double*) calloc (maxnd+1, sizeof(double));
  double * startingfreq = (double*) calloc (maxnd+1, sizeof(double));

  int n1d=n1*2;
  int n2d=n2*2;

  /*
    Contributions of mutations to AFV before the change in population size.
  */

  //obtain AFV at equilibrium

  eqf_using_matrix_inversion(n1d,s,egf_0,FVW);
  FVW[1]+=1.0;//cumulative vector includes 1st mutant
  //perform change in population size N1->N2
  setupmatrix(egf_0,n1d,n2d,s,H);
  tmiterate(egf_0,FVW,1, n1d, n2d,1,FVW);
  setupmatrix(egf_0,n2d,n2d,s,H);
  tmiterate(egf_0,FVW,t-1, n2d, n2d,1,FVW);

  //dumpvector(FVW,0,n2,"FVW");
  /*
    Contributions of mutations to AFV after the change in population size.
  */
  setupmatrix(egf_0,n2d,n2d,s,H);
  startingfreq[1]=1.0;
  tmiterate(egf_0,startingfreq, t-1, n2d, n2d,0,FVX);
  FVX[1]+=1.0;//cumulative vector includes 1st mutant
  //dumpvector(FVX,0,n2,"FVX");

  /*
    averaging of vectors FVW (N1) and FVX(N2)
  */
  vector_average(n1,n2,FVW,FVX,FV);

  for(i = 0; i < maxnd+1; i++){
    free(egf_0[i]);
  }

  free(egf_0);
  free(FVX);
  free(FVW);
  free(startingfreq);
}
////////////////////////////////////////
//print given matrix
////////////////////////////////////////

dumpmatrix(double **m, int s1,
	   int rows, int s2, int cols, char *s){
  int i, j;
  printf("\n%s\n", s);
  for (i = s1; i<=rows; i++)
    {
      for (j = s2; j<=cols; j++)
	{
	  printf(" %2.3f", m[i][j]);
	}
      printf("\n");
    }
  printf("\n");
}
////////////////////////////////////////
//print given vector
////////////////////////////////////////
dumpvector(double *v, int min, int max, char *s){
  int i, control = 0;
  //   printf("\n");
  printf("%s", s); printf("\n");
  for (i = min; i<=max; i++)
    {
      printf("%3d %lf ", i, v[i]);
      control++;
      if (control==5)
	{
	  control = 0;
	  printf("\n");
	}
    }
  printf("\n");
}

/*
////////////////////////////////////////
MAIN START
////////////////////////////////////////
*/

int main(int argc, char **argv)
{

  //initialise input variables and use default values
  int N1=100,n=10,sampleS=1000,sampleN=1000;
  int N2=N1;
  int t=N1;
  double f0=0.9,S=0;
  char outfile[maxnd+1]="sfs.out";
  long seed=time (NULL) * getpid();

  //other variables
  int option_index = 0;
  int c=0;
  static int verbose_flag;

  const gsl_rng_type * T;
  T = gsl_rng_taus;
  gsl_rng  *rgen =gsl_rng_alloc(T);
  gsl_rng_set(rgen,seed);

  int * discrete0 = (int*) calloc (maxnd+1, sizeof(int));
  int * discrete1 = (int*) calloc (maxnd+1, sizeof(int));

  double * FV0 = (double*) calloc (maxnd+1, sizeof(double));
  double * FVS = (double*) calloc (maxnd+1, sizeof(double));

  //parse input
  while (1)
    {

      static struct option long_options[] =
	{
	  {"h",no_argument,  &verbose_flag, 1},
	  {"N1",     required_argument,0, 'a'},
	  {"N2",  required_argument,0, 'b'},
	  {"n",  required_argument, 0, 'c'},
	  {"f0",  required_argument, 0, 'd'},
	  {"s",    required_argument, 0, 'e'},
	  {"LS",    required_argument, 0, 'f'},
	  {"LN",    required_argument, 0, 'g'},
	  {"t",    required_argument, 0, 'h'},
	  {"o",    required_argument, 0, 'i'},
	  {"seed",    required_argument, 0, 'j'},
	  {NULL, 0, 0, 0}
	};
      /* getopt_long stores the option index here. */
      c = getopt_long_only (argc, argv, "",
			    long_options, &option_index);

      /* Detect the end of the options. */
      if (c == -1)
	break;

      switch (c)
	{
	case 0:
	  /* If this option set a flag, do nothing else now. */
	  if (long_options[option_index].flag != 0)
	    {
	      printf("\
option\tdefault_value\tdescription\n\
N1\t100\tsize N1\n\
N2\t100\tsize N2\n\
t\t100\ttime since size change\n\
n\t10\tallele sample size\n\
f0\t0.9\t1-f0 is proportional to mutation rate\n\
s\t0\tselection coeffient\n\
LS\t1000\tno. neutral sites\n\
LN\t1000\tno. selected sites\n\
o\tsfs.out\toutput file\n\
seed\ttime*pid\tset seed for random generator\n\
");
	    }
	  break;
	  printf ("option %s", long_options[option_index].name);
	  if (optarg)
	    printf (" with arg %s", optarg);
	  printf ("\n");
	  break;

	case 'a':
	  N1=atoi(optarg);
	  break;
	case 'b':
	  N2=atoi(optarg);
	  break;
	case 'c':
	  n=atoi(optarg);
	  break;
	case 'd':
	  f0=atof(optarg);
	  break;
	case 'e':
	  S=atof(optarg);
	  break;
	case 'f':
	  sampleS=atoi(optarg);
	  break;
	case 'g':
	  sampleN=atoi(optarg);
	  break;
	case 'h':
	  t=atoi(optarg);
	  break;
	case 'i':
	  strncpy(outfile,optarg,maxnd+1);
	  break;
	case 'j':
	  seed=atoi(optarg);
	  gsl_rng_set(rgen,seed);
	  break;
	case '?':
	  abort();
	  break;
	default:
	  abort ();
	}
    }

  if (optind < argc)
    {
      printf ("non-option ARGV-elements: ");
      while (optind < argc)
	printf ("%s ", argv[optind++]);
      putchar ('\n');
    }

  if (N1>1000||N2>1000){printf ("Use N1,N2 <1000\n"); abort();}
  /*
    calculate allele frequency vectors for neutral (FV0)
    and selected sites (FVS)
  */
  if (N1==N2){
    calculate_FV_one_epoch(N1,0.0,FV0);
    calculate_FV_one_epoch(N1,S,FVS);
  }else{
    calculate_FV_two_epoch(N1,N2,t,0.0,FV0);
    calculate_FV_two_epoch(N1,N2,t,S,FVS);
  }

  /*
    scale FVS to include the frequency of sites at which mutations
    have been eliminated by selection or have not experienced a mutation
  */
  egf_scaling_s(N2,FV0,FVS, f0);
  egf_scaling_f0(N2,FVS,f0);

  /*
    scale FV0 to account for sites that have never experienced a mutation.
  */
  egf_scaling_s(N2,FV0, FV0,f0);
  egf_scaling_f0(N2,FV0,f0);

  /*sampling
    sample sampleN neutral sites from FV0
    sample sampleS selected sites from FVS
  */

  binomial_sampling(N2,n,sampleN,FV0,discrete0,rgen);
  binomial_sampling(N2,n,sampleS, FVS,discrete1,rgen);

  output_sfs_to_file(n,discrete0,discrete1,outfile,seed);

  //free vectors
  free(FV0);
  free(FVS);

  free(discrete0);
  free(discrete1);

  gsl_rng_free (rgen);

  return 0;
}

/*
////////////////////////////////////////
END OF MAIN
////////////////////////////////////////
*/

max_stars_repo_name: kousathanas/simsfs

max_issues_repo_name: kousathanas/simsfs

max_forks_repo_name: kousathanas/simsfs


text:
/***************************************************************************
 * data_cf.h is part of Math Graphic Library
 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru>       *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 3 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#ifndef _MGL_DATAC_CF_H_
#define _MGL_DATAC_CF_H_
//-----------------------------------------------------------------------------
#include "mgl2/abstract.h"
//-----------------------------------------------------------------------------
#if MGL_HAVE_GSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#else
struct gsl_vector;
struct gsl_matrix;
#endif
//-----------------------------------------------------------------------------
#ifdef __cplusplus
class mglDataC;
typedef mglDataC* HADT;
extern "C" {
#else
typedef void *HADT;
#endif
/// Get integer power of x
cmdual MGL_EXPORT_CONST mgl_ipowc(mdual x,int n);
cmdual MGL_EXPORT mgl_ipowc_(mdual *x,int *n);
/// Get complex number from string. Parse (%g,%g), {%g,%g} and [%g,%g] if adv!=0.
cmdual MGL_EXPORT mgl_atoc(const char *s, int adv);
/// Get exp(i*a)
cmdual MGL_EXPORT_CONST mgl_expi(mdual a);

/// Create HMDT object
HADT MGL_EXPORT mgl_create_datac();
uintptr_t MGL_EXPORT mgl_create_datac_();
/// Create HMDT object with specified sizes
HADT MGL_EXPORT mgl_create_datac_size(long nx, long ny, long nz);
uintptr_t MGL_EXPORT mgl_create_datac_size_(int *nx, int *ny, int *nz);
/// Create HMDT object with data from file
HADT MGL_EXPORT mgl_create_datac_file(const char *fname);
uintptr_t MGL_EXPORT mgl_create_datac_file_(const char *fname, int len);
/// Delete HMDT object
void MGL_EXPORT mgl_delete_datac(HADT dat);
void MGL_EXPORT mgl_delete_datac_(uintptr_t *dat);

/// Rearange data dimensions
void MGL_EXPORT mgl_datac_rearrange(HADT dat, long mx,long my,long mz);
void MGL_EXPORT mgl_datac_rearrange_(uintptr_t *dat, int *mx, int *my, int *mz);
/// Link external data array (don't delete it at exit)
void MGL_EXPORT mgl_datac_link(HADT dat, mdual *A,long mx,long my,long mz);
void MGL_EXPORT mgl_datac_link_(uintptr_t *d, mdual *A, int *nx,int *ny,int *nz);
/// Allocate memory and copy the data from the (float *) array
void MGL_EXPORT mgl_datac_set_float(HADT dat, const float *A,long mx,long my,long mz);
void MGL_EXPORT mgl_datac_set_float_(uintptr_t *dat, const float *A,int *NX,int *NY,int *NZ);
/// Allocate memory and copy the data from the (double *) array
void MGL_EXPORT mgl_datac_set_double(HADT dat, const double *A,long mx,long my,long mz);
void MGL_EXPORT mgl_datac_set_double_(uintptr_t *dat, const double *A,int *NX,int *NY,int *NZ);
/// Allocate memory and copy the data from the (dual *) array
void MGL_EXPORT mgl_datac_set_complex(HADT dat, const mdual *A,long mx,long my,long mz);
void MGL_EXPORT mgl_datac_set_complex_(uintptr_t *d, const mdual *A,int *NX,int *NY,int *NZ);
/// Import data from abstract type
void MGL_EXPORT mgl_datac_set(HADT dat, HCDT a);
void MGL_EXPORT mgl_datac_set_(uintptr_t *dat, uintptr_t *a);
/// Allocate memory and copy the data from the gsl_vector
void MGL_EXPORT mgl_datac_set_vector(HADT dat, gsl_vector *v);
/// Allocate memory and copy the data from the gsl_matrix
void MGL_EXPORT mgl_datac_set_matrix(HADT dat, gsl_matrix *m);
/// Set value of data element [i,j,k]
void MGL_EXPORT mgl_datac_set_value(HADT dat, mdual v, long i, long j, long k);
void MGL_EXPORT mgl_datac_set_value_(uintptr_t *d, mdual *v, int *i, int *j, int *k);
/// Get value of data element [i,j,k]
cmdual MGL_EXPORT mgl_datac_get_value(HCDT dat, long i, long j, long k);
cmdual MGL_EXPORT mgl_datac_get_value_(uintptr_t *d, int *i, int *j, int *k);
/// Allocate memory and scanf the data from the string
void MGL_EXPORT mgl_datac_set_values(HADT dat, const char *val, long nx, long ny, long nz);
void MGL_EXPORT mgl_datac_set_values_(uintptr_t *d, const char *val, int *nx, int *ny, int *nz, int l);

/// Get array as solution of tridiagonal matrix solution a[i]*x[i-1]+b[i]*x[i]+c[i]*x[i+1]=d[i]
/** String \a how may contain:
 * 'x', 'y', 'z' for solving along x-,y-,z-directions, or
 * 'h' for solving along hexagonal direction at x-y plain (need nx=ny),
 * 'c' for using periodical boundary conditions,
 * 'd' for diffraction/diffuse calculation.
 * NOTE: It work for flat data model only (i.e. for a[i,j]==a[i+nx*j]) */
HADT MGL_EXPORT mgl_datac_tridmat(HCDT A, HCDT B, HCDT C, HCDT D, const char *how);
uintptr_t MGL_EXPORT mgl_datac_tridmat_(uintptr_t *A, uintptr_t *B, uintptr_t *C, uintptr_t *D, const char *how, int);

/// Returns pointer to internal data array
MGL_EXPORT mdual *mgl_datac_data(HADT dat);
/// Returns pointer to data element [i,j,k]
MGL_EXPORT mdual *mgl_datac_value(HADT dat, long i,long j,long k);

/// Set the data from HCDT objects for real and imaginary parts
void MGL_EXPORT mgl_datac_set_ri(HADT dat, HCDT re, HCDT im);
void MGL_EXPORT mgl_datac_set_ri_(uintptr_t *dat, uintptr_t *re, uintptr_t *im);
/// Set the data from HCDT objects as amplitude and phase of complex data
void MGL_EXPORT mgl_datac_set_ap(HADT dat, HCDT abs, HCDT phi);
void MGL_EXPORT mgl_datac_set_ap_(uintptr_t *dat, uintptr_t *abs, uintptr_t *phi);

/// Read data from tab-separated text file with auto determining size
int MGL_EXPORT mgl_datac_read(HADT dat, const char *fname);
int MGL_EXPORT mgl_datac_read_(uintptr_t *d, const char *fname,int l);
/// Read data from text file with size specified at beginning of the file
int MGL_EXPORT mgl_datac_read_mat(HADT dat, const char *fname, long dim);
int MGL_EXPORT mgl_datac_read_mat_(uintptr_t *dat, const char *fname, int *dim, int);
/// Read data from text file with specifeid size
int MGL_EXPORT mgl_datac_read_dim(HADT dat, const char *fname,long mx,long my,long mz);
int MGL_EXPORT mgl_datac_read_dim_(uintptr_t *dat, const char *fname,int *mx,int *my,int *mz,int);
/// Read data from tab-separated text files with auto determining size which filenames are result of sprintf(fname,templ,t) where t=from:step:to
int MGL_EXPORT mgl_datac_read_range(HADT d, const char *templ, double from, double to, double step, int as_slice);
int MGL_EXPORT mgl_datac_read_range_(uintptr_t *d, const char *fname, mreal *from, mreal *to, mreal *step, int *as_slice,int l);
/// Read data from tab-separated text files with auto determining size which filenames are satisfied to template (like "t_*.dat")
int MGL_EXPORT mgl_datac_read_all(HADT dat, const char *templ, int as_slice);
int MGL_EXPORT mgl_datac_read_all_(uintptr_t *d, const char *fname, int *as_slice,int l);
/// Save whole data array (for ns=-1) or only ns-th slice to text file
void MGL_EXPORT mgl_datac_save(HCDT dat, const char *fname,long ns);
void MGL_EXPORT mgl_datac_save_(uintptr_t *dat, const char *fname,int *ns,int);

/// Read data array from HDF file (parse HDF4 and HDF5 files)
int MGL_EXPORT mgl_datac_read_hdf(HADT d,const char *fname,const char *data);
int MGL_EXPORT mgl_datac_read_hdf_(uintptr_t *d, const char *fname, const char *data,int l,int n);
/// Save data to HDF file
void MGL_EXPORT mgl_datac_save_hdf(HCDT d,const char *fname,const char *data,int rewrite);
void MGL_EXPORT mgl_datac_save_hdf_(uintptr_t *d, const char *fname, const char *data, int *rewrite,int l,int n);

/// Create or recreate the array with specified size and fill it by zero
void MGL_EXPORT mgl_datac_create(HADT dat, long nx,long ny,long nz);
void MGL_EXPORT mgl_datac_create_(uintptr_t *dat, int *nx,int *ny,int *nz);
/// Transpose dimensions of the data (generalization of Transpose)
void MGL_EXPORT mgl_datac_transpose(HADT dat, const char *dim);
void MGL_EXPORT mgl_datac_transpose_(uintptr_t *dat, const char *dim,int);

/// Get sub-array of the data with given fixed indexes
HADT MGL_EXPORT mgl_datac_subdata(HCDT dat, long xx,long yy,long zz);
uintptr_t MGL_EXPORT mgl_datac_subdata_(uintptr_t *dat, int *xx,int *yy,int *zz);
/// Get sub-array of the data with given fixed indexes (like indirect access)
HADT MGL_EXPORT mgl_datac_subdata_ext(HCDT dat, HCDT xx, HCDT yy, HCDT zz);
uintptr_t MGL_EXPORT mgl_datac_subdata_ext_(uintptr_t *dat, uintptr_t *xx,uintptr_t *yy,uintptr_t *zz);
/// Get column (or slice) of the data filled by formulas of named columns
HADT MGL_EXPORT mgl_datac_column(HCDT dat, const char *eq);
uintptr_t MGL_EXPORT mgl_datac_column_(uintptr_t *dat, const char *eq,int l);
/// Get trace of the data array
HADT MGL_EXPORT mgl_datac_trace(HCDT d);
uintptr_t MGL_EXPORT mgl_datac_trace_(uintptr_t *d);
/// Resize the data to new sizes
HADT MGL_EXPORT mgl_datac_resize(HCDT dat, long mx,long my,long mz);
uintptr_t MGL_EXPORT mgl_datac_resize_(uintptr_t *dat, int *mx,int *my,int *mz);
/// Resize the data to new sizes of box [x1,x2]*[y1,y2]*[z1,z2]
HADT MGL_EXPORT mgl_datac_resize_box(HCDT dat, long mx,long my,long mz,mreal x1,mreal x2,mreal y1,mreal y2,mreal z1,mreal z2);
uintptr_t MGL_EXPORT mgl_datac_resize_box_(uintptr_t *dat, int *mx,int *my,int *mz,mreal *x1,mreal *x2,mreal *y1,mreal *y2,mreal *z1,mreal *z2);
/// Get momentum (1D-array) of data along direction 'dir'. String looks like "x1" for median in x-direction, "x2" for width in x-dir and so on.
HADT MGL_EXPORT mgl_datac_momentum(HCDT dat, char dir, const char *how);
uintptr_t MGL_EXPORT mgl_datac_momentum_(uintptr_t *dat, char *dir, const char *how, int,int);
/// Get array which values is result of interpolation this for coordinates from other arrays
HADT MGL_EXPORT mgl_datac_evaluate(HCDT dat, HCDT idat, HCDT jdat, HCDT kdat, int norm);
uintptr_t MGL_EXPORT mgl_datac_evaluate_(uintptr_t *dat, uintptr_t *idat, uintptr_t *jdat, uintptr_t *kdat, int *norm);
/// Get array which is result of summation in given direction or directions
HADT MGL_EXPORT mgl_datac_sum(HCDT dat, const char *dir);
uintptr_t MGL_EXPORT mgl_datac_sum_(uintptr_t *dat, const char *dir,int);
/// Get the data which is direct multiplication (like, d[i,j] = this[i]*a[j] and so on)
HADT MGL_EXPORT mgl_datac_combine(HCDT dat1, HCDT dat2);
uintptr_t MGL_EXPORT mgl_datac_combine_(uintptr_t *dat1, uintptr_t *dat2);
/// Get data from sections ids, separated by value val along specified direction.
/** If section id is negative then reverse order is used (i.e. -1 give last section). */
HADT MGL_EXPORT mgl_datac_section(HCDT dat, HCDT ids, char dir, mreal val);
uintptr_t MGL_EXPORT mgl_datac_section_(uintptr_t *d, uintptr_t *ids, const char *dir, mreal *val,int);
/// Get data from section id, separated by value val along specified direction.
/** If section id is negative then reverse order is used (i.e. -1 give last section). */
HADT MGL_EXPORT mgl_datac_section_val(HCDT dat, long id, char dir, mreal val);
uintptr_t MGL_EXPORT mgl_datac_section_val_(uintptr_t *d, int *id, const char *dir, mreal *val,int);

/// Equidistantly fill the data to range [x1,x2] in direction dir
void MGL_EXPORT mgl_datac_fill(HADT dat, mdual x1,mdual x2,char dir);
void MGL_EXPORT mgl_datac_fill_(uintptr_t *dat, mdual *x1,mdual *x2,const char *dir,int);
/// Modify the data by specified formula assuming x,y,z in range [r1,r2]
void MGL_EXPORT mgl_datac_fill_eq(HMGL gr, HADT dat, const char *eq, HCDT vdat, HCDT wdat,const char *opt);
void MGL_EXPORT mgl_datac_fill_eq_(uintptr_t *gr, uintptr_t *dat, const char *eq, uintptr_t *vdat, uintptr_t *wdat,const char *opt, int, int);
/// Fill dat by interpolated values of vdat parametrically depended on xdat for x in range [x1,x2] using global spline
void MGL_EXPORT mgl_datac_refill_gs(HADT dat, HCDT xdat, HCDT vdat, mreal x1, mreal x2, long sl);
void MGL_EXPORT mgl_datac_refill_gs_(uintptr_t *dat, uintptr_t *xdat, uintptr_t *vdat, mreal *x1, mreal *x2, long *sl);
/// Fill dat by interpolated values of vdat parametrically depended on xdat for x in range [x1,x2]
void MGL_EXPORT mgl_datac_refill_x(HADT dat, HCDT xdat, HCDT vdat, mreal x1, mreal x2, long sl);
void MGL_EXPORT mgl_datac_refill_x_(uintptr_t *dat, uintptr_t *xdat, uintptr_t *vdat, mreal *x1, mreal *x2, long *sl);
/// Fill dat by interpolated values of vdat parametrically depended on xdat,ydat for x,y in range [x1,x2]*[y1,y2]
void MGL_EXPORT mgl_datac_refill_xy(HADT dat, HCDT xdat, HCDT ydat, HCDT vdat, mreal x1, mreal x2, mreal y1, mreal y2, long sl);
void MGL_EXPORT mgl_datac_refill_xy_(uintptr_t *dat, uintptr_t *xdat, uintptr_t *ydat, uintptr_t *vdat, mreal *x1, mreal *x2, mreal *y1, mreal *y2, long *sl);
/// Fill dat by interpolated values of vdat parametrically depended on xdat,ydat,zdat for x,y,z in range [x1,x2]*[y1,y2]*[z1,z2]
void MGL_EXPORT mgl_datac_refill_xyz(HADT dat, HCDT xdat, HCDT ydat, HCDT zdat, HCDT vdat, mreal x1, mreal x2, mreal y1, mreal y2, mreal z1, mreal z2);
void MGL_EXPORT mgl_datac_refill_xyz_(uintptr_t *dat, uintptr_t *xdat, uintptr_t *ydat, uintptr_t *zdat, uintptr_t *vdat, mreal *x1, mreal *x2, mreal *y1, mreal *y2, mreal *z1, mreal *z2);
/// Fill dat by interpolated values of vdat parametrically depended on xdat,ydat,zdat for x,y,z in axis range
void MGL_EXPORT mgl_datac_refill_gr(HMGL gr, HADT dat, HCDT xdat, HCDT ydat, HCDT zdat, HCDT vdat, long sl, const char *opt);
void MGL_EXPORT mgl_datac_refill_gr_(uintptr_t *gr, uintptr_t *dat, uintptr_t *xdat, uintptr_t *ydat, uintptr_t *zdat, uintptr_t *vdat, long *sl, const char *opt,int);

/// Modify the data by specified formula
void MGL_EXPORT mgl_datac_modify(HADT dat, const char *eq,long dim);
void MGL_EXPORT mgl_datac_modify_(uintptr_t *dat, const char *eq,int *dim,int);
/// Modify the data by specified formula
void MGL_EXPORT mgl_datac_modify_vw(HADT dat, const char *eq,HCDT vdat,HCDT wdat);
void MGL_EXPORT mgl_datac_modify_vw_(uintptr_t *dat, const char *eq, uintptr_t *vdat, uintptr_t *wdat,int);

/// Limit the data to be inside [-v,v], keeping the original sign
void MGL_EXPORT mgl_datac_limit(HADT dat, mreal v);
void MGL_EXPORT mgl_datac_limit_(uintptr_t *dat, mreal *v);

/// Put value to data element(s)
void MGL_EXPORT mgl_datac_put_val(HADT dat, mdual val, long i, long j, long k);
void MGL_EXPORT mgl_datac_put_val_(uintptr_t *dat, mdual *val, int *i, int *j, int *k);
/// Put array to data element(s)
void MGL_EXPORT mgl_datac_put_dat(HADT dat, HCDT val, long i, long j, long k);
void MGL_EXPORT mgl_datac_put_dat_(uintptr_t *dat, uintptr_t *val, int *i, int *j, int *k);

/// Reduce size of the data
void MGL_EXPORT mgl_datac_squeeze(HADT dat, long rx,long ry,long rz,long smooth);
void MGL_EXPORT mgl_datac_squeeze_(uintptr_t *dat, int *rx,int *ry,int *rz,int *smooth);
/// Extend data dimensions
void MGL_EXPORT mgl_datac_extend(HADT dat, long n1, long n2);
void MGL_EXPORT mgl_datac_extend_(uintptr_t *dat, int *n1, int *n2);
/// Insert data rows/columns/slices
void MGL_EXPORT mgl_datac_insert(HADT dat, char dir, long at, long num);
void MGL_EXPORT mgl_datac_insert_(uintptr_t *dat, const char *dir, int *at, int *num, int);
/// Delete data rows/columns/slices
void MGL_EXPORT mgl_datac_delete(HADT dat, char dir, long at, long num);
void MGL_EXPORT mgl_datac_delete_(uintptr_t *dat, const char *dir, int *at, int *num, int);
/// Joind another data array
void MGL_EXPORT mgl_datac_join(HADT dat, HCDT d);
void MGL_EXPORT mgl_datac_join_(uintptr_t *dat, uintptr_t *d);

/// Smooth the data on specified direction or directions
/** String \a dir may contain:
 *  ‘x’, ‘y’, ‘z’ for 1st, 2nd or 3d dimension;
 *  ‘dN’ for linear averaging over N points;
 *  ‘3’ for linear averaging over 3 points;
 *  ‘5’ for linear averaging over 5 points.
 *  By default quadratic averaging over 5 points is used. */
void MGL_EXPORT mgl_datac_smooth(HADT d, const char *dirs, mreal delta);
void MGL_EXPORT mgl_datac_smooth_(uintptr_t *dat, const char *dirs, mreal *delta,int);
/// Cumulative summation the data in given direction or directions
void MGL_EXPORT mgl_datac_cumsum(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_cumsum_(uintptr_t *dat, const char *dir,int);
/// Integrate (cumulative summation) the data in given direction or directions
void MGL_EXPORT mgl_datac_integral(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_integral_(uintptr_t *dat, const char *dir,int);
/// Differentiate the data in given direction or directions
void MGL_EXPORT mgl_datac_diff(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_diff_(uintptr_t *dat, const char *dir,int);
/// Differentiate the parametrically specified data along direction v1 with v2,v3=const (v3 can be NULL)
void MGL_EXPORT mgl_datac_diff_par(HADT dat, HCDT v1, HCDT v2, HCDT v3);
void MGL_EXPORT mgl_datac_diff_par_(uintptr_t *dat, uintptr_t *v1, uintptr_t *v2, uintptr_t *v3);
/// Double-differentiate (like Laplace operator) the data in given direction
void MGL_EXPORT mgl_datac_diff2(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_diff2_(uintptr_t *dat, const char *dir,int);
/// Swap left and right part of the data in given direction (useful for Fourier spectrum)
void MGL_EXPORT mgl_datac_swap(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_swap_(uintptr_t *dat, const char *dir,int);
/// Roll data along direction dir by num slices
void MGL_EXPORT mgl_datac_roll(HADT dat, char dir, long num);
void MGL_EXPORT mgl_datac_roll_(uintptr_t *dat, const char *dir, int *num, int);
/// Mirror the data in given direction (useful for Fourier spectrum)
void MGL_EXPORT mgl_datac_mirror(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_mirror_(uintptr_t *dat, const char *dir,int);
/// Crop the data
void MGL_EXPORT mgl_datac_crop(HADT dat, long n1, long n2, char dir);
void MGL_EXPORT mgl_datac_crop_(uintptr_t *dat, int *n1, int *n2, const char *dir,int);
/// Crop the data to be most optimal for FFT (i.e. to closest value of 2^n*3^m*5^l)
void MGL_EXPORT mgl_datac_crop_opt(HADT dat, const char *how);
void MGL_EXPORT mgl_datac_crop_opt_(uintptr_t *dat, const char *how,int);

/// Multiply the data by other one for each element
void MGL_EXPORT mgl_datac_mul_dat(HADT dat, HCDT d);
void MGL_EXPORT mgl_datac_mul_dat_(uintptr_t *dat, uintptr_t *d);
/// Divide the data by other one for each element
void MGL_EXPORT mgl_datac_div_dat(HADT dat, HCDT d);
void MGL_EXPORT mgl_datac_div_dat_(uintptr_t *dat, uintptr_t *d);
/// Add the other data
void MGL_EXPORT mgl_datac_add_dat(HADT dat, HCDT d);
void MGL_EXPORT mgl_datac_add_dat_(uintptr_t *dat, uintptr_t *d);
/// Subtract the other data
void MGL_EXPORT mgl_datac_sub_dat(HADT dat, HCDT d);
void MGL_EXPORT mgl_datac_sub_dat_(uintptr_t *dat, uintptr_t *d);
/// Multiply each element by the number
void MGL_EXPORT mgl_datac_mul_num(HADT dat, mdual d);
void MGL_EXPORT mgl_datac_mul_num_(uintptr_t *dat, mdual *d);
/// Divide each element by the number
void MGL_EXPORT mgl_datac_div_num(HADT dat, mdual d);
void MGL_EXPORT mgl_datac_div_num_(uintptr_t *dat, mdual *d);
/// Add the number
void MGL_EXPORT mgl_datac_add_num(HADT dat, mdual d);
void MGL_EXPORT mgl_datac_add_num_(uintptr_t *dat, mdual *d);
/// Subtract the number
void MGL_EXPORT mgl_datac_sub_num(HADT dat, mdual d);
void MGL_EXPORT mgl_datac_sub_num_(uintptr_t *dat, mdual *d);

/// Apply Hankel transform
void MGL_EXPORT mgl_datac_hankel(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_hankel_(uintptr_t *dat, const char *dir,int);
/// Apply Sin-Fourier transform
void MGL_EXPORT mgl_datac_sinfft(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_sinfft_(uintptr_t *dat, const char *dir,int);
/// Apply Cos-Fourier transform
void MGL_EXPORT mgl_datac_cosfft(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_cosfft_(uintptr_t *dat, const char *dir,int);
/// Apply Fourier transform
void MGL_EXPORT mgl_datac_fft(HADT dat, const char *dir);
void MGL_EXPORT mgl_datac_fft_(uintptr_t *dat, const char *dir,int);
/// Find correlation between 2 data arrays
HADT MGL_EXPORT mgl_datac_correl(HCDT dat1, HCDT dat2, const char *dir);
uintptr_t MGL_EXPORT mgl_datac_correl_(uintptr_t *dat1, uintptr_t *dat2, const char *dir,int);
/// Calculate one step of diffraction by finite-difference method with parameter q
void MGL_EXPORT mgl_datac_diffr(HADT dat, const char *how, mreal q);
void MGL_EXPORT mgl_datac_diffr_(uintptr_t *d, const char *how, double q,int l);
/// Apply wavelet transform
/** Parameter \a dir may contain:
 * ‘x‘,‘y‘,‘z‘ for directions,
 * ‘d‘ for daubechies, ‘D‘ for centered daubechies,
 * ‘h‘ for haar, ‘H‘ for centered haar,
 * ‘b‘ for bspline, ‘B‘ for centered bspline,
 * ‘i‘ for applying inverse transform. */
void MGL_EXPORT mgl_datac_wavelet(HADT dat, const char *how, int k);
void MGL_EXPORT mgl_datac_wavelet_(uintptr_t *d, const char *dir, int *k,int);
/// Set as the data envelop
void MGL_EXPORT mgl_datac_envelop(HADT dat, char dir);
void MGL_EXPORT mgl_datac_envelop_(uintptr_t *dat, const char *dir, int);

/// Get real part of data values
HMDT MGL_EXPORT mgl_datac_real(HCDT dat);
uintptr_t MGL_EXPORT mgl_datac_real_(uintptr_t *dat);
/// Get imaginary part of data values
HMDT MGL_EXPORT mgl_datac_imag(HCDT dat);
uintptr_t MGL_EXPORT mgl_datac_imag_(uintptr_t *dat);
/// Get absolute value of data values, i.e. |u|
HMDT MGL_EXPORT mgl_datac_abs(HCDT dat);
uintptr_t MGL_EXPORT mgl_datac_abs_(uintptr_t *dat);
/// Get argument of data values
HMDT MGL_EXPORT mgl_datac_arg(HCDT dat);
uintptr_t MGL_EXPORT mgl_datac_arg_(uintptr_t *dat);
/// Get square of absolute value of data values, i.e. |u|^2
HMDT MGL_EXPORT mgl_datac_norm(HCDT dat);
uintptr_t MGL_EXPORT mgl_datac_norm_(uintptr_t *dat);

/// Interpolate by linear function the data to given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
cmdual MGL_EXPORT mgl_datac_linear(HCDT d, mreal x,mreal y,mreal z);
cmdual MGL_EXPORT mgl_datac_linear_(uintptr_t *d, mreal *x,mreal *y,mreal *z);
/// Interpolate by linear function the data and return its derivatives at given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
cmdual MGL_EXPORT mgl_datac_linear_ext(HCDT d, mreal x,mreal y,mreal z, mdual *dx,mdual *dy,mdual *dz);
cmdual MGL_EXPORT mgl_datac_linear_ext_(uintptr_t *d, mreal *x,mreal *y,mreal *z, mdual *dx,mdual *dy,mdual *dz);
/// Interpolate by cubic spline the data to given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
cmdual MGL_EXPORT mgl_datac_spline(HCDT dat, mreal x,mreal y,mreal z);
cmdual MGL_EXPORT mgl_datac_spline_(uintptr_t *dat, mreal *x,mreal *y,mreal *z);
/// Interpolate by cubic spline the data and return its derivatives at given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
cmdual MGL_EXPORT mgl_datac_spline_ext(HCDT dat, mreal x,mreal y,mreal z, mdual *dx,mdual *dy,mdual *dz);
cmdual MGL_EXPORT mgl_datac_spline_ext_(uintptr_t *dat, mreal *x,mreal *y,mreal *z, mdual *dx,mdual *dy,mdual *dz);
/// Prepare coefficients for global spline interpolation
HADT MGL_EXPORT mgl_gsplinec_init(HCDT x, HCDT v);
uintptr_t MGL_EXPORT mgl_gspline_init_(uintptr_t *x, uintptr_t *v);
/// Evaluate global spline (and its derivatives d1, d2 if not NULL) using prepared coefficients \a coef
cmdual MGL_EXPORT mgl_gsplinec(HCDT coef, mreal dx, mdual *d1, mdual *d2);
cmdual MGL_EXPORT mgl_gsplinec_(uintptr_t *c, mreal *dx, mdual *d1, mdual *d2);

/// Find roots for set of nonlinear equations defined by textual formulas
HADT MGL_EXPORT mgl_find_roots_txt_c(const char *func, const char *vars, HCDT ini);
uintptr_t MGL_EXPORT mgl_find_roots_txt_c_(const char *func, const char *vars, uintptr_t *ini,int,int);

#ifdef __cplusplus
}
#endif
//-----------------------------------------------------------------------------
#endif

max_stars_repo_name: angelamsj/cruise-control

max_issues_repo_name: angelamsj/cruise-control

max_forks_repo_name: angelamsj/cruise-control