text:
// MIT License Copyright (c) 2020 Jarrett Wendt

#pragma once

#include "Macros.h"
#include "Memory.h"
#include "Util.h"

#include <algorithm>			// std::min/max
#include <compare>				// std::strong_ordering
#include <functional>
#include <initializer_list>		// std::initializer/list
#include <iterator>				// std::random_access_iterator

#include <gsl/gsl>				// gsl::owner

// Some macros to save on typing.
// These get #undef-ed at the end of the .inl file.
// Also pretty convenient in case I ever decide I want to change/rename the template args.
// The only downside is that intellisense thinks my methods are unimplemented.
#define TEMPLATE template<typename T, Concept::ReserveStrategy ReserveStrategy>
#define ARRAY Array<T, ReserveStrategy>

namespace Library
{
	// TODO: Array<bool> can be specialized as a bit-array to save memory.

	/**
	 * A contiguous-memory dynamically resizing array.
	 * Internally uses std::vector and provides more convenience methods to it.
	 * Name chosen to disambiguate between vector the data container and vector the mathematical concept.
	 *
	 * @param <T>					The type for this container to store.
	 * @param <ReserveStrategy>		A callable type for determining a new capacity given the current memory usage.
	 */
	template<typename T, Concept::ReserveStrategy ReserveStrategy = Util::DefaultReserveStrategy>
	class Array final
	{
	public:
		using value_type = T;
		using reference = T&;
		using const_reference = const T&;
		using pointer = T*;
		using const_pointer = const T*;
		using size_type = size_t;
		using difference_type = ptrdiff_t;
		using reserve_strategy = ReserveStrategy;

	private:
		using internal_size_type = uint32_t;

		T* array{ nullptr };
		internal_size_type size{ 0 };
		internal_size_type capacity{ 0 };

	public:
#pragma region Special Members
		/**
		 * default ctor
		 * creates an empty container owning no memory
		 */
		Array() = default;

		/**
		 * explicit ctor
		 *
		 * @param capacity		how many elements worth of space to reserve
		 */
		explicit Array(size_type capacity);

		/**
		 * explicit ctor
		 *
		 * @param count			how many elements of space to reserve and construct
		 * @param prototype		prototypical value to be copy-constructed into each element
		 */
		Array(size_type count, const T& prototype);

		/**
		 * explicit ctor
		 * Fills the container with values read through the iterators.
		 *
		 * @param <It>			the iterator type
		 * @param first			the beginning iterator, inclusive
		 * @param last			the last iterator, exclusive
		 */
		template<std::forward_iterator It>
		Array(It first, It last);

		/**
		 * explicit ctor
		 * Fills the container with values read through the iterators.
		 *
		 * @param <It>			the iterator type
		 * @param first			the beginning iterator, inclusive
		 * @param last			the last iterator, exclusive
		 */
		template<std::random_access_iterator It>
		Array(It first, It last);

		/**
		 * typecast ctor which accepts any range
		 * will at most resize once
		 *
		 * @param range			the range to construct from
		 */
		template<Concept::RangeOf<T> Range>
		Array(const Range& range);

		/**
		 * typecast ctor which accepts any range
		 * will at most resize once
		 *
		 * @param range			the range to construct from
		 */
		template<Concept::RangeOf<T> Range>
		Array& operator=(const Range& range);

		/**
		 * initializer_list ctor
		 *
		 * @param list		the initializer_list to construct this container from
		 */
		Array(std::initializer_list<T> list);

		/**
		 * initializer_list operator=
		 *
		 * @param list		the initializer_list to assign this container to
		 * @return			this container after assignment
		 */
		Array& operator=(std::initializer_list<T> list);

		/**
		 * copy ctor
		 *
		 * @param other		the container to copy
		 */
		Array(const Array& other);

		/**
		 * move ctor
		 *
		 * @param other		the container to move, will be made empty after this operation
		 */
		Array(Array&& other) noexcept;

		/**
		 * copy operator=
		 *
		 * @param other		the container to move, will be made empty after this operation
		 * @returns			this container after assignment
		 */
		Array& operator=(const Array& other);

		/**
		 * move operator=
		 *
		 * @param other		the container to move, will be made empty after this operation
		 * @returns			this container after assignment
		 */
		Array& operator=(Array&& other) noexcept;

		/**
		 * templated copy ctor accepting containers with other ReserveStrategies
		 * this container will keep its original reserve strategy
		 *
		 * @param other		the container to move, will be made empty after this operation
		 * @returns			this container after assignment
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		Array(const Array<T, OtherReserveStrategy>& other);

		/**
		 * templated move ctor accepting containers with other ReserveStrategies
		 * this container will keep its original reserve strategy
		 *
		 * @param other		the container to move, will be made empty after this operation
		 * @returns			this container after assignment
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		Array(Array<T, OtherReserveStrategy>&& other) noexcept;

		/**
		 * templated copy operator= accepting containers with other ReserveStrategies
		 * this container will keep its original reserve strategy
		 *
		 * @param other		the container to move, will be made empty after this operation
		 * @returns			this container after assignment
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		Array& operator=(const Array<T, OtherReserveStrategy>& other);

		/**
		 * templated move operator= accepting containers with other ReserveStrategies
		 * this container will keep its original reserve strategy
		 *
		 * @param other		the container to move, will be made empty after this operation
		 * @returns			this container after assignment
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		Array& operator=(Array<T, OtherReserveStrategy>&& other) noexcept;

		/**
		 * dtor
		 * frees all associated memory
		 */
		~Array();
#pragma endregion

#pragma region iterator
		class iterator final
		{
			friend class Array;
			friend class const_iterator;

		private:
			T* data{ nullptr };
#ifdef _DEBUG
			Array* owner{ nullptr };
#endif

		public:
			using iterator_category = std::random_access_iterator_tag;
			using value_type = T;
			using difference_type = ptrdiff_t;
			using pointer = T*;
			using reference = T&;

		private:
			iterator(T* data, [[maybe_unused]] Array* owner) noexcept :
				data(data)
#ifdef _DEBUG
				,owner(owner)
#endif
			{}

		public:
			/**
			 * explicit constructor
			 *
			 * @param index		index for this iterator to start at
			 * @param owner		the container this iterator references
			 */
			iterator(const size_type index, Array& owner) noexcept :
				iterator(owner.array + index, &owner) {}

			SPECIAL_MEMBERS(iterator, default)

			/**
			 * Difference of two iterators.
			 *
			 * @param left		lhs of the operator
			 * @param right		rhs of the operator
			 * @returns			The difference between the indices of the two iterators.
			 *
			 * @asserts			The iterators are initialized and they belong to the same container.
			 */
			[[nodiscard]] friend difference_type operator-(const iterator left, const iterator right) noexcept
			{
				// No need to AssertInitialize anything because it will be called by operator +=
				// AssertInitialize is therefore implicitly called on right since we're comparing it to left's owner
				assertm(left.owner == right.owner, "iterators must belong to the same container");
				return left.data - right.data;
			}

			/**
			 * Will not go beyond end().
			 *
			 * @param i		How much to increment the iterator by.
			 * @returns		A new iterator at that index.
			 *
			 * @asserts		this iterator is initialized
			 */
			[[nodiscard]] iterator operator+(difference_type i) const noexcept;

			/**
			 * @returns		The value of the index the iterator is at.
			 *
			 * @asserts		this iterator is initialized and not at end().
			 */
			[[nodiscard]] reference operator*() const;

			RANDOM_ITER_OPS(iterator)

		private:
			/**
			 * @asserts		this iterator is initialized.
			 */
			constexpr void AssertInitialized() const noexcept;
		};

		class const_iterator final
		{
			friend class Array;
			CONST_RANDOM_ACCESS_ITERATOR(const_iterator, iterator)

			/**
			 * explicit constructor
			 *
			 * @param index		index for this iterator to start at
			 * @param owner		the container this iterator references
			 */
			const_iterator(const size_t index, const Array& owner) noexcept :
				it{ index, const_cast<Array&>(owner) } {};
		};

		BEGIN_END(iterator, const_iterator, Array)
		MAKE_REVERSE_BEGIN_END
#pragma endregion

#pragma region Properties
		/**
		 * @returns		true if the container is empty, false otherwise
		 */
		[[nodiscard]] constexpr bool IsEmpty() const noexcept;

		/**
		 * @returns		true if the container is full, false otherwise
		 */
		[[nodiscard]] constexpr bool IsFull() const noexcept;

		/**
		 * @returns		How many elements are in the container.
		 */
		[[nodiscard]] constexpr size_type Size() const noexcept;

		/**
		 * @returns		The number of elements this container can hold without resizing.
		 */
		[[nodiscard]] constexpr size_type Capacity() const noexcept;

		/**
		 * @param it	an iterator
		 * @returns		the index represented by the iterator
		 */
		[[nodiscard]] constexpr size_type IndexOf(const_iterator it) const noexcept;
#pragma endregion

#pragma region Element Access
		/**
		 * Returns a reference to the given index with bounds checking.
		 *
		 * @param index		position of the element to return
		 * @returns			reference to the requested element
		 *
		 * @throws std::out_of_range	if index >= Size()
		 */
		[[nodiscard]] T& At(size_type index);

		/**
		 * Returns a reference to the given index with bounds checking.
		 *
		 * @param index		position of the element to return
		 * @returns			reference to the requested element
		 *
		 * @throws std::out_of_range	if index >= Size()
		 */
		[[nodiscard]] const T& At(size_type index) const;

		/**
		 * Returns a reference to the given index with bounds checking.
		 *
		 * @param index		position of the element to return
		 * @returns			reference to the requested element
		 *
		 * @throws std::out_of_range	if index >= Size()
		 */
		[[nodiscard]] T& operator[](size_type index);

		/**
		 * Returns a reference to the given index with bounds checking.
		 *
		 * @param index		position of the element to return
		 * @returns			reference to the requested element
		 *
		 * @throws std::out_of_range	if index >= Size()
		 */
		[[nodiscard]] const T& operator[](size_type index) const;

		/**
		 * @returns		reference to the first element.
		 *
		 * @throws std::out_of_range	if the container is empty
		 */
		[[nodiscard]] T& Front();

		/**
		 * @returns		reference to the first element.
		 *
		 * @throws std::out_of_range	if the container is empty
		 */
		[[nodiscard]] const T& Front() const;

		/**
		 * @returns		reference to the last element.
		 *
		 * @throws std::out_of_range	if the container is empty
		 */
		[[nodiscard]] T& Back();

		/**
		 * @returns		reference to the last element.
		 *
		 * @throws std::out_of_range	if the container is empty
		 */
		[[nodiscard]] const T& Back() const;
#pragma endregion

#pragma region Insert
		/**
		 * @param index		where the new element will be constructed
		 * @param t			value to insert
		 */
		void Insert(size_t index, const T& t);

		/**
		 * @param index		where the new element will be constructed
		 * @param t			value to insert
		 */
		void Insert(size_t index, T&& t);

		/**
		 * @param index		where the new element will be constructed
		 * @param count		how many copies of the value to insert
		 * @param t			value to insert
		 * @returns			iterator at the position of the last element inserted
		 */
		iterator Insert(size_t index, size_t count, const T& t);

		/**
		 * @param index		where the new element will be constructed
		 * @param list		initializer_list of values to insert
		 * @returns			iterator at the position of the last element inserted
		 */
		iterator Insert(size_t index, std::initializer_list<T> list);

		/**
		 * Insert elements through the values read from the iterators in the range [first, last).
		 *
		 * @param index		where the first new element will be constructed
		 * @param first		beginning iterator, inclusive
		 * @param last		ending iterator, exclusive
		 * @returns			iterator at the position of the last element inserted
		 *
		 * @throws			std::invalid_argument if the iterators belong to this container
		 */
		template<std::forward_iterator It>
		iterator Insert(size_t index, It first, It last);

		/**
		 * @param it		position of where to insert
		 * @param t			prototypical value to insert
		 *
		 * @throws std::invalid_argument	if the iterator does not belong to this container
		 */
		void Insert(const_iterator it, const T& t);

		/**
		 * @param it		position of where to insert
		 * @param t			prototypical value to insert
		 *
		 * @throws std::invalid_argument	if the iterator does not belong to this container
		 */
		void Insert(const_iterator it, T&& t);

		/**
		 * @param it		initial position of where to insert
		 * @param count		how many copies of the prototype to insert
		 * @param t			prototypical value to insert
		 * @returns			iterator at the last element inserted
		 *
		 * @throws std::invalid_argument	if the iterator does not belong to this container
		 */
		iterator Insert(const_iterator it, size_t count, const T& t);

		/**
		 * @param it		initial position of where to insert
		 * @param list		initializer_list of values to insert
		 * @returns			iterator at the last element inserted
		 *
		 * @throws std::invalid_argument	if the iterator does not belong to this container
		 */
		iterator Insert(const_iterator it, std::initializer_list<T> list);

		/**
		 * Insert elements through the values read from the iterators in the range [first, last).
		 *
		 * @param it		initial position of where to insert
		 * @param first		beginning iterator, inclusive
		 * @param last		ending iterator, exclusive
		 * @returns			iterator at the position of the last element inserted
		 *
		 * @throws			std::invalid_argument if the iterators do not belong to this container
		 */
		template<std::forward_iterator It>
		iterator Insert(const_iterator it, It first, It last);

		/**
		 * Inserts a new element into the container directly at index.
		 *
		 * @param index		Position in the container to put this element.
		 * @param args		Arguments to forward to the constructor of the element.
		 *
		 * @throws std::out_of_range	if index > Size().
		 */
		template<typename... Args>
		reference Emplace(size_t index, Args&&... args);

		/**
		 * Inserts a new element into the container directly at index.
		 * The passed iterator is guaranteed to be valid and pointing at the position of the newly emplaced element after the operation.
		 *
		 * @param it		Position in the container to put this element.
		 * @param args		Arguments to forward to the constructor of the element.
		 *
		 * @throws std::out_of_range	if index > Size().
		 */
		template<typename... Args>
		reference Emplace(const_iterator it, Args&&... args);

		/**
		 * Appends a new element to the end of the container.
		 * O(1) most cases.
		 * O(n) if the operation requires a Resize().
		 *
		 * @param args		Arguments to forward to the constructor of the element.
		 */
		template<typename... Args>
		reference EmplaceBack(Args&&... args);

		/**
		 * Appends a new element to the front of the container.
		 * O(1) most cases.
		 * O(n) if the operation requires a Resize().
		 *
		 * @param args		Arguments to forward to the constructor of the element.
		 */
		template<typename... Args>
		reference EmplaceFront(Args&&... args);

		/**
		 * Appends the element to the end of the container.
		 * O(1) most cases.
		 * O(n) if the operation requires a Resize().
		 *
		 * @param t		The value of the element to append.
		 */
		void PushBack(const T& t);

		/**
		 * Appends the element to the end of the container.
		 * O(n) where n is Size().
		 *
		 * @param t		The value of the element to append.
		 */
		void PushBack(T&& t);

		/**
		 * Appends the element to the front of the container.
		 * O(n) where n is Size().
		 *
		 * @param t		The value of the element to prepend.
		 */
		void PushFront(const T& t);

		/**
		 * Appends the element to the front of the container.
		 * O(n) where n is Size().
		 *
		 * @param t		The value of the element to prepend.
		 */
		void PushFront(T&& t);

		/**
		 * Insert elements through the values read from the iterators in the range [first, last).
		 * Only resizes once (this is most efficient if It is a random access iterator).
		 *
		 * @param it		initial position of where to insert
		 * @param first		beginning iterator, inclusive
		 * @param last		ending iterator, exclusive
		 * @returns			iterator at the position of the last element inserted
		 */
		template<std::forward_iterator It>
		void PushBack(It first, It last);

		/**
		 * Insert elements through the values read from the iterators in the range [first, last).
		 * Only resizes once (this is most efficient if It is a random access iterator).
		 *
		 * @param it		initial position of where to insert
		 * @param first		beginning iterator, inclusive
		 * @param last		ending iterator, exclusive
		 * @returns			iterator at the position of the last element inserted
		 */
		template<std::forward_iterator It>
		void PushFront(It first, It last);

		/**
		 * Copies the passed prototypical value into the remaining empty slots.
		 * O(n) where n is Capacity() - Size().
		 *
		 * @param prototype		The prottype to fill the container with.
		 */
		void Fill(const T& prototype);
#pragma endregion

#pragma region Remove
		/**
		 * Removes the first instance of the passed value from the container.
		 * Will only return false if the value did not exist in the container.
		 * A return value of true does not imply that no duplicate values do not still exist in the container.
		 * All elements after the removed one get shifted to the left in memory.
		 * O(n)
		 *
		 * @param t		The value to remove.
		 * @returns		Whether or not a removal operation was performed.
		 */
		bool Remove(const T& t);

		/**
		 * Removes the first element matching the passed predicate.
		 * Will only return false if the value did not exist in the container.
		 * A return value of true does not imply that no duplicate values do not still exist in the container.
		 * All elements after the removed one get shifted to the left in memory.
		 * O(n)
		 *
		 * @param predicate		The predicate to query with.
		 * @returns				Whether or not a removal was performed.
		 */
		template<std::predicate<T> Predicate>
		bool Remove(Predicate predicate);

		/**
		 * Removes the index specified and shift values to the left over.
		 * O(n) where n is Size() - index;
		 *
		 * @param index		The index to remove.
		 */
		void RemoveAt(size_t index);

		/**
		 * Removes all elements matching the specified value.
		 * O(n)
		 *
		 * @param t		The value to remove all instances of.
		 * @returns		How many removals were performed.
		 */
		size_t RemoveAll(const T& t);

		/**
		 * Removes all elements matching the passed predicate.
		 * O(n)
		 *
		 * @param predicate		The predicate to query with.
		 * @returns				How many elements were removed.
		 */
		template<std::predicate<T> Predicate>
		size_t RemoveAll(Predicate predicate);

		template<std::predicate<T> Predicate>
		friend size_t erase_if(Array& arr, const Predicate predicate)
		{
			return arr.RemoveAll(predicate);
		}

		/**
		 * Removes the range [first, last) elements between the two passed indices.
		 * If last is beyond Size() then all elements from first to the end will be removed safely.
		 * Does nothing if the indices are both beyond Size().
		 *
		 * @param first		index of the first value to remove, inclusive.
		 * @param last		index of the last value to remove, exclusive.
		 * @returns			The number of elements removed, equal to first - last unless last > Size().
		 *
		 * @throws std::invalid_argument	if first > last.
		 */
		size_t Remove(size_t first, size_t last);

		/**
		 * Removes the range [first, last) elements between the two passed iterators.
		 * O(n) where n = Size() - first
		 *
		 * @param first		iterator at the first element to remove, inclusive
		 * @param last		iterator at the last element to remove, exclusive
		 * @returns			The number of elements removed.
		 */
		size_t Remove(const_iterator first, const_iterator last);

		/**
		 * O(1)
		 * Removes the last element of the container.
		 * Does nothing if the container is empty.
		 */
		void PopBack();

		/**
		 * O(n) where n = Size() - 1 before the operation.
		 * Removes the first element of the container.
		 * Does nothing if the container is empty.
		 */
		void PopFront();

		/**
		 * Erases all elements from the container by calling the destructors on each of them.
		 * Afterwards Size() == 0.
		 */
		void Clear() noexcept;

		/**
		 * Erases all memory held by this container.
		 */
		void Empty() noexcept;
#pragma endregion

#pragma region Query
		/**
		 * Does an O(n) search for the first occurrence of the passed element.
		 *
		 * @param t				the element to query for.
		 * @returns				the index of the first element found, or Size() if it was not found.
		 */
		size_type IndexOf(const T& t) const;

		/**
		 * Does an O(n) search for the first element matching the passed predicate.
		 *
		 * @param predicate		the predicate to query with.
		 * @returns				the index of the first element found, or Size() if it was not found.
		 */
		template<std::predicate<T> Predicate>
		size_type IndexOf(Predicate predicate) const;
#pragma endregion

#pragma region Memory
		/**
		 * Warning! Use at your own risk! Manipulating this memory can invalidate the state of the Array.
		 *
		 * @returns		pointer to the underlying data. nullptr if the container is empty.
		 */
		[[nodiscard]] T* Data() noexcept;

		/**
		 * Warning! Use at your own risk! Manipulating this memory can invalidate the state of the Array.
		 *
		 * @returns		pointer to the underlying data. nullptr if the container is empty.
		 */
		[[nodiscard]] const T* Data() const noexcept;

		/**
		 * Gets the underlying data and empties this container without deallocating any memory.
		 * Guaranteed to cause a memory leak if the caller discards this return value.
		 * Caller assumes all responsibility for managing this data.
		 *
		 * @returns		pointer, size, capacity
		 */
		[[nodiscard]]
		std::tuple<gsl::owner<T*>, size_t, size_t> TakeData() noexcept;

		/**
		 * Force this container to take on the given memory.
		 * Will free all currently associated memory and destruct all objects.
		 * The container will manage the given data as if it had allocated it itself,
		 * possibly reallocating it in the event of an operation that would do such a thing.
		 *
		 * @param data	pointer, size, capacity
		 */
		void SetData(std::tuple<T*, size_t, size_t> data) noexcept;

		/**
		 * Force this container to take on the given memory.
		 * Will free all currently associated memory and destruct all objects.
		 * The container will manage the given data as if it had allocated it itself,
		 * possibly reallocating it in the event of an operation that would do such a thing.
		 *
		 * @param array		pointer to array of T
		 * @param size		how many elements array contains
		 * @param capacity	how many elements array can store
		 */
		void SetData(T* array, size_t size, size_t capacity) noexcept;

		/**
		 * O(n) where n is the current size of the container.
		 * Does nothing if newCapacity <= Capacity().
		 *
		 * @param newCapacity	The new capacity for the container.
		 *
		 * @throws std::length_error	if newCapacity > MaxSize()
		 * @throws std::bad_alloc		if memory allocation fails.
		 */
		void Reserve(size_t newCapacity);

		/**
		 * O(n) where n = newSize - Size(), or n = newSize if the operation requires a Reserve().
		 * If the current Size() is greater than newSize, the container is reduced to its first newSize elements.
		 *
		 * @param newSize		The new size for the container.
		 * @param prototype		The value to initialize each new element with.
		 */
		void Resize(size_t newSize, const T& prototype = T());

		/**
		 * Reduces memory usage such that Size() == Capacity().
		 * O(n) where n = Size().
		 * Does nothing if IsFull()
		 */
		void ShrinkToFit();

		/**
		 * Reduces memory usage such that Capacity() == std::max(count, Size()).
		 * O(n) where n = Size().
		 * Does nothing if Capacity() already == std::max(count, Size()).
		 */
		void ShrinkToFit(size_t count);

		/**
		 * TODO: It might be that the STL expects swap to be lowercase and a non-member function (so a friend function will do).
		 * O(1)
		 *
		 * @param <OtherReserveStrategy>	template parameter allowing swapping Arrays with different ReserveStrategies.
		 * @param other						container to exchange contents with.
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		void Swap(Array<T, OtherReserveStrategy>& other) noexcept;

		/**
		 * Reverses the order of the elements in the container.
		 * O(n)
		 */
		void Reverse() noexcept;
#pragma endregion

#pragma region Operators
		/**
		 * templated operator== to allow for comparing containers with different reserve strategies
		 *
		 * @param <OtherReserveStrategy>	reserve strategy for the other container
		 * @param left						lhs container
		 * @param right						rhs container
		 * @returns							whether or not the two containers are equal.
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		[[nodiscard]] friend bool operator==(const Array& left, const Array<T, OtherReserveStrategy>& right)
		{
			return &left == reinterpret_cast<const Array*>(&right) || left.Size() == right.Size() && std::equal(left.begin(), left.end(), right.begin(), right.end());
		}

		/**
		 * templated operator== to allow for comparing containers with different reserve strategies
		 *
		 * @param <OtherReserveStrategy>	reserve strategy for the other container
		 * @param left						lhs container
		 * @param right						rhs container
		 * @returns							whether or not the two containers are equal.
		 */
		template<Concept::ReserveStrategy OtherReserveStrategy>
		[[nodiscard]] friend bool operator!=(const Array& left, const Array<T, OtherReserveStrategy>& right)
		{
			return !operator==(left, right);
		}

		friend std::ostream& operator<<(std::ostream& stream, const ARRAY& array) noexcept
		{
			Util::StreamTo(stream, array.begin(), array.end());
			return stream;
		}
#pragma endregion

#pragma region Helpers
	private:
		/**
		 * assumes Capacity() is sufficient
		 *
		 * @param first		the beginning iterator to read from (inclusive)
		 * @param last		the ending iterator to read from (exclusive)
		 */
		template<std::forward_iterator It>
		void PushBackNoResize(It first, It last);

		/**
		 * Helper for methods that accept iterators.
		 *
		 * @param it		The iterator who's owner to compare against this.
		 *
		 * @asserts			owner == this
		 */
		void AssertOwner(iterator it) const;

		/**
		 * @param it		The iterator who's owner to compare against this.
		 *
		 * @asserts			owner == this
		 */
		void AssertOwner(const_iterator it) const;

		/**
		 * Calls the ReserveStrategy template argument and sanity checks its return value.
		 * @returns		A safe new Capacity() for this container.
		 */
		[[nodiscard]] size_type InvokeReserveStrategy() const noexcept;

		/**
		 * Calls destructor on all elements.
		 */
		void DestructAll() noexcept;

		/**
		 * Calls the destructor on all elements in range [first, last).
		 *
		 * @param first		Index of first element to be deconstructed, inclusive.
		 * @param last		Index of last element to be decondstructed, exclusive.
		 */
		void DestructAll(const size_type first, const size_type last) noexcept;

		/**
		 * Shifts the internal array to the right starting at the passed index.
		 * Moves memory as if calling the destructor on elements before they're overwritten and after they're moved.
		 *
		 * @param startIndex	The index to start the shifts at.
		 * @param shiftAmount	How much to shift by.
		 */
		void ShiftRight(size_type startIndex = 0, size_type shiftAmount = 1) noexcept;

		/**
		 * Shifts the internal array to the left starting at the passed index.
		 * Moves memory as if calling the destructor on elements before they're overwritten and after they're moved.
		 *
		 * @param startIndex	The index to start the shifts at.
		 * @param shiftAmount	How much to shift by.
		 */
		void ShiftLeft(size_type startIndex = 0, size_type shiftAmount = 1) noexcept;

		/**
		 * Helper called from constructors.
		 * Assumes size is 0 and capacity is sufficient.
		 *
		 * @param list		list to copy elements from.
		 */
		void Copy(std::initializer_list<T> list);

		/**
		 * Helper called from constructors.
		 * Assumes size is 0 and capacity is sufficient.
		 *
		 * @param list		list to copy elements from.
		 */
		template<typename OtherReserveStrategy>
		void Copy(const Array<T, OtherReserveStrategy>& other);

		/**
		 * Helper called from copy assignment operators.
		 *
		 * @param other		The array being copied.
		 * @returns			this after the copy.
		 */
		template<typename OtherReserveStrategy>
		[[nodiscard]] Array& CopyAssign(const Array<T, OtherReserveStrategy>& other);

		/**
		 * Helper called from move assignment operators.
		 *
		 * @param other		The array being moved.
		 * @returns			this after the move.
		 */
		template<typename OtherReserveStrategy>
		[[nodiscard]] Array& MoveAssign(Array<T, OtherReserveStrategy>&& other);

		/**
		 * Helper called from constructors.
		 * Sets this Array's state after std::move.
		 */
		void SetPostMoveState();
#pragma endregion
	};
}

#include "Array.inl"

ext: h

max_stars_repo_name: JarrettWendt/FIEAEngine

max_issues_repo_name: JarrettWendt/FIEAEngine

max_forks_repo_name: JarrettWendt/FIEAEngine


text:
/* multifit/fsolver.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU 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 General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multifit_nlin.h>

gsl_multifit_fsolver *
gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T,
                            size_t n, size_t p)
{
  int status;

  gsl_multifit_fsolver * s;

  if (n < p)
    {
      GSL_ERROR_VAL ("insufficient data points, n < p", GSL_EINVAL, 0);
    }

  s = (gsl_multifit_fsolver *) malloc (sizeof (gsl_multifit_fsolver));

  if (s == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for multifit solver struct",
                        GSL_ENOMEM, 0);
    }

  s->x = gsl_vector_calloc (p);

  if (s->x == 0)
    {
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
    }

  s->f = gsl_vector_calloc (n);

  if (s->f == 0)
    {
      gsl_vector_free (s->x);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0);
    }

  s->dx = gsl_vector_calloc (p);

  if (s->dx == 0)
    {
      gsl_vector_free (s->x);
      gsl_vector_free (s->f);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
    }

  s->state = malloc (T->size);

  if (s->state == 0)
    {
      gsl_vector_free (s->dx);
      gsl_vector_free (s->x);
      gsl_vector_free (s->f);
      free (s);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for multifit solver state",
                        GSL_ENOMEM, 0);
    }

  s->type = T ;

  status = (s->type->alloc)(s->state, n, p);

  if (status != GSL_SUCCESS)
    {
      (s->type->free)(s->state);
      free (s->state);
      gsl_vector_free (s->dx);
      gsl_vector_free (s->x);
      gsl_vector_free (s->f);
      free (s);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to set solver", status, 0);
    }

  s->function = NULL;

  return s;
}

int
gsl_multifit_fsolver_set (gsl_multifit_fsolver * s,
                          gsl_multifit_function * f,
                          const gsl_vector * x)
{
  if (s->f->size != f->n)
    {
      GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
    }

  if (s->x->size != x->size)
    {
      GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
    }

  s->function = f;
  gsl_vector_memcpy(s->x,x);

  return (s->type->set) (s->state, s->function, s->x, s->f, s->dx);
}

int
gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * s)
{
  return (s->type->iterate) (s->state, s->function, s->x, s->f, s->dx);
}

void
gsl_multifit_fsolver_free (gsl_multifit_fsolver * s)
{
  RETURN_IF_NULL (s);
  (s->type->free) (s->state);
  free (s->state);
  gsl_vector_free (s->dx);
  gsl_vector_free (s->x);
  gsl_vector_free (s->f);
  free (s);
}

const char *
gsl_multifit_fsolver_name (const gsl_multifit_fsolver * s)
{
  return s->type->name;
}

gsl_vector *
gsl_multifit_fsolver_position (const gsl_multifit_fsolver * s)
{
  return s->x;
}

ext: c


text:
/* movstat/sumacc.c
 *
 * Moving window sum accumulator
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU 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 General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_movstat.h>

typedef double ringbuf_type_t;

#include "ringbuf.c"

typedef struct
{
  double sum;     /* current window sum */
  ringbuf *rbuf;  /* ring buffer storing current window */
} sumacc_state_t;

static size_t
sumacc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(sumacc_state_t);
  size += ringbuf_size(n);

  return size;
}

static int
sumacc_init(const size_t n, void * vstate)
{
  sumacc_state_t * state = (sumacc_state_t *) vstate;

  state->sum = 0.0;

  state->rbuf = (ringbuf *) ((unsigned char *) vstate + sizeof(sumacc_state_t));
  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
sumacc_insert(const double x, void * vstate)
{
  sumacc_state_t * state = (sumacc_state_t *) vstate;

  if (ringbuf_is_full(state->rbuf))
    {
      /* subtract oldest element from sum */
      state->sum -= ringbuf_peek_back(state->rbuf);
    }

  /* add new element to sum and ring buffer */
  state->sum += x;
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
sumacc_delete(void * vstate)
{
  sumacc_state_t * state = (sumacc_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    {
      state->sum -= ringbuf_peek_back(state->rbuf);
      ringbuf_pop_back(state->rbuf);
    }

  return GSL_SUCCESS;
}

static int
sumacc_get(void * params, double * result, const void * vstate)
{
  const sumacc_state_t * state = (const sumacc_state_t *) vstate;
  (void) params;
  *result = state->sum;
  return GSL_SUCCESS;
}

static const gsl_movstat_accum sum_accum_type =
{
  sumacc_size,
  sumacc_init,
  sumacc_insert,
  sumacc_delete,
  sumacc_get
};

const gsl_movstat_accum *gsl_movstat_accum_sum = &sum_accum_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
static const char help[] =
"Solves time-dependent nonlinear ice sheet problem in 2D:\n"
"(*)    H_t + div q = m\n"
"where q = (q^x,q^y) is the nonsliding shallow ice approximation (SIA) flux,\n"
"       q = - Gamma H^{n+2} |grad s|^{n-1} grad s\n"
"always subject to the constraint\n"
"       H(t,x,y) >= 0.\n"
"In these equations  H(t,x,y)  is ice thickness,  b(x,y)  is bed elevation,\n"
"s(t,x,y) = H(t,x,y) + b(x,y)  is surface elevation,  and  m(x,y)  is the\n"
"climatic mass balance.  Constants are  n > 1  and  Gamma = 2 A (rho g)^n / (n+2).\n"
"The domain is square,  Omega = [0,L] x [0,L],  with periodic boundary conditions.\n"
"Equation (*) is semi-discretized in space by a Q1 structured-grid FVE method\n"
"(Bueler, 2016).  The resulting method-of-lines ODE in time is in the form\n"
"      F(t,H,H_t) = G(t,H)\n"
"For this ODE we define three callbacks:\n"
"      FormIFunction()   evaluates F(H,H_t) = H_t + div q\n"
"      FormRHSFunction() evaluates G(t,H) = m\n"
"      FormIJacobian()   evaluates (shift) dF/dH_t + dF/dH.\n"
"Options -snes_fd_color works well but is slower than the default Jacobian\n"
"by a factor of two or so.  Requires SNESVI types (i.e.\n"
"-snes_type vinewton{rsls|ssls}) because of constraint.  With current PETSc\n"
"design, explicit TS types do not work.\n\n";

/* try:

./ice -da_refine 3                  # only meaningful at this res and higher

./ice -snes_atol 1.0e-14            # CRITICAL; some lower bound on SNES norm is
                                    #   critical to allowing convergence near steady state

./ice                               # DEFAULT uses analytical jacobian
./ice -snes_fd_color
./ice -snes_type vinewtonrsls       # DEFAULT
./ice -snes_type vinewtonssls       # slightly more robust?
(other -snes_types like newtonls not allowed because don't support bounds)

./ice -ts_type arkimex              # DEFAULT
./ice -ts_type beuler               # robust
./ice -ts_type cn                   # mis-behaves on long time steps
./ice -ts_type cn -ts_theta_adapt                         # good
./ice -ts_type theta -ts_theta_adapt                      # good
./ice -ts_type bdf -ts_bdf_adapt -ts_bdf_order 2|3|4|5|6  # good

# time-stepping control options
-ts_adapt_type basic                # DEFAULT
-ts_adapt_basic_clip 0.5,1.2        # DEFAULT and recommended;
                                      doesn't lengthen too much, but allows
                                      significantly shorter, in response to estimate of
                                      local truncation error (vs default: 0.1,10.0)
-ts_max_snes_failures -1            # recommended: do retry solve
-ts_adapt_scale_solve_failed 0.9    # recommended: try a slightly-easier problem
-ts_max_reject 50                   # recommended?:  keep trying if lte is too big

./ice -ice_dtlimits                 # info on comparison to explicit
./ice -ts_adapt_monitor             # info on adapt

# PC possibilities
-ksp_type gmres                                               # DEFAULT
-pc_type ilu                                                  # DEFAULT
-pc_type gamg -pc_gamg_threshold 0.0 -pc_gamg_agg_nsmooths 1  # these are GAMG defaults
-pc_type gamg -pc_gamg_threshold 0.2 -pc_gamg_agg_nsmooths 1  # a little faster?
-pc_type lu
-pc_type asm -sub_pc_type lu
-pc_type asm -sub_pc_type ilu
-pc_type mg
-pc_type mg -pc_mg_levels 4 -mg_levels_ksp_monitor

# shows nontriviality converging ice caps on mountains:
mpiexec -n 2 ./ice -da_refine 5 -ts_monitor_solution draw -snes_converged_reason -ice_tf 10000.0 -ice_dtinit 100.0 -ts_max_snes_failures -1 -ts_adapt_scale_solve_failed 0.9

# start with short time step and it will find good time scale
./ice -snes_converged_reason -da_refine 4 -ice_dtinit 0.1

# recovery from convergence failures works!  (failure here triggered by n=4):
./ice -da_refine 3 -ice_n 4 -ts_max_snes_failures -1 -snes_converged_reason

verif with dome=1, halfar=2:
for TEST in 1 2; do
    for N in 2 3 4 5 6; do
        mpiexec -n 2 ./ice -ice_monitor 0 -ice_verif $TEST -ice_eps 0.0 -ice_dtinit 50.0 -ice_tf 2000.0 -da_refine $N -ts_type beuler
    done
done

actual test B from Bueler et al 2005:
./ice -ice_verif 2 -ice_eps 0 -ice_dtinit 100 -ice_tf 25000 -ice_L 2200e3 -da_refine $N

recommended "new PISM":
mpiexec -n N ./ice -da_refine M \
   -snes_type vinewtonrsls \
   -ts_type arkimex \    #(OR -ts_type bdf -ts_bdf_adapt -ts_bdf_order 4)
   -ts_adapt_type basic -ts_adapt_basic_clip 0.5,1.2 \
   -ts_max_snes_failures -1 -ts_adapt_scale_solve_failed 0.9 \
   -pc_type gamg

run to generate final-time result in file ice_192_50000.dat:
mpiexec -n 4 ./ice -da_refine 6 -pc_type mg -pc_mg_levels 4 -ice_dtlimits -ice_tf 50000 -ice_dtinit 1.0 -ts_max_snes_failures -1 -ts_adapt_scale_solve_failed 0.9 -ice_dump
*/

#include <petsc.h>
#include "icecmb.h"

// context is entirely grid-independent info
typedef struct {
    double    secpera,// number of seconds in a year
              L,      // spatial domain is [0,L] x [0,L]
              tf,     // final time; time domain is [0,tf]
              dtinit, // user-requested initial time step
              dtmax,  // set TS maximum time step
              g,      // acceleration of gravity
              rho_ice,// ice density
              n_ice,  // Glen exponent for SIA flux term
              A_ice,  // ice softness
              Gamma,  // coefficient for SIA flux term
              D0,     // representative value of diffusivity (used in regularizing D)
              eps,    // regularization parameter for D
              delta,  // dimensionless regularization for slope in SIA formulas
              lambda, // amount of upwinding; lambda=0 is none and lambda=1 is "full"
              locmaxD,// maximum of diffusivity from last residual evaluation
              dtexplicitsum;// running sum of explicit dt limit
    int       verif;  // 0 = not verification, 1 = dome, 2 = Halfar (1983)
    PetscBool monitor,// use -ice_monitor
              monitor_dt_limits,// also monitor time step limits for explicit schemes
              dump;   // dump state (H,b) at final time
    CMBModel  *cmb;// defined in cmbmodel.h
} AppCtx;

#include "iceverif.h"

extern PetscErrorCode SetFromOptionsAppCtx(AppCtx*);
extern PetscErrorCode IceMonitor(TS, int, double, Vec, void*);
extern PetscErrorCode ExplicitLimitsMonitor(TS, int, double, Vec, void*);
extern PetscErrorCode FormBedLocal(DMDALocalInfo*, int, double**, AppCtx*);
extern PetscErrorCode FormBounds(SNES,Vec,Vec);
extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo*, double,
                          double**, double**, double**, AppCtx*);
extern PetscErrorCode FormIJacobianLocal(DMDALocalInfo*, double,
                          double**, double**, double, Mat, Mat, AppCtx *user);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, double,
                          double**, double**, AppCtx*);

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  DM             da;
  TS             ts;
  SNES           snes;   // no need to destroy (owned by TS)
  TSAdapt        adapt;
  Vec            H;
  AppCtx         user;
  CMBModel       cmb;
  DMDALocalInfo  info;
  double         dx,dy,**aH;

  PetscInitialize(&argc,&argv,(char*)0,help);

  ierr = SetFromOptionsAppCtx(&user); CHKERRQ(ierr);
  ierr = SetFromOptions_CMBModel(&cmb,user.secpera);
  user.cmb = &cmb;

  // this DMDA is the cell-centered grid
  ierr = DMDACreate2d(PETSC_COMM_WORLD,
                      DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,
                      DMDA_STENCIL_BOX,
                      3,3,PETSC_DECIDE,PETSC_DECIDE,
                      1, 1,        // dof=1, stencilwidth=1
                      NULL,NULL,&da);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);  // this must be called BEFORE SetUniformCoordinates
  ierr = DMSetApplicationContext(da, &user);CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(da, 0.0, user.L, 0.0, user.L, 0.0,1.0); CHKERRQ(ierr);

  // initialize the TS and configure its time-axis, including its TSAdapt
  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetType(ts,TSARKIMEX); CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,user.tf); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,user.dtinit); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSGetAdapt(ts,&adapt); CHKERRQ(ierr);
  ierr = TSAdaptSetType(adapt,TSADAPTBASIC); CHKERRQ(ierr);
  ierr = TSAdaptSetClip(adapt,0.5,1.2); CHKERRQ(ierr);
  ierr = TSAdaptSetStepLimits(adapt,0.0,user.dtmax); CHKERRQ(ierr);
  if (user.monitor) {
      ierr = TSMonitorSet(ts,IceMonitor,&user,NULL); CHKERRQ(ierr);
  }
  if (user.monitor_dt_limits) {
      ierr = TSMonitorSet(ts,ExplicitLimitsMonitor,&user,NULL); CHKERRQ(ierr);
  }

  // set methods to compute parts of ODE system, for the DMDA grid, for TS callback
  ierr = TSSetDM(ts,da); CHKERRQ(ierr);
  ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,
           (DMDATSIFunctionLocal)FormIFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDATSSetIJacobianLocal(da,
           (DMDATSIJacobianLocal)FormIJacobianLocal,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);

  // configure the SNES and DM to solve a NCP/VI at each step
  ierr = TSGetSNES(ts,&snes); CHKERRQ(ierr);
  ierr = SNESSetType(snes,SNESVINEWTONRSLS); CHKERRQ(ierr);
  ierr = SNESVISetComputeVariableBounds(snes,&FormBounds); CHKERRQ(ierr);

  // should be done setting up TS
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);

  // report on space-time grid
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  dx = user.L / (double)(info.mx);
  dy = user.L / (double)(info.my);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
     "solving on domain [0,L] x [0,L] (L=%.3f km) and time interval [0,tf] (tf=%.3f a)\n"
     "grid: %d x %d points, spacing dx=%.3f km x dy=%.3f km, dtinit=%.3f a\n",
     user.L/1000.0,user.tf/user.secpera,
     info.mx,info.my,dx/1000.0,dy/1000.0,user.dtinit/user.secpera);

  // set up initial condition on fine grid
  ierr = DMCreateGlobalVector(da,&H);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)H,"H"); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,H,&aH); CHKERRQ(ierr);
  if (user.verif == 1) {
      ierr = DomeThicknessLocal(&info,aH,&user); CHKERRQ(ierr);
  } else if (user.verif == 2) {
      double t0;
      ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
      ierr = HalfarThicknessLocal(&info,t0,aH,&user); CHKERRQ(ierr);
  } else {
      // fill H according to chop-scale-CMB
      ierr = FormBedLocal(&info,0,aH,&user); CHKERRQ(ierr);  // H(x,y) <- b(x,y)
      ierr = ChopScaleInitialHLocal_CMBModel(&cmb,&info,aH,aH); CHKERRQ(ierr);
  }
  ierr = DMDAVecRestoreArray(da,H,&aH); CHKERRQ(ierr);

  // solve
  ierr = TSSolve(ts,H); CHKERRQ(ierr);

  // time-stepping summary if -ice_dtlimits
  if (user.monitor_dt_limits) {
      int count;
      ierr = TSGetStepNumber(ts,&count); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,
          "average dt %.5f a, average dtexplicit %.5f a\n",
          (user.tf/user.secpera)/(double)count,
          (user.dtexplicitsum/user.secpera)/(double)count); CHKERRQ(ierr);
  }

  // dump state if requested
  if (user.dump) {
      char           filename[1024];
      PetscViewer    viewer;
      Vec            b;
      double         **ab;
      ierr = VecDuplicate(H,&b);CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject)b,"b"); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(da,b,&ab); CHKERRQ(ierr);
      ierr = FormBedLocal(&info,0,ab,&user); CHKERRQ(ierr);
      ierr = DMDAVecRestoreArray(da,b,&ab); CHKERRQ(ierr);
      ierr = sprintf(filename,"ice_%d_%d.dat",info.mx,(int)(user.tf/user.secpera));
      ierr = PetscPrintf(PETSC_COMM_WORLD,"writing PETSC binary file %s ...\n",filename); CHKERRQ(ierr);
      ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
      ierr = VecView(b,viewer); CHKERRQ(ierr);
      ierr = VecView(H,viewer); CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
      VecDestroy(&b);
  }

  // compute error in verification case
  if (user.verif > 0) {
      Vec Hexact;
      double infnorm, onenorm;
      ierr = VecDuplicate(H,&Hexact); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(da,Hexact,&aH); CHKERRQ(ierr);
      if (user.verif == 1) {
          ierr = DomeThicknessLocal(&info,aH,&user); CHKERRQ(ierr);
      } else if (user.verif == 2) {
          double tf;
          ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
          ierr = HalfarThicknessLocal(&info,tf,aH,&user); CHKERRQ(ierr);
      } else {
          SETERRQ(PETSC_COMM_WORLD,3,"invalid user.verif ... how did I get here?\n");
      }
      ierr = DMDAVecRestoreArray(da,Hexact,&aH); CHKERRQ(ierr);
      ierr = VecAXPY(H,-1.0,Hexact); CHKERRQ(ierr);    // H <- H + (-1.0) Hexact
      VecDestroy(&Hexact);
      ierr = VecNorm(H,NORM_INFINITY,&infnorm); CHKERRQ(ierr);
      ierr = VecNorm(H,NORM_1,&onenorm); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,
          "errors on verif %d: |H-Hexact|_inf = %.3f, |H-Hexact|_average = %.3f\n",
          user.verif,infnorm,onenorm/(double)(info.mx*info.my)); CHKERRQ(ierr);
  }

  // clean up
  VecDestroy(&H);  TSDestroy(&ts);  DMDestroy(&da);
  return PetscFinalize();
}

PetscErrorCode SetFromOptionsAppCtx(AppCtx *user) {
  PetscErrorCode ierr;
  PetscBool      set;

  user->secpera= 31556926.0;  // number of seconds in a year
  user->L      = 1800.0e3;    // m; note  domeL=750.0e3 is radius of verification ice sheet
  user->tf     = 100.0 * user->secpera;  // default to 100 years
  user->dtinit = 10.0 * user->secpera;   // default to 10 year as initial step
  user->dtmax  = 1.0e6 * user->secpera;   // default to million years; huge
  user->g      = 9.81;        // m/s^2
  user->rho_ice= 910.0;       // kg/m^3
  user->n_ice  = 3.0;
  user->A_ice  = 3.1689e-24;  // 1/(Pa^3 s); EISMINT I value
  user->D0     = 1.0;         // m^2 / s
  user->eps    = 0.001;
  user->delta  = 1.0e-4;
  user->lambda = 0.25;
  user->dtexplicitsum = 0.0;
  user->verif  = 0;
  user->monitor = PETSC_TRUE;
  user->monitor_dt_limits = PETSC_FALSE;
  user->dump   = PETSC_FALSE;
  user->cmb    = NULL;

  PetscFunctionBeginUser;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"ice_","options to ice","");CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-A", "set value of ice softness A in units Pa-3 s-1",
      "ice.c",user->A_ice,&user->A_ice,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-D0", "representative value of diffusivity (used in regularizing D) in units m2 s-1",
      "ice.c",user->D0,&user->D0,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-delta", "dimensionless regularization for slope in SIA formulas",
      "ice.c",user->delta,&user->delta,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-dt_init", "initial time step; input units are years",
      "ice.c",user->dtinit,&user->dtinit,&set);CHKERRQ(ierr);
  if (set)   user->dtinit *= user->secpera;
  ierr = PetscOptionsReal(
      "-dt_max", "maximum time step; input units are years",
      "ice.c",user->dtmax,&user->dtmax,&set);CHKERRQ(ierr);
  if (set)   user->dtmax *= user->secpera;
  ierr = PetscOptionsBool(
      "-dump", "save final state (H, b)",
      "ice.c",user->dump,&user->dump,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-eps", "dimensionless regularization for diffusivity D",
      "ice.c",user->eps,&user->eps,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-L", "side length of domain in meters",
      "ice.c",user->L,&user->L,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-lambda", "amount of upwinding; lambda=0 is none and lambda=1 is full",
      "ice.c",user->lambda,&user->lambda,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool(
      "-monitor", "use the ice monitor which shows ice sheet volume and area",
      "ice.c",user->monitor,&user->monitor,&set);CHKERRQ(ierr);
  if (!set)   user->monitor = PETSC_TRUE;
  ierr = PetscOptionsBool(
      "-monitor_dt_limits", "monitor the time-step limits which would apply to an explicit scheme",
      "ice.c",user->monitor_dt_limits,&user->monitor_dt_limits,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-n", "value of Glen exponent n",
      "ice.c",user->n_ice,&user->n_ice,NULL);CHKERRQ(ierr);
  if (user->n_ice <= 1.0) {
      SETERRQ1(PETSC_COMM_WORLD,1,
          "ERROR: n = %f not allowed ... n > 1 is required\n",user->n_ice);
  }
  ierr = PetscOptionsReal(
      "-rho", "ice density in units kg m3",
      "ice.c",user->rho_ice,&user->rho_ice,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-tf", "final time in seconds; input units are years",
      "ice.c",user->tf,&user->tf,&set);CHKERRQ(ierr);
  if (set)   user->tf *= user->secpera;
  ierr = PetscOptionsInt(
      "-verif","1 = dome exact solution; 2 = halfar exact solution",
      "ice.c",user->verif,&(user->verif),&set);CHKERRQ(ierr);
  if ((set) && ((user->verif < 0) || (user->verif > 2))) {
      SETERRQ1(PETSC_COMM_WORLD,2,
          "ERROR: verif = %d not allowed ... 0 <= verif <= 2 is required\n",
          user->verif);
  }
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  // derived constant computed after other ice properties are set
  user->Gamma = 2.0 * PetscPowReal(user->rho_ice*user->g,user->n_ice)
                    * user->A_ice / (user->n_ice+2.0);

  PetscFunctionReturn(0);
}

// this basic monitor gives current time, volume, area
PetscErrorCode IceMonitor(TS ts, int step, double time, Vec H, void *ctx) {
    PetscErrorCode ierr;
    AppCtx         *user = (AppCtx*)ctx;
    double         lvol = 0.0, vol, larea = 0.0, area, darea, **aH;
    int            j, k;
    MPI_Comm       com;
    DM             da;
    DMDALocalInfo  info;

    PetscFunctionBeginUser;
    ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    darea = user->L * user->L / (double)(info.mx * info.my);
    ierr = DMDAVecGetArrayRead(da,H,&aH); CHKERRQ(ierr);
    for (k = info.ys; k < info.ys + info.ym; k++) {
        for (j = info.xs; j < info.xs + info.xm; j++) {
            if (aH[k][j] > 1.0) {  // for volume/area its helpful to not count tinys
                larea += darea;
                lvol += aH[k][j];
            }
        }
    }
    ierr = DMDAVecRestoreArrayRead(da,H,&aH); CHKERRQ(ierr);
    lvol *= darea;
    ierr = PetscObjectGetComm((PetscObject)(da),&com); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&lvol,&vol,1,MPI_DOUBLE,MPI_SUM,com); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&larea,&area,1,MPI_DOUBLE,MPI_SUM,com); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
        "%3d: time %.3f a,  volume %.1f 10^3 km^3,  area %.1f 10^3 km^2\n",
        step,time/user->secpera,vol/1.0e12,area/1.0e9); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

// this monitor reports on max diffusivity, velocity, explicit time-step limits (for comparison)
PetscErrorCode ExplicitLimitsMonitor(TS ts, int step, double time, Vec H, void *ctx) {
    PetscErrorCode ierr;
    AppCtx         *user = (AppCtx*)ctx;
    double         maxD, dd, dtD=PETSC_INFINITY;
    MPI_Comm       com;
    DM             da;
    DMDALocalInfo  info;

    PetscFunctionBeginUser;
    if (time <= 0.0) {
        PetscFunctionReturn(0);
    }
    // globalize maxD
    ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
    ierr = PetscObjectGetComm((PetscObject)(da),&com); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&(user->locmaxD),&maxD,1,MPI_DOUBLE,MPI_MAX,com); CHKERRQ(ierr);
    // compute explicit limits
    if (maxD <= 0.0) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,"    [NO -ice_monitor_dt_limits output because maxD is zero]\n"); CHKERRQ(ierr);
    } else {
        ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
        dd = PetscMin(user->L / (double)(info.mx), user->L / (double)(info.my));
        dtD = dd * dd / (4.0*maxD);
        ierr = PetscPrintf(PETSC_COMM_WORLD,"    max D_SIA %.3f m2s-1,  dtD %.3e a\n",
            maxD,dtD/user->secpera); CHKERRQ(ierr);
        user->dtexplicitsum += dtD;
    }
    PetscFunctionReturn(0);
}

PetscErrorCode FormBedLocal(DMDALocalInfo *info, int stencilwidth, double **ab, AppCtx *user) {
  int          j,k,r,s;
  const double dx = user->L / (double)(info->mx),
               dy = user->L / (double)(info->my),
               Z = PETSC_PI / user->L;
  double       x, y, b;
  // vaguely-random frequencies and coeffs generated by fiddling; see randbed.py
  const int    nc = 4,
               jc[4] = {1, 3, 6, 8},
               kc[4] = {1, 3, 4, 7};
  const double scalec = 750.0,
               C[4][4] = { { 2.00000000,  0.33000000, -0.55020034,  0.54495520},
                           { 0.50000000,  0.45014486,  0.60551833, -0.52250644},
                           { 0.93812068,  0.32638429, -0.24654812,  0.33887052},
                           { 0.17592361, -0.35496741,  0.22694547, -0.05280704} };
  PetscFunctionBeginUser;
  // go through owned portion of grid and compute  b(x,y)
  for (k = info->ys-stencilwidth; k < info->ys + info->ym+stencilwidth; k++) {
      y = k * dy;
      for (j = info->xs-stencilwidth; j < info->xs + info->xm+stencilwidth; j++) {
          x = j * dx;
          // b(x,y) is sum of a few sines
          b = 0.0;
          for (r = 0; r < nc; r++) {
              for (s = 0; s < nc; s++) {
                  b += C[r][s] * sin(jc[r] * Z * x) * sin(kc[s] * Z * y);
              }
          }
          ab[k][j] = scalec * b;
      }
  }
  PetscFunctionReturn(0);
}

//  for call-back: tell SNESVI (variational inequality) that we want
//    0.0 <= H < +infinity
PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu) {
  PetscErrorCode ierr;
  PetscFunctionBeginUser;
  ierr = VecSet(Xl,0.0); CHKERRQ(ierr);
  ierr = VecSet(Xu,PETSC_INFINITY); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// value of gradient at a point
typedef struct {
    double x,y;
} Grad;

/* We factor the SIA flux as
    q = - H^{n+2} sigma(|grad s|) grad s
where sigma is the slope-dependent part
    sigma(z) = Gamma z^{n-1}.
Also
    D = H^{n+2} sigma(|grad s|)
so that q = - D grad s.  */
static double sigma(Grad gH, Grad gb, const AppCtx *user) {
    const double n = user->n_ice;
    if (n > 1.0) {
        const double sx = gH.x + gb.x,
                     sy = gH.y + gb.y,
                     slopesqr = sx * sx + sy * sy + user->delta * user->delta;
        return user->Gamma * PetscPowReal(slopesqr,(n-1.0)/2);
    } else {
        return user->Gamma;
    }
}

/* Regularized derivative of sigma with respect to nodal value H_l:
   d sigma / dl = (d sigma / d gH.x) * (d gH.x / dl) + (d sigma / d gH.y) * (d gH.y / dl),
but
   d sigma / d gH.x = Gamma * (n-1) * |gH + gb|^{n-3.0} * (gH.x + gb.x)
   d sigma / d gH.y = Gamma * (n-1) * |gH + gb|^{n-3.0} * (gH.y + gb.y)
However, power n-3.0 can be negative, which generates NaN in areas where the
surface gradient gH+gb is zero or tiny, so we add delta^2 to |grad s|^2. */
static double DsigmaDl(Grad gH, Grad gb, Grad dgHdl, const AppCtx *user) {
    const double n = user->n_ice;
    if (n > 1.0) {
        const double sx = gH.x + gb.x,
                     sy = gH.y + gb.y,
                     slopesqr = sx * sx + sy * sy + user->delta * user->delta,
                     tmp = user->Gamma * (n-1) * PetscPowReal(slopesqr,(n-3.0)/2);
        return tmp * sx * dgHdl.x + tmp * sy * dgHdl.y;
    } else {
        return 0.0;
    }
}

/* Pseudo-velocity from bed slope:  W = - sigma * grad b. */
static Grad W(double sigma, Grad gb) {
    Grad W;
    W.x = - sigma * gb.x;
    W.y = - sigma * gb.y;
    return W;
}

/* Derivative of pseudo-velocity W with respect to nodal value H_l. */
static Grad DWDl(double dsigmadl, Grad gb) {
    Grad dWdl;
    dWdl.x = - dsigmadl * gb.x;
    dWdl.y = - dsigmadl * gb.y;
    return dWdl;
}

/* DCS = diffusivity from the continuation scheme:
     D(eps) = (1-eps) sigma H^{n+2} + eps D_0
so D(1)=D_0 and D(0)=sigma H^{n+2}. */
static double DCS(double sigma, double H, const AppCtx *user) {
  return (1.0 - user->eps) * sigma * PetscPowReal(PetscAbsReal(H),user->n_ice+2.0)
         + user->eps * user->D0;
}

/* Derivative of diffusivity D = DCS with respect to nodal value H_l.  Since
  D = (1-eps) sigma H^{n+2} + eps D0
it follows that
  d D / dl = (1-eps) [ (d sigma / dl) H^{n+2} + sigma (n+2) H^{n+1} (d H / dl) ]
           = (1-eps) H^{n+1} [ (d sigma / dl) H + sigma (n+2) (d H / dl) ]    */
static double DDCSDl(double sigma, double dsigmadl, double H, double dHdl,
                     const AppCtx *user) {
    const double Hpow = PetscPowReal(PetscAbsReal(H),user->n_ice+1.0);
    return (1.0 - user->eps) * Hpow * ( dsigmadl * H + sigma * (user->n_ice+2.0) * dHdl );
}

/* Flux component from the non-sliding SIA on a general bed. */
PetscErrorCode SIAflux(Grad gH, Grad gb, double H, double Hup, PetscBool xdir,
                       double *D, double *q, const AppCtx *user) {
  const double mysig = sigma(gH,gb,user),
               myD   = DCS(mysig,H,user);
  const Grad   myW   = W(mysig,gb);
  PetscFunctionBeginUser;
  if (D) {
      *D = myD;
  }
  if (xdir && q) {
      *q = - myD * gH.x + myW.x * PetscPowReal(PetscAbsReal(Hup),user->n_ice+2.0);
  } else {
      *q = - myD * gH.y + myW.y * PetscPowReal(PetscAbsReal(Hup),user->n_ice+2.0);
  }
  PetscFunctionReturn(0);
}

static double DSIAfluxDl(Grad gH, Grad gb, Grad dgHdl,
                         double H, double dHdl, double Hup, double dHupdl,
                         PetscBool xdir, const AppCtx *user) {
    const double Huppow   = PetscPowReal(PetscAbsReal(Hup),user->n_ice+1.0),
                 dHuppow  = (user->n_ice+2.0) * Huppow * dHupdl,
                 mysig    = sigma(gH,gb,user),
                 myD      = DCS(mysig,H,user),
                 dsigmadl = DsigmaDl(gH,gb,dgHdl,user),
                 dDdl     = DDCSDl(mysig,dsigmadl,H,dHdl,user);
    const Grad   myW      = W(mysig,gb),
                 dWdl     = DWDl(dsigmadl,gb);
    if (xdir)
        return - dDdl * gH.x - myD * dgHdl.x + dWdl.x * Huppow * Hup + myW.x * dHuppow;
    else
        return - dDdl * gH.y - myD * dgHdl.y + dWdl.y * Huppow * Hup + myW.y * dHuppow;
}

// gradients of weights for Q^1 interpolant
static const double gx[4] = {-1.0,  1.0, 1.0, -1.0},
                    gy[4] = {-1.0, -1.0, 1.0,  1.0};

static double fieldatpt(double xi, double eta, double f[4]) {
  // weights for Q^1 interpolant
  double x[4] = { 1.0-xi,      xi,  xi, 1.0-xi},
         y[4] = {1.0-eta, 1.0-eta, eta,    eta};
  return   x[0] * y[0] * f[0] + x[1] * y[1] * f[1]
         + x[2] * y[2] * f[2] + x[3] * y[3] * f[3];
}

static double fieldatptArray(int u, int v, double xi, double eta, double **f) {
  double ff[4] = {f[v][u], f[v][u+1], f[v+1][u+1], f[v+1][u]};
  return fieldatpt(xi,eta,ff);
}

static double dfieldatpt(int l, double xi, double eta) {
  const double x[4] = { 1.0-xi,      xi,  xi, 1.0-xi},
               y[4] = {1.0-eta, 1.0-eta, eta,    eta};
  return x[l] * y[l];
}

static Grad gradfatpt(double xi, double eta, double dx, double dy, double f[4]) {
  Grad gradf;
  double x[4] = { 1.0-xi,      xi,  xi, 1.0-xi},
         y[4] = {1.0-eta, 1.0-eta, eta,    eta};
  gradf.x =   gx[0] * y[0] * f[0] + gx[1] * y[1] * f[1]
            + gx[2] * y[2] * f[2] + gx[3] * y[3] * f[3];
  gradf.y =    x[0] *gy[0] * f[0] +  x[1] *gy[1] * f[1]
            +  x[2] *gy[2] * f[2] +  x[3] *gy[3] * f[3];
  gradf.x /= dx;
  gradf.y /= dy;
  return gradf;
}

static Grad gradfatptArray(int u, int v, double xi, double eta, double dx, double dy,
                           double **f) {
  double ff[4] = {f[v][u], f[v][u+1], f[v+1][u+1], f[v+1][u]};
  return gradfatpt(xi,eta,dx,dy,ff);
}

static Grad dgradfatpt(int l, double xi, double eta, double dx, double dy) {
  Grad dgradfdl;
  const double x[4] = { 1.0-xi,      xi,  xi, 1.0-xi},
               y[4] = {1.0-eta, 1.0-eta, eta,    eta},
               gx[4] = {-1.0,  1.0, 1.0, -1.0},
               gy[4] = {-1.0, -1.0, 1.0,  1.0};
  dgradfdl.x = gx[l] *  y[l] / dx;
  dgradfdl.y =  x[l] * gy[l] / dy;
  return dgradfdl;
}

// indexing of the 8 quadrature points along the boundary of the control volume in M*
// point s=0,...,7 is in element (j,k) = (j+je[s],k+ke[s])
static const int  je[8] = {0,  0, -1, -1, -1, -1,  0,  0},
                  ke[8] = {0,  0,  0,  0, -1, -1, -1, -1},
                  ce[8] = {0,  3,  1,  0,  2,  1,  3,  2};

// direction of flux at 4 points in each element
static const PetscBool xdire[4] = {PETSC_TRUE, PETSC_FALSE, PETSC_TRUE, PETSC_FALSE};

// local (element-wise) coords of quadrature points for M*
static const double locx[4] = {  0.5, 0.75,  0.5, 0.25},
                    locy[4] = { 0.25,  0.5, 0.75,  0.5};

/* FormIFunctionLocal  =  IFunction call-back by TS using DMDA info.

Evaluates residual FF on local process patch:
   FF_{j,k} = \int_{\partial V_{j,k}} \mathbf{q} \cdot \mathbf{n}
              - m_{j,k} \Delta x \Delta y
where V_{j,k} is the control volume centered at (x_j,y_k).

Regarding indexing locations along the boundary of the control volume where
flux is evaluated, this figure shows four elements and one control volume
centered at (x_j,y_k).  The boundary of the control volume has 8 points,
numbered s=0,...,7:
   -------------------
  |         |         |
  |    ..2..|..1..    |
  |   3:    |    :0   |
k |--------- ---------|
  |   4:    |    :7   |
  |    ..5..|..6..    |
  |         |         |
   -------------------
            j

Regarding flux-component indexing on the element indexed by (j,k) node,
the value  (aqquad[c])[k][j] for c=0,1,2,3 is an x-component at "*" and
a y-component at "%"; note (x_j,y_k) is lower-left corner:
   -------------------
  |         :         |
  |         *2        |
  |    3    :    1    |
  |....%.... ....%....|
  |         :         |
  |         *0        |
  |         :         |
  @-------------------
(j,k)
*/
PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, double t,
                                  double **aH, double **aHdot, double **FF,
                                  AppCtx *user) {
  PetscErrorCode  ierr;
  const double    dx = user->L / (double)(info->mx),
                  dy = user->L / (double)(info->my);
  // coefficients of quadrature evaluations along the boundary of the control volume in M*
  const double    coeff[8] = {dy/2, dx/2, dx/2, -dy/2, -dy/2, -dx/2, -dx/2, dy/2};
  const PetscBool upwind = (user->lambda > 0.0);
  const double    upmin = (1.0 - user->lambda) * 0.5,
                  upmax = (1.0 + user->lambda) * 0.5;
  int             c, j, k, s;
  double          H, Hup, lxup, lyup, **aqquad[4], **ab, DSIA_ckj, qSIA_ckj;
  Grad            gH, gb;
  Vec             qquad[4], b;

  PetscFunctionBeginUser;

#if 0
  // optionally check admissibility
  for (k = info->ys; k < info->ys + info->ym; k++) {
      for (j = info->xs; j < info->xs + info->xm; j++) {
          if (aH[k][j] < 0.0) {
              SETERRQ3(PETSC_COMM_WORLD,1,"ERROR: non-admissible H[k][j] = %.3e < 0.0 detected at j,k = %d,%d ... stopping\n",aH[k][j],j,k);
          }
      }
  }
#endif

  user->locmaxD = 0.0;
  ierr = DMGetLocalVector(info->da, &b); CHKERRQ(ierr);
  if (user->verif > 0) {
      ierr = VecSet(b,0.0); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(info->da,b,&ab); CHKERRQ(ierr);
  } else {
      ierr = DMDAVecGetArray(info->da,b,&ab); CHKERRQ(ierr);
      ierr = FormBedLocal(info,1,ab,user); CHKERRQ(ierr);  // get stencil width
  }
  for (c = 0; c < 4; c++) {
      ierr = DMGetLocalVector(info->da, &(qquad[c])); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(info->da,qquad[c],&(aqquad[c])); CHKERRQ(ierr);
  }

  // loop over locally-owned elements, including ghosts, to get fluxes q at
  // c = 0,1,2,3 points in element;  note start at (xs-1,ys-1)
  for (k = info->ys-1; k < info->ys + info->ym; k++) {
      for (j = info->xs-1; j < info->xs + info->xm; j++) {
          for (c=0; c<4; c++) {
              H  = fieldatptArray(j,k,locx[c],locy[c],aH);
              gH = gradfatptArray(j,k,locx[c],locy[c],dx,dy,aH);
              gb = gradfatptArray(j,k,locx[c],locy[c],dx,dy,ab);
              if (upwind) {
                  if (xdire[c] == PETSC_TRUE) {
                      lxup = (gb.x <= 0.0) ? upmin : upmax;
                      lyup = locy[c];
                  } else {
                      lxup = locx[c];
                      lyup = (gb.y <= 0.0) ? upmin : upmax;
                  }
                  Hup = fieldatptArray(j,k,lxup,lyup,aH);
              } else
                  Hup = H;
              ierr = SIAflux(gH,gb,H,Hup,xdire[c],
                             &DSIA_ckj,&qSIA_ckj,user); CHKERRQ(ierr);
              aqquad[c][k][j] = qSIA_ckj;
              if (user->monitor_dt_limits) {
                  user->locmaxD = PetscMax(user->locmaxD,DSIA_ckj);
              }
          }
      }
  }

  // loop over nodes, not including ghosts, to get function F(t,H,H') from quadature over
  // s = 0,1,...,7 points on boundary of control volume (rectangle) around node
  for (k=info->ys; k<info->ys+info->ym; k++) {
      for (j=info->xs; j<info->xs+info->xm; j++) {
          FF[k][j] = aHdot[k][j];
          // now add integral over control volume boundary using two
          // quadrature points on each side
          for (s=0; s<8; s++)
              FF[k][j] += coeff[s] * aqquad[ce[s]][k+ke[s]][j+je[s]] / (dx * dy);
      }
  }

  for (c = 0; c < 4; c++) {
      ierr = DMDAVecRestoreArray(info->da,qquad[c],&(aqquad[c])); CHKERRQ(ierr);
      ierr = DMRestoreLocalVector(info->da, &(qquad[c])); CHKERRQ(ierr);
  }
  ierr = DMDAVecRestoreArray(info->da,b,&ab); CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(info->da, &b); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// use j,k for x,y directions in loops and MatSetValuesStencil
typedef struct {
  PetscInt foo,k,j,bar;
} MyStencil;

PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info, double t,
                                  double **aH, double **aHdot, double shift,
                                  Mat J, Mat P, AppCtx *user) {
  PetscErrorCode  ierr;
  const double    dx = user->L / (double)(info->mx),
                  dy = user->L / (double)(info->my);
  const double    coeff[8] = {dy/2, dx/2, dx/2, -dy/2, -dy/2, -dx/2, -dx/2, dy/2};
  const PetscBool upwind = (user->lambda > 0.0);
  const double    upmin = (1.0 - user->lambda) * 0.5,
                  upmax = (1.0 + user->lambda) * 0.5;
  int             j, k, c, l, s, u, v;
  double          H, Hup, **aDqDlquad[16], **ab, val[33], DqSIADl_clkj;
  Grad            gH, gb;
  Vec             DqDlquad[16], b;
  MyStencil       col[33],row;

  PetscFunctionBeginUser;
  ierr = MatZeroEntries(P); CHKERRQ(ierr);  // because using ADD_VALUES below

  ierr = DMGetLocalVector(info->da, &b); CHKERRQ(ierr);
  if (user->verif > 0) {
      ierr = VecSet(b,0.0); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(info->da,b,&ab); CHKERRQ(ierr);
  } else {
      ierr = DMDAVecGetArray(info->da,b,&ab); CHKERRQ(ierr);
      ierr = FormBedLocal(info,1,ab,user); CHKERRQ(ierr);  // get stencil width
  }
  for (c = 0; c < 16; c++) {
      ierr = DMGetLocalVector(info->da, &(DqDlquad[c])); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(info->da,DqDlquad[c],&(aDqDlquad[c])); CHKERRQ(ierr);
  }

  // loop over locally-owned elements, including ghosts, to get DfluxDl for
  // l=0,1,2,3 at c=0,1,2,3 points in element;  note start at (xs-1,ys-1)
  for (k = info->ys-1; k < info->ys + info->ym; k++) {
      for (j = info->xs-1; j < info->xs + info->xm; j++) {
          for (c=0; c<4; c++) {
              double lxup = locx[c], lyup = locy[c];
              H  = fieldatptArray(j,k,locx[c],locy[c],aH);
              gH = gradfatptArray(j,k,locx[c],locy[c],dx,dy,aH);
              gb = gradfatptArray(j,k,locx[c],locy[c],dx,dy,ab);
              if (upwind) {
                  if (xdire[c] == PETSC_TRUE) {
                      lxup = (gb.x <= 0.0) ? upmin : upmax;
                      lyup = locy[c];
                  } else {
                      lxup = locx[c];
                      lyup = (gb.y <= 0.0) ? upmin : upmax;
                  }
                  Hup = fieldatptArray(j,k,lxup,lyup,aH);
              } else {
                  Hup = H;
              }
              for (l=0; l<4; l++) {
                  Grad   dgHdl;
                  double dHdl, dHupdl;
                  dgHdl  = dgradfatpt(l,locx[c],locy[c],dx,dy);
                  dHdl   = dfieldatpt(l,locx[c],locy[c]);
                  dHupdl = (upwind) ? dfieldatpt(l,lxup,lyup) : dHdl;
                  DqSIADl_clkj = DSIAfluxDl(gH,gb,dgHdl,H,dHdl,Hup,dHupdl,xdire[c],user);
                  aDqDlquad[4*c+l][k][j] = DqSIADl_clkj;
              }
          }
      }
  }

  // loop over nodes, not including ghosts, to get derivative of residual
  // with respect to nodal value
  for (k=info->ys; k<info->ys+info->ym; k++) {
      row.k = k;
      for (j=info->xs; j<info->xs+info->xm; j++) {
          row.j = j;
          for (s=0; s<8; s++) {
              u = j + je[s];
              v = k + ke[s];
              for (l=0; l<4; l++) {
                  const int djfroml[4] = { 0,  1,  1,  0},
                            dkfroml[4] = { 0,  0,  1,  1};
                  col[4*s+l].j = u + djfroml[l];
                  col[4*s+l].k = v + dkfroml[l];
                  val[4*s+l]   = coeff[s] * aDqDlquad[4*ce[s]+l][v][u] / (dx * dy);
              }
          }
          col[32].j = j;
          col[32].k = k;
          val[32]   = shift;
          ierr = MatSetValuesStencil(P,1,(MatStencil*)&row,
                                     33,(MatStencil*)col,val,ADD_VALUES);CHKERRQ(ierr);
      }
  }

  for (c = 0; c < 16; c++) {
      ierr = DMDAVecRestoreArray(info->da,DqDlquad[c],&(aDqDlquad[c])); CHKERRQ(ierr);
      ierr = DMRestoreLocalVector(info->da, &(DqDlquad[c])); CHKERRQ(ierr);
  }
  ierr = DMDAVecRestoreArray(info->da,b,&ab); CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(info->da, &b); CHKERRQ(ierr);

  ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if (J != P) {
    ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  }
  //ierr = MatView(J,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info, double t, double **aH,
                                    double **GG, AppCtx *user) {
  PetscErrorCode  ierr;
  const double    dx = user->L / (double)(info->mx),
                  dy = user->L / (double)(info->my);
  int             j, k;
  Vec             b;
  double          **ab, y, x, m;

  PetscFunctionBeginUser;
  ierr = DMGetLocalVector(info->da, &b); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(info->da,b,&ab); CHKERRQ(ierr);
  ierr = FormBedLocal(info,0,ab,user); CHKERRQ(ierr);  // stencil width NOT needed
  for (k=info->ys; k<info->ys+info->ym; k++) {
      y = k * dy;
      for (j=info->xs; j<info->xs+info->xm; j++) {
          x = j * dx;
          if (user->verif == 1) {
              m = DomeCMB(x,y,user);
          } else if (user->verif == 2) {
              m = 0.0;
          } else {
              m = M_CMBModel(user->cmb,ab[k][j] + aH[k][j]);  // s = b + H is surface elevation
          }
          GG[k][j] = m;
      }
  }
  ierr = DMDAVecRestoreArray(info->da,b,&ab); CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(info->da, &b); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}


text:
/*
GENETIC - A simple genetic algorithm.

Copyright 2014, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
	list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation
	and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file evolution.c
 * \brief Source file to define evolution functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2014 Javier Burguete Tolosa. All rights reserved.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
//#include "config.h"
#include "bits.h"
#include "sort.h"
#include "entity.h"
#include "population.h"
#include "mutation.h"
#include "reproduction.h"
#include "adaptation.h"
#include "selection.h"
#include "evolution.h"

#define DEBUG_EVOLUTION 0       ///< Macro to debug the evolution functions.

/**
 * Function to sort the survivals of the evaluated population
 */
void
evolution_sort (Population * population)        ///< Population.
{
  unsigned int i, j, index[population->nentities];
  double objective[population->nsurvival];
  Entity entity[population->nsurvival];
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_sort: start\n");
#endif
  index_new (population->objective, index, population->nentities);
  for (i = 0; i < population->nsurvival; ++i)
    {
      j = index[i];
      objective[i] = population->objective[j];
      entity_new (entity + i, population->genome_nbytes, i);
      memcpy (entity[i].genome, population->entity[j].genome,
              population->genome_nbytes);
    }
  for (i = 0; i < population->nsurvival; ++i)
    {
      population->objective[i] = objective[i];
      memcpy (population->entity[i].genome, entity[i].genome,
              population->genome_nbytes);
      entity_free (entity + i);
    }
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_sort: end\n");
#endif
}

/**
 * Funtion to apply the mutation evolution.
 */
void
evolution_mutation (Population * population,    ///< Population.
                    gsl_rng * rng)      ///< GSL random numbers generator.
{
  unsigned int i;
  Entity *mother, *son;
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_mutation: start\n");
#endif
  for (i = population->mutation_max; i > population->mutation_min;)
    {
#if DEBUG_EVOLUTION
      fprintf (stderr, "evolution_mutation: selection\n");
#endif
      selection_mutation (population, &mother, rng);
      son = population->entity + --i;
#if DEBUG_EVOLUTION
      fprintf (stderr, "evolution_mutation: mutation in %u\n", i);
#endif
      mutation (population, mother, son, rng);
    }
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_mutation: end\n");
#endif
}

/**
 * Funtion to apply the reproduction evolution.
 */
void
evolution_reproduction (Population * population,        ///< Population.
                        gsl_rng * rng)  ///< GSL random numbers generator.
{
  unsigned int i;
  Entity *mother, *father, *son;
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_reproduction: start\n");
#endif
  for (i = population->reproduction_max; i > population->reproduction_min;)
    {
#if DEBUG_EVOLUTION
      fprintf (stderr, "evolution_reproduction: selection\n");
#endif
      selection_reproduction (population, &mother, &father, rng);
      son = population->entity + --i;
#if DEBUG_EVOLUTION
      fprintf (stderr, "evolution_reproduction: reproduction in %u\n", i);
#endif
      reproduction (mother, father, son, population->genome_nbits, rng);
    }
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_reproduction: end\n");
#endif
}

/**
 * Funtion to apply the adaptation evolution.
 */
void
evolution_adaptation (Population * population,  ///< Population.
                      gsl_rng * rng)    ///< GSL random numbers generator.
{
  unsigned int i;
  Entity *mother, *son;
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_adaptation: start\n");
#endif
  for (i = population->adaptation_max; i > population->adaptation_min;)
    {
#if DEBUG_EVOLUTION
      fprintf (stderr, "evolution_adaptation: selection\n");
#endif
      selection_adaptation (population, &mother, rng);
      son = population->entity + --i;
#if DEBUG_EVOLUTION
      fprintf (stderr, "evolution_adaptation: adaptation in %u\n", i);
#endif
      adaptation (population, mother, son, rng);
    }
#if DEBUG_EVOLUTION
  fprintf (stderr, "evolution_adaptation: end\n");
#endif
}

max_stars_repo_name: jburguete/genetic

max_issues_repo_name: jburguete/genetic

max_forks_repo_name: jburguete/genetic


text:
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz, Conor Williams, Oliver Schonrock
//
// 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 <units/bits/external/hacks.h>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <tuple>
#include <type_traits>
#include <gsl/gsl-lite.hpp>

namespace units::detail {

template<typename T>
[[nodiscard]] constexpr T abs(T v) noexcept
{
  return v < 0 ? -v : v;
}

// the following functions enable gcd and related computations on ratios
// with exponents. They avoid overflow. Further information here:
// https://github.com/mpusz/units/issues/62#issuecomment-588152833

// Computes (a * b) mod m relies on unsigned integer arithmetic, should not
// overflow
[[nodiscard]] constexpr std::uint64_t mulmod(std::uint64_t a, std::uint64_t b, std::uint64_t m)
{
  std::uint64_t res = 0;

  if (b >= m) {
    if (m > UINT64_MAX / 2u) {
      b -= m;
    } else {
      b %= m;
    }
  }

  while (a != 0) {
    if (a & 1) {
      if (b >= m - res) {
        res -= m;
      }
      res += b;
    }
    a >>= 1;

    std::uint64_t temp_b = b;
    if (b >= m - b) {
      temp_b -= m;
    }
    b += temp_b;
  }

  return res;
}

// Calculates (a ^ e) mod m , should not overflow.
[[nodiscard]] constexpr std::uint64_t modpow(std::uint64_t a, std::uint64_t e, std::uint64_t m)
{
  a %= m;
  std::uint64_t result = 1;

  while (e > 0) {
    if (e & 1) {
      result = mulmod(result, a, m);
    }
    a = mulmod(a, a, m);
    e >>= 1;
  }
  return result;
}

// gcd(a * 10 ^ e, b), should not overflow
[[nodiscard]] constexpr std::intmax_t gcdpow(std::intmax_t a, std::intmax_t e, std::intmax_t b) noexcept
{
  assert(a > 0);
  assert(e >= 0);
  assert(b > 0);

  // gcd(i, j) = gcd(j, i mod j) for j != 0 Euclid;
  //
  // gcd(a 10^e, b) = gcd(b, a 10^e mod b)
  //
  // (a 10^e) mod b -> [ (a mod b) (10^e mod b) ] mod b

  return std::gcd(
      b, static_cast<std::intmax_t>(mulmod(static_cast<std::uint64_t>(a % b),
                                           modpow(10, static_cast<std::uint64_t>(e), static_cast<std::uint64_t>(b)),
                                           static_cast<std::uint64_t>(b))));
}

constexpr void cwap(std::intmax_t& lhs, std::intmax_t& rhs)
{
  std::intmax_t tmp = lhs;
  lhs = rhs;
  rhs = tmp;
}

// Computes the rational gcd of n1/d1 x 10^e1 and n2/d2 x 10^e2
[[nodiscard]] constexpr auto gcd_frac(std::intmax_t n1, std::intmax_t d1, std::intmax_t e1, std::intmax_t n2, std::intmax_t d2,
                                      std::intmax_t e2) noexcept
{
  // Short cut for equal ratios
  if (n1 == n2 && d1 == d2 && e1 == e2) {
    return std::array{n1, d1, e1};
  }

  if (e2 > e1) {
    detail::cwap(n1, n2);
    detail::cwap(d1, d2);
    detail::cwap(e1, e2);
  }

  std::intmax_t exp = e2;  // minimum

  // gcd(a/b,c/d) = gcd(a⋅d, c⋅b) / b⋅d

  assert(std::numeric_limits<std::intmax_t>::max() / n1 > d2);
  assert(std::numeric_limits<std::intmax_t>::max() / n2 > d1);

  std::intmax_t num = detail::gcdpow(n1 * d2, e1 - e2, n2 * d1);

  assert(std::numeric_limits<std::intmax_t>::max() / d1 > d2);

  std::intmax_t den = d1 * d2;

  std::intmax_t gcd = std::gcd(num, den);

  return std::array{num / gcd, den / gcd, exp};
}

constexpr void normalize(std::intmax_t& num, std::intmax_t& den, std::intmax_t& exp)
{
  if(num == 0) {
    den = 1;
    exp = 0;
    return;
  }

  std::intmax_t gcd = std::gcd(num, den);
  num = num * (den < 0 ? -1 : 1) / gcd;
  den = detail::abs(den) / gcd;

  while (num % 10 == 0) {
    num /= 10;
    ++exp;
  }
  while (den % 10 == 0) {
    den /= 10;
    --exp;
  }
}

[[nodiscard]] constexpr std::intmax_t safe_multiply(std::intmax_t lhs, std::intmax_t rhs)
{
  constexpr std::intmax_t c = std::uintmax_t(1) << (sizeof(std::intmax_t) * 4);

  const std::intmax_t a0 = detail::abs(lhs) % c;
  const std::intmax_t a1 = detail::abs(lhs) / c;
  const std::intmax_t b0 = detail::abs(rhs) % c;
  const std::intmax_t b1 = detail::abs(rhs) / c;

  gsl_Expects(a1 == 0 || b1 == 0);                               // overflow in multiplication
  gsl_Expects(a0 * b1 + b0 * a1 < (c >> 1));                     // overflow in multiplication
  gsl_Expects(b0 * a0 <= INTMAX_MAX);                            // overflow in multiplication
  gsl_Expects((a0 * b1 + b0 * a1) * c <= INTMAX_MAX - b0 * a0);  // overflow in multiplication

  return lhs * rhs;
}

}  // namespace units::detail

max_stars_repo_name: bebuch/units

max_issues_repo_name: bebuch/units

max_forks_repo_name: bebuch/units


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "parmt_utils.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"

/*!
 * @brief Computes the [mrows x 6] forward modeling matrix for synthesizing
 *        a waveform.
 *
 * @param[in] sacGxx   data structure describing the Green's functions to be
 *                     scaled by the mxx Green's functions
 * @param[in] sacGyy   data structure describing the Green's functions to be
 *                     scaled by the myy Green's functions
 * @param[in] sacGzz   data structure describing the Green's functions to be
 *                     scaled by the mzz Green's functions
 * @param[in] sacGxy   data structure describing the Green's functions to be
 *                     scaled by the mxy Green's functions
 * @param[in] sacGxz   data structure describing the Green's functions to be
 *                     scaled by the mxz Green's functions
 * @param[in] sacGyz   data structure describing the Green's functions to be
 *                     scaled by the myz Green's functions
 *
 * @param[out] mrows   number of rows in forward modeling matrix
 * @param[out] ierr    0 indicates success
 *
 * @result the [npts x 6] column major forward modeling matirx packed
 *         \f$ \{G_{xx}, G_{yy}, G_{zz}, G_{xy}, G_{xz}, G_{yz} \f$.
 *         this must be freed with ISCL's memory_free64f.

 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
double *parmt_utils_sacGrns2GrnsMat(const struct sacData_struct sacGxx,
                                    const struct sacData_struct sacGyy,
                                    const struct sacData_struct sacGzz,
                                    const struct sacData_struct sacGxy,
                                    const struct sacData_struct sacGxz,
                                    const struct sacData_struct sacGyz,
                                    int *mrows, int *ierr)
{
    const char *fcnm = "parmt_utils_sacGrns2GrnsMat\0";
    double *G;
    int npts, nptsGxx, nptsGyy, nptsGzz, nptsGxy, nptsGxz, nptsGyz;
    G = NULL;
    *mrows = 0;
    // Verify the inputs
    *ierr  = sacio_getIntegerHeader(SAC_INT_NPTS, sacGxx.header, &nptsGxx);
    *ierr += sacio_getIntegerHeader(SAC_INT_NPTS, sacGyy.header, &nptsGyy);
    *ierr += sacio_getIntegerHeader(SAC_INT_NPTS, sacGzz.header, &nptsGzz);
    *ierr += sacio_getIntegerHeader(SAC_INT_NPTS, sacGxy.header, &nptsGxy);
    *ierr += sacio_getIntegerHeader(SAC_INT_NPTS, sacGxz.header, &nptsGxz);
    *ierr += sacio_getIntegerHeader(SAC_INT_NPTS, sacGyz.header, &nptsGyz);
    if (*ierr != 0 || nptsGxx != nptsGyy ||
        nptsGxx != nptsGzz || nptsGxx != nptsGxy ||
        nptsGxx != nptsGxz || nptsGxx != nptsGyz)
    {
        if (*ierr != 0){printf("%s: Error reading headers\n", fcnm);}
        printf("%s: Inconsistent Green's functions trace lenths\n", fcnm);
        return G;
    }
    if (nptsGxx < 1)
    {
        printf("%s: No points in Green's functions\n", fcnm);
        return G;
    }
    if (nptsGxx != sacGxx.npts || nptsGxx != sacGyy.npts ||
        nptsGxx != sacGzz.npts || nptsGxx != sacGxy.npts ||
        nptsGxx != sacGxz.npts || nptsGxx != sacGyz.npts)
    {
        printf("%s: header npts differs from data npts\n", fcnm);
        return G;
    }
    if (sacGxx.data == NULL || sacGyy.data == NULL || sacGzz.data == NULL ||
        sacGxy.data == NULL || sacGxz.data == NULL || sacGyz.data == NULL)
    {
        printf("%s: null grns fns detected\n", fcnm);
        return G;
    }
    // Finally set the column major modeling matrix
    npts = nptsGxx;
    *mrows = npts;
    G = memory_calloc64f(npts*6);
    cblas_dcopy(npts, sacGxx.data, 1, &G[0*npts], 1);
    cblas_dcopy(npts, sacGyy.data, 1, &G[1*npts], 1);
    cblas_dcopy(npts, sacGzz.data, 1, &G[2*npts], 1);
    cblas_dcopy(npts, sacGxy.data, 1, &G[3*npts], 1);
    cblas_dcopy(npts, sacGxz.data, 1, &G[4*npts], 1);
    cblas_dcopy(npts, sacGyz.data, 1, &G[5*npts], 1);
    return G;
}
//============================================================================//
/*!
 * @brief Computes a synthetic seismogram from Green's functions
 *
 */
int parmt_utils_sacGrnsToEst(const struct sacData_struct data,
                             const struct sacData_struct sacGxx,
                             const struct sacData_struct sacGyy,
                             const struct sacData_struct sacGzz,
                             const struct sacData_struct sacGxy,
                             const struct sacData_struct sacGxz,
                             const struct sacData_struct sacGyz,
                             const double *__restrict__ mt,
                             const bool lrescale,
                             struct sacData_struct *est)
{
    const char *fcnm = "parmt_utils_sacGrnsToEst\0";
    double *G;
    int ierr, npts, nptsGxx, nptsGyy, nptsGzz, nptsGxy, nptsGxz, nptsGyz;
    const int ldg = 6;
    double xnum, xden, xscal;
    // Release memory on the target
    sacio_free(est);
    // Build the modeling matrix
    G = parmt_utils_sacGrns2GrnsMat(sacGxx, sacGyy, sacGzz,
                                    sacGxy, sacGxz, sacGyz,
                                    &npts, &ierr);
    if (ierr != 0)
    {
        printf("%s: Failed to make forward modeling matrix\n", fcnm);
        memory_free64f(&G);
    }
    // Copy the sacGxx to the estimate
    sacio_copy(data, est);
    est->header.isynth = 501; // set a number
    // Set the modeling matrix
    G = memory_calloc64f(npts*ldg);
    cblas_dcopy(npts, sacGxx.data, 1, &G[0*npts], 1);
    cblas_dcopy(npts, sacGyy.data, 1, &G[1*npts], 1);
    cblas_dcopy(npts, sacGzz.data, 1, &G[2*npts], 1);
    cblas_dcopy(npts, sacGxy.data, 1, &G[3*npts], 1);
    cblas_dcopy(npts, sacGxz.data, 1, &G[4*npts], 1);
    cblas_dcopy(npts, sacGyz.data, 1, &G[5*npts], 1);
    // Compute Gm = u to compute the estimate
    cblas_dgemv(CblasColMajor, CblasNoTrans,
                npts, 6, 1.0, G, npts, mt, 1, 0.0, est->data, 1);
    memory_free64f(&G);
    // Potentially rescale the synthetic
    if (lrescale)
    {
        xnum = array_maxAbs64f(npts, data.data, &ierr);
        xden = array_maxAbs64f(npts, est->data, &ierr);
        xscal = xnum/xden;
        cblas_dscal(npts, xscal, est->data, 1);
    }
    return 0;
}


text:
#include <math.h>
#include <stdio.h>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_vector.h>

#include <argtable2.h>

#include "qdm.h"

void
qdm_mspline_vector(
    gsl_vector *result,
    const double tau,
    const size_t spline_df,
    const gsl_vector *knots
)
{
  size_t bin = qdm_vector_search(knots, tau);

  gsl_vector_set(result, 0, 0);

  for (size_t m = 0; m < result->size - 1; m++) {
    double v;

    if (bin < m) {
      v = 0;
    } else if (bin - spline_df + 1 > m) {
      v = 0;
    } else if (bin == m) {
      double n = 3 * gsl_sf_pow_int(tau - gsl_vector_get(knots, m), 2);
      double d1 = gsl_vector_get(knots, m + 1) - gsl_vector_get(knots, m);
      double d2 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m);
      double d3 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);

      v = n / (d1 * d2 * d3);
    } else if (bin == m + 1) {
      double i1 = 1 / (gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1));

      double i2n = gsl_sf_pow_int(tau - gsl_vector_get(knots, m), 2);
      double i2d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double i2d2 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m);
      double i2d3 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i2 = i2n / (i2d1 * i2d2 * i2d3);

      double i3n = gsl_sf_pow_int(gsl_vector_get(knots, m + 3) - tau, 2);
      double i3d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double i3d2 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 1);
      double i3d3 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i3 = i3n / (i3d1 * i3d2 * i3d3);

      v = 3 * (i1 - i2 - i3);
    } else {
      double n = 3 * gsl_sf_pow_int(gsl_vector_get(knots, m + 3) - tau, 2);
      double d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double d2 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 1);
      double d3 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 2);
      double d = d1 * d2 * d3;

      v = n / d;
    }

    gsl_vector_set(result, m + 1, v);
  }
}

void
qdm_mspline_matrix(gsl_matrix *result, const size_t spline_df, const gsl_vector *x, const gsl_vector *knots)
{
  for (size_t i = 0; i < result->size1; i++) {
    gsl_vector_view row = gsl_matrix_row(result, i);
    qdm_mspline_vector(&row.vector, gsl_vector_get(x, i), spline_df, knots);
  }
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#pragma once
#include "halley/maths/vector3.h"
#include <gsl/span>

namespace Halley
{
	class AudioListenerData;
	class AudioChannelData;

	class AudioPosition
	{
	public:
		struct SpatialSource {
			SpatialSource();
			SpatialSource(Vector2f pos, float referenceDistance = 200.0f, float maxDistance = 400.0f);
			SpatialSource(Vector3f pos, float referenceDistance = 200.0f, float maxDistance = 400.0f);

			Vector3f pos;
			float referenceDistance;
			float maxDistance;
		};

		AudioPosition();

		static AudioPosition makeUI(float pan = 0.0f); // -1.0f = left, 1.0f = right
		static AudioPosition makePositional(Vector2f pos, float referenceDistance = 200.0f, float maxDistance = 400.0f);
		static AudioPosition makePositional(Vector3f pos, float referenceDistance = 200.0f, float maxDistance = 400.0f);
		static AudioPosition makePositional(std::vector<SpatialSource> sources);
		static AudioPosition makeFixed();

		void setMix(size_t srcChannels, gsl::span<const AudioChannelData> dstChannels, gsl::span<float, 16> dst, float gain, const AudioListenerData& listener) const;
		void setPosition(Vector3f position);

	private:
		std::vector<SpatialSource> sources;
		float pan = 0;
		bool isUI = false;
		bool isPannable = false;

		void setMixFixed(size_t srcChannels, gsl::span<const AudioChannelData> dstChannels, gsl::span<float, 16> dst, float gain, const AudioListenerData& listener) const;
		void setMixUI(gsl::span<const AudioChannelData> dstChannels, gsl::span<float, 16> dst, float gain, const AudioListenerData& listener) const;
		void setMixPositional(gsl::span<const AudioChannelData> dstChannels, gsl::span<float, 16> dst, float gain, const AudioListenerData& listener) const;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley