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

#include <Hyperion/error/GenericDomain.h>
#include <Hyperion/error/StatusCode.h>
#include <cstring>
#include <gsl/gsl>

namespace hyperion::error {

	class PosixDomain;

	using PosixStatusCode = StatusCode<PosixDomain>;
	using PosixErrorCode = ErrorCode<PosixDomain>;

	/// @brief `PosixDomain` is the `StatusCodeDomain` that covers status codes covering a
	/// platform's specific implementation of `errno` values in addition to those __strictly__
	/// required by POSIX (those represented by `Errno`).
	/// @ingroup error
	/// @headerfile "Hyperion/error/PosixDomain.h"
	class [[nodiscard("A StatusCodeDomain should always be used")]] PosixDomain {
	  public:
		/// @brief The value type of `PosixDomain` status codes is `i64`
		/// @ingroup error
		using value_type = i64;

		static const constexpr char(&UUID)[num_chars_in_uuid] // NOLINT
			= "4a6a9b0f-c335-473e-bc42-d23974a25bb0";

		static constexpr u64 ID = parse_uuid_from_string(UUID);

		/// @brief Constructs a `PosixDomain` with the default UUID
		/// @ingroup error
		constexpr PosixDomain() noexcept = default;
		/// @brief Constructs a `PosixDomain` with a user-specific UUID
		///
		/// @note When using a custom UUID __**ALL**__ instances of `PosixDomain` in the program
		/// should be constructed with the same custom UUID, otherwise equality comparison between
		/// other domains and `PosixDomain` instances could give erroneous results, and equality
		/// comparison between different `PosixDomain` instances will give erroneous results.
		/// As a result, this constructor should only be used when you specifically require a custom
		/// UUID and **YOU KNOW WHAT YOU ARE DOING™**
		///
		/// @param uuid - The UUID to use for `PosixDomain`
		/// @ingroup error
		explicit constexpr PosixDomain(u64 uuid) noexcept : m_uuid(uuid) {
		}
		/// @brief Constructs a `PosixDomain` with a user-specific UUID
		///
		/// @note When using a custom UUID __**ALL**__ instances of `PosixDomain` in the program
		/// should be constructed with the same custom UUID, otherwise equality comparison between
		/// other domains and `PosixDomain` instances could give erroneous results, and equality
		/// comparison between different `PosixDomain` instances will give erroneous results.
		/// As a result, this constructor should only be used when you specifically require a custom
		/// UUID and **YOU KNOW WHAT YOU ARE DOING™**
		///
		/// @param uuid - The UUID to use for `PosixDomain`
		/// @ingroup error
		template<UUIDString UUID>
		explicit constexpr PosixDomain(UUID && uuid) noexcept // NOLINT (forwarding reference)
			: m_uuid(parse_uuid_from_string(std::forward<UUID>(uuid))) {
		}
		/// @brief Copy-Constructor
		/// @ingroup error
		constexpr PosixDomain(const PosixDomain&) noexcept = default;
		/// @brief Move-Constructor
		/// @ingroup error
		constexpr PosixDomain(PosixDomain &&) noexcept = default;
		/// @brief Destructor
		/// @ingroup error
		constexpr ~PosixDomain() noexcept = default;

		/// @brief Returns the UUID of the domain
		///
		/// @return the domain UUID
		/// @ingroup error
		[[nodiscard]] constexpr auto id() const noexcept->u64 {
			return m_uuid;
		}

		/// @brief Returns the name of the domain
		///
		/// @return the domain name
		/// @ingroup error
		[[nodiscard]] constexpr auto name() const noexcept->std::string_view { // NOLINT
			return "POSIX domain";
		}

		/// @brief Returns the textual message associated with the given status code
		///
		/// @param code - The status code to get the message for
		///
		/// @return the message associated with the code
		/// @ingroup error
		[[nodiscard]] auto message(value_type code) // NOLINT
			const noexcept->std::string {

			return as_string(code);
		}

		/// @brief Returns the textual message associated with the given status code
		///
		/// @param code - The status code to get the message for
		///
		/// @return the message associated with the code
		/// @ingroup error
		[[nodiscard]] auto message(const PosixStatusCode& code) // NOLINT
			const noexcept->std::string {

			return as_string(code.code());
		}

		/// @brief Returns whether the given status code represents an error
		///
		/// @param code - The status code to check
		///
		/// @return `true` if the code represents an error, otherwise `false`
		/// @ingroup error
		[[nodiscard]] constexpr auto is_error(const PosixStatusCode& code) // NOLINT
			const noexcept->bool {
			return code.code() != 0;
		}

		/// @brief Returns whether the given status code represents success
		///
		/// @param code - The status code to check
		///
		/// @return `true` if the code represents success, otherwise `false`
		/// @ingroup error
		[[nodiscard]] constexpr auto is_success(const PosixStatusCode& code) // NOLINT
			const noexcept->bool {
			return code.code() == 0;
		}

		/// @brief Returns whether the given status codes are semantically equivalent
		///
		/// Checks if the given codes are semantically equivalent. For most `StatusCodeDomain`s,
		/// this usually means checking the codes for equality after being converted to
		/// `GenericStatusCode`s.
		///
		/// @tparam Domain - The `StatusCodeDomain` of the second status code
		/// @param lhs - The first status code to compare
		/// @param rhs - The second status code to compare
		/// @return `true` if the codes are semantically equivalent, `false` otherwise
		/// @ingroup error
		template<typename Domain>
		[[nodiscard]] constexpr auto are_equivalent(const PosixStatusCode& lhs,
													const StatusCode<Domain>& rhs)
			const noexcept->bool {
			if constexpr(ConvertibleToGenericStatusCode<StatusCode<Domain>>) {
				return as_generic_code(lhs) == rhs.as_generic_code();
			}
			else if(rhs.domain() == *this) {
				const auto lhs_code = lhs.code();
				const auto rhs_code = rhs.code();
				return lhs_code == rhs_code && lhs_code != -1 && rhs_code != -1;
			}
			else {
				return false;
			}
		}

		/// @brief Converts the given status code to a `GenericStatusCode`
		///
		/// This will convert the given code to its semantically equivalent counterpart in the
		/// `GenericDomain`.
		///
		/// @param code - The status code to convert to a `GenericStatusCode`
		/// @return The given status code as a `GenericStatusCode`
		/// @note Not all status code values are convertible to the `GenericDomain`, even
		/// from domains fully compatible with `GenericDomain` and that satisfy
		/// `ConvertibleToGenericStatusCode`. In this case, they will map to `Errno::Unknown`.
		/// Codes of value `Errno::Unknown` will never compare as semantically equivalent.
		/// @ingroup error
		[[nodiscard]] constexpr auto as_generic_code(const PosixStatusCode& code) // NOLINT
			const noexcept->GenericStatusCode {
			return make_status_code(to_generic_code(code.code()));
		}

		/// @brief Returns the value indicating success for this domain
		///
		/// @return The domain's success value
		/// @ingroup error
		[[nodiscard]] static inline constexpr auto success_value() noexcept->value_type {
			return 0;
		}

		/// @brief Returns the most recent value of `errno`
		///
		/// @return the value of `errno` at the time of the call
		/// @ingroup error
		[[nodiscard]] static inline auto get_last_error() noexcept->value_type {
			return errno;
		}

		/// @brief Domain equality comparison operator
		///
		/// @tparam Domain - The type of the second domain to compare
		///
		/// @param lhs - The left-hand domain to compare
		/// @param rhs - The right-hand domain to compare
		///
		/// @return Whether the two domains are equal
		/// @ingroup error
		template<typename Domain>
		friend constexpr auto operator==(const PosixDomain& lhs, const Domain& rhs) noexcept->bool {
			return lhs.id() == rhs.id();
		}

		/// @brief Domain inequality comparison operator
		///
		/// @tparam Domain - The type of the second domain to compare
		///
		/// @param lhs - The left-hand domain to compare
		/// @param rhs - The right-hand domain to compare
		///
		/// @return Whether the two domains are __not__ equal
		/// @ingroup error
		template<typename Domain>
		friend constexpr auto operator!=(const PosixDomain& lhs, const Domain& rhs) noexcept->bool {
			return lhs.id() != rhs.id();
		}

		/// @brief Copy-assignment operator
		/// @ingroup error
		constexpr auto operator=(const PosixDomain&) noexcept->PosixDomain& = default;
		/// @brief Move-assignment operator
		/// @ingroup error
		constexpr auto operator=(PosixDomain&&) noexcept->PosixDomain& = default;

	  private:
		u64 m_uuid = ID;

		[[nodiscard]] static inline auto as_string(value_type code) noexcept->std::string {
			char buffer[1024]; // NOLINT
#if HYPERION_PLATFORM_WINDOWS
			strerror_s(buffer, 1024, gsl::narrow_cast<i32>(code)); // NOLINT
#elif defined(__gnu_linux__) && !defined(__ANDROID__)
			char* message = strerror_r(gsl::narrow_cast<i32>(code), buffer, 1024); // NOLINT
			if(message != nullptr) {
				strncpy(buffer, message, 1024); // NOLINT
				buffer[1023] = 0;				// NOLINT
			}
#else
			str_error_r(code, buffer, 1024); // NOLINT
#endif
			usize length = strlen(buffer);		// NOLINT
			return std::string(buffer, length); // NOLINT
		}

		[[nodiscard]] static inline constexpr auto to_generic_code(
			value_type code) noexcept->Errno {
			switch(code) {
				case 0: return Errno::Success;
				case EAFNOSUPPORT: return Errno::AddressFamilyNotSupported;
				case EADDRINUSE: return Errno::AddressInUse;
				case EADDRNOTAVAIL: return Errno::AddressNotAvailable;
				case EISCONN: return Errno::AlreadyConnected;
				case E2BIG: return Errno::ArgumentListTooLong;
				case EDOM: return Errno::ArgumentOutOfDomain;
				case EFAULT: return Errno::BadAddress;
				case EBADF: return Errno::BadFileDescriptor;
				case EBADMSG: return Errno::BadMessage;
				case EPIPE: return Errno::BrokenPipe;
				case ECONNABORTED: return Errno::ConnectionAborted;
				case EALREADY: return Errno::ConnectionAlreadyInProgress;
				case ECONNREFUSED: return Errno::ConnectionRefused;
				case ECONNRESET: return Errno::ConnectionReset;
				case EXDEV: return Errno::CrossDeviceLink;
				case EDESTADDRREQ: return Errno::DestinationAddressRequired;
				case EBUSY: return Errno::DeviceOrResourceBusy;
				case ENOTEMPTY: return Errno::DirectoryNotEmpty;
				case ENOEXEC: return Errno::ExecutableFormatError;
				case EEXIST: return Errno::FileExists;
				case EFBIG: return Errno::FileTooLarge;
				case ENAMETOOLONG: return Errno::FilenameTooLong;
				case ENOSYS: return Errno::FunctionNotSupported;
				case EHOSTUNREACH: return Errno::HostUnreachable;
				case EIDRM: return Errno::IdentifierRemoved;
				case EILSEQ: return Errno::IllegalByteSequence;
				case ENOTTY: return Errno::InappropriateIOControlOperation;
				case EINTR: return Errno::Interrupted;
				case EINVAL: return Errno::InvalidArgument;
				case ESPIPE: return Errno::InvalidSeek;
				case EIO: return Errno::IOError;
				case EISDIR: return Errno::IsADirectory;
				case EMSGSIZE: return Errno::MessageSize;
				case ENETDOWN: return Errno::NetworkDown;
				case ENETRESET: return Errno::NetworkReset;
				case ENETUNREACH: return Errno::NetworkUnreachable;
				case ENOBUFS: return Errno::NoBufferSpace;
				case ECHILD: return Errno::NoChildProcess;
				case ENOLINK: return Errno::NoLink;
				case ENOLCK: return Errno::NoLockAvailable;
				case ENODATA: return Errno::NoMessageAvailable;
				case ENOMSG: return Errno::NoMessage;
				case ENOPROTOOPT: return Errno::NoProtocolOption;
				case ENOSPC: return Errno::NoSpaceOnDevice;
				case ENOSR: return Errno::NoStreamResources;
				case ENXIO: return Errno::NoSuchDeviceOrAddress;
				case ENODEV: return Errno::NoSuchDevice;
				case ENOENT: return Errno::NoSuchFileOrDirectory;
				case ESRCH: return Errno::NoSuchProcess;
				case ENOTDIR: return Errno::NotADirectory;
				case ENOTSOCK: return Errno::NotASocket;
				case ENOSTR: return Errno::NotAStream;
				case ENOTCONN: return Errno::NotConnected;
				case ENOMEM: return Errno::NotEnoughMemory;
#if ENOTSUP != EOPNOTSUPP
				case ENOTSUP: return Errno::NotSupported;
#endif
				case ECANCELED: return Errno::OperationCanceled;
				case EINPROGRESS: return Errno::OperationInProgress;
				case EPERM: return Errno::OperationNotPermitted;
				case EOPNOTSUPP: return Errno::OperationNotSupported;
#if EWOULDBLOCK != EAGAIN
				case EWOULDBLOCK: return Errno::OperationWouldBlock;
#endif
				case EOWNERDEAD: return Errno::OwnerDead;
				case EACCES: return Errno::PermissionDenied;
				case EPROTO: return Errno::ProtocolError;
				case EPROTONOSUPPORT: return Errno::ProtocolNotSupported;
				case EROFS: return Errno::ReadOnlyFileSystem;
				case EDEADLK: return Errno::ResourceDeadlockWouldOccur;
				case EAGAIN: return Errno::ResourceUnavailableTryAgain;
				case ERANGE: return Errno::ResultOutOfRange;
				case ENOTRECOVERABLE: return Errno::StateNotRecoverable;
				case ETIME: return Errno::StreamTimeout;
				case ETXTBSY: return Errno::TextFileBusy;
				case ETIMEDOUT: return Errno::TimedOut;
				case ENFILE: return Errno::TooManyFilesOpenInSystem;
				case EMFILE: return Errno::TooManyFilesOpen;
				case EMLINK: return Errno::TooManyLinks;
				case ELOOP: return Errno::TooManySymbolicLinkLevels;
				case EOVERFLOW: return Errno::ValueTooLarge;
				case EPROTOTYPE: return Errno::WrongProtocolType;
				default: return Errno::Unknown;
			}
		}
	};
} // namespace hyperion::error

/// @brief Specialize `make_status_code_domain` for `PosixDomain` and `u64`.
/// Creates a `PosixDomain` with a custom UUID.
///
/// @note When using a custom UUID __**ALL**__ instances of `PosixDomain` in the program
/// should be constructed with the same custom UUID, otherwise equality comparison between
/// other domains and `PosixDomain` instances could give erroneous results, and equality
/// comparison between different `PosixDomain` instances will give erroneous results.
/// As a result, this constructor should only be used when you specifically require a custom
/// UUID and **YOU KNOW WHAT YOU ARE DOING™**
///
/// @param uuid - The UUID to use for `PosixDomain`
///
/// @return a `PosixDomain`
/// @ingroup error
template<>
[[nodiscard]] inline constexpr auto
// NOLINTNEXTLINE(readability-inconsistent-declaration-parameter-name)
make_status_code_domain<hyperion::error::PosixDomain, hyperion::u64>(hyperion::u64&& uuid) noexcept
	-> hyperion::error::PosixDomain {
	return hyperion::error::PosixDomain(uuid);
}

/// @brief Specialize `make_status_code_domain` for `PosixDomain` and a `UUIDString`.
/// Creates a `PosixDomain` with a custom UUID.
///
/// @note When using a custom UUID __**ALL**__ instances of `PosixDomain` in the program
/// should be constructed with the same custom UUID, otherwise equality comparison between
/// other domains and `PosixDomain` instances could give erroneous results, and equality
/// comparison between different `PosixDomain` instances will give erroneous results.
/// As a result, this constructor should only be used when you specifically require a custom
/// UUID and **YOU KNOW WHAT YOU ARE DOING™**
///
/// @param uuid - The UUID to use for `PosixDomain`
///
/// @return a `PosixDomain`
/// @ingroup error
template<>
[[nodiscard]] inline constexpr auto
	// NOLINTNEXTLINE(readability-inconsistent-declaration-parameter-name)
	make_status_code_domain<hyperion::error::PosixDomain,
							const char (&)[hyperion::error::num_chars_in_uuid]> // NOLINT
	(const char (&uuid)[hyperion::error::num_chars_in_uuid]) noexcept			// NOLINT
	-> hyperion::error::PosixDomain {
	return hyperion::error::PosixDomain(uuid);
}

/// @brief Specialize `make_status_code_domain` for `PosixDomain` and an MS-style `UUIDString`.
/// Creates a `PosixDomain` with a custom UUID.
///
/// @note When using a custom UUID __**ALL**__ instances of `PosixDomain` in the program
/// should be constructed with the same custom UUID, otherwise equality comparison between
/// other domains and `PosixDomain` instances could give erroneous results, and equality
/// comparison between different `PosixDomain` instances will give erroneous results.
/// As a result, this constructor should only be used when you specifically require a custom
/// UUID and **YOU KNOW WHAT YOU ARE DOING™**
///
/// @param uuid - The UUID to use for `PosixDomain`
///
/// @return a `PosixDomain`
/// @ingroup error
template<>
[[nodiscard]] inline constexpr auto
	// NOLINTNEXTLINE(readability-inconsistent-declaration-parameter-name)
	make_status_code_domain<hyperion::error::PosixDomain,
							const char (&)[hyperion::error::num_chars_in_ms_uuid]> // NOLINT
	(const char (&uuid)[hyperion::error::num_chars_in_ms_uuid]) noexcept		   // NOLINT
	-> hyperion::error::PosixDomain {
	return hyperion::error::PosixDomain(uuid);
}

/// @brief Specialize `make_status_code_domain` for `PosixDomain` with no arguments. Creates
/// a `PosixDomain` with the default UUID.
///
/// @return a `PosixDomain`
/// @ingroup error
template<>
inline constexpr auto
make_status_code_domain<hyperion::error::PosixDomain>() noexcept -> hyperion::error::PosixDomain {
	return {};
}

namespace hyperion::error {
	static_assert(StatusCodeDomain<PosixDomain>);
} // namespace hyperion::error


text:
#ifndef __MATHINT_H_
#define __MATHINT_H_

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <immintrin.h>
#include <stdbool.h>
#if (defined __AVX__ || defined __SSE4_2__)
#include "math_x86.h"
#endif
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>

#define pi     3.14159265358979323846264338327950288
#define MAX(a, b) ((a > b) ? a : b)
#define MIN(a, b) ((a > b) ? b : a)
#define EULER  0.577215664901532860606512090082402431042
#define EPSILON __DBL_EPSILON__
#define STACK_SIZE 100000  /*FIXME: This must be adpative.*/

double IncompBesselK0_Simpson(double tol, int *cnt, double bes_a, double bes_b, int der);
double computeK0(double a, double b);
double computeINCBK0(double a, double b, int der);

typedef struct {
  double a;
  double b;
} gsl_params;
double bessel_f(double x, void * p);
double bessel_f_der(double x, void * p);
double call_gsl_bessel_integrator(double a, double b,
				  gsl_integration_workspace *w, int der);

void inline set_zero_4(double *x)
{
  x[0] = 0; x[1] = 0; x[2] = 0; x[3] = 0;
}

#endif

ext: h


text:
#ifndef TTT_GFX_TEXT_H
#define TTT_GFX_TEXT_H

#include <SDL_render.h>
#include <SDL_ttf.h>

#include <gsl/pointers>
#include <gsl/util>

#include <memory>
#include <string_view>
#include <string>
#include <stdexcept>
#include <cstdint>

namespace ttt::gfx
{
	class TextError final : public std::runtime_error
	{
	public:
		TextError(const std::string &msg) noexcept;
	};

	inline TextError::TextError(const std::string &msg) noexcept
		: runtime_error{msg}
	{
	}

	// A very simple class for rendering text.
	class Text final
	{
	public:
		Text() = default;
		Text(std::string_view str, gsl::not_null<TTF_Font *> font, const SDL_Colour &colour, bool shadow, gsl::not_null<SDL_Renderer *> renderer);

		void setText(std::string_view str, gsl::not_null<TTF_Font *> font, const SDL_Colour &colour, bool shadow, gsl::not_null<SDL_Renderer *> renderer);
		void clear() noexcept;

		void render(gsl::not_null<SDL_Renderer *> renderer, const SDL_FPoint &pos, double angle = 0.0, const SDL_FPoint *centre = nullptr, SDL_RendererFlip flip = SDL_FLIP_NONE) const noexcept;

		[[nodiscard]] SDL_Point getSize() const noexcept { return size; }

	private:
		using TexturePointer = std::unique_ptr<SDL_Texture, decltype(&SDL_DestroyTexture)>;

		TexturePointer primary{nullptr, SDL_DestroyTexture};
		TexturePointer secondary{nullptr, SDL_DestroyTexture};
		SDL_Point size{0, 0};

		[[nodiscard]] TexturePointer createTexture(std::string_view str, gsl::not_null<TTF_Font *> font, const SDL_Colour &colour, gsl::not_null<SDL_Renderer *> renderer);
		[[nodiscard]] SDL_Colour createShadowColour(const SDL_Colour &colour) const noexcept;
	};

	inline SDL_Colour Text::createShadowColour(const SDL_Colour &colour) const noexcept
	{
		return {
			gsl::narrow_cast<std::uint8_t>(colour.r % 5),
			gsl::narrow_cast<std::uint8_t>(colour.g % 5),
			gsl::narrow_cast<std::uint8_t>(colour.b % 5),
			colour.a
		};
	}
}

#endif

max_stars_repo_name: itsArtem/TicTacToe

max_issues_repo_name: itsArtem/TicTacToe

max_forks_repo_name: itsArtem/TicTacToe


text:
/* $Id$      */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2013                                               */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  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 2 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., 675 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*;Correspondence about this software should be addressed as follows: */
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/
#ifndef OBITRMFIT_H
#define OBITRMFIT_H

#include "Obit.h"
#include "ObitErr.h"
#include "ObitImage.h"
#include "ObitBeamShape.h"
#include "ObitThread.h"
#include "ObitInfoList.h"
#ifdef HAVE_GSL
#include <gsl/gsl_multifit_nlin.h>
#endif /* HAVE_GSL */

/*-------- Obit:  Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitRMFit.h
 *
 * ObitRMFit Class for fitting rotatioin measures (RM) to Stokes Q and U image pixels
 *
 * This class does least squares fitting of RM and EVPA at 0 lambda^2
 * Either an image cube or a set of single plane images at arbitrary
 * lambda squares may be fitted.
 * The result is an image cube of RM and  EVPA at 0 lambda^2 as the planes.
 * The function ObitRMFitEval will evaluate this fit and return an image
 * with the flux densities at the desired frequencies.
 *
 * \section ObitRMFitaccess Creators and Destructors
 * An ObitRMFit will usually be created using ObitRMFitCreate which allows
 * specifying a name for the object as well as other information.
 *
 * A copy of a pointer to an ObitRMFit should always be made using the
 * #ObitRMFitRef function which updates the reference count in the object.
 * Then whenever freeing an ObitRMFit or changing a pointer, the function
 * #ObitRMFitUnref will decrement the reference count and destroy the object
 * when the reference count hits 0.
 * There is no explicit destructor.
 */

/*--------------Class definitions-------------------------------------*/
/** ObitRMFit Class structure. */
typedef struct {
#include "ObitRMFitDef.h"   /* this class definition */
} ObitRMFit;

/*----------------- Macroes ---------------------------*/
/**
 * Macro to unreference (and possibly destroy) an ObitRMFit
 * returns a ObitRMFit*.
 * in = object to unreference
 */
#define ObitRMFitUnref(in) ObitUnref (in)

/**
 * Macro to reference (update reference count) an ObitRMFit.
 * returns a ObitRMFit*.
 * in = object to reference
 */
#define ObitRMFitRef(in) ObitRef (in)

/**
 * Macro to determine if an object is the member of this or a
 * derived class.
 * Returns TRUE if a member, else FALSE
 * in = object to reference
 */
#define ObitRMFitIsA(in) ObitIsA (in, ObitRMFitGetClass())

/*---------------Public functions---------------------------*/
/** Public: Class initializer. */
void ObitRMFitClassInit (void);

/** Public: Default Constructor. */
ObitRMFit* newObitRMFit (gchar* name);

/** Public: Create/initialize ObitRMFit structures */
ObitRMFit* ObitRMFitCreate (gchar* name, olong nterm);
/** Typedef for definition of class pointer structure */
typedef ObitRMFit* (*ObitRMFitCreateFP) (gchar* name,
						     olong nterm);

/** Public: ClassInfo pointer */
gconstpointer ObitRMFitGetClass (void);

/** Public: Copy (deep) constructor. */
ObitRMFit* ObitRMFitCopy  (ObitRMFit *in, ObitRMFit *out, ObitErr *err);

/** Public: Copy structure. */
void ObitRMFitClone (ObitRMFit *in, ObitRMFit *out, ObitErr *err);

/** Public: Fit spectrum to an image cube */
void ObitRMFitCube (ObitRMFit* in, ObitImage *inQImage, ObitImage *inUImage,
		    ObitImage *outImage, ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef void(*ObitRMFitCubeFP) (ObitRMFit* in,
				ObitImage *inQImage, ObitImage *inUImage,
				ObitImage *outImage, ObitErr *err);

/** Public: Fit spectrum to an array of images */
void ObitRMFitImArr (ObitRMFit* in, olong nimage,
		     ObitImage **imQArr, ObitImage **imUArr,
		     ObitImage *outImage, ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef void(*ObitRMFitImArrFP) (ObitRMFit* in, olong nimage,
				 ObitImage **imQArr, ObitImage **imUArr,
				 ObitImage *outImage, ObitErr *err);

/** Public: Fit single spectrum */
ofloat* ObitRMFitSingle (olong nlamb2, olong nterm, odouble refLamb2, odouble *lamb2,
			 ofloat *qflux, ofloat *qsigma, ofloat *uflux, ofloat *usigma,
			 ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef ofloat*(*ObitRMFitSingleFP) (olong nlamb2, olong nterm, odouble refLamb2, odouble *lamb2,
				     ofloat *qflux, ofloat *qsigma, ofloat *uflux, ofloat *usigma,
				     ObitErr *err);

/** Public: Make fitting arg structure */
gpointer ObitRMFitMakeArg (olong nlamb2, olong nterm,
			   odouble refLamb2, odouble *lamb2,
			   ofloat **out, ObitErr *err);

/** Public: Fit single RM using arg */
void ObitRMFitSingleArg (gpointer arg,
			 ofloat *qflux, ofloat *qsigma,
			 ofloat *uflux, ofloat *usigma,
			 ofloat *out);

/** Public: Kill fitting arg structure */
void ObitRMFitKillArg (gpointer arg);
/*----------- ClassInfo Structure -----------------------------------*/
/**
 * ClassInfo Structure.
 * Contains class name, a pointer to any parent class
 * (NULL if none) and function pointers.
 */
typedef struct  {
#include "ObitRMFitClassDef.h"
} ObitRMFitClassInfo;

#endif /* OBITFRMFIT_H */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit


text:
/* multifit_nlinear/svd.c
 *
 * Copyright (C) 2016 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.
 */

/*
 * This module handles the solution of the linear least squares
 * system:
 *
 * [     J      ] dx = - [ f ]
 * [ sqrt(mu)*D ]        [ 0 ]
 *
 * using an SVD approach. The system above is transformed to "standard form"
 * via:
 *
 * J~ = J D^{-1}
 * dx~ = D dx
 *
 * so that
 *
 * [     J~     ] dx~ = - [ f ]
 * [ sqrt(mu)*I ]         [ 0 ]
 *
 * can be solved with a standard SVD method, and then dx is recovered
 * from dx~ via: dx = D^{-1} dx~
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

typedef struct
{
  size_t n;                  /* number of residuals */
  size_t p;                  /* number of parameters */
  gsl_matrix *U;             /* U factor of J, n-by-p */
  gsl_matrix *V;             /* V factor of J, p-by-p */
  gsl_vector *S;             /* singular values, size p */
  gsl_vector *workp;         /* workspace, length p */
  double mu;                 /* LM parameter */
} svd_state_t;

static int svd_init(const void * vtrust_state, void * vstate);
static int svd_presolve(const double mu, const void * vtrust_state, void * vstate);
static int svd_solve(const gsl_vector * f, gsl_vector *x,
                     const void * vtrust_state, void *vstate);
static int svd_rcond(double * rcond, void * vstate);

static void *
svd_alloc (const size_t n, const size_t p)
{
  svd_state_t *state;

  (void)n;

  state = calloc(1, sizeof(svd_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate svd state", GSL_ENOMEM);
    }

  state->U = gsl_matrix_alloc(n, p);
  if (state->U == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for U", GSL_ENOMEM);
    }

  state->V = gsl_matrix_alloc(p, p);
  if (state->V == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for V", GSL_ENOMEM);
    }

  state->S = gsl_vector_alloc(p);
  if (state->S == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for S",
                      GSL_ENOMEM);
    }

  state->workp = gsl_vector_alloc(p);
  if (state->workp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workp",
                      GSL_ENOMEM);
    }

  state->mu = 0.0;
  state->n = n;
  state->p = p;

  return state;
}

static void
svd_free(void *vstate)
{
  svd_state_t *state = (svd_state_t *) vstate;

  if (state->U)
    gsl_matrix_free(state->U);

  if (state->V)
    gsl_matrix_free(state->V);

  if (state->S)
    gsl_vector_free(state->S);

  if (state->workp)
    gsl_vector_free(state->workp);

  free(state);
}

/* compute svd of J */
static int
svd_init(const void * vtrust_state, void * vstate)
{
  int status;
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  svd_state_t *state = (svd_state_t *) vstate;
  size_t i;

  gsl_matrix_set_zero(state->U);

  /* compute U = J D^{-1} */
  for (i = 0; i < state->p; ++i)
    {
      gsl_vector_const_view Ji = gsl_matrix_const_column(trust_state->J, i);
      gsl_vector_view ui = gsl_matrix_column(state->U, i);
      double di = gsl_vector_get(trust_state->diag, i);

      gsl_blas_daxpy(1.0 / di, &Ji.vector, &ui.vector);
    }

  status = gsl_linalg_SV_decomp(state->U, state->V, state->S, state->workp);

  return status;
}

static int
svd_presolve(const double mu, const void * vtrust_state, void * vstate)
{
  svd_state_t *state = (svd_state_t *) vstate;

  state->mu = mu;

  (void)vtrust_state;

  return GSL_SUCCESS;
}

static int
svd_solve(const gsl_vector * f, gsl_vector *x,
          const void * vtrust_state, void *vstate)
{
  int status = GSL_SUCCESS;
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  svd_state_t *state = (svd_state_t *) vstate;
  const size_t p = state->p;
  const double tol = GSL_DBL_EPSILON;
  const double s0 = gsl_vector_get(state->S, 0);
  size_t j;

  /* compute workp = - U^T f */
  gsl_blas_dgemv(CblasTrans, -1.0, state->U, f, 0.0, state->workp);

  /*
   * compute:
   *
   * workp = sum_i s_i / (s_i^2 + mu) (-u_i^T f)
   */

  if (state->mu == 0.0)
    {
      /*
       * compute Gauss-Newton direction by solving
       * J x = -f
       */

      for (j = 0; j < p; ++j)
        {
          double sj = gsl_vector_get(state->S, j);
          double *ptr = gsl_vector_ptr(state->workp, j);
          double alpha;

          if (sj <= tol * s0)
            alpha = 0.0;
          else
            alpha = 1.0 / sj;

          *ptr *= alpha;
        }
    }
  else
    {
      /*
       * solve:
       *
       * [  J D^{-1}  ] (D x) = -[ f ]
       * [ sqrt(mu) I ]          [ 0 ]
       *
       * using SVD factorization of J D^{-1}
       */

      for (j = 0; j < p; ++j)
        {
          double sj = gsl_vector_get(state->S, j);
          double *ptr = gsl_vector_ptr(state->workp, j);

          *ptr *= sj / (sj*sj + state->mu);
        }
    }

  /* compute: x = V * workp */
  gsl_blas_dgemv(CblasNoTrans, 1.0, state->V, state->workp, 0.0, x);

  /* compute D^{-1} x */
  gsl_vector_div(x, trust_state->diag);

  return status;
}

static int
svd_rcond(double * rcond, void * vstate)
{
  int status = GSL_SUCCESS;
  svd_state_t *state = (svd_state_t *) vstate;
  double smax = gsl_vector_get(state->S, 0);
  double smin = gsl_vector_get(state->S, state->p - 1);

  *rcond = smin / smax;

  return status;
}

static const gsl_multifit_nlinear_solver svd_type =
{
  "svd",
  svd_alloc,
  svd_init,
  svd_presolve,
  svd_solve,
  svd_rcond,
  svd_free
};

const gsl_multifit_nlinear_solver *gsl_multifit_nlinear_solver_svd = &svd_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include "Cesium3DTilesReader/Library.h"

#include <Cesium3DTiles/Tileset.h>
#include <CesiumJsonReader/ExtensionReaderContext.h>

#include <gsl/span>

#include <functional>
#include <memory>
#include <optional>
#include <string>
#include <unordered_map>
#include <vector>

namespace Cesium3DTilesReader {

/**
 * @brief The result of reading a tileset with
 * {@link TilesetReader::readTileset}.
 */
struct CESIUM3DTILESREADER_API TilesetReaderResult {
  /**
   * @brief The read tileset, or std::nullopt if the tileset could not be read.
   */
  std::optional<Cesium3DTiles::Tileset> tileset;

  /**
   * @brief Errors, if any, that occurred during the load process.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warnings, if any, that occurred during the load process.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief Options for how to read a tileset.
 */
struct CESIUM3DTILESREADER_API ReadTilesetOptions {};

/**
 * @brief Reads tilesets.
 */
class CESIUM3DTILESREADER_API TilesetReader {
public:
  /**
   * @brief Constructs a new instance.
   */
  TilesetReader();

  /**
   * @brief Gets the context used to control how extensions are loaded from a
   * tileset.
   */
  CesiumJsonReader::ExtensionReaderContext& getExtensions();

  /**
   * @brief Gets the context used to control how extensions are loaded from a
   * tileset.
   */
  const CesiumJsonReader::ExtensionReaderContext& getExtensions() const;

  /**
   * @brief Reads a tileset.
   *
   * @param data The buffer from which to read the tileset.
   * @param options Options for how to read the tileset.
   * @return The result of reading the tileset.
   */
  TilesetReaderResult readTileset(
      const gsl::span<const std::byte>& data,
      const ReadTilesetOptions& options = ReadTilesetOptions()) const;

private:
  CesiumJsonReader::ExtensionReaderContext _context;
};

} // namespace Cesium3DTilesReader


text:
#ifndef __TEST_UTILS_H__
#define __TEST_UTILS_H__

#define THEN_WHEN(x)
#define THEN_CHECK(x)

#include <initializer_list>
#include <memory>
#include <string>
#include <vector>

#include <gsl/gsl>
#include <gsl/string_span>

#include "base-utils/yaml.h"
#include "config/path.h"
#include "config/pattern.h"
#include "plugins/pluginUtils.h"

#include "executorStub.h"

namespace execHelper {
namespace config {
class SettingsNode;
} // namespace config

namespace test {
namespace baseUtils {
class ConfigFileWriter;
} // namespace baseUtils
} // namespace test
} // namespace execHelper

namespace std {
template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& stream) {
    for(const auto& element : stream) {
        os << element << ", ";
    }
    return os;
}

template <typename T, typename U>
std::ostream& operator<<(std::ostream& os, const std::pair<T, U>& stream) {
    os << stream.first << ": " << stream.second;
    return os;
}
} // namespace std

namespace execHelper {
namespace test {
namespace utils {
using Patterns = std::vector<config::Pattern>;

using Arguments = std::vector<std::string>;
struct MainVariables {
    int argc;
    std::unique_ptr<char*[]> argv;

    explicit MainVariables(const Arguments& arguments);
};

template <typename T> void appendVectors(T& appendTo, const T& appendFrom) {
    appendTo.insert(std::end(appendTo), std::begin(appendFrom),
                    std::end(appendFrom));
}

baseUtils::YamlWriter toYaml(const config::SettingsNode& settings,
                             const config::Patterns& patterns) noexcept;
void writeSettingsFile(
    gsl::not_null<baseUtils::ConfigFileWriter*> configFileWriter,
    const config::SettingsNode& settings,
    const config::Patterns& patterns) noexcept;
std::string convertToConfig(const Patterns& patterns) noexcept;
std::string
convertToConfig(const config::SettingsNode& rootSettings,
                const std::string& prepend = std::string()) noexcept;
std::string
convertToConfig(const config::SettingsNode& settings, const Patterns& patterns,
                const std::string& prepend = std::string()) noexcept;

std::string convertToConfig(std::string key, std::string value,
                            const std::string& prepend = std::string());
std::string convertToConfig(const std::string& key,
                            const std::initializer_list<std::string>& values,
                            const std::string& prepend = std::string());
std::string convertToConfig(const std::string& key,
                            const std::vector<std::string>& values,
                            const std::string& prepend = std::string());

std::string
convertToConfig(const std::initializer_list<std::string>& keys,
                const std::string& value,
                const std::string& prepend = std::string()) noexcept;
std::string
convertToConfig(const std::initializer_list<std::string>& keys,
                const std::initializer_list<std::string>& values,
                const std::string& prepend = std::string()) noexcept;
std::string
convertToConfig(const std::initializer_list<std::string>& keys,
                const std::vector<std::string>& values,
                const std::string& prepend = std::string()) noexcept;
std::string
convertToConfig(const std::vector<std::string>& keys,
                const std::initializer_list<std::string>& values,
                const std::string& prepend = std::string()) noexcept;
std::string
convertToConfig(const std::vector<std::string>& keys,
                const std::vector<std::string>& values,
                const std::string& prepend = std::string()) noexcept;
std::string basename(const std::string& file);

config::PatternCombinations createPatternCombination(
    const std::initializer_list<config::PatternKey>& keys,
    const std::initializer_list<config::PatternValue>& values) noexcept;

config::PatternCombinations
createPatternCombination(const config::PatternKeys& keys,
                         const config::PatternValues& values) noexcept;
plugins::PatternPermutator
makePatternPermutator(const config::Patterns& patterns) noexcept;
core::test::ExecutorStub::TaskQueue
getExpectedTasks(const core::Task& task,
                 const config::Patterns patterns) noexcept;
core::test::ExecutorStub::TaskQueue
getExpectedTasks(const core::test::ExecutorStub::TaskQueue& tasks,
                 const config::Patterns patterns) noexcept;

std::string toString(const config::SettingsNode& settings,
                     unsigned int nbOfTabs = 0) noexcept;

std::string inheritWorkingDirKey() noexcept;
Patterns getPredefinedPatterns() noexcept;
} // namespace utils
} // namespace test
} // namespace execHelper

#endif /* __TEST_UTILS_H__ */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
/*
 *     Copyright 2016 Couchbase, Inc
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include <memcached/dockey.h>
#include <gsl/gsl>
#include <limits>
#include <memory>
#include <string>
#include <type_traits>

#include "ep_types.h"

class SerialisedDocKey;

/**
 * StoredDocKey is a container for key data
 *
 * Internally an n byte key is stored in a n + sizeof(CollectionID) std::string.
 *  a) We zero terminate so that data() is safe for printing as a c-string.
 *  b) We store the CollectionID in byte 0:3. This is because StoredDocKey
 *    typically ends up being written to disk and the CollectionID forms part
 *    of the on-disk key. Accounting and for for the CollectionID means storage
 *    components don't have to create a new buffer into which they can layout
 *    CollectionID and key data.
 */
class StoredDocKey : public DocKeyInterface<StoredDocKey> {
public:
    /**
     * Construct empty - required for some std containers
     */
    StoredDocKey() : keydata() {
    }

    /**
     * Create a StoredDocKey from a DocKey
     * @param key DocKey that is to be copied-in
     */
    StoredDocKey(const DocKey& key) {
        if (key.getEncoding() == DocKeyEncodesCollectionId::Yes) {
            keydata.resize(key.size());
            std::copy(key.data(), key.data() + key.size(), keydata.begin());
        } else {
            // This key is for the default collection (which has a fixed value)
            keydata.resize(key.size() + 1);
            keydata[0] = DefaultCollectionLeb128Encoded;
            std::copy(key.data(), key.data() + key.size(), keydata.begin() + 1);
        }
    }

    /**
     * Create a StoredDocKey from a std::string (test code uses this)
     * @param key std::string to be copied-in
     * @param cid the CollectionID that the key applies to (and will be encoded
     *        into the stored data)
     */
    StoredDocKey(const std::string& key, CollectionID cid);

    const uint8_t* data() const {
        return reinterpret_cast<const uint8_t*>(keydata.data());
    }

    size_t size() const {
        return keydata.size();
    }

    DocNamespace getDocNamespace() const {
        return getCollectionID();
    }

    CollectionID getCollectionID() const;

    DocKeyEncodesCollectionId getEncoding() const {
        return DocKeyEncodesCollectionId::Yes;
    }

    /**
     * @return a DocKey that views this StoredDocKey but without any
     * collection-ID prefix.
     */
    DocKey makeDocKeyWithoutCollectionID() const;

    /**
     * Intended for debug use only
     * @returns cid:key
     */
    std::string to_string() const;

    /**
     * For tests only
     * @returns the 'key' part of the StoredDocKey
     */
    const char* c_str() const;

    int compare(const StoredDocKey& rhs) const {
        return keydata.compare(rhs.keydata);
    }

    bool operator==(const StoredDocKey& rhs) const {
        return keydata == rhs.keydata;
    }

    bool operator!=(const StoredDocKey& rhs) const {
        return !(*this == rhs);
    }

    bool operator<(const StoredDocKey& rhs) const {
        return keydata < rhs.keydata;
    }

    operator DocKey() const {
        return {keydata, DocKeyEncodesCollectionId::Yes};
    }

protected:
    std::string keydata;
};

std::ostream& operator<<(std::ostream& os, const StoredDocKey& key);

static_assert(sizeof(CollectionID) == sizeof(uint32_t),
              "StoredDocKey: CollectionID has changed size");

/**
 * A hash function for StoredDocKey so they can be used in std::map and friends.
 */
namespace std {
template <>
struct hash<StoredDocKey> {
    std::size_t operator()(const StoredDocKey& key) const {
        return key.hash();
    }
};
}

class MutationLogEntryV2;
class StoredValue;

/**
 * SerialisedDocKey maintains the key data in an allocation that is not owned by
 * the class. The class is essentially immutable, providing a "view" onto the
 * larger block.
 *
 * For example where a StoredDocKey needs to exist as part of a bigger block of
 * data, SerialisedDocKey is the class to use.
 *
 * A limited number of classes are friends and only those classes can construct
 * a SerialisedDocKey.
 */
class SerialisedDocKey : public DocKeyInterface<SerialisedDocKey> {
public:
    /**
     * The copy constructor is deleted due to the bytes living outside of the
     * object.
     */
    SerialisedDocKey(const SerialisedDocKey& obj) = delete;

    const uint8_t* data() const {
        return bytes;
    }

    size_t size() const {
        return length;
    }

    DocNamespace getDocNamespace() const {
        return getCollectionID();
    }

    CollectionID getCollectionID() const;

    DocKeyEncodesCollectionId getEncoding() const {
        return DocKeyEncodesCollectionId::Yes;
    }

    bool operator==(const DocKey& rhs) const;

    /**
     * Return how many bytes are (or need to be) allocated to this object
     */
    size_t getObjectSize() const {
        return getObjectSize(length);
    }

    /**
     * Return how many bytes are needed to store the DocKey
     * @param key a DocKey that needs to be stored in a SerialisedDocKey
     */
    static size_t getObjectSize(const DocKey key) {
        return getObjectSize(key.size());
    }

    /**
     * Create a SerialisedDocKey and return a unique_ptr to the object.
     * Note that the allocation is bigger than sizeof(SerialisedDocKey)
     * @param key a DocKey to be stored as a SerialisedDocKey
     */
    struct SerialisedDocKeyDelete {
        void operator()(SerialisedDocKey* p) {
            p->~SerialisedDocKey();
            delete[] reinterpret_cast<uint8_t*>(p);
        }
    };

    operator DocKey() const {
        return {bytes, length, DocKeyEncodesCollectionId::Yes};
    }

    /**
     * make a SerialisedDocKey and return a unique_ptr to it - this is used
     * in test code only.
     */
    static std::unique_ptr<SerialisedDocKey, SerialisedDocKeyDelete> make(
            const StoredDocKey& key) {
        std::unique_ptr<SerialisedDocKey, SerialisedDocKeyDelete> rval(
                reinterpret_cast<SerialisedDocKey*>(
                        new uint8_t[getObjectSize(key)]));
        new (rval.get()) SerialisedDocKey(key);
        return rval;
    }

protected:
    /**
     * These following classes are "white-listed". They know how to allocate
     * and construct this object so are allowed access to the constructor.
     */
    friend class MutationLogEntryV2;
    friend class MutationLogEntryV3;
    friend class StoredValue;

    SerialisedDocKey() : length(0), bytes() {
    }

    /**
     * Create a SerialisedDocKey from a DocKey. Protected constructor as
     * this must be used by friends who know how to pre-allocate the object
     * storage
     * @param key a DocKey to be copied in
     */
    SerialisedDocKey(const DocKey& key)
        : length(gsl::narrow_cast<uint8_t>(key.size())) {
        if (key.getEncoding() == DocKeyEncodesCollectionId::Yes) {
            std::copy(key.data(), key.data() + key.size(), bytes);
        } else {
            // This key is for the default collection
            bytes[0] = DefaultCollectionLeb128Encoded;
            std::copy(key.data(), key.data() + key.size(), bytes + 1);
            length++;
        }
    }

    /**
     * Create a SerialisedDocKey from a byte_buffer that has no collection data
     * and requires the caller to state the collection-ID
     * This is used by MutationLogEntryV1/V2 to V3 upgrades
     */
    SerialisedDocKey(cb::const_byte_buffer key, CollectionID cid);

    /**
     * Create a SerialisedDocKey from a byte_buffer that has collection data
     */
    SerialisedDocKey(cb::const_byte_buffer key)
        : length(gsl::narrow_cast<uint8_t>(key.size())) {
        std::copy(key.begin(), key.end(), bytes);
    }

    /**
     * Returns the size in bytes of this object - fixed size plus the variable
     * length for the bytes making up the key.
     */
    static size_t getObjectSize(size_t len) {
        return sizeof(SerialisedDocKey) +
               (len - sizeof(SerialisedDocKey().bytes));
    }

    uint8_t length{0};
    uint8_t bytes[1];
};

std::ostream& operator<<(std::ostream& os, const SerialisedDocKey& key);

static_assert(std::is_standard_layout<SerialisedDocKey>::value,
              "SeralisedDocKey: must satisfy is_standard_layout");


text:
/*	$Id$ */
/*
 * Copyright (c) 2014, 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#ifdef __linux__
#include <bsd/stdlib.h> /* arc4random() */
#endif
#include <string.h>

#include <glib.h>
#include <gtk/gtk.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_histogram.h>
#include <kplot.h>

#include "extern.h"

enum	kmltype {
	KML_INIT = 0,
	KML_KML,
	KML_DOCUMENT,
	KML_FOLDER,
	KML_PLACEMARK,
	KML_POINT,
	KML_COORDINATES,
	KML_DESCRIPTION,
	KML__MAX
};

enum	kmlkey {
	KMLKEY_MEAN,
	KMLKEY_MEAN_PERCENT,
	KMLKEY_STDDEV,
	KMLKEY_POPULATION,
	KMLKEY__MAX
};

struct	kmlsave {
	FILE		*f;
	struct sim	*cursim;
	size_t		 curisland;
	gchar		*buf;
	gsize		 bufsz;
	gsize		 bufmax;
	int		 buffering;
};

struct	kmlparse {
	enum kmltype	 elem; /* current element */
	struct kmlplace	*cur; /* currently-parsed kmlplace */
	gchar		*buf; /* current parse text buffer */
	gsize	 	 bufsz; /* size in parse buffer */
	gsize	 	 bufmax; /* maximum sized buffer */
	gchar		*altbuf;
	gsize	 	 altbufsz;
	gsize	 	 altbufmax;
	GList		*places; /* parsed places */
	gchar		*ign; /* element we're currently ignoring */
	size_t		 ignstack; /* stack of "ign" while ignoring */
#define	KML_STACKSZ	 128
	enum kmltype	 stack[128];
	size_t		 stackpos;
};

static	const char *const kmlkeys[KMLKEY__MAX] = {
	"mean", /* KMLKEY_MEAN */
	"meanpct", /* KMLKEY_MEAN_PERCENT */
	"stddev", /* KMLKEY_STDDEV */
	"population", /* KMLKEY_POPULATION */
};

static	const char *const kmltypes[KML__MAX] = {
	NULL, /* KML_INIT */
	"kml", /* KML_KML */
	"Document", /* KML_DOCUMENT */
	"Folder", /* KML_FOLDER */
	"Placemark", /* KML_PLACEMARK */
	"Point", /* KML_POINT */
	"coordinates", /* KML_COORDINATES */
	"description", /* KML_DESCRIPTION */
};

static	const double DEG_TO_RAD = 0.017453292519943295769236907684886;
static	const double EARTH_RADIUS_IN_METERS = 6372797.560856;

static double
kml_dist(const struct kmlplace *from, const struct kmlplace *to)
{
	double	 latitudeArc, longitudeArc,
		 latitudeH, lontitudeH, tmp;

	latitudeArc = (from->lat - to->lat) * DEG_TO_RAD;
	longitudeArc = (from->lng - to->lng) * DEG_TO_RAD;
	latitudeH = sin(latitudeArc * 0.5);
	latitudeH *= latitudeH;
	lontitudeH = sin(longitudeArc * 0.5);
	lontitudeH *= lontitudeH;
	tmp = cos(from->lat * DEG_TO_RAD) *
		cos(to->lat * DEG_TO_RAD);

	return(EARTH_RADIUS_IN_METERS * 2.0 *
		asin(sqrt(latitudeH + tmp * lontitudeH)));
}

static void
kml_append(gchar **buf, gsize *bufsz,
	gsize *bufmax, const char *text, gsize sz)
{

	/* XXX: no check for overflow... */
	if (*bufsz + sz + 1 > *bufmax) {
		*bufmax = *bufsz + sz + 1024;
		*buf = g_realloc(*buf, *bufmax);
	}

	memcpy(*buf + *bufsz, text, sz);
	*bufsz += sz;
	(*buf)[*bufsz] = '\0';
	g_assert(*bufsz <= *bufmax);
}

static enum kmltype
kml_lookup(const gchar *name)
{
	enum kmltype	i;

	for (i = 0; i < KML__MAX; i++) {
		if (NULL == kmltypes[i])
			continue;
		if (0 == g_strcmp0(kmltypes[i], name))
			break;
	}

	return(i);
}

static void
kmlparse_free(gpointer dat)
{
	struct kmlplace	*place = dat;

	if (NULL == place)
		return;
	free(place);
}

void
kml_free(struct kml *kml)
{

	if (NULL == kml)
		return;

	g_list_free_full(kml->kmls, kmlparse_free);

	if (NULL != kml->file)
		g_mapped_file_unref(kml->file);
}

/*
 * Try to find the string "@@population=NN@@", which stipulates the
 * population for this particular island.
 * If we don't find it, just return--we'll use the default.
 * If we find a bad population, raise an error.
 */
static int
kml_placemark(const gchar *buf, struct kmlplace *place)
{
	const gchar	*cp;
	gchar		*ep;
	gsize		 keysz;
	gchar		 nbuf[22];

	keysz = strlen("@@population=");
	while (NULL != (cp = strstr(buf, "@@population="))) {
		buf = cp + keysz;
		if (NULL == (cp = strstr(buf, "@@")))
			break;
		if ((gsize)(cp - buf) >= sizeof(nbuf) - 1)
			return(0);
		memcpy(nbuf, buf, cp - buf);
		nbuf[cp - buf] = '\0';
		place->pop = g_ascii_strtoull(buf, &ep, 10);
		if (ERANGE == errno || EINVAL == errno || ep == buf)
			return(0);
		break;
	}

	return(1);
}

static void
kml_elem_end(GMarkupParseContext *ctx,
	const gchar *name, gpointer dat, GError **er)
{
	struct kmlparse	 *p = dat;
	enum kmltype	  t;
	gchar		**set;

	if (NULL != p->ign) {
		g_assert(p->ignstack > 0);
		if (0 == g_strcmp0(p->ign, name))
			p->ignstack--;
		if (p->ignstack > 0)
			return;
		g_free(p->ign);
		p->ign = NULL;
		return;
	}

	t = kml_lookup(name);
	g_assert(p->stackpos > 0);
	g_assert(t == p->stack[p->stackpos - 1]);
	p->stackpos--;

	switch (p->stack[p->stackpos]) {
	case (KML_PLACEMARK):
		/*
		 * if we're ending a Placemark, first check whether
		 * we've listed a population somewhere in any of the
		 * nested text segments.
		 * Also make sure that we have some coordinates (the
		 * default value was 360 for both, which of course isn't
		 * a valid coordinate).
		 */
		g_assert(NULL != p->cur);
		if ( ! kml_placemark(p->altbuf, p->cur)) {
			*er = g_error_new_literal
				(G_MARKUP_ERROR,
				 G_MARKUP_ERROR_INVALID_CONTENT,
				 "Cannot parse population");
			break;
		} else if (p->cur->lat > 180 || p->cur->lng > 180) {
			*er = g_error_new_literal
				(G_MARKUP_ERROR,
				 G_MARKUP_ERROR_INVALID_CONTENT,
				 "No coordinates for placemark");
			break;
		}
		p->places = g_list_append(p->places, p->cur);
		p->cur = NULL;
		break;
	case (KML_COORDINATES):
		/*
		 * Parse coordinates from the mix.
		 * Coordinates are longitude,latitude,altitude.
		 * Parse quickly and just make sure they're valid.
		 */
		set = g_strsplit(p->buf, ",", 3);
		p->cur->lng = g_ascii_strtod(set[0], NULL);
		if (ERANGE == errno)
			*er = g_error_new_literal
				(G_MARKUP_ERROR,
				 G_MARKUP_ERROR_INVALID_CONTENT,
				 "Cannot parse longitude");
		else if (p->cur->lng > 180.0 || p->cur->lng < -180.0)
			*er = g_error_new_literal
				(G_MARKUP_ERROR,
				 G_MARKUP_ERROR_INVALID_CONTENT,
				 "Invalid longitude");
		p->cur->lat = g_ascii_strtod(set[1], NULL);
		if (ERANGE == errno)
			*er = g_error_new_literal
				(G_MARKUP_ERROR,
				 G_MARKUP_ERROR_INVALID_CONTENT,
				 "Cannot parse latitude");
		else if (p->cur->lat > 90.0 || p->cur->lat < -90.0)
			*er = g_error_new_literal
				(G_MARKUP_ERROR,
				 G_MARKUP_ERROR_INVALID_CONTENT,
				 "Invalid latitude");
		g_strfreev(set);
		break;
	default:
		break;
	}
}

static void
kml_elem_start(GMarkupParseContext *ctx,
	const gchar *name, const gchar **attrn,
	const gchar **attrv, gpointer dat, GError **er)
{
	enum kmltype	 t;
	struct kmlparse	*p = dat;

	if (NULL != p->ign) {
		g_assert(p->ignstack > 0);
		if (0 == g_strcmp0(p->ign, name))
			p->ignstack++;
		return;
	}

	if (KML__MAX == (t = kml_lookup(name))) {
		assert(0 == p->ignstack);
		p->ign = g_strdup(name);
		p->ignstack = 1;
		return;
	}

	p->stack[p->stackpos++] = t;
	p->bufsz = 0;

	switch (p->stack[p->stackpos - 1]) {
	case (KML_PLACEMARK):
		/*
		 * If we're starting a Placemark, initialise ourselves
		 * with bad coorindates and a default population.
		 */
		if (NULL == p->cur) {
			p->cur = g_malloc0(sizeof(struct kmlplace));
			p->cur->lat = p->cur->lng = 360;
			p->cur->pop = 2;
			p->altbufsz = 0;
			break;
		}
		*er = g_error_new_literal
			(G_MARKUP_ERROR,
			 G_MARKUP_ERROR_INVALID_CONTENT,
			 "Nested placemarks not allowed.");
		break;
	default:
		break;
	}
}

static void
kml_text(GMarkupParseContext *ctx,
	const gchar *txt, gsize sz, gpointer dat, GError **er)
{
	struct kmlparse	*p = dat;

	if (NULL != p->cur)
		kml_append(&p->altbuf, &p->altbufsz,
			&p->altbufmax, txt, sz);

	/* No collection w/o element or while ignoring. */
	if (NULL != p->ign || 0 == p->stackpos)
		return;

	/* Each element for which we're going to collect text. */
	switch (p->stack[p->stackpos - 1]) {
	case (KML_COORDINATES):
		break;
	default:
		return;
	}

	kml_append(&p->buf, &p->bufsz, &p->bufmax, txt, sz);
}

static void
kml_error(GMarkupParseContext *ctx, GError *er, gpointer dat)
{

	g_warning("%s", er->message);
}

struct kml *
kml_torus(size_t islands, size_t islanders)
{
	struct kml	*kml;
	struct kmlplace	*p;
	size_t		 i;

	kml = g_malloc0(sizeof(struct kml));

	for (i = 0; i < islands; i++) {
		p = g_malloc0(sizeof(struct kmlplace));
		p->pop = islanders;
		p->lng = 360 * i / (double)islands - 180.0;
		p->lat = 0.0;
		kml->kmls = g_list_append(kml->kmls, p);
	}

	return(kml);
}

struct kml *
kml_rand(size_t islands, size_t islanders)
{
	struct kml	*kml;
	struct kmlplace	*p;
	size_t		 i;

	kml = g_malloc0(sizeof(struct kml));

	for (i = 0; i < islands; i++) {
		p = g_malloc0(sizeof(struct kmlplace));
		p->pop = islanders;
		p->lng = 360 * arc4random() / (double)UINT32_MAX - 180.0;
		p->lat = 180 * arc4random() / (double)UINT32_MAX - 90.0;
		kml->kmls = g_list_append(kml->kmls, p);
	}

	return(kml);
}

struct kml *
kml_parse(const gchar *file, GError **er)
{
	GMarkupParseContext	*ctx;
	GMarkupParser	  	 parse;
	GMappedFile		*f;
	int			 rc;
	struct kmlparse		 data;
	struct kml		*kml;

	if (NULL != er)
		*er = NULL;

	memset(&parse, 0, sizeof(GMarkupParser));
	memset(&data, 0, sizeof(struct kmlparse));

	data.elem = KML_INIT;

	parse.start_element = kml_elem_start;
	parse.end_element = kml_elem_end;
	parse.text = kml_text;
	parse.error = kml_error;

	if (NULL == (f = g_mapped_file_new(file, FALSE, er)))
		return(NULL);

	ctx = g_markup_parse_context_new
		(&parse, 0, &data, NULL);
	g_assert(NULL != ctx);
	rc = g_markup_parse_context_parse
		(ctx, g_mapped_file_get_contents(f),
		 g_mapped_file_get_length(f), er);

	g_markup_parse_context_free(ctx);
	g_free(data.buf);
	g_free(data.altbuf);
	g_free(data.ign);
	kmlparse_free(data.cur);

	if (0 == rc) {
		g_list_free_full(data.places, kmlparse_free);
		g_mapped_file_unref(f);
		return(NULL);
	} else if (NULL == data.places) {
		g_mapped_file_unref(f);
		return(NULL);
	}

	g_assert(NULL != data.places);
	kml = g_malloc0(sizeof(struct kml));
	kml->file = f;
	kml->kmls = data.places;
	return(kml);
}

double **
kml_migration_twonearest(GList *list, enum maptop map)
{
	double		**p;
	double		  dist, min;
	size_t		  i, j, len, minj, min2j;
	struct kmlplace	 *pl1, *pl2;

	len = (size_t)g_list_length(list);
	g_assert(len > 0);
	if (1 == len) {
		g_debug("Request for two nearest falling "
			"back to nearest");
		return(kml_migration_nearest(list, map));
	}

	p = g_malloc0_n(len, sizeof(double *));

	/*
	 * Special case the torus, which has a well-defined layout
	 * between nodes, such that the next ("right") and previous
	 * ("left") islands, wrapping around, get the migrants.
	 */
	if (MAPTOP_TORUS == map) {
		for (i = 0; i < len; i++) {
			p[i] = g_malloc0_n(len, sizeof(double));
			if (i < len - 1)
				p[i][i + 1] = 0.5;
			else
				p[i][0] = 0.5;
			if (i > 0)
				p[i][i - 1] = 0.5;
			else
				p[i][len - 1] = 0.5;
		}
		return(p);
	}

	for (i = 0; i < len; i++) {
		pl1 = g_list_nth_data(list, i);
		p[i] = g_malloc0_n(len, sizeof(double));
		for (minj = j = 0, min = DBL_MAX; j < len; j++) {
			if (i == j)
				continue;
			pl2 = g_list_nth_data(list, j);
			if ((dist = kml_dist(pl1, pl2)) < min) {
				min = dist;
				minj = j;
			}
		}
		p[i][minj] = 0.5;
		for (min2j = j = 0, min = DBL_MAX; j < len; j++) {
			if (i == j || j == minj)
				continue;
			pl2 = g_list_nth_data(list, j);
			if ((dist = kml_dist(pl1, pl2)) < min) {
				min = dist;
				min2j = j;
			}
		}
		p[i][min2j] = 0.5;
		g_assert(min2j != minj);
		g_assert(i != min2j);
	}

	return(p);
}

double **
kml_migration_nearest(GList *list, enum maptop map)
{
	double		**p;
	double		  dist, min;
	size_t		  i, j, len, minj;
	struct kmlplace	 *pl1, *pl2;

	len = (size_t)g_list_length(list);
	g_assert(len > 0);
	p = g_malloc0_n(len, sizeof(double *));

	/*
	 * Special case the torus, which has a well-defined layout
	 * between nodes, such that the next ("right") island, wrapping
	 * around, gets the migrant.
	 */
	if (MAPTOP_TORUS == map) {
		for (i = 0; i < len; i++) {
			p[i] = g_malloc0_n(len, sizeof(double));
			p[i][(i + 1) % len] = 1.0;
		}
		return(p);
	}

	for (i = 0; i < len; i++) {
		pl1 = g_list_nth_data(list, i);
		p[i] = g_malloc0_n(len, sizeof(double));
		for (minj = j = 0, min = DBL_MAX; j < len; j++) {
			if (i == j)
				continue;
			pl2 = g_list_nth_data(list, j);
			if ((dist = kml_dist(pl1, pl2)) < min) {
				min = dist;
				minj = j;
			}
		}
		p[i][minj] = 1.0;
	}

	return(p);
}

double **
kml_migration_distance(GList *list, enum maptop map)
{
	double		**p;
	double		  dist, sum;
	size_t		  i, j, len;
	struct kmlplace	 *pl1, *pl2;

	len = (size_t)g_list_length(list);
	g_assert(len > 0);
	p = g_malloc0_n(len, sizeof(double *));
	for (i = 0; i < len; i++) {
		pl1 = g_list_nth_data(list, i);
		p[i] = g_malloc0_n(len, sizeof(double));
		for (sum = 0.0, j = 0; j < len; j++) {
			if (i == j) {
				p[i][j] = 0.0;
				continue;
			}
			pl2 = g_list_nth_data(list, j);
			dist = kml_dist(pl1, pl2);
			p[i][j] = 1.0 / (dist * dist);
			sum += p[i][j];
		}
		for (j = 0; j < len; j++)
			if (i != j)
				p[i][j] = p[i][j] / sum;
	}

	return(p);
}

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate


text:
#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_errno.h>

void inc_gamma_error_msg(const char* fn_name, const char* reason, const char* file, int line, int gsl_errno) {
	fprintf(stderr, "Error in GSL %s\n", fn_name);
	fprintf(stderr, "GSL Error %i: %s.\n", gsl_errno, gsl_strerror(gsl_errno));
	fprintf(stderr, "%s in %s:%i.\n", reason, file, line);
}

void gamma_inc_lower_error_handler(const char* reason, const char* file, int line, int gsl_errno) {
	if (gsl_errno != GSL_EUNDRFLW) {
		inc_gamma_error_msg("gamma_inc_P while computing lower incomplete gamma ratio.", reason, file, line, gsl_errno);
		abort();
	}
}

void gamma_inc_upper_error_handler(const char* reason, const char* file, int line, int gsl_errno) {
	if (gsl_errno != GSL_EUNDRFLW) {
		inc_gamma_error_msg("gamma_inc while computing upper incomplete gamma ratio.", reason, file, line, gsl_errno);
		abort();
	}
}

#define INC_GAMMA_NO_UNDERFLOW(fn, handler) do { \
		gsl_error_handler_t* old_handler = gsl_set_error_handler(& handler ); \
		fn; \
		gsl_set_error_handler(old_handler); \
	} while (0)

void inc_gamma_check_value(long double val, const char* fn_name) {
	if (isnan(val) || !isfinite(val)) {
		fprintf(stderr, "Error while computing %s.\n", fn_name);
		if (isnan(val)) {
			fprintf(stderr, "NaN detected. ");
		} else {
			fprintf(stderr, "Infinity detected. ");
		}
		fprintf(stderr, "This is likely the result of an overflow.\n");
		abort();
	}
}

/*
	Computes the upper Gamma(x+y,a)/Gamma(x) for x = [1 ... x_max]
*/
void inc_gamma_upper_ratio_ary(int x_max, double y, double a, double* result) {

	// x = 1 case
	size_t idx = 0;
	int x = 1;
	long double ln_gamma_xp1 = 0;
	// With fn(x) = Gamma(x+y, a)/Gamma(x)
	long double fn_x;
	INC_GAMMA_NO_UNDERFLOW(fn_x = gsl_sf_gamma_inc(x+y, a), gamma_inc_upper_error_handler);
	result[idx] = fn_x;

	for (;x < x_max;) {
		ln_gamma_xp1 += logl((long double) x);
		// Using the recurrence relation Gamma(s+1,a) = s Gamma(s,a) + a^s exp(-a)
		// We have Gamma((x+1)+y)/Gamma(x+1) = (x+y) ( Gamma(x+y)/Gamma(x) )/x + a^(x+y) exp(-a)/Gamma(x+1)
		// -> fn(x+1) = (x+y) fn(x)/x + exp((x+y) ln a - a - ln Gamma(x+1))
		long double fn_xp1 = (x+y) * fn_x/x + expl((x+y)*logl(a) - a - ln_gamma_xp1);
		result[idx+1] = (double)fn_xp1;
		fn_x = fn_xp1;
		idx += 1;
		x += 1;
	}
	// Make use of the fact that NaNs and infinities will propagate to the last value computed.
	inc_gamma_check_value(result[x_max-1], "upper incomplete gamma ratio");
}

/*
	Computes the lower gamma(x+y,a)/Gamma(x) for x = [1 ... x_max]
*/
void inc_gamma_lower_ratio_ary(int x_max, double y, double a, double* result) {

	// x = x_max case
	size_t idx = x_max - 1;
	int x = x_max;
	long double ln_gamma_x = gsl_sf_lngamma(x);
	// With fn(x) = gamma(x+y, a)/Gamma(x)
	long double fn_x;
	//INC_GAMMA_NO_UNDERFLOW(fn_x = expl((x+y)*logl(a) - a - logl(x+y) - ln_gamma_x)*gsl_sf_hyperg_1F1(1, x+y+1, a), hyperg_error_handler);
	INC_GAMMA_NO_UNDERFLOW(fn_x = gsl_sf_gamma_inc_P(x+y,a) * expl(gsl_sf_lngamma(x+y) - ln_gamma_x), gamma_inc_lower_error_handler);
	result[idx] = fn_x;

	for (;x>1;) {
		x -= 1;
		idx -= 1;
		ln_gamma_x -= logl((long double) x);
		long double fn_xp1 = fn_x;
		// We have gamma(s,a) = (gamma(s+1, a) + a^s exp(-a))/s
		// So gamma(x+y,a)/Gamma(x) = ( x gamma((x+1)+y, a)/Gamma(x+1) + a^(x+y) exp(-a)/Gamma(x) ) / s
		// -> fn(x) = ( x fn(x+1) + exp( (x+y) ln a - a - ln Gamma(x)! ) ) / (x+y)
		fn_x = ( x*fn_xp1 + expl( (x+y)*logl(a) - a - ln_gamma_x ) ) / (x+y);
		result[idx] = (double)fn_x;
	}
	// Make use of the fact that NaNs and infinities will propagate to the last value computed.
	inc_gamma_check_value(result[0], "lower incomplete gamma ratio");
}

max_stars_repo_name: bsafdi/NPTFit

max_issues_repo_name: bsafdi/NPTFit

max_forks_repo_name: bsafdi/NPTFit


text:
#ifndef AMICI_SUNDIALS_MATRIX_WRAPPER_H
#define AMICI_SUNDIALS_MATRIX_WRAPPER_H

#include <sunmatrix/sunmatrix_band.h>   // SUNMatrix_Band
#include <sunmatrix/sunmatrix_dense.h>  // SUNMatrix_Dense
#include <sunmatrix/sunmatrix_sparse.h> // SUNMatrix_Sparse

#include <gsl/gsl-lite.hpp>

#include <vector>

#include "amici/vector.h"

namespace amici {

/**
 * @brief A RAII wrapper for SUNMatrix structs.
 *
 * This can create dense, sparse, or banded matrices using the respective
 * constructor.
 */
class SUNMatrixWrapper {
  public:
    SUNMatrixWrapper() = default;

    /**
     * @brief Create sparse matrix. See SUNSparseMatrix in sunmatrix_sparse.h
     * @param M Number of rows
     * @param N Number of columns
     * @param NNZ Number of nonzeros
     * @param sparsetype Sparse type
     */
    SUNMatrixWrapper(sunindextype M, sunindextype N, sunindextype NNZ,
                     int sparsetype);

    /**
     * @brief Create dense matrix. See SUNDenseMatrix in sunmatrix_dense.h
     * @param M Number of rows
     * @param N Number of columns
     */
    SUNMatrixWrapper(sunindextype M, sunindextype N);

    /**
     * @brief Create banded matrix. See SUNBandMatrix in sunmatrix_band.h
     * @param M Number of rows and columns
     * @param ubw Upper bandwidth
     * @param lbw Lower bandwidth
     */
    SUNMatrixWrapper(sunindextype M, sunindextype ubw, sunindextype lbw);

    /**
     * @brief Create sparse matrix from dense or banded matrix. See
     * SUNSparseFromDenseMatrix and SUNSparseFromBandMatrix in
     * sunmatrix_sparse.h
     * @param A Wrapper for dense matrix
     * @param droptol tolerance for dropping entries
     * @param sparsetype Sparse type
     */
    SUNMatrixWrapper(const SUNMatrixWrapper &A, realtype droptol,
                     int sparsetype);

    /**
     * @brief Wrap existing SUNMatrix
     * @param mat
     */
    explicit SUNMatrixWrapper(SUNMatrix mat);

    ~SUNMatrixWrapper();

    /**
     * @brief Copy constructor
     * @param other
     */
    SUNMatrixWrapper(const SUNMatrixWrapper &other);

    /**
     * @brief Move constructor
     * @param other
     */
    SUNMatrixWrapper(SUNMatrixWrapper &&other);

    /**
     * @brief Copy assignment
     * @param other
     * @return
     */
    SUNMatrixWrapper &operator=(const SUNMatrixWrapper &other);

    /**
     * @brief Move assignment
     * @param other
     * @return
     */
    SUNMatrixWrapper &operator=(SUNMatrixWrapper &&other);

    /**
     * @brief Reallocate space for sparse matrix according to specified nnz
     * @param nnz new number of nonzero entries
     */
    void reallocate(sunindextype nnz);

    /**
     * @brief Reallocate space for sparse matrix to used space according to last entry in indexptrs
     */
    void realloc();

    /**
     * @brief Get the wrapped SUNMatrix
     * @return raw SunMatrix object
     * @note Even though the returned matrix_ pointer is const qualified, matrix_->content will not be const.
     * This is a shortcoming in the underlying C library, which we cannot address and it is not intended that
     * any of those values are modified externally. If matrix_->content is manipulated,
     * cpp:meth:SUNMatrixWrapper:`refresh` needs to be called.
     */
    SUNMatrix get() const;

    /**
     * @brief Get the number of rows
     * @return number of rows
     */
    sunindextype rows() const;

    /**
     * @brief Get the number of columns
     * @return number of columns
     */
    sunindextype columns() const;

    /**
     * @brief Get the number of specified non-zero elements (sparse matrices only)
     * @note value will be 0 before indexptrs are set.
     * @return number of nonzero entries
     */
    sunindextype num_nonzeros() const;

    /**
     * @brief Get the number of indexptrs that can be specified (sparse matrices only)
     * @return number of indexptrs
     */
    sunindextype num_indexptrs() const;

    /**
     * @brief Get the number of allocated data elements
     * @return number of allocated entries
     */
    sunindextype capacity() const;

    /**
     * @brief Get  raw data of a sparse matrix
     * @return pointer to first data entry
     */
    realtype *data();

    /**
     * @brief Get const raw data of a sparse matrix
     * @return pointer to first data entry
     */
    const realtype *data() const;

    /**
     * @brief Get data of a sparse matrix
     * @param idx data index
     * @return idx-th data entry
     */
    realtype get_data(sunindextype idx) const;

    /**
     * @brief Get data entry for a dense matrix
     * @param irow row
     * @param icol col
     * @return A(irow,icol)
     */
    realtype get_data(sunindextype irow, sunindextype icol) const;

    /**
     * @brief Set data entry for a sparse matrix
     * @param idx data index
     * @param data data for idx-th entry
     */
    void set_data(sunindextype idx, realtype data);

    /**
     * @brief Set data entry for a dense matrix
     * @param irow row
     * @param icol col
     * @param data data for idx-th entry
     */
    void set_data(sunindextype irow, sunindextype icol, realtype data);

    /**
     * @brief Get the index value of a sparse matrix
     * @param idx data index
     * @return row (CSC) or column (CSR) for idx-th data entry
     */
    sunindextype get_indexval(sunindextype idx) const;

    /**
     * @brief Set the index value of a sparse matrix
     * @param idx data index
     * @param val row (CSC) or column (CSR) for idx-th data entry
     */
    void set_indexval(sunindextype idx, sunindextype val);

    /**
     * @brief Set the index values of a sparse matrix
     * @param vals rows (CSC) or columns (CSR) for data entries
     */
    void set_indexvals(const gsl::span<const sunindextype> vals);

    /**
     * @brief Get the index pointer of a sparse matrix
     * @param ptr_idx pointer index
     * @return index where the ptr_idx-th column (CSC) or row (CSR) starts
     */
    sunindextype get_indexptr(sunindextype ptr_idx) const;

    /**
     * @brief Set the index pointer of a sparse matrix
     * @param ptr_idx pointer index
     * @param ptr data-index where the ptr_idx-th column (CSC) or row (CSR) starts
     */
    void set_indexptr(sunindextype ptr_idx, sunindextype ptr);

    /**
     * @brief Set the index pointers of a sparse matrix
     * @param ptrs starting data-indices where the columns (CSC) or rows (CSR) start
     */
    void set_indexptrs(const gsl::span<const sunindextype> ptrs);

    /**
     * @brief Get the type of sparse matrix
     * @return matrix type
     */
    int sparsetype() const;

    /**
     * @brief multiply with a scalar (in-place)
     * @param a scalar value to multiply matrix
     */
    void scale(realtype a);

    /**
     * @brief N_Vector interface for multiply
     * @param c output vector, may already contain values
     * @param b multiplication vector
     * @param alpha scalar coefficient for matrix
     */
    void multiply(N_Vector c, const_N_Vector b, realtype alpha = 1.0) const;

    /**
     * @brief Perform matrix vector multiplication c += alpha * A*b
     * @param c output vector, may already contain values
     * @param b multiplication vector
     * @param alpha scalar coefficient
     */
    void multiply(gsl::span<realtype> c, gsl::span<const realtype> b,
                  const realtype alpha = 1.0) const;

    /**
     * @brief Perform reordered matrix vector multiplication c += A[:,cols]*b
     * @param c output vector, may already contain values
     * @param b multiplication vector
     * @param cols int vector for column reordering
     * @param transpose bool transpose A before multiplication
     */
    void multiply(N_Vector c,
                  const_N_Vector b,
                  gsl::span <const int> cols,
                  bool transpose) const;

    /**
     * @brief Perform reordered matrix vector multiplication c += A[:,cols]*b
     * @param c output vector, may already contain values
     * @param b multiplication vector
     * @param cols int vector for column reordering
     * @param transpose bool transpose A before multiplication
     */
    void multiply(gsl::span<realtype> c,
                  gsl::span<const realtype> b,
                  gsl::span <const int> cols,
                  bool transpose) const;

    /**
     * @brief Perform matrix matrix multiplication C = A * B for sparse A, B, C
     * @param C output matrix,
     * @param B multiplication matrix
     * @note will overwrite existing data, indexptrs, indexvals for C, but will use preallocated space for these vars
     */
    void sparse_multiply(SUNMatrixWrapper &C,
                         const SUNMatrixWrapper &B) const;

    /**
     * @brief Perform sparse matrix matrix addition C = alpha * A +  beta * B
     * @param A addition matrix
     * @param alpha scalar A
     * @param B addition matrix
     * @param beta scalar B
     * @note will overwrite existing data, indexptrs, indexvals for C, but will use preallocated space for these vars
     */
    void sparse_add(const SUNMatrixWrapper &A, realtype alpha,
                    const SUNMatrixWrapper &B, realtype beta);

    /**
     * @brief Perform matrix-matrix addition A = sum(mats(0)...mats(len(mats)))
     * @param mats vector of sparse matrices
     * @note will overwrite existing data, indexptrs, indexvals for A, but will use preallocated space for these vars
     */
    void sparse_sum(const std::vector<SUNMatrixWrapper> &mats);

    /**
     * @brief Compute x = x + beta * A(:,k), where x is a dense vector and A(:,k) is sparse, and update
     * the sparsity pattern for C(:,j) if applicable
     *
     * This function currently has two purposes:
     *   - perform parts of sparse matrix-matrix multiplication C(:,j)=A(:,k)*B(k,j)
     *    enabled by passing beta=B(k,j), x=C(:,j), C=C, w=sparsity of C(:,j) from B(k,0...j-1), nnz=nnz(C(:,0...j-1)
     *   - add the k-th column of the sparse matrix A multiplied by beta to the dense vector x.
     *    enabled by passing beta=*, x=x, C=nullptr, w=nullptr, nnz=*
     *
     * @param k column index
     * @param beta scaling factor
     * @param w index workspace, (w[i]<mark) indicates non-zeroness of C(i,j) (dimension: m),
     * if this is a nullptr, sparsity pattern of C will not be updated (if applicable).
     * @param x dense output vector (dimension: m)
     * @param mark marker for w to indicate nonzero pattern
     * @param C sparse output matrix, if this is a nullptr, sparsity pattern of C will not be updated
     * @param nnz number of nonzeros that were already written to C
     * @return updated number of nonzeros in C
     */
    sunindextype scatter(const sunindextype k, const realtype beta,
                         sunindextype *w, gsl::span<realtype> x,
                         const sunindextype mark,
                         SUNMatrixWrapper *C, sunindextype nnz) const;

    /**
     * @brief Compute transpose A' of sparse matrix A and writes it to the matrix C = alpha * A'
     *
     * @param C output matrix (sparse or dense)
     * @param alpha scalar multiplier
     * @param blocksize blocksize for transposition. For full matrix transpose set to ncols/nrows
     */
    void transpose(SUNMatrixWrapper &C, const realtype alpha,
                   sunindextype blocksize) const;

    /**
     * @brief Writes a sparse matrix A to a dense matrix D.
     *
     * @param D dense output matrix
     */
    void to_dense(SUNMatrixWrapper &D) const;

    /**
     * @brief Writes the diagonal of sparse matrix A to a dense vector v.
     *
     * @param v dense outut vector
     */
    void to_diag(N_Vector v) const;

    /**
     * @brief Set to 0.0, for sparse matrices also resets indexptr/indexvals
     */
    void zero();

    /**
     * @brief Get matrix id
     * @return SUNMatrix_ID
     */
    SUNMatrix_ID matrix_id() const {return id_;};

    /**
     * @brief Update internal cache, needs to be called after external manipulation of matrix_->content
     */
    void refresh();

  private:

    /**
     * @brief SUNMatrix to which all methods are applied
     */
    SUNMatrix matrix_ {nullptr};

    /**
     * @brief cache for SUNMatrixGetId(matrix_)
     */
    SUNMatrix_ID id_ {SUNMATRIX_CUSTOM};

    /**
     * @brief cache for SUNMatrixGetId(matrix_)
     */
    int sparsetype_ {CSC_MAT};

    /**
     * @brief cache for SM_INDEXPTRS_S(matrix_)[SM_NP_S(matrix_)]
     */
    sunindextype num_nonzeros_ {0};
    /**
     * @brief cache for SM_NNZ_S(matrix_)
     */
    sunindextype capacity_ {0};

    /**
     * @brief cache for SM_DATA_S(matrix_)
     */
    realtype *data_ {nullptr};
    /**
     * @brief cache for SM_INDEXPTRS_S(matrix_)
     */
    sunindextype *indexptrs_ {nullptr};
    /**
     * @brief cache for SM_INDEXVALS_S(matrix_)
     */
    sunindextype *indexvals_ {nullptr};

    /**
     * @brief cache for SM_ROWS_X(matrix_)
     */
    sunindextype num_rows_ {0};
    /**
     * @brief cache for SM_COLUMS_X(matrix_)
     */
    sunindextype num_columns_ {0};
    /**
     * @brief cache for SM_NP_S(matrix_)
     */
    sunindextype num_indexptrs_ {0};

    /**
     * @brief call update_ptrs & update_size
     */
    void finish_init();
    /**
     * @brief update data_, indexptrs_, indexvals_ if applicable
     */
    void update_ptrs();
    /**
     * @brief update num_rows_, num_columns_, num_indexptrs if applicable
     */
    void update_size();
    /**
     * @brief indicator whether this wrapper allocated matrix_ and is responsible for deallocation
     */
    bool ownmat = true;
};

} // namespace amici

namespace gsl {
/**
 * @brief Create span from SUNMatrix
 * @param m SUNMatrix
 * @return Created span
 */
inline span<realtype> make_span(SUNMatrix m)
{
    switch (SUNMatGetID(m)) {
    case SUNMATRIX_DENSE:
        return span<realtype>(SM_DATA_D(m), SM_LDATA_D(m));
    case SUNMATRIX_SPARSE:
        return span<realtype>(SM_DATA_S(m), SM_NNZ_S(m));
    default:
        throw amici::AmiException("Unimplemented SUNMatrix type for make_span");
    }
}
} // namespace gsl

#endif // AMICI_SUNDIALS_MATRIX_WRAPPER_H

max_stars_repo_name: PaulJonasJost/AMICI

max_issues_repo_name: PaulJonasJost/AMICI

max_forks_repo_name: PaulJonasJost/AMICI


text:
#include <stdio.h>

#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_math.h>
#include "timer.c"
#include "timer.h"
#include "integrand.h"

double dipole_approx (double r);

double gaussian (double *x, int dim);

int main (void)
{
    double res, err;

    size_t dim = 6;
    double x1[] = { 0., 0., 0., 0., 0., 0., };
    double xu[] = { 1., 1., 1., 1., 1., 1., };
    double distmin = 1.001;
    double distmax = 4.;
    double dist;
    int numpoints = 20;
    double nt = (distmax - distmin) / (numpoints - 1);
    double vegas[20], dipole[20], distance[20];
//vegas integration
    gsl_rng *r = gsl_rng_alloc (gsl_rng_taus2);
    unsigned long seed = 1UL;

    gsl_rng_set (r, seed);

    size_t calls = 1000000;

    dist = distmin;

    gsl_monte_function G = { &g, dim, &dist };

    gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (dim);

    gsl_monte_vegas_init (sv);

    timer_start ();

    for (int i = 0; i < numpoints; i++){
        gsl_monte_vegas_integrate (&G, x1, xu, dim, calls / 5, r, sv, &res,
            &err);
        do
        {
            gsl_monte_vegas_integrate (&G, x1, xu, dim, calls, r, sv, &res,
                &err);
	    fflush(stdout);
        }
        while (fabs (gsl_monte_vegas_chisq (sv) - 1.0) > 0.2);
        dist += nt;
        vegas[i] = res;
    distance[i] = dist;
    dipole[i] = -2. / pow (dist, 3.);
    }
    timer_stop();
    //double vegastime = timer_stop();
     gsl_monte_vegas_free (sv);

//Monte carlo integration
    double sum;
    double x[6];

    long i, j, nn;
    nn = 1000000;

    timer_start ();

    double monte[20];

    dist = distmin;
    for (j = 0; j < numpoints; j++)
    {
        sum = 0.;
        for (i = 0; i < nn; i++)
        {
            for (int k = 0; k < (int) dim; k++)
            {
                x[k] = gsl_rng_uniform (r);
            }
            sum += g (x, dim, &dist);
        }
        res = sum/nn;
	fflush(stdout);
        dist += nt;
	monte[j] = res;
    }
    timer_stop ();
    //double montetime = timer_stop();
    //printf("Vegas time: %f Monte time: %f", vegastime, montetime);

    gsl_rng_free(r);

    double monteerr = 0.0;

    for (int jj = 0; jj < numpoints; jj++)
    {
    monteerr += fabs(monte[jj] - vegas[jj]);
    }
//prints the output of each integration
   printf("#    Dist               Vegas           Monte       Dipolapprox\n");
    for( int l = 0; l < numpoints; l++)
    {
    double dd = distance[l];
    double vv = fabs(vegas[l]);
    double mm = fabs(monte[l]);
    double di = fabs(dipole[l]);
    printf("   %.6f         %.6f       %.6f      %.6f\n", dd, vv, mm, di);
    }

    return 0;
}