text:
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 Microsoft Corporation. All rights reserved.
//
// This code is licensed under the MIT License (MIT).
//
// 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.
//
///////////////////////////////////////////////////////////////////////////////

#ifndef GSL_STRING_SPAN_H
#define GSL_STRING_SPAN_H

#include <gsl/gsl_assert> // for Ensures, Expects
#include <gsl/gsl_util>   // for narrow_cast
#include <gsl/span>       // for operator!=, operator==, dynamic_extent
#include <gsl/pointers>   // for not_null

#include <algorithm> // for equal, lexicographical_compare
#include <array>     // for array
#include <cstddef>   // for ptrdiff_t, size_t, nullptr_t
#include <cstdint>   // for PTRDIFF_MAX
#include <cstring>
#include <string>      // for basic_string, allocator, char_traits
#include <type_traits> // for declval, is_convertible, enable_if_t, add_...

#ifdef _MSC_VER
#pragma warning(push)

// Turn MSVC /analyze rules that generate too much noise. TODO: fix in the tool.
#pragma warning(disable : 26446) // TODO: bug in parser - attributes and templates
#pragma warning(disable : 26481) // TODO: suppress does not work inside templates sometimes

#if _MSC_VER < 1910
#pragma push_macro("constexpr")
#define constexpr /*constexpr*/

#endif // _MSC_VER < 1910
#endif // _MSC_VER

namespace gsl
{
//
// czstring and wzstring
//
// These are "tag" typedefs for C-style strings (i.e. null-terminated character arrays)
// that allow static analysis to help find bugs.
//
// There are no additional features/semantics that we can find a way to add inside the
// type system for these types that will not either incur significant runtime costs or
// (sometimes needlessly) break existing programs when introduced.
//

template <typename CharT, std::ptrdiff_t Extent = dynamic_extent>
using basic_zstring = CharT*;

template <std::ptrdiff_t Extent = dynamic_extent>
using czstring = basic_zstring<const char, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cwzstring = basic_zstring<const wchar_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cu16zstring = basic_zstring<const char16_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cu32zstring = basic_zstring<const char32_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using zstring = basic_zstring<char, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using wzstring = basic_zstring<wchar_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using u16zstring = basic_zstring<char16_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using u32zstring = basic_zstring<char32_t, Extent>;

namespace details
{
    template <class CharT>
    std::ptrdiff_t string_length(const CharT* str, std::ptrdiff_t n)
    {
        if (str == nullptr || n <= 0) return 0;

        const span<const CharT> str_span{str, n};

        std::ptrdiff_t len = 0;
        while (len < n && str_span[len]) len++;

        return len;
    }
} // namespace details

//
// ensure_sentinel()
//
// Provides a way to obtain an span from a contiguous sequence
// that ends with a (non-inclusive) sentinel value.
//
// Will fail-fast if sentinel cannot be found before max elements are examined.
//
template <typename T, const T Sentinel>
span<T, dynamic_extent> ensure_sentinel(T* seq, std::ptrdiff_t max = PTRDIFF_MAX)
{
    Ensures(seq != nullptr);

    GSL_SUPPRESS(f.23) // NO-FORMAT: attribute // TODO: false positive // TODO: suppress does not work
    auto cur = seq;
    Ensures(cur != nullptr); // workaround for removing the warning

    GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute // TODO: suppress does not work
    while ((cur - seq) < max && *cur != Sentinel) ++cur;
    Ensures(*cur == Sentinel);
    return {seq, cur - seq};
}

//
// ensure_z - creates a span for a zero terminated strings.
// Will fail fast if a null-terminator cannot be found before
// the limit of size_type.
//
template <typename CharT>
span<CharT, dynamic_extent> ensure_z(CharT* const& sz, std::ptrdiff_t max = PTRDIFF_MAX)
{
    return ensure_sentinel<CharT, CharT(0)>(sz, max);
}

template <typename CharT, std::size_t N>
span<CharT, dynamic_extent> ensure_z(CharT (&sz)[N])
{
    return ensure_z(&sz[0], narrow_cast<std::ptrdiff_t>(N));
}

template <class Cont>
span<typename std::remove_pointer<typename Cont::pointer>::type, dynamic_extent>
ensure_z(Cont& cont)
{
    return ensure_z(cont.data(), narrow_cast<std::ptrdiff_t>(cont.size()));
}

template <typename CharT, std::ptrdiff_t>
class basic_string_span;

namespace details {
    template <typename T>
    struct is_basic_string_span_oracle : std::false_type
    {
    };

    template <typename CharT, std::ptrdiff_t Extent>
    struct is_basic_string_span_oracle<basic_string_span<CharT, Extent>> : std::true_type
    {
    };

    template <typename T>
    struct is_basic_string_span : is_basic_string_span_oracle<std::remove_cv_t<T>>
    {
    };
} // namespace details

//
// string_span and relatives
//
template <typename CharT, std::ptrdiff_t Extent = dynamic_extent>
class basic_string_span
{
public:
    using element_type = CharT;
    using pointer = std::add_pointer_t<element_type>;
    using reference = std::add_lvalue_reference_t<element_type>;
    using const_reference = std::add_lvalue_reference_t<std::add_const_t<element_type>>;
    using impl_type = span<element_type, Extent>;

    using index_type = typename impl_type::index_type;
    using iterator = typename impl_type::iterator;
    using const_iterator = typename impl_type::const_iterator;
    using reverse_iterator = typename impl_type::reverse_iterator;
    using const_reverse_iterator = typename impl_type::const_reverse_iterator;

    // default (empty)
    constexpr basic_string_span() noexcept = default;

    // copy
    constexpr basic_string_span(const basic_string_span& other) noexcept = default;

    // assign
    constexpr basic_string_span& operator=(const basic_string_span& other) noexcept = default;

    constexpr basic_string_span(pointer ptr, index_type length) : span_(ptr, length) {}
    constexpr basic_string_span(pointer firstElem, pointer lastElem) : span_(firstElem, lastElem) {}

    // From static arrays - if 0-terminated, remove 0 from the view
    // All other containers allow 0s within the length, so we do not remove them
    template <std::size_t N>
    constexpr basic_string_span(element_type (&arr)[N]) : span_(remove_z(arr))
    {}

    template <std::size_t N, class ArrayElementType = std::remove_const_t<element_type>>
    constexpr basic_string_span(std::array<ArrayElementType, N>& arr) noexcept : span_(arr)
    {}

    template <std::size_t N, class ArrayElementType = std::remove_const_t<element_type>>
    constexpr basic_string_span(const std::array<ArrayElementType, N>& arr) noexcept : span_(arr)
    {}

    // Container signature should work for basic_string after C++17 version exists
    template <class Traits, class Allocator>
    // GSL_SUPPRESS(bounds.4) // NO-FORMAT: attribute // TODO: parser bug
    constexpr basic_string_span(std::basic_string<element_type, Traits, Allocator>& str)
        : span_(&str[0], narrow_cast<std::ptrdiff_t>(str.length()))
    {}

    template <class Traits, class Allocator>
    constexpr basic_string_span(const std::basic_string<element_type, Traits, Allocator>& str)
        : span_(&str[0], str.length())
    {}

    // from containers. Containers must have a pointer type and data() function signatures
    template <class Container,
              class = std::enable_if_t<
                  !details::is_basic_string_span<Container>::value &&
                  std::is_convertible<typename Container::pointer, pointer>::value &&
                  std::is_convertible<typename Container::pointer,
                                      decltype(std::declval<Container>().data())>::value>>
    constexpr basic_string_span(Container& cont) : span_(cont)
    {}

    template <class Container,
              class = std::enable_if_t<
                  !details::is_basic_string_span<Container>::value &&
                  std::is_convertible<typename Container::pointer, pointer>::value &&
                  std::is_convertible<typename Container::pointer,
                                      decltype(std::declval<Container>().data())>::value>>
    constexpr basic_string_span(const Container& cont) : span_(cont)
    {}

    // from string_span
    template <
        class OtherValueType, std::ptrdiff_t OtherExtent,
        class = std::enable_if_t<std::is_convertible<
            typename basic_string_span<OtherValueType, OtherExtent>::impl_type, impl_type>::value>>
    constexpr basic_string_span(basic_string_span<OtherValueType, OtherExtent> other)
        : span_(other.data(), other.length())
    {}

    template <index_type Count>
    constexpr basic_string_span<element_type, Count> first() const
    {
        return {span_.template first<Count>()};
    }

    constexpr basic_string_span<element_type, dynamic_extent> first(index_type count) const
    {
        return {span_.first(count)};
    }

    template <index_type Count>
    constexpr basic_string_span<element_type, Count> last() const
    {
        return {span_.template last<Count>()};
    }

    constexpr basic_string_span<element_type, dynamic_extent> last(index_type count) const
    {
        return {span_.last(count)};
    }

    template <index_type Offset, index_type Count>
    constexpr basic_string_span<element_type, Count> subspan() const
    {
        return {span_.template subspan<Offset, Count>()};
    }

    constexpr basic_string_span<element_type, dynamic_extent>
    subspan(index_type offset, index_type count = dynamic_extent) const
    {
        return {span_.subspan(offset, count)};
    }

    constexpr reference operator[](index_type idx) const { return span_[idx]; }
    constexpr reference operator()(index_type idx) const { return span_[idx]; }

    constexpr pointer data() const { return span_.data(); }

    constexpr index_type length() const noexcept { return span_.size(); }
    constexpr index_type size() const noexcept { return span_.size(); }
    constexpr index_type size_bytes() const noexcept { return span_.size_bytes(); }
    constexpr index_type length_bytes() const noexcept { return span_.length_bytes(); }
    constexpr bool empty() const noexcept { return size() == 0; }

    constexpr iterator begin() const noexcept { return span_.begin(); }
    constexpr iterator end() const noexcept { return span_.end(); }

    constexpr const_iterator cbegin() const noexcept { return span_.cbegin(); }
    constexpr const_iterator cend() const noexcept { return span_.cend(); }

    constexpr reverse_iterator rbegin() const noexcept { return span_.rbegin(); }
    constexpr reverse_iterator rend() const noexcept { return span_.rend(); }

    constexpr const_reverse_iterator crbegin() const noexcept { return span_.crbegin(); }
    constexpr const_reverse_iterator crend() const noexcept { return span_.crend(); }

private:
    static impl_type remove_z(pointer const& sz, std::ptrdiff_t max)
    {
        return {sz, details::string_length(sz, max)};
    }

    template <std::size_t N>
    static impl_type remove_z(element_type (&sz)[N])
    {
        return remove_z(&sz[0], narrow_cast<std::ptrdiff_t>(N));
    }

    impl_type span_;
};

template <std::ptrdiff_t Extent = dynamic_extent>
using string_span = basic_string_span<char, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cstring_span = basic_string_span<const char, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using wstring_span = basic_string_span<wchar_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cwstring_span = basic_string_span<const wchar_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using u16string_span = basic_string_span<char16_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cu16string_span = basic_string_span<const char16_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using u32string_span = basic_string_span<char32_t, Extent>;

template <std::ptrdiff_t Extent = dynamic_extent>
using cu32string_span = basic_string_span<const char32_t, Extent>;

//
// to_string() allow (explicit) conversions from string_span to string
//

template <typename CharT, std::ptrdiff_t Extent>
std::basic_string<typename std::remove_const<CharT>::type>
to_string(basic_string_span<CharT, Extent> view)
{
    return {view.data(), narrow_cast<std::size_t>(view.length())};
}

template <typename CharT, typename Traits = typename std::char_traits<CharT>,
          typename Allocator = std::allocator<CharT>, typename gCharT, std::ptrdiff_t Extent>
std::basic_string<CharT, Traits, Allocator> to_basic_string(basic_string_span<gCharT, Extent> view)
{
    return {view.data(), narrow_cast<std::size_t>(view.length())};
}

template <class ElementType, std::ptrdiff_t Extent>
basic_string_span<const byte, details::calculate_byte_size<ElementType, Extent>::value>
as_bytes(basic_string_span<ElementType, Extent> s) noexcept
{
    GSL_SUPPRESS(type.1) // NO-FORMAT: attribute
    return {reinterpret_cast<const byte*>(s.data()), s.size_bytes()};
}

template <class ElementType, std::ptrdiff_t Extent,
          class = std::enable_if_t<!std::is_const<ElementType>::value>>
basic_string_span<byte, details::calculate_byte_size<ElementType, Extent>::value>
as_writeable_bytes(basic_string_span<ElementType, Extent> s) noexcept
{
    GSL_SUPPRESS(type.1) // NO-FORMAT: attribute
    return {reinterpret_cast<byte*>(s.data()), s.size_bytes()};
}

// zero-terminated string span, used to convert
// zero-terminated spans to legacy strings
template <typename CharT, std::ptrdiff_t Extent = dynamic_extent>
class basic_zstring_span {
public:
    using value_type = CharT;
    using const_value_type = std::add_const_t<CharT>;

    using pointer = std::add_pointer_t<value_type>;
    using const_pointer = std::add_pointer_t<const_value_type>;

    using zstring_type = basic_zstring<value_type, Extent>;
    using const_zstring_type = basic_zstring<const_value_type, Extent>;

    using impl_type = span<value_type, Extent>;
    using string_span_type = basic_string_span<value_type, Extent>;

    constexpr basic_zstring_span(impl_type s) : span_(s)
    {
        // expects a zero-terminated span
        Expects(s[s.size() - 1] == '\0');
    }

    // copy
    constexpr basic_zstring_span(const basic_zstring_span& other) = default;

    // move
    constexpr basic_zstring_span(basic_zstring_span&& other) = default;

    // assign
    constexpr basic_zstring_span& operator=(const basic_zstring_span& other) = default;

    // move assign
    constexpr basic_zstring_span& operator=(basic_zstring_span&& other) = default;

    constexpr bool empty() const noexcept { return span_.size() == 0; }

    constexpr string_span_type as_string_span() const noexcept
    {
        const auto sz = span_.size();
        return {span_.data(), sz > 1 ? sz - 1 : 0};
    }
    constexpr string_span_type ensure_z() const { return gsl::ensure_z(span_); }

    constexpr const_zstring_type assume_z() const noexcept { return span_.data(); }

private:
    impl_type span_;
};

template <std::ptrdiff_t Max = dynamic_extent>
using zstring_span = basic_zstring_span<char, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using wzstring_span = basic_zstring_span<wchar_t, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using u16zstring_span = basic_zstring_span<char16_t, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using u32zstring_span = basic_zstring_span<char32_t, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using czstring_span = basic_zstring_span<const char, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using cwzstring_span = basic_zstring_span<const wchar_t, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using cu16zstring_span = basic_zstring_span<const char16_t, Max>;

template <std::ptrdiff_t Max = dynamic_extent>
using cu32zstring_span = basic_zstring_span<const char32_t, Max>;

// operator ==
template <class CharT, std::ptrdiff_t Extent, class T,
          class = std::enable_if_t<
              details::is_basic_string_span<T>::value ||
              std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>>>::value>>
bool operator==(const gsl::basic_string_span<CharT, Extent>& one, const T& other)
{
    const gsl::basic_string_span<std::add_const_t<CharT>> tmp(other);
    return std::equal(one.begin(), one.end(), tmp.begin(), tmp.end());
}

template <class CharT, std::ptrdiff_t Extent, class T,
          class = std::enable_if_t<
              !details::is_basic_string_span<T>::value &&
              std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>>>::value>>
bool operator==(const T& one, const gsl::basic_string_span<CharT, Extent>& other)
{
    const gsl::basic_string_span<std::add_const_t<CharT>> tmp(one);
    return std::equal(tmp.begin(), tmp.end(), other.begin(), other.end());
}

// operator !=
template <typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
          typename = std::enable_if_t<std::is_convertible<
              T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value>>
bool operator!=(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return !(one == other);
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename = std::enable_if_t<
        std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value &&
        !gsl::details::is_basic_string_span<T>::value>>
bool operator!=(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return !(one == other);
}

// operator<
template <typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
          typename = std::enable_if_t<std::is_convertible<
              T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value>>
bool operator<(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    const gsl::basic_string_span<std::add_const_t<CharT>, Extent> tmp(other);
    return std::lexicographical_compare(one.begin(), one.end(), tmp.begin(), tmp.end());
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename = std::enable_if_t<
        std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value &&
        !gsl::details::is_basic_string_span<T>::value>>
bool operator<(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    gsl::basic_string_span<std::add_const_t<CharT>, Extent> tmp(one);
    return std::lexicographical_compare(tmp.begin(), tmp.end(), other.begin(), other.end());
}

#ifndef _MSC_VER

// VS treats temp and const containers as convertible to basic_string_span,
// so the cases below are already covered by the previous operators

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator<(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    gsl::basic_string_span<std::add_const_t<CharT>, Extent> tmp(other);
    return std::lexicographical_compare(one.begin(), one.end(), tmp.begin(), tmp.end());
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator<(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    gsl::basic_string_span<std::add_const_t<CharT>, Extent> tmp(one);
    return std::lexicographical_compare(tmp.begin(), tmp.end(), other.begin(), other.end());
}
#endif

// operator <=
template <typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
          typename = std::enable_if_t<std::is_convertible<
              T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value>>
bool operator<=(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return !(other < one);
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename = std::enable_if_t<
        std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value &&
        !gsl::details::is_basic_string_span<T>::value>>
bool operator<=(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return !(other < one);
}

#ifndef _MSC_VER

// VS treats temp and const containers as convertible to basic_string_span,
// so the cases below are already covered by the previous operators

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator<=(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return !(other < one);
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator<=(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return !(other < one);
}
#endif

// operator>
template <typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
          typename = std::enable_if_t<std::is_convertible<
              T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value>>
bool operator>(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return other < one;
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename = std::enable_if_t<
        std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value &&
        !gsl::details::is_basic_string_span<T>::value>>
bool operator>(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return other < one;
}

#ifndef _MSC_VER

// VS treats temp and const containers as convertible to basic_string_span,
// so the cases below are already covered by the previous operators

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator>(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return other < one;
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator>(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return other < one;
}
#endif

// operator >=
template <typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
          typename = std::enable_if_t<std::is_convertible<
              T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value>>
bool operator>=(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return !(one < other);
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename = std::enable_if_t<
        std::is_convertible<T, gsl::basic_string_span<std::add_const_t<CharT>, Extent>>::value &&
        !gsl::details::is_basic_string_span<T>::value>>
bool operator>=(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return !(one < other);
}

#ifndef _MSC_VER

// VS treats temp and const containers as convertible to basic_string_span,
// so the cases below are already covered by the previous operators

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator>=(gsl::basic_string_span<CharT, Extent> one, const T& other)
{
    return !(one < other);
}

template <
    typename CharT, std::ptrdiff_t Extent = gsl::dynamic_extent, typename T,
    typename DataType = typename T::value_type,
    typename = std::enable_if_t<
        !gsl::details::is_span<T>::value && !gsl::details::is_basic_string_span<T>::value &&
        std::is_convertible<DataType*, CharT*>::value &&
        std::is_same<std::decay_t<decltype(std::declval<T>().size(), *std::declval<T>().data())>,
                     DataType>::value>>
bool operator>=(const T& one, gsl::basic_string_span<CharT, Extent> other)
{
    return !(one < other);
}
#endif
} // namespace gsl

#ifdef _MSC_VER
#pragma warning(pop)

#if _MSC_VER < 1910
#undef constexpr
#pragma pop_macro("constexpr")

#endif // _MSC_VER < 1910
#endif // _MSC_VER

#endif // GSL_STRING_SPAN_H

ext: h

max_stars_repo_name: Redchards/CVRP

max_issues_repo_name: Redchards/CVRP

max_forks_repo_name: Redchards/CVRP


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "DirectionalLight.h"

namespace Library
{
	class ProxyModel;
}

namespace Rendering
{
	class BlinnPhongMaterial;

	class BlinnPhongDemo final : public Library::DrawableGameComponent
	{
	public:
		BlinnPhongDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		BlinnPhongDemo(const BlinnPhongDemo&) = delete;
		BlinnPhongDemo(BlinnPhongDemo&&) = default;
		BlinnPhongDemo& operator=(const BlinnPhongDemo&) = default;
		BlinnPhongDemo& operator=(BlinnPhongDemo&&) = default;
		~BlinnPhongDemo();

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);
		void ToggleAnimation();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float DirectionalLightIntensity() const;
		void SetDirectionalLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightDirection() const;
		void RotateDirectionalLight(const DirectX::XMFLOAT2& amount);

		float SpecularIntensity() const;
		void SetSpecularIntensity(float intensity);

		float SpecularPower() const;
		void SetSpecularPower(float power);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<BlinnPhongMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		Library::DirectionalLight mDirectionalLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		float mModelRotationAngle{ 0.0f };
		bool mAnimationEnabled{ true };
		bool mUpdateMaterial{ true };
	};
}


text:
/*
 * Copyright (c) 2016-2020 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * 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

#ifndef basis_40d2b174_3729_49b6_a8d6_5b7ac75a2483_h
#define basis_40d2b174_3729_49b6_a8d6_5b7ac75a2483_h

#include <gslib/std.h>
#include <gslib/string.h>
#include <gslib/tree.h>
#include <ariel/widget.h>

namespace ui_editor {

using gs::rect;
using gs::rectf;
using gs::point;
using gs::pointf;
using gs::string;
using gs::list;
using gs::unordered_map;
using gs::gchar;
using gs::uint;
using gs::real;
using gs::real32;
using gs::ariel::painter;
using gs::ariel::wsys_manager;
using gs::ariel::unikey;
typedef gs::ariel::widget core_widget;

struct ui_proc_node;
typedef unordered_map<string, string> ui_prop_map;
typedef unordered_map<string, ui_proc_node> ui_proc_map;

struct ui_proc_node
{
    string              args;
    string              injection;
};

// class ui_editor_configs
// {
//
// };

struct ui_macro
{
    string              name;
    string              definition;

public:
    ui_macro(const string& n): name(n) {}
    ui_macro(const string& n, const string& d): name(n), definition(d) {}
};

typedef list<string> ui_includes;
typedef list<ui_macro> ui_macros;

class ui_src_entry
{
public:
    void add_include(const string& str);
    void add_include_group(const string& strgroup);
    void add_macro(const string& name) { _macro_list.push_back(ui_macro(name)); }
    void add_macro(const string& name, const string& def) { _macro_list.push_back(ui_macro(name, def)); }
    void set_preproc(const string& preproc) { _preproc_injections = preproc; }
    void set_postproc(const string& postproc) { _postproc_injections = postproc; }
    void set_output(const string& str) { _output_name = str; }
    void set_export(const string& str) { _export_name = str; }
    const ui_includes& get_includes() const { return _include_list; }
    const ui_macros& get_macros() const { return _macro_list; }
    const string& get_preproc() const { return _preproc_injections; }
    const string& get_postproc() const { return _postproc_injections; }
    const string& get_output() const { return _output_name; }
    const string& get_export() const { return _export_name; }

protected:
    ui_includes         _include_list;
    ui_macros           _macro_list;
    string              _preproc_injections;
    string              _postproc_injections;
    string              _output_name;       /* for output file name */
    string              _export_name;       /* for export class name */
};

class ui_editor_context
{
public:
    typedef unordered_map<string, int> var_index_map;

public:
    ui_src_entry& get_src_entry() { return _src_entry; }
    const ui_src_entry& const_src_entry() const { return _src_entry; }
    int acquire_variable_index(const string& vartype);
    int query_variable_index(const string& vartype) const;
    int acquire_variable_name(string& name, const string& vartype);

protected:
    ui_src_entry        _src_entry;
    var_index_map       _index_map;
};

class ui_node
{
public:
    rect                position;
    string              type;
    string              name;
    string              variable;
    bool                is_instance;
    ui_prop_map         prop_map;
    ui_proc_map         proc_map;
    core_widget*        assoc_widget;

public:
    ui_node(): is_instance(false), assoc_widget(nullptr) {}
};

typedef gs::_treenode_cpy_wrapper<ui_node> ui_node_wrapper;
typedef gs::tree<ui_node, ui_node_wrapper> ui_tree;

typedef core_widget* (__stdcall *pfn_create_widget)(wsys_manager* wsys, core_widget* parent, ui_node& uinode);

class ui_transparent_layer:
    public core_widget
{
public:
    ui_transparent_layer(wsys_manager* m);
    void install_reflection_widget(core_widget* refw);
    core_widget* get_reflection_widget() const { return _reflection_widget; }

public:
    virtual ~ui_transparent_layer();
    virtual void on_reflection_close();
    virtual void on_reflection_show(bool b);
    virtual void on_reflection_enable(bool b);
    virtual void on_reflection_move(const rect& rc);
    virtual void on_reflection_char(uint um, uint ch) {}
    virtual void on_reflection_caret() {}
    virtual void on_reflection_capture(bool b) {}
    virtual void on_reflection_focus(bool b) {}
    virtual void on_reflection_scroll(const point& pt, real32 scr, bool vert) {}
    virtual void on_reflection_accelerator(unikey key, uint mask) {}

protected:
    core_widget*        _reflection_widget;

protected:
    void sync_size_with_reflection();
};

class ui_edit_control;

class ui_editor_layer:
    public ui_transparent_layer
{
public:
    enum size_orient
    {
        so_none,
        so_n,
        so_s,
        so_e,
        so_w,
        so_ne,
        so_nw,
        so_se,
        so_sw,
    };

public:
    ui_editor_layer(wsys_manager* m);
    void initialize(ui_node& uinode);
    void install_reflection_widget(core_widget* refw);
    ui_node* get_ui_node() const { return _assoc_ui_node; }

public:
    virtual void on_reflection_press(uint um, unikey uk, const point& pt);
    virtual void on_reflection_click(uint um, unikey uk, const point& pt);
    virtual void on_reflection_hover(uint um, const point& pt);
    virtual void on_reflection_leave(uint um, const point& pt);
    virtual void on_reflection_keydown(uint um, unikey uk);
    virtual void on_reflection_keyup(uint um, unikey uk);

protected:
    ui_node*            _assoc_ui_node;
    ui_edit_control*    _edit_ctl;

protected:
    void refresh_cursor(const point& pt);
    size_orient get_size_orient(const point& pt, int accuracy = 3) const;
    void try_install_edit_control(const point& pt);
    void try_uninstall_edit_control();
    void refresh_edit_control(uint um, const point& pt);
};

extern bool ui_create_widgets(wsys_manager* wsys, core_widget* root_layer, ui_tree& uicore);
//extern bool ui_;

}

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/**
 * @file tracker.h
 * @brief Beamforming and speaker tracking with a spherical microphone array.
 * @author John McDonough
 */
#ifndef TRACKER_H
#define TRACKER_H

#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_complex_math.h>

#include "common/refcount.h"
#include "common/jexception.h"

#include "stream/stream.h"
#include "beamformer/spectralinfoarray.h"
#include "modulated/modulated.h"
#include "aec/aec.h"

// ----- definition of class 'BaseDecomposition' -----
//
class BaseDecomposition : public Countable {
public:
  class SubbandEntry {
  public:
    SubbandEntry(unsigned subbandX, const gsl_complex& Bkl);

    unsigned subbandX() const { return subbandX_; }
    gsl_complex bkl()   const { return bkl_; }

    void* operator new(size_t sz) { return memoryManager().new_elem(); }
    void  operator delete(void* e) { memoryManager().delete_elem(e); }

    static MemoryManager<SubbandEntry>& memoryManager();

  private:
    const unsigned			subbandX_;
    const gsl_complex			bkl_;
  };

  class GreaterThan {
  public:
    bool operator()(SubbandEntry* sbX1, SubbandEntry* sbX2) {
      return (gsl_complex_abs(sbX1->bkl()) > gsl_complex_abs(sbX2->bkl()));
    }
  };

  class Iterator;

  class SubbandList : public Countable {
  public:
    SubbandList(const gsl_vector_complex* bkl, unsigned useSubbandsN = 0);
    ~SubbandList();

    unsigned useSubbandsN() const { return useSubbandsN_; }
    SubbandEntry** subbands() { return subbands_; }

    friend class Iterator;

  private:
    const unsigned			subbandsN_;
    const unsigned			useSubbandsN_;
    SubbandEntry**			subbands_;
  };

  typedef refcountable_ptr<SubbandList> SubbandListPtr;

  class Iterator {
  public:
    Iterator(const SubbandListPtr& subbandList)
      : subbandX_(0), useSubbandsN_(subbandList->useSubbandsN()), subbands_(subbandList->subbands()) { }
    void operator++(int) { subbandX_++; }
    bool more()          { return subbandX_ < useSubbandsN_; }
    const SubbandEntry& operator*() { return *(subbands_[subbandX_]); }

  private:
    unsigned				subbandX_;
    const unsigned			useSubbandsN_;
    SubbandEntry**			subbands_;
  };

 public:
  BaseDecomposition(unsigned orderN, unsigned subbandsN, double a, double sampleRate, unsigned useSubbandsN = 0, bool spatial = false);
  ~BaseDecomposition();

  unsigned orderN()         const { return orderN_;         }
  unsigned modesN()         const { return modesN_;         }
  unsigned subbandsN2()     const { return subbandsN2_;     }
  unsigned subbandsN()      const { return subbandsN_;      }
  unsigned useSubbandsN()   const { return useSubbandsN_;   }
  unsigned subbandLengthN() const { return subbandLengthN_; }

  SubbandListPtr subbandList()  { return subbandList_;    }

  virtual void reset();

  static gsl_complex harmonic(int order, int degree, double theta, double phi);
  gsl_complex harmonic(int order, int degree, unsigned channelX) const;
  static gsl_complex harmonic_deriv_polar_angle(int order, int degree, double theta, double phi);
  static gsl_complex harmonic_deriv_azimuth(int order, int degree, double theta, double phi);
  static gsl_complex modal_coefficient(unsigned order, double ka);
  gsl_complex modal_coefficient(unsigned order, unsigned subbandX) const;
  virtual void transform(const gsl_vector_complex* initial, gsl_vector_complex* transformed) { }
  virtual void estimate_Bkl(double theta, double phi, const gsl_vector_complex* snapshot, unsigned subbandX) = 0;
  virtual void calculate_gkl(double theta, double phi, unsigned subbandX) = 0;
  virtual const gsl_matrix_complex* linearize(gsl_vector* xk, int frame_no) = 0;
  virtual const gsl_vector_complex* predicted_observation(gsl_vector* xk, int frame_no) = 0;

#ifdef ENABLE_LEGACY_BTK_API
  static gsl_complex harmonicDerivPolarAngle(int order, int degree, double theta, double phi){
    return harmonic_deriv_polar_angle(order, degree, theta, phi);
  }
  static gsl_complex harmonicDerivAzimuth(int order, int degree, double theta, double phi){
    return harmonic_deriv_azimuth(order, degree, theta, phi);
  }
  static gsl_complex modalCoefficient(unsigned order, double ka){
    return modal_coefficient(order, ka);
  }
#endif

  static const unsigned				StateN_;
  static const unsigned				ChannelsN_;
  static const gsl_complex 			ComplexZero_;
  static const gsl_complex 			ComplexOne_;
  static const double				SpeedOfSound_;

 protected:
  static gsl_complex calc_in_(int n);
  void set_eigenmike_geometry_();
  gsl_complex calc_Gnm_(unsigned subbandX, int n, int m, double theta, double phi);
  gsl_complex calc_dGnm_dtheta_(unsigned subbandX, int n, int m, double theta, double phi);
  static double calc_Pnm_(int order, int degree, double theta);
  static double calc_dPnm_dtheta_(int n, int m, double theta);
  static double calc_normalization_(int order, int degree);

  const unsigned				orderN_;
  const unsigned				modesN_;
  const unsigned				subbandsN_;
  const unsigned				subbandsN2_;
  const unsigned				useSubbandsN_;
  const unsigned				subbandLengthN_;
  const double					samplerate_;
  const double					a_;
  gsl_vector_complex**				bn_;
  gsl_vector*					theta_s_;
  gsl_vector*					phi_s_;
  gsl_vector_complex**				spherical_component_;
  gsl_vector_complex*				bkl_;
  gsl_vector_complex*				dbkl_dtheta_;
  gsl_vector_complex*				dbkl_dphi_;
  gsl_vector_complex**				gkl_;
  gsl_vector_complex**				dgkl_dtheta_;
  gsl_vector_complex**				dgkl_dphi_;
  gsl_vector_complex*				vkl_;
  gsl_matrix_complex*				Hbar_k_;
  gsl_vector_complex*				yhat_k_;

  SubbandListPtr				subbandList_;
};

typedef refcountable_ptr<BaseDecomposition> BaseDecompositionPtr;

// ----- definition of class 'ModalDecomposition' -----
//
class ModalDecomposition : public BaseDecomposition {
public:
  ModalDecomposition(unsigned orderN, unsigned subbandsN, double a, double sampleRate, unsigned useSubbandsN = 0);
  ~ModalDecomposition() { }

  virtual void estimate_Bkl(double theta, double phi, const gsl_vector_complex* snapshot, unsigned subbandX);
  virtual void transform(const gsl_vector_complex* initial, gsl_vector_complex* transformed);
  virtual const gsl_matrix_complex* linearize(gsl_vector* xk, int frame_no);
  virtual const gsl_vector_complex* predicted_observation(gsl_vector* xk, int frame_no);
  virtual void calculate_gkl(double theta, double phi, unsigned subbandX);
};

typedef Inherit<ModalDecomposition, BaseDecompositionPtr>	ModalDecompositionPtr;

// ----- definition of class 'SpatialDecomposition' -----
//
class SpatialDecomposition : public BaseDecomposition {
public:
  SpatialDecomposition(unsigned orderN, unsigned subbandsN, double a, double sampleRate, unsigned useSubbandsN = 0);
  ~SpatialDecomposition() { }

  virtual void estimate_Bkl(double theta, double phi, const gsl_vector_complex* snapshot, unsigned subbandX);
  virtual const gsl_matrix_complex* linearize(gsl_vector* xk, int frame_no);
  virtual const gsl_vector_complex* predicted_observation(gsl_vector* xk, int frame_no);
  virtual void calculate_gkl(double theta, double phi, unsigned subbandX);
};

typedef Inherit<SpatialDecomposition, BaseDecompositionPtr>	SpatialDecompositionPtr;

// ----- definition of class 'BaseSphericalArrayTracker' -----
//
class BaseSphericalArrayTracker : public VectorFloatFeatureStream {
  typedef list<VectorComplexFeatureStreamPtr>	ChannelList_;
  typedef ChannelList_::iterator		ChannelIterator_;
  typedef BaseDecomposition::SubbandEntry	SubbandEntry;
  typedef BaseDecomposition::SubbandListPtr	SubbandListPtr;
  typedef BaseDecomposition::Iterator		Iterator;

public:
  BaseSphericalArrayTracker(BaseDecompositionPtr& baseDecomposition, double sigma2_u = 10.0, double sigma2_v = 10.0, double sigma2_init = 10.0,
			    unsigned maxLocalN = 1, const String& nm = "BaseSphericalArrayTracker");
  ~BaseSphericalArrayTracker();

  virtual const gsl_vector_float* next(int frame_no = -5) = 0;

  virtual void reset();

  unsigned chanN() const { return channelList_.size(); }
  void set_channel(VectorComplexFeatureStreamPtr& chan);
  void set_V(const gsl_matrix_complex* Vk, unsigned subbandX);
  void next_speaker();
  void set_initial_position(double theta, double phi);

#ifdef ENABLE_LEGACY_BTK_API
  void setChannel(VectorComplexFeatureStreamPtr& chan){ set_channel(chan); }
  void setV(const gsl_matrix_complex* Vk, unsigned subbandX){ set_V(Vk, subbandX); }
  void nextSpeaker(){ next_speaker(); }
  void setInitialPosition(double theta, double phi){ set_initial_position(theta, phi); }
#endif

protected:
  static void printMatrix_(const gsl_matrix_complex* mat);
  static void printMatrix_(const gsl_matrix* mat);
  static void printVector_(const gsl_vector_complex* vec);
  static void printVector_(const gsl_vector* vec);

  static double calc_givens_rotation_(double v1, double v2, double& c, double& s);
  static void apply_givens_rotation_(double v1, double v2, double c, double s, double& v1p, double& v2p);

  void alloc_image_();
  void update_(const gsl_matrix_complex* Hbar_k, const gsl_vector_complex* yhat_k, const SubbandListPtr& subbandList);
  void lower_triangularize_();
  void copy_position_();
  void check_physical_constraints_();
  double calc_residual_();
  void realify_(const gsl_matrix_complex* Hbar_k, const gsl_vector_complex* yhat_k);
  void realify_residual_();

  static const unsigned 			StateN_;
  static const gsl_complex 			ComplexZero_;
  static const gsl_complex 			ComplexOne_;
  static const double      			Epsilon_;
  static const double				Tolerance_;

  bool						firstFrame_;
  const unsigned				subbandsN_;
  const unsigned				subbandsN2_;
  const unsigned				useSubbandsN_;
  const unsigned 				modesN_;
  const unsigned 				subbandLengthN_;
  const unsigned 				observationN_;
  const unsigned 				maxLocalN_;
  bool						is_end_;
  const double 					sigma_init_;

  SnapShotArrayPtr				snapshot_array_;
  BaseDecompositionPtr				base_decomposition_;
  ChannelList_					channelList_;

  // these quantities are stored as Cholesky factors
  gsl_matrix*					U_;
  gsl_matrix*					V_;
  gsl_matrix*					K_k_k1_;

  // work space for state estimate update
  gsl_matrix*					prearray_;
  gsl_vector_complex*				vk_;
  gsl_matrix*					Hbar_k_;
  gsl_vector*					yhat_k_;
  gsl_vector*					correction_;
  gsl_vector*					position_;
  gsl_vector*					eta_i_;
  gsl_vector*					delta_;
  gsl_vector_complex*				residual_;
  gsl_vector*					residual_real_;
  gsl_vector*					scratch_;
};

typedef Inherit<BaseSphericalArrayTracker, VectorFloatFeatureStreamPtr> BaseSphericalArrayTrackerPtr;

// ----- definition of class 'ModalSphericalArrayTracker' -----
//
class ModalSphericalArrayTracker : public BaseSphericalArrayTracker {
  typedef list<VectorComplexFeatureStreamPtr>	ChannelList_;
  typedef ChannelList_::iterator		ChannelIterator_;

public:
  ModalSphericalArrayTracker(ModalDecompositionPtr& modalDecomposition, double sigma2_u = 10.0, double sigma2_v = 10.0, double sigma2_init = 10.0,
			     unsigned maxLocalN = 1, const String& nm = "ModalSphericalArrayTracker");
  ~ModalSphericalArrayTracker() { }

  const gsl_vector_float* next(int frame_no = -5);
};

typedef Inherit<ModalSphericalArrayTracker, VectorFloatFeatureStreamPtr> ModalSphericalArrayTrackerPtr;

// ----- definition of class 'SpatialSphericalArrayTracker' -----
//
class SpatialSphericalArrayTracker : public BaseSphericalArrayTracker {
  typedef list<VectorComplexFeatureStreamPtr>	ChannelList_;
  typedef ChannelList_::iterator		ChannelIterator_;

public:
 SpatialSphericalArrayTracker(SpatialDecompositionPtr& spatialDecomposition, double sigma2_u = 10.0, double sigma2_v = 10.0, double sigma2_init = 10.0,
			       unsigned maxLocalN = 1, const String& nm = "SpatialSphericalArrayTracker");
  ~SpatialSphericalArrayTracker() { }

  const gsl_vector_float* next(int frame_no= -5);
};

typedef Inherit<SpatialSphericalArrayTracker, VectorFloatFeatureStreamPtr> SpatialSphericalArrayTrackerPtr;

// ----- definition of class 'PlaneWaveSimulator' -----
//
class PlaneWaveSimulator : public VectorComplexFeatureStream {
public:
  PlaneWaveSimulator(const VectorComplexFeatureStreamPtr& source, ModalDecompositionPtr& modalDecomposition,
		     unsigned channelX, double theta, double phi, const String& nm = "Plane Wave Simulator");
  ~PlaneWaveSimulator();

  const gsl_vector_complex* next(int frame_no = -5);

  virtual void reset();

private:
  static const gsl_complex			ComplexZero_;

  const unsigned				subbandsN_;
  const unsigned				subbandsN2_;
  const unsigned				channelX_;
  const double					theta_;
  const double					phi_;

  VectorComplexFeatureStreamPtr			source_;
  ModalDecompositionPtr				modalDecomposition_;
  gsl_vector_complex*				subbandCoefficients_;
};

typedef Inherit<PlaneWaveSimulator, VectorComplexFeatureStreamPtr> PlaneWaveSimulatorPtr;

#ifdef HAVE_GSL_V11X

int gsl_matrix_complex_add( gsl_matrix_complex *am, gsl_matrix_complex *bm )
{
  if( am->size1 != bm->size1 ){
    fprintf(stderr,"Dimension error \n",am->size1, bm->size1 );
    return 0;
  }

  if( am->size2 != bm->size2 ){
    fprintf(stderr,"Dimension error \n",am->size2, bm->size2 );
    return 0;
  }

  for(size_t i=0;i<am->size1;i++){
    for(size_t j=0;j<am->size2;j++){
      gsl_complex val = gsl_complex_add( gsl_matrix_complex_get(am,i,j), gsl_matrix_complex_get(bm,i,j) );
      gsl_matrix_complex_set(am,i,j,val);
    }
  }

  return 1;
}

int gsl_vector_complex_add( gsl_vector_complex *av, gsl_vector_complex *bv )
{
  if( av->size != bv->size ){
    fprintf(stderr,"Dimension error \n",av->size, bv->size);
    return 0;
  }

  for(size_t i=0;i<av->size;i++){
    gsl_complex val = gsl_complex_add( gsl_vector_complex_get(av,i), gsl_vector_complex_get(bv,i) );
    gsl_vector_complex_set(av,i,val);
  }

  return 1;
}

int gsl_vector_complex_sub( const gsl_vector_complex *av, const gsl_vector_complex *bv )
{
  if( av->size != bv->size ){
    fprintf(stderr,"Dimension error \n",av->size, bv->size);
    return 0;
  }

  for(size_t i=0;i<av->size;i++){
    gsl_complex val = gsl_complex_sub( gsl_vector_complex_get(av,i), gsl_vector_complex_get(bv,i) );
    gsl_vector_complex_set((gsl_vector_complex*)av,i,val);
  }

  return 1;
}

#endif /* End of HAVE_GSL_V11X */

#endif


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, D. Strubbe

 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, 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include <ctype.h>
#include <math.h>
#include <gsl/gsl_complex_math.h>

#include "symbols.h"
#include "gsl_userdef.h"

/* The symbol table: a chain of `struct symrec'.  */
symrec *sym_table = (symrec *)0;

void str_tolower(char *in)
{
  for(; *in; in++)
    *in = (char)tolower(*in);
}

void sym_mark_table_used ()
{
  symrec *ptr;

  for (ptr = sym_table; ptr != (symrec *) 0;
       ptr = (symrec *)ptr->next)
  {
    ptr->used = 1;
  }
}

symrec *putsym (const char *sym_name, symrec_type sym_type)
{
  symrec *ptr;
  ptr = (symrec *)malloc(sizeof(symrec));

  /* names are always lowercase */
  ptr->name = strdup(sym_name);
  str_tolower(ptr->name);

  ptr->def  = 0;
  ptr->used = 0;
  ptr->type = sym_type;
  GSL_SET_COMPLEX(&ptr->value.c, 0, 0); /* set value to 0 even if fctn.  */
  ptr->next = (struct symrec *)sym_table;
  sym_table = ptr;
  return ptr;
}

symrec *getsym (const char *sym_name)
{
  symrec *ptr;
  for (ptr = sym_table; ptr != (symrec *) 0;
       ptr = (symrec *)ptr->next)
    if (strcasecmp(ptr->name,sym_name) == 0){
      ptr->used = 1;
      return ptr;
    }
  return (symrec *) 0;
}

int rmsym (const char *sym_name)
{
  symrec *ptr, *prev;
  for (prev = (symrec *) 0, ptr = sym_table; ptr != (symrec *) 0;
       prev = ptr, ptr = ptr->next)
    if (strcasecmp(ptr->name,sym_name) == 0){
      if(prev == (symrec *) 0)
	sym_table = ptr->next;
      else
	prev->next = ptr->next;
      free(ptr->name);
      free(ptr);

      return 1;
    }

  return 0;
}

struct init_fntc{
  char *fname;
  int  nargs;
  gsl_complex (*fnctptr)();
};

void sym_notdef (symrec *sym)
{
  fprintf(stderr, "Parser error: symbol '%s' used before being defined.\n", sym->name);
  exit(1);
}

void sym_redef (symrec *sym)
{
  fprintf(stderr, "Parser warning: redefining symbol, previous value ");
  sym_print(stderr, sym);
  fprintf(stderr, "\n");
}

void sym_wrong_arg (symrec *sym)
{
  if(sym->type == S_BLOCK) {
    fprintf(stderr, "Parser error: block name '%s' used in variable context.\n", sym->name);
  } else if(sym->type == S_STR) {
    fprintf(stderr, "Parser error: string variable '%s' used in expression context.\n", sym->name);
  } else {
    fprintf(stderr, "Parser error: function '%s' requires %d argument(s).\n", sym->name, sym->nargs);
  }
  exit(1);
}

static struct init_fntc arith_fncts[] = {
  {"sqrt",   1, (gsl_complex (*)()) &gsl_complex_sqrt},
  {"exp",    1, (gsl_complex (*)()) &gsl_complex_exp},
  {"ln",     1, (gsl_complex (*)()) &gsl_complex_log},
  {"log",    1, (gsl_complex (*)()) &gsl_complex_log},
  {"log10",  1, (gsl_complex (*)()) &gsl_complex_log10},
  {"logb",   2, (gsl_complex (*)()) &gsl_complex_log_b}, /* takes two arguments logb(z, b) = log_b(z) */

  {"arg",    1, (gsl_complex (*)()) &gsl_complex_carg},
  {"abs",    1, (gsl_complex (*)()) &gsl_complex_cabs},
  {"abs2",   1, (gsl_complex (*)()) &gsl_complex_cabs2},
  {"logabs", 1, (gsl_complex (*)()) &gsl_complex_clogabs},

  {"conjg",  1, (gsl_complex (*)()) &gsl_complex_conjugate},
  {"inv",    1, (gsl_complex (*)()) &gsl_complex_inverse},

  {"sin",    1, (gsl_complex (*)()) &gsl_complex_sin},
  {"cos",    1, (gsl_complex (*)()) &gsl_complex_cos},
  {"tan",    1, (gsl_complex (*)()) &gsl_complex_tan},
  {"sec",    1, (gsl_complex (*)()) &gsl_complex_sec},
  {"csc",    1, (gsl_complex (*)()) &gsl_complex_csc},
  {"cot",    1, (gsl_complex (*)()) &gsl_complex_cot},

  {"asin",   1, (gsl_complex (*)()) &gsl_complex_arcsin},
  {"acos",   1, (gsl_complex (*)()) &gsl_complex_arccos},
  {"atan",   1, (gsl_complex (*)()) &gsl_complex_arctan},
  {"atan2",  2, (gsl_complex (*)()) &gsl_complex_arctan2}, /* takes two arguments atan2(y,x) = atan(y/x) */
  {"asec",   1, (gsl_complex (*)()) &gsl_complex_arcsec},
  {"acsc",   1, (gsl_complex (*)()) &gsl_complex_arccsc},
  {"acot",   1, (gsl_complex (*)()) &gsl_complex_arccot},

  {"sinh",   1, (gsl_complex (*)()) &gsl_complex_sinh},
  {"cosh",   1, (gsl_complex (*)()) &gsl_complex_cosh},
  {"tanh",   1, (gsl_complex (*)()) &gsl_complex_tanh},
  {"sech",   1, (gsl_complex (*)()) &gsl_complex_sech},
  {"csch",   1, (gsl_complex (*)()) &gsl_complex_csch},
  {"coth",   1, (gsl_complex (*)()) &gsl_complex_coth},

  {"asinh",  1, (gsl_complex (*)()) &gsl_complex_arcsinh},
  {"acosh",  1, (gsl_complex (*)()) &gsl_complex_arccosh},
  {"atanh",  1, (gsl_complex (*)()) &gsl_complex_arctanh},
  {"asech",  1, (gsl_complex (*)()) &gsl_complex_arcsech},
  {"acsch",  1, (gsl_complex (*)()) &gsl_complex_arccsch},
  {"acoth",  1, (gsl_complex (*)()) &gsl_complex_arccoth},

/* user-defined step function. this is not available in GSL,
   but we use GSL namespacing and macros here. */
  {"step",   1, (gsl_complex (*)()) &gsl_complex_step_real},

/* Minimum and maximum of two arguments (comparing real parts) */
  {"min",    2, (gsl_complex (*)()) &gsl_complex_min_real},
  {"max",    2, (gsl_complex (*)()) &gsl_complex_max_real},

  {"erf",    1, (gsl_complex (*)()) &gsl_complex_erf},

  {"realpart", 1, (gsl_complex (*)()) &gsl_complex_realpart},
  {"imagpart", 1, (gsl_complex (*)()) &gsl_complex_imagpart},
  {"round",   1, (gsl_complex (*)()) &gsl_complex_round},
  {"floor",   1, (gsl_complex (*)()) &gsl_complex_floor},
  {"ceiling", 1, (gsl_complex (*)()) &gsl_complex_ceiling},

  {"rand",    0, (gsl_complex (*)()) &gsl_complex_rand},

  {0, 0, 0}
};

struct init_cnst{
	char *fname;
	double re;
	double im;
};

static struct init_cnst arith_cnts[] = {
	{"pi",    M_PI, 0},
	{"e",      M_E, 0},
	{"i",        0, 1},
	{"true",     1, 0},
	{"yes",      1, 0},
	{"false",    0, 0},
	{"no",       0, 0},
	{0,          0, 0}
};

char *reserved_symbols[] = {
  "x", "y", "z", "r", "w", "t", 0
};

void sym_init_table ()  /* puts arithmetic functions in table. */
{
  int i;
  symrec *ptr;
  for (i = 0; arith_fncts[i].fname != 0; i++){
    ptr = putsym (arith_fncts[i].fname, S_FNCT);
    ptr->def = 1;
    ptr->used = 1;
    ptr->nargs = arith_fncts[i].nargs;
    ptr->value.fnctptr = arith_fncts[i].fnctptr;
  }

  /* now the constants */
  for (i = 0; arith_cnts[i].fname != 0; i++){
    ptr = putsym(arith_cnts[i].fname, S_CMPLX);
    ptr->def = 1;
    ptr->used = 1;
    GSL_SET_COMPLEX(&ptr->value.c, arith_cnts[i].re, arith_cnts[i].im);
  }
}

void sym_end_table()
{
  symrec *ptr, *ptr2;

  for (ptr = sym_table; ptr != NULL;){
    free(ptr->name);
    switch(ptr->type){
    case S_STR:
      free(ptr->value.str);
      break;
    case S_BLOCK:
      if(ptr->value.block->n > 0){
	free(ptr->value.block->lines);
      }
      free(ptr->value.block);
      break;
    case S_CMPLX:
    case S_FNCT:
      break;
    }
    ptr2 = ptr->next;
    free(ptr);
    ptr = ptr2;
  }

  sym_table = NULL;
}

/* this function is defined in src/basic/varinfo_low.c */
int varinfo_variable_exists(const char * var_name);

void sym_output_table(int only_unused, int mpiv_node)
{
  FILE *f;
  symrec *ptr;
  int any_unused = 0;

  if(mpiv_node != 0) {
    return;
  }

  if(only_unused) {
    f = stderr;
  } else {
    f = stdout;
  }

  for(ptr = sym_table; ptr != NULL; ptr = ptr->next){
    if(only_unused && ptr->used == 1) continue;
    if(only_unused && varinfo_variable_exists(ptr->name)) continue;
    if(any_unused == 0) {
      fprintf(f, "\nParser warning: possible mistakes in input file.\n");
      fprintf(f, "List of variable assignments not used by parser:\n");
      any_unused = 1;
    }

    sym_print(f, ptr);
  }
  if(any_unused == 1) {
    fprintf(f, "\n");
  }
}

void sym_print(FILE *f, const symrec *ptr)
{
  fprintf(f, "%s", ptr->name);
  switch(ptr->type){
  case S_CMPLX:
    if(fabs(GSL_IMAG(ptr->value.c)) < 1.0e-14){
      fprintf(f, " = %f\n", GSL_REAL(ptr->value.c));
    } else {
      fprintf(f, " = (%f,%f)\n", GSL_REAL(ptr->value.c), GSL_IMAG(ptr->value.c));
    }
    break;
  case S_STR:
    fprintf(f, " = \"%s\"\n", ptr->value.str);
    break;
  case S_BLOCK:
    fprintf(f, "%s\n", " <= BLOCK");
    break;
  case S_FNCT:
    fprintf(f, "%s\n", " <= FUNCTION");
    break;
  }
}

ext: c

max_stars_repo_name: shunsuke-sato/octopus

max_issues_repo_name: shunsuke-sato/octopus

max_forks_repo_name: shunsuke-sato/octopus


text:
/*
* Copyright 2019 Virtru Corporation
*
* SPDX - License Identifier: BSD-3-Clause-Clear
*
*/
//
//  TDF SDK
//
//  Created by Sujan Reddy on 2019/04/07.
//

#ifndef VIRTRU_BYTES_H
#define VIRTRU_BYTES_H

#include <gsl/span>
#include <array>

/// This provides a simple wrapper around gsl to work with bytes.
namespace virtru::crypto {

    /// Constants
    constexpr auto kKeyLength = 32;
    constexpr auto kIVLength = 16;

   ///
   /// Define the bytes(Read only).
   ///
   template <std::ptrdiff_t Extent = gsl :: dynamic_extent>
   using BytesT = gsl::span <const gsl::byte, Extent>;
   using Bytes = BytesT<>;

   ///
   /// Wrapper around std::array
   ///
   template <std::size_t N>
   using ByteArray = std::array <gsl :: byte, N>;
   using WrappedKey = ByteArray<kKeyLength>;
   using IV = ByteArray<kIVLength>;

    ///
    /// Define writeable bytes
    ///
    template <std::ptrdiff_t Extent = gsl::dynamic_extent>
    using WriteableBytesT = gsl::span <gsl::byte, Extent>;
    using WriteableBytes = WriteableBytesT<>;

    template <std :: ptrdiff_t Extent >
    constexpr auto toBytes(BytesT < Extent > data) noexcept {
        return data;
    }

    template <std::ptrdiff_t Extent>
    constexpr auto toWriteableBytes(WriteableBytesT <Extent> data) noexcept {
        return data;
    }

    /// Adding constness
    template <std::ptrdiff_t Extent>
    constexpr auto toBytes (gsl::span <gsl::byte, Extent> data) noexcept {
        return BytesT <Extent> {data};
    }

    /// Any containers inherited from basic_string<ElementType> to bytes(Read only).
    template <typename ElementType,
              std::ptrdiff_t Extent,
              typename = std::enable_if_t <std::has_unique_object_representations_v <ElementType>> >
    auto toBytes (gsl::span <ElementType, Extent> data) noexcept {
        return as_bytes(data);
    }

    /// Any containers inherited from basic_string<ElementType> to writable bytes.
    template <typename ElementType,
              std::ptrdiff_t Extent,
              typename = std::enable_if_t <std::has_unique_object_representations_v <ElementType>> >
    auto toWriteableBytes(gsl::span <ElementType, Extent> data) noexcept {
        return as_writeable_bytes(data);
    }

    /// Form std::array to bytes(Read only).
    template <typename ElementType,
              std::size_t N,
              typename = std::enable_if_t <std::has_unique_object_representations_v <ElementType>> >
    constexpr auto toBytes(const std::array <ElementType, N> & data) noexcept {
        return toBytes(gsl::span <const ElementType, N > { data });
    }

    /// Form std::array to writable bytes.
    template <typename ElementType,
              std::size_t N,
              typename = std::enable_if_t <std::has_unique_object_representations_v <ElementType>> >
    constexpr auto toWriteableBytes(std::array <ElementType, N> & data) noexcept {
        return toWriteableBytes(gsl::span <ElementType, N> { data });
    }

    // Form Plain old C array to bytes(Read only).
    template <typename ElementType,
              std::size_t N,
              typename = std::enable_if_t <std::has_unique_object_representations_v <ElementType>> >
    constexpr auto toBytes(ElementType (& arr)[N]) noexcept {
        return toBytes(gsl::make_span(arr));
    }

    // Form Plain old C array to writable bytes.
    template <typename ElementType,
              std::size_t N,
              typename = std::enable_if_t <std::has_unique_object_representations_v <ElementType>> >
    constexpr auto toWriteableBytes(ElementType(&arr)[N]) noexcept {
        return toWriteableBytes(gsl::make_span(arr));
    }

    // From container(std::string/std::vector<char>) to bytes(Read only).
    template <class Cont,
              typename = std::enable_if_t <std::has_unique_object_representations_v <std::remove_reference_t <
                                           decltype(* std::declval <const Cont&> ().data ())>> > >
    constexpr auto toBytes(const Cont & cont) noexcept {
        return toBytes(gsl::make_span(cont));
    }

    // From container(std::string/std::vector<char>) to to writable bytes.
    template <class Cont,
              typename = std::enable_if_t <std::has_unique_object_representations_v <std::remove_reference_t <
                                           decltype(* std::declval <Cont &> ().data ())>> > >
    constexpr auto toWriteableBytes(Cont & cont) noexcept {
        return toWriteableBytes(gsl::make_span(cont));
    }

    ///
    /// Simple conversions.
    ///
    template <typename T>
    using ConvertibleToBytes = decltype(toBytes(std::declval <const T&> ()));

    template <typename T, std::ptrdiff_t Extent = gsl::dynamic_extent>
    inline constexpr bool isBytes = std::is_convertible_v <const T&, BytesT <Extent>>;

    template <typename T>
    using ExplicitlyConvertibleToBytes = std::enable_if_t <! isBytes <T>, ConvertibleToBytes <T>>;

    template <typename T, typename = ConvertibleToBytes <T>>
    auto toDynamicBytes (const T& t) {
        return Bytes { toBytes(t) };
    }

    template <typename T, typename = ExplicitlyConvertibleToBytes <T>>
    auto toDynamicBytesExplicitly (const T& t) {
        return toDynamicBytes(t);
    }

    inline unsigned char* toUchar(gsl::byte* p) {
        return reinterpret_cast <unsigned char *> (p);
    }

    inline const unsigned char *toUchar (const gsl::byte *p) {
        return reinterpret_cast <const unsigned char *> (p);
    }

    template <std::ptrdiff_t Extent>
    auto toUchar(BytesT<Extent> src) {
        return gsl::span <const unsigned char, Extent> {reinterpret_cast <const unsigned char *> (src.data ()), src.size()};
    }

    inline char* toChar(gsl::byte *p) {
        return reinterpret_cast <char*> (p);
    }

    inline const char* toChar(const gsl::byte* p) {
        return reinterpret_cast <const char *> (p);
    }

    template <std::ptrdiff_t Extent>
    auto toChar(BytesT <Extent> src) {
        return gsl :: span <const char, Extent> { reinterpret_cast <const char*> (src.data()), src.size ()};
    }

    template <std::ptrdiff_t Extent>
    auto toChar(WriteableBytesT <Extent> src) {
        return gsl::span <const char, Extent> { reinterpret_cast <const char *> (src.data()), src.size ()};
    }

    inline auto finalizeSize(WriteableBytes& buffer, const int& size) {
        return gsl::finally( [&buffer, &size ] { buffer = buffer.first(size); });
    }
} // namespace virtru::crypto

#endif //VIRTRU_BYTES_H

max_stars_repo_name: opentdf/client-cpp

max_issues_repo_name: opentdf/client-cpp

max_forks_repo_name: opentdf/client-cpp


text:
// ********************  Functions for segSiteHMM

#include "stdio.h"
#include "math.h"
#include "segSites.h"
#include "../hmm/popGenTools.h"
#include "../hmm/numerical.h"
#include "assert.h"
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>

//core function for probs
double my_h(double q, void * p){
  struct my_f_params * params = (struct my_f_params *) p;
  double beta = params->beta;
  int i = params->i;
  int n = params->n;

   return ((1.0 - exp(-2.0 * beta * (1.0 - q))) / (1.0 - exp(-2.0 * beta))) \
     * (2.0 / (q * (1.0 - q))) * gsl_ran_binomial_pdf(i, q, n);
}

/* integral of f over allele freq -- returns Prob(Segsite | mu, beta) */
double my_H(gsl_vector *pVector, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_f_params fParams;
  gsl_function f;
  gsl_integration_workspace * w;
  double result,error;

  if (gsl_vector_get(pVector, 0) == 0){
    return gsl_vector_get(pVector, 1) / params->i;
  }
  else{
    w = gsl_integration_workspace_alloc (100000);
    fParams.beta = gsl_vector_get(pVector, 0);
    fParams.i = params->i;
    fParams.n = params->n;
    f.function = &(my_h);
    f.params = &fParams;
    gsl_integration_qags(&f,0.0,1.0,1e-5,1e-5, 100000,w, &result,&error);
    gsl_integration_workspace_free(w);
    return result * gsl_vector_get(pVector,1);
  }
}

/* siteProbMatrix- fills up a matrix of log siteProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (0 to max -1), freq 0 is for monomorphic sites */
void siteProbMatrix(gsl_matrix *probs, gsl_vector *pVector,  int maxSampleSize, gsl_vector *sampleSizeVector){
  struct my_snpProb_params FParams;
  int i, j;
  double tot, prob;

  //check & init probs
  assert(probs->size1 == maxSampleSize + 1 && probs->size2 == maxSampleSize);
  gsl_matrix_set_zero(probs);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //set tot to zero for p{monomorphic}
      tot = 0;
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//calc prob
	prob = 0;
	FParams.i = j;
	FParams.n = i;
	//need to impose constraint that p{SegSite = 0 | theta < 0}
	if(gsl_vector_get(pVector,1) < 0){
	  prob = 0;
	}
	else{
	  prob = my_H(pVector, &FParams);
	}
	gsl_matrix_set(probs, i, j, log(prob));
	tot += prob;
      }
      //set p{monomorphic}
      gsl_matrix_set(probs,i, 0, log(1.0 - tot));
    }
  }
}

double weightedLikLookSites(gsl_vector *pVector, void * p){
  struct my_jointLik_params * params = (struct my_jointLik_params *) p;
  gsl_matrix *probs;
  double lik;
  int i;

  //make prob matrix (log space)
  probs = gsl_matrix_alloc(params->maxSampleSize + 1, params->maxSampleSize);
  siteProbMatrix(probs, pVector,  params->maxSampleSize, params->sampleSizeVector);

  //tally likelihood
  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)  * gsl_vector_get(params->weights,i);
  }
  gsl_matrix_free(probs);
  return -lik;
}
//  weightedLikLookSites2State -- for use in 2 state HMM
double weightedLikLookSites2State(gsl_vector *pVector, void * p){
  struct my_jointLik_params * params = (struct my_jointLik_params *) p;
  gsl_matrix *probs, *probsN;
  gsl_vector *npVector;
  double lik;
  int i;

  //setup neutral params
  npVector = gsl_vector_alloc(2);
  gsl_vector_set(npVector,0,0);
  gsl_vector_set(npVector,1,gsl_vector_get(pVector,1));

  //alloc & make probs and probsN matrix (log space)
  probs = gsl_matrix_alloc(params->maxSampleSize + 1, params->maxSampleSize);
  siteProbMatrix(probs, pVector,  params->maxSampleSize, params->sampleSizeVector);
  probsN = gsl_matrix_alloc(params->maxSampleSize + 1, params->maxSampleSize);
  siteProbMatrix(probsN,npVector,  params->maxSampleSize, params->sampleSizeVector);

  //tally likelihood
  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)  * gsl_vector_get(params->weights,i);
    lik += gsl_matrix_get(probsN,params->data[i].n,params->data[i].i)  * (1.0 - gsl_vector_get(params->weights,i));
  }
  gsl_matrix_free(probs);
  gsl_matrix_free(probsN);
  gsl_vector_free(npVector);
  return -lik;
}

gsl_vector *jointMLEst(double *lik, void * p){
  size_t np = 2;
  struct my_jointLik_params * params = (struct my_jointLik_params *) p;
  const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
  gsl_multimin_fminimizer *s = NULL;
  gsl_multimin_function minex_func;
  gsl_vector *ss, *x, *results;

  size_t iter = 0, i;
  int status;
  double size;

  /* Initial vertex size vector */
  ss = gsl_vector_alloc (np);
  //initialize results vector
  results = gsl_vector_alloc (np);
  /* Set all step sizes to 1 */
  gsl_vector_set_all (ss, 0.01);

  /* Starting point */
  x = gsl_vector_alloc (np);
  gsl_vector_set (x, 0, gsl_vector_get(params->startingPoint,0));
  gsl_vector_set (x, 1, gsl_vector_get(params->startingPoint,1));

  /* Initialize method and iterate */
  minex_func.f = &weightedLikLookSites2State;
  minex_func.n = np;
  minex_func.params = params;

  s = gsl_multimin_fminimizer_alloc (T, np);
  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);
      if (status)
	break;
      size = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, 1e-3);
    }
  while (status == GSL_CONTINUE && iter < 500);
  for (i = 0; i < np; i++){
    gsl_vector_set(results,i,gsl_vector_get (s->x, i));
  }
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);
  return(results);
}

//  weightedLikLookSitesNState -- for use in n state HMM
// pVector now consists of n parameters, not just two
// pVector = [theta_hat, beta_hat_1, beta_hat_2,
double weightedLikLookSitesNState(gsl_vector *pVector, void * p){
  struct my_jointLik_params * params = (struct my_jointLik_params *) p;
  gsl_matrix *probs;
  gsl_vector *currentVector;
  double lik;
  int i,j;

  //setup neutral params
  currentVector = gsl_vector_alloc(params->nstates);
  gsl_vector_set(currentVector,0,0);
  gsl_vector_set(currentVector,1,gsl_vector_get(pVector,0));

  //make prob matrix (log space)
  probs = gsl_matrix_alloc(params->maxSampleSize + 1, params->maxSampleSize);
  siteProbMatrix(probs, currentVector,  params->maxSampleSize, params->sampleSizeVector);

  //tally likelihood due to neutral
  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)  * gsl_matrix_get(params->posts,0,i);
  }
  //go through states tallying lik
  for (i = 1; i < params->nstates; i++){
    gsl_vector_set(currentVector,0,gsl_vector_get(pVector,i));
    gsl_vector_set(currentVector,1,gsl_vector_get(pVector,0));
    siteProbMatrix(probs,currentVector,  params->maxSampleSize, params->sampleSizeVector);
    for(j = 0; j < params->snpNumber; j++){
      lik += gsl_matrix_get(probs,params->data[j].n,params->data[j].i)  * gsl_matrix_get(params->posts,i,j);
    }
  }
  gsl_matrix_free(probs);
  gsl_vector_free(currentVector);
  return -lik;
}

//N state version
void jointMLEstNState(gsl_vector *results, double *lik, void * p){
  struct my_jointLik_params * params = (struct my_jointLik_params *) p;
  size_t np = params->nstates;
  const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
  gsl_multimin_fminimizer *s = NULL;
  gsl_multimin_function minex_func;
  gsl_vector *ss, *x;

  size_t iter = 0, i;
  int status,j;
  double size;

  //check results vector
  assert(results->size == np);

  /* Initial vertex size vector */
  ss = gsl_vector_alloc (np);

  /* Set all step sizes to 1 */
  gsl_vector_set_all (ss, 0.01);

  /* Starting point */
  x = gsl_vector_alloc (np);
  for(j = 0; j < np; j++){
    gsl_vector_set (x, j, gsl_vector_get(params->startingPoint,j));
  }

  /* Initialize method and iterate */
  minex_func.f = &weightedLikLookSitesNState;
  minex_func.n = np;
  minex_func.params = params;

  s = gsl_multimin_fminimizer_alloc (T, np);
  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);
      if (status)
	break;
      size = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, 1e-5);
    }
  while (status == GSL_CONTINUE && iter < 500);
  for (i = 0; i < np; i++){
    gsl_vector_set(results,i,gsl_vector_get (s->x, i));
  }
  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);
}

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
#include "jfftw_complex_Plan.h"
#include <fftw.h>

/*
 * Class:     jfftw_complex_Plan
 * Method:    createPlan
 * Signature: (III)V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_Plan_createPlan( JNIEnv *env, jobject obj, jint n, jint dir, jint flags )
{
	jclass clazz;
	jfieldID id;
	jbyteArray arr;
	unsigned char* carr;

	if( sizeof( jdouble ) != sizeof( fftw_real ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/RuntimeException" ), "jdouble and fftw_real are incompatible" );
		return;
	}

	clazz = (*env)->GetObjectClass( env, obj );
	id    = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	arr   = (*env)->NewByteArray( env, sizeof( fftw_plan ) );
	carr  = (*env)->GetByteArrayElements( env, arr, 0 );

	(*env)->MonitorEnter( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	*(fftw_plan*)carr = fftw_create_plan( n, dir, flags );

	(*env)->MonitorExit( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->SetObjectField( env, obj, id, arr );
}
/*
 * Class:     jfftw_complex_Plan
 * Method:    createPlanSpecific
 * Signature: (III[DI[DI)V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_Plan_createPlanSpecific( JNIEnv *env, jobject obj, jint n, jint dir, jint flags, jdoubleArray in, jint idist, jdoubleArray out, jint odist )
{
	jclass clazz;
	jfieldID id;
	jbyteArray arr;
	unsigned char* carr;
	double *cin, *cout;

	if( sizeof( jdouble ) != sizeof( fftw_real ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/RuntimeException" ), "jdouble and fftw_real are incompatible" );
		return;
	}

	clazz = (*env)->GetObjectClass( env, obj );
	id    = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	arr   = (*env)->NewByteArray( env, sizeof( fftw_plan ) );
	carr  = (*env)->GetByteArrayElements( env, arr, 0 );
	cin   = (*env)->GetDoubleArrayElements( env, in, 0 );
	cout  = (*env)->GetDoubleArrayElements( env, out, 0 );

	(*env)->MonitorEnter( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	*(fftw_plan*)carr = fftw_create_plan_specific( n, dir, flags, (fftw_complex*)cin, idist, (fftw_complex*)cout, odist );

	(*env)->MonitorExit( env, (*env)->FindClass( env, "jfftw/Plan" ) );

	(*env)->ReleaseDoubleArrayElements( env, in, cin, 0 );
	(*env)->ReleaseDoubleArrayElements( env, out, cout, 0 );
	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->SetObjectField( env, obj, id, arr );
}
/*
 * Class:     jfftw_complex_Plan
 * Method:    destroyPlan
 * Signature: ()V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_Plan_destroyPlan( JNIEnv* env, jobject obj )
{
	jclass clazz = (*env)->GetObjectClass( env, obj );
	jfieldID id = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	jbyteArray arr = (jbyteArray)(*env)->GetObjectField( env, obj, id );
	unsigned char* carr = (*env)->GetByteArrayElements( env, arr, 0 );

	fftw_destroy_plan( *(fftw_plan*)carr );

	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->SetObjectField( env, obj, id, NULL );
}
/*
 * Class:     jfftw_complex_Plan
 * Method:    transform
 * Signature: ([D)[D
 */
JNIEXPORT jdoubleArray JNICALL Java_jfftw_complex_Plan_transform___3D( JNIEnv* env, jobject obj, jdoubleArray in )
{
	jdouble *cin, *cout;
	jdoubleArray out;
	int i;

	jclass clazz = (*env)->GetObjectClass( env, obj );
	jfieldID id = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	jbyteArray arr = (jbyteArray)(*env)->GetObjectField( env, obj, id );
	unsigned char* carr = (*env)->GetByteArrayElements( env, arr, 0 );
	fftw_plan plan = *(fftw_plan*)carr;
	if( plan->n * 2 != (*env)->GetArrayLength( env, in ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IndexOutOfBoundsException" ), "the Plan was created for a different length" );
		(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
		return NULL;
	}

	cin = (*env)->GetDoubleArrayElements( env, in, 0 );

	if( ! plan->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorEnter( env, obj );
	}

	if( plan->flags & FFTW_IN_PLACE )
	{
		out = in;

		fftw_one( plan, (fftw_complex*)cin, NULL );
	}
	else
	{
		out = (*env)->NewDoubleArray( env, plan->n * 2 );
		cout = (*env)->GetDoubleArrayElements( env, out, 0 );

		fftw_one( plan, (fftw_complex*)cin, (fftw_complex*)cout );

		(*env)->ReleaseDoubleArrayElements( env, out, cout, 0 );
	}

	if( ! plan->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorEnter( env, obj );
	}

	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->ReleaseDoubleArrayElements( env, in, cin, 0 );
	return out;
}
/*
 * Class:     jfftw_complex_Plan
 * Method:    transform
 * Signature: (I[DII[DII)V
 */
JNIEXPORT void JNICALL Java_jfftw_complex_Plan_transform__I_3DII_3DII( JNIEnv *env, jobject obj, jint howmany, jdoubleArray in, jint istride, jint idist, jdoubleArray out, jint ostride, jint odist )
{
	jdouble *cin, *cout;
	int i;

	jclass clazz = (*env)->GetObjectClass( env, obj );
	jfieldID id = (*env)->GetFieldID( env, clazz, "plan", "[B" );
	jbyteArray arr = (jbyteArray)(*env)->GetObjectField( env, obj, id );
	unsigned char* carr = (*env)->GetByteArrayElements( env, arr, 0 );
	fftw_plan plan = *(fftw_plan*)carr;
	if( (howmany - 1) * idist * 2 + plan->n * istride * 2 != (*env)->GetArrayLength( env, in ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IndexOutOfBoundsException" ), "the Plan was created for a different length (in)" );
		(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
		return;
	}
	if( (howmany - 1) * odist * 2 + plan->n * ostride * 2 != (*env)->GetArrayLength( env, out ) )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IndexOutOfBoundsException" ), "the Plan was created for a different length (out)" );
		(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
		return;
	}

	cin = (*env)->GetDoubleArrayElements( env, in, 0 );
	cout = (*env)->GetDoubleArrayElements( env, out, 0 );

	if( ! plan->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorEnter( env, obj );
	}
	fftw( plan, howmany, (fftw_complex*)cin, istride, idist, (fftw_complex*)cout, ostride, odist );
	if( ! plan->flags & FFTW_THREADSAFE )
	{
		// synchronization
		(*env)->MonitorExit( env, obj );
	}

	(*env)->ReleaseByteArrayElements( env, arr, carr, 0 );
	(*env)->ReleaseDoubleArrayElements( env, in, cin, 0 );
	(*env)->ReleaseDoubleArrayElements( env, out, cout, 0 );
}

max_stars_repo_name: CSCSI/Triana

max_issues_repo_name: CSCSI/Triana

max_forks_repo_name: CSCSI/Triana


text:
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include "../globals.h"
#include "galaxy_cluster.h"

#define NTABLE 256
#define GSL_SPLINE gsl_interp_cspline

#define Zero 1
#define Infinity 1e25

struct HaloProperties Halo = { 0 };

void setup_GalaxyCluster()
{
    Problem.Boxsize[0] = Problem.Boxsize[1] = Problem.Boxsize[2] = 1000;  // kpc
    Halo.Rho0 = 1e-26;  // g/cm^3
    Halo.Beta = 2.0 / 3.0;
    Halo.Rcore = 20;    // kpc
    Halo.R_Sample = Problem.Boxsize[0];

    sprintf ( Problem.Name, "IC_GalaxyCluster" );

    Problem.Rho_Max = Halo.Rho0;

    Density_Func_Ptr = &GalaxyCluster_Density;
    U_Func_Ptr = &GalaxyCluster_U;
    Velocity_Func_Ptr = &GalaxyCluster_Velocity;
}

float betamodel ( const float r )
{
    float rho0 = Halo.Rho0;
    float beta = Halo.Beta;
    float rc = Halo.Rcore;

    return rho0 * pow ( 1 + p2 ( r / rc ), -1.5 * beta );
}

float GalaxyCluster_Density ( const int ipart , const double bias )
{
    const float x = P[ipart].Pos[0] - 0.5 * Problem.Boxsize[0];
    const float y = P[ipart].Pos[1] - 0.5 * Problem.Boxsize[1];;
    const float z = P[ipart].Pos[2] - 0.5 * Problem.Boxsize[2];;

    const float r = sqrt ( x * x + y * y + z * z );

    return betamodel ( r );
}

void GalaxyCluster_Velocity ( const int ipart, float out[3] )
{
    out[0] = 0.0;
    out[1] = 0.0;
    out[2] = 0.0;
}

float GalaxyCluster_U ( const int ipart )
{
    return 0;
}

max_stars_repo_name: elehcim/WVTICs

max_issues_repo_name: elehcim/WVTICs

max_forks_repo_name: elehcim/WVTICs


text:
#pragma once

#include <cstddef> // std::byte
#include <map>
#include <memory> // std::align
#include <memory_resource>
#include <vector>

#include <gsl-lite/gsl-lite.hpp>

#include <cuda/runtime_api.hpp>

namespace thrustshift {

namespace pmr {

//! Unified CUDA Memory Resource. Buffers cannet be deallocated partially. Alignment is ignored.
class managed_resource_type : public std::pmr::memory_resource {
	void* do_allocate(std::size_t bytes,
	                  [[maybe_unused]] std::size_t alignment) override {
		auto region = cuda::memory::managed::detail::allocate(bytes);
		return region.get();
	}

	void do_deallocate(void* p,
	                   [[maybe_unused]] std::size_t bytes,
	                   [[maybe_unused]] std::size_t alignment) override {
		cuda::memory::managed::detail::free(p);
	}

	bool do_is_equal(
	    const std::pmr::memory_resource& other) const noexcept override {
		return this == &other;
	}
};

//! Device CUDA Memory Resource. Buffers cannot be deallocated partially. Alignment is ignored.
class device_resource_type : public std::pmr::memory_resource {
	void* do_allocate(std::size_t bytes,
	                  [[maybe_unused]] std::size_t alignment) override {
		auto region = cuda::memory::device::detail::allocate(bytes);
		return region.get();
	}

	void do_deallocate(void* p,
	                   [[maybe_unused]] std::size_t bytes,
	                   [[maybe_unused]] std::size_t alignment) override {
		cuda::memory::device::free(p);
	}

	bool do_is_equal(
	    const std::pmr::memory_resource& other) const noexcept override {
		return this == &other;
	}
};

class host_resource_type : public std::pmr::memory_resource {
	void* do_allocate(std::size_t bytes,
	                  [[maybe_unused]] std::size_t alignment) override {
		return malloc(bytes);
	}

	void do_deallocate(void* p,
	                   [[maybe_unused]] std::size_t bytes,
	                   [[maybe_unused]] std::size_t alignment) override {
		free(p);
	}

	bool do_is_equal(
	    const std::pmr::memory_resource& other) const noexcept override {
		return this == &other;
	}
};

namespace detail {

struct page_id_type {
	size_t bytes;
	size_t alignment;
};

inline bool operator<(const page_id_type& a, const page_id_type& b) {
	if (a.bytes == b.bytes) {
		return a.alignment < b.alignment;
	}
	return a.bytes < b.bytes;
}

} // namespace detail

/*! \brief Allocates only if requested buffer size is currently not in the free pool.
 *
 *  This pool is useful for functions, which require temporary GPU memory. The
 *  host can allocate memory via this pool, launch the kernel with the corresponding
 *  pointer and exit the function without deallocating the memory because the
 *  host is not aware about the runtime of the GPU kernel and when the kernel
 *  may read the temporary required memory. This pool should not be used if the
 *  byte size of the allocations differ often becaues every time an allocation
 *  is called it actually results into a new allocation if not the exact same
 *  size was allocated previously.
 *
 *  This class is useful if functions delegate the memory resource to other functions, which
 *  are maybe called iteratively and if the function requires different buffers of
 *  the same size.
 *
 *  \note std::pmr compatible
 */
template <class Upstream>
class delayed_pool_type : public std::pmr::memory_resource {
   private:
	struct page_item_type {
		void* ptr;
		bool allocated;
	};

   public:
	delayed_pool_type() = default;

	~delayed_pool_type() noexcept {
		for (auto& [page_id, book_page] : book_) {
			for (auto& page_item : book_page) {
				res_.deallocate(
				    page_item.ptr, page_id.bytes, page_id.alignment);
			}
		}
	}

	auto& get_book() const {
		return book_;
	}

   private:
	void* do_allocate(size_t bytes, size_t alignment) override {
		if (bytes == 0) {
			return nullptr;
		}
		const detail::page_id_type page_id({bytes, alignment});
		if (auto it = book_.find(page_id); it != book_.end()) {
			for (auto& page_item : it->second) {
				if (!page_item.allocated) {
					page_item.allocated = true;
					return page_item.ptr;
				}
			}
			// no page item was free
			void* ptr = res_.allocate(bytes, alignment);
			it->second.push_back({ptr, true});
			return ptr;
		}
		// the required size was never allocated before
		void* ptr = res_.allocate(bytes, alignment);
		book_[page_id] = {{ptr, true}};
		return ptr;
	}

	void do_deallocate(void* ptr,
	                   size_t bytes,
	                   size_t alignment) noexcept override {
		const detail::page_id_type page_id({bytes, alignment});
		if (ptr == nullptr || bytes == 0) {
			return;
		}
		for (auto& page_item : book_[page_id]) {
			if (page_item.ptr == ptr && page_item.allocated) {
				page_item.allocated = false;
				return;
			}
		}
		// The pointer which should be deallocated was never allocated before
		std::terminate();
	}

	bool do_is_equal(
	    const std::pmr::memory_resource& other) const noexcept override {
		return this == &other;
	}

	Upstream res_;
	// size in bytes -> (ptr, alignment)
	std::map<detail::page_id_type, std::vector<page_item_type>> book_;
};

/*! \brief Allocates only if there is no buffer in the free pool which is of the required size or larger.
 *
 *  This pool searches for the smallest buffer in the pool which has the required alignment and is
 *  of equal or larger byte size than the required buffer and returns that buffer if do_allocate is called.
 *  Morevover, this pool has the same 'delayed' properties than the `delayed_pool_type`.
 *
 *  \sa delayed_pool_type
 *  \note std::pmr compatible
 */
template <class Upstream>
class delayed_fragmenting_pool_type : public std::pmr::memory_resource {
   private:
	struct page_item_type {
		void* ptr;
		bool allocated;
	};

   public:
	delayed_fragmenting_pool_type() = default;

	~delayed_fragmenting_pool_type() noexcept {
		for (auto& [page_id, book_page] : book_) {
			for (auto& page_item : book_page) {
				res_.deallocate(
				    page_item.ptr, page_id.bytes, page_id.alignment);
			}
		}
	}

	auto& get_book() const {
		return book_;
	}

   private:
	void* do_allocate(size_t bytes, size_t alignment) override {
		if (bytes == 0) {
			return nullptr;
		}
		const detail::page_id_type page_id({bytes, alignment});
		for (auto& [k, v] : book_) {
			if (k.alignment == alignment && k.bytes >= bytes) {
				for (auto& p : v) {
					if (!p.allocated) {
						p.allocated = true;
						return p.ptr;
					}
				}
			}
		}
		// no page item was free
		if (auto it = book_.find(page_id); it != book_.end()) {
			void* ptr = res_.allocate(bytes, alignment);
			it->second.push_back({ptr, true});
			return ptr;
		}
		// the required size was never allocated before
		void* ptr = res_.allocate(bytes, alignment);
		book_[page_id] = {{ptr, true}};
		return ptr;
	}

	void do_deallocate(void* ptr,
	                   size_t bytes,
	                   size_t alignment) noexcept override {
		const detail::page_id_type page_id({bytes, alignment});

		if (ptr == nullptr || bytes == 0) {
			return;
		}

		for (auto& [k, v] : book_) {
			if (k.alignment == alignment && k.bytes >= bytes) {
				for (auto& p : v) {
					// no comparison based on the bytes because the buffer can be larger
					// than the amount of requested bytes.
					if (p.ptr == ptr && p.allocated) {
						p.allocated = false;
						return;
					}
				}
			}
		}

		// The pointer which should be deallocated was never allocated before
		std::terminate();
	}

	bool do_is_equal(
	    const std::pmr::memory_resource& other) const noexcept override {
		return this == &other;
	}

	Upstream res_;
	// size in bytes -> (ptr, alignment)
	std::map<detail::page_id_type, std::vector<page_item_type>> book_;
};

/*! \brief Wraps a memory resource around an existing buffer of fixed size.
 *
 *  This class takes a pointer and a size of an existing buffer and wraps a memory
 *  resource around that buffer. The resource can provide parts of that buffer by
 *  calling `do_allocate` and simply returns subsequently aligned parts of the buffer.
 *  If the buffer size is exceeded `gsl_FailFast()` is called.
 *
 *  \note std::pmr compatible
 */
class wrapping_resource_type : public std::pmr::memory_resource {

   public:
	wrapping_resource_type(void* ptr, size_t size_in_bytes)
	    : size_in_bytes_(size_in_bytes), ptr_(ptr) {
		static_assert(sizeof(std::byte) == 1);
	}

   private:
	void* do_allocate(size_t bytes, size_t alignment) noexcept override {
		if (bytes == 0) {
			return nullptr;
		}
		void* ptr = std::align(alignment, bytes, ptr_, size_in_bytes_);
		gsl_Expects(ptr); // not nullptr
		gsl_Expects(
		    ptr ==
		    ptr_); // this is how I understand the function documentation of std::align
		gsl_Expects(size_in_bytes_ >= bytes);
		size_in_bytes_ -= bytes;
		ptr_ = reinterpret_cast<std::byte*>(ptr_) + bytes;
		return ptr;
	}

	void do_deallocate([[maybe_unused]] void* ptr,
	                   [[maybe_unused]] size_t bytes,
	                   [[maybe_unused]] size_t alignment) noexcept override {
	}

	bool do_is_equal(
	    const std::pmr::memory_resource& other) const noexcept override {
		return this == &other;
	}

	size_t size_in_bytes_;
	void* ptr_;
};

static managed_resource_type default_resource;

} // namespace pmr

} // namespace thrustshift

max_stars_repo_name: codecircuit/thrustshift

max_issues_repo_name: codecircuit/thrustshift

max_forks_repo_name: codecircuit/thrustshift