text:
// Copyright (C) 2022. Huawei Technologies Co., Ltd. All rights reserved.

// 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.

#ifndef COMMON_H
#define COMMON_H

#include <cmath>
#include <optional>
#include <yato/array_view.h>

//#define CHECK_ASSERT

#include <training/system/Errors.h>
#include <training/system/Name.h>
#include <training/system/TypeHalf.h>
#include <training/system/Types.h>

#include <map>
#include <utility>

#if defined(_BLAS) && !defined(_BLAS_ENHANCE)
extern "C"
{
#include <cblas.h>
}
#else
#ifndef OPENBLAS_CONST
#define OPENBLAS_CONST const
#endif

typedef enum CBLAS_TRANSPOSE
{
    CblasNoTrans = 111,
    CblasTrans = 112,
    CblasConjTrans = 113,
    CblasConjNoTrans = 114
} CBLAS_TRANSPOSE;

typedef enum CBLAS_UPLO
{
    CblasUpper = 121,
    CblasLower = 122
} CBLAS_UPLO;
#endif

#ifdef CHECK_ASSERT
#define CHECK_NEAR ASSERT_NEAR
#else
#define CHECK_NEAR EXPECT_NEAR
#endif

#define RAUL_E 2.71828182845904523536        // e
#define RAUL_LOG2E 1.44269504088896340736    // log2(e)
#define RAUL_LOG10E 0.434294481903251827651  // log10(e)
#define RAUL_LN2 0.693147180559945309417     // ln(2)
#define RAUL_LN10 2.30258509299404568402     // ln(10)
#define RAUL_PI 3.14159265358979323846       // pi
#define RAUL_PI_2 1.57079632679489661923     // pi/2
#define RAUL_PI_4 0.785398163397448309616    // pi/4
#define RAUL_1_PI 0.318309886183790671538    // 1/pi
#define RAUL_2_PI 0.636619772367581343076    // 2/pi
#define RAUL_2_SQRTPI 1.12837916709551257390 // 2/sqrt(pi)
#define RAUL_SQRT2_PI 0.79788456080286535588 // sqrt(2/pi)
#define RAUL_SQRT2 1.41421356237309504880    // sqrt(2)
#define RAUL_SQRT1_2 0.707106781186547524401 // 1/sqrt(2)

#define GELU_CONST 0.044715

namespace raul
{

enum class Limit : int
{
    Left = 0,
    Middle = 1,
    Right = 2
};

enum class Dimension : int
{
    Default = -1,
    Batch = 0,
    Depth = 1,
    Height = 2,
    Width = 3
};

#if defined(_MSC_VER)
#define INLINE __forceinline
#else
#define INLINE __attribute__((always_inline))
#endif

template<typename Type>
class TensorImpl;
typedef TensorImpl<dtype> Tensor;
typedef TensorImpl<half> TensorFP16;

#if defined(ANDROID)
#define TOMMTYPE(var) static_cast<typename MM::type>(var)
#else
#define TOMMTYPE(var) castHelper<typename MM::type>::cast(var)
#endif

using shape = yato::dimensionality<4U, size_t>;
} // raul namespace

namespace raul
{

enum class NetworkMode
{
    Train = 0,
    Test = 1,
    TrainCheckpointed = 2
};

enum class CompressionMode
{
    NONE = -1,
    FP16 = 0,
    INT8 = 1
};

enum class CalculationMode
{
    DETERMINISTIC = 0,
#if defined(_OPENMP)
    FAST = 1,
#endif
};

/**
 * @brief Hardware target platform
 *
 */
enum class ExecutionTarget
{
    CPU = 0,
    CPUFP16 = 1
};

/**
 * @brief Hardware target platform per layer
 *
 * \note Might override execution target for workflow, useful for mixed precision
 */
enum class LayerExecutionTarget
{
    Default = -1, // use same as ExecutionTarget
    CPU = 0,      // from this point enums should be aligned with ExecutionTarget (due to LayerExecutionTarget = static_cast<ExecutionTarget>(enum))
    CPUFP16 = 1
};

/**
 * @brief Memory allocation mode
 */
enum class AllocationMode
{
    STANDARD,
    POOL
};

enum class DeclarationType
{
    Tensor = 0,
    Shape = 1,
    //    Alias = 2
};

class OpenclInitializer;

class Common
{
  public:

    // generate vector of random index permutation of [0..n-1]
    static void generate_permutation(size_t n, std::vector<size_t>& ind_vector, unsigned int seed = 0);

    /*
     * [cols x rows]
     * A[k x m]
     * B[n x k]
     * C[n x m]
     * https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemm
     * C = alpha * A * B + beta * C
     * bOffset - in elements (not bytes)
     */
    static void gemm(OPENBLAS_CONST CBLAS_TRANSPOSE transA,
                     OPENBLAS_CONST CBLAS_TRANSPOSE transB,
                     size_t m,
                     size_t n,
                     size_t k,
                     OPENBLAS_CONST dtype alpha,
                     OPENBLAS_CONST dtype* a,
                     OPENBLAS_CONST dtype* b,
                     OPENBLAS_CONST dtype beta,
                     dtype* c);

    static void gemm(OPENBLAS_CONST CBLAS_TRANSPOSE transA,
                     OPENBLAS_CONST CBLAS_TRANSPOSE transB,
                     size_t m,
                     size_t n,
                     size_t k,
                     OPENBLAS_CONST dtype alpha,
                     OPENBLAS_CONST half* a,
                     OPENBLAS_CONST half* b,
                     OPENBLAS_CONST dtype beta,
                     half* c);

    /**
     * @brief : Basic Linear Algebra Subroutine y = y + ax
     *
     *  \f[
     *      \vec{y} = \vec{y} + \alpha * \vec{x},
     *  \f]
     *
     * @param n The number of elements in vectors x and y.
     * @param sa The scalar alpha.
     * @param sx The vector x of length n. Specified as: a one-dimensional array of (at least) length \f$ 1+(n-1)|incx| \f$.
     * @param incx The stride for vector x. Specified as: an integer. It can have any value.
     * @param sy The vector y of length n. Specified as: a one-dimensional array of (at least) length \f$ 1+(n-1)|incy| \f$.
     * @param incy The stride for vector y.
     * @param xOffset The offset for vector x.
     * @param yOffset The offset for vector y.
     * @return The vector y, containing the results of the computation.
     */
    static void axpy(size_t n, OPENBLAS_CONST dtype sa, OPENBLAS_CONST dtype* sx, size_t incx, dtype* sy, size_t incy, size_t xOffset = 0, size_t yOffset = 0);
    static void axpy(size_t n, OPENBLAS_CONST dtype sa, OPENBLAS_CONST half* sx, size_t incx, half* sy, size_t incy, size_t xOffset = 0, size_t yOffset = 0);

    /**
     * @brief : Basic Linear Algebra Subroutine y = ax + by
     *
     *  \f[
     *      \vec{y} = \alpha \vec{x} + \beta \vec{y},
     *  \f]
     *
     * @param n The number of elements in vectors x and y.
     * @param alpha The scalar alpha.
     * @param x The vector x of length n. Specified as: a one-dimensional array of (at least) length \f$ 1+(n-1)|incx| \f$.
     * @param incx The stride for vector x. Specified as: an integer. It can have any value.
     * @param beta The scalar beta.
     * @param y The vector y of length n. Specified as: a one-dimensional array of (at least) length \f$ 1+(n-1)|incy| \f$.
     * @param incy The stride for vector y.
     * @param xOffset The offset for vector x.
     * @param yOffset The offset for vector y.
     * @return The vector y, containing the results of the computation.
     */
    static int axpby(OPENBLAS_CONST size_t n,
                     OPENBLAS_CONST dtype alpha,
                     OPENBLAS_CONST dtype* x,
                     OPENBLAS_CONST size_t incx,
                     OPENBLAS_CONST dtype beta,
                     dtype* y,
                     OPENBLAS_CONST size_t incy,
                     size_t xOffset,
                     size_t yOffset);

    static int axpby(OPENBLAS_CONST size_t n,
                     OPENBLAS_CONST dtype alpha,
                     OPENBLAS_CONST half* x,
                     OPENBLAS_CONST size_t incx,
                     OPENBLAS_CONST dtype beta,
                     half* y,
                     OPENBLAS_CONST size_t incy,
                     size_t xOffset,
                     size_t yOffset);
    /**
     * @brief : Basic Linear Algebra Subroutine y = alpha * a * x + beta * y
     *
     * Vector by vector element wise multiplication
     *
     *  \f[
     *      \vec{y} = \alpha \vec{a} \vec{x} + \beta \vec{y},
     *  \f]
     *
     * @param n The number of elements in vectors x and y.
     * @param alpha The scalar alpha.
     * @param a The vector of length n.
     * @param x The vector x of length n. Specified as: a one-dimensional array of (at least) length \f$ 1+(n-1)|incx| \f$.
     * @param incx The stride for vector x. Specified as: an integer. It can have any value.
     * @param beta The scalar beta.
     * @param y The vector y of length n. Specified as: a one-dimensional array of (at least) length \f$ 1+(n-1)|incy| \f$.
     * @param incy The stride for vector y.
     */

    static void hadamard(OPENBLAS_CONST size_t n,
                         OPENBLAS_CONST dtype alpha,
                         OPENBLAS_CONST dtype* a,
                         OPENBLAS_CONST dtype* x,
                         OPENBLAS_CONST size_t incx,
                         OPENBLAS_CONST dtype beta,
                         dtype* y,
                         OPENBLAS_CONST size_t incy);

    static dtype dot(size_t n, OPENBLAS_CONST dtype* sx, size_t incx, OPENBLAS_CONST dtype* sy, size_t incy);

    static void scal(size_t n, OPENBLAS_CONST dtype sa, dtype* sx, size_t incx);

    static void transpose(Tensor& tensor, size_t cols);
    static void transpose(TensorFP16& tensor, size_t cols);

    /*
     * memory for dst should be allocated externaly
     */
    static void addPadding1D(const dtype* src, dtype* dst, size_t srcChannels, size_t srcSize, size_t dstSize, bool reversedOrder = false);
    template<typename T>
    static void addPadding2D(const T* src, T* dst, size_t srcChannels, size_t srcWidth, size_t srcHeight, size_t dstWidth, size_t dstHeight)
    {
        if ((dstWidth >= srcWidth) && (dstHeight >= srcHeight))
        {
            size_t padWidth = dstWidth - srcWidth;
            size_t padHeight = dstHeight - srcHeight;

            size_t leftPad = padWidth / 2;
            // size_t rightPad = padWidth - leftPad;
            size_t topPad = padHeight / 2;
            size_t bottomPad = padHeight - topPad;

            for (size_t d = 0; d < srcChannels; ++d)
            {
                // top
                for (size_t y = 0; y < topPad; ++y)
                {
                    for (size_t x = 0; x < dstWidth; ++x)
                    {
                        dst[d * dstWidth * dstHeight + dstWidth * y + x] = static_cast<T>(0.0_dt);
                    }
                }

                for (size_t y = topPad; y < topPad + srcHeight; ++y)
                {
                    // left
                    for (size_t x = 0; x < leftPad; ++x)
                    {
                        dst[d * dstWidth * dstHeight + dstWidth * y + x] = static_cast<T>(0.0_dt);
                    }

                    // src
                    for (size_t x = leftPad; x < leftPad + srcWidth; ++x)
                    {
                        dst[d * dstWidth * dstHeight + dstWidth * y + x] = src[d * srcWidth * srcHeight + srcWidth * (y - topPad) + x - leftPad];
                    }

                    // right
                    for (size_t x = leftPad + srcWidth; x < dstWidth; ++x)
                    {
                        dst[d * dstWidth * dstHeight + dstWidth * y + x] = static_cast<T>(0.0_dt);
                    }
                }

                // bottom
                for (size_t y = dstHeight - bottomPad; y < dstHeight; ++y)
                {
                    for (size_t x = 0; x < dstWidth; ++x)
                    {
                        dst[d * dstWidth * dstHeight + dstWidth * y + x] = static_cast<T>(0.0_dt);
                    }
                }
            }
        }
    }

    /*
     * memory for dst should be allocated externaly
     */
    static void removePadding1D(const dtype* src, dtype* dst, size_t srcChannels, size_t srcSize, size_t dstSize, bool reversedOrder = false, bool overwrite = true);
    template<typename T>
    static void removePadding2D(const T* src, T* dst, size_t srcChannels, size_t srcWidth, size_t srcHeight, size_t dstWidth, size_t dstHeight, bool overwrite = true)
    {
        if ((dstWidth <= srcWidth) && (dstHeight <= srcHeight))
        {
            size_t padWidth = srcWidth - dstWidth;
            size_t padHeight = srcHeight - dstHeight;

            size_t leftPad = padWidth / 2;
            // size_t rightPad = padWidth - leftPad;
            size_t topPad = padHeight / 2;
            // size_t bottomPad = padHeight - topPad;

            if (overwrite)
            {
                for (size_t d = 0; d < srcChannels; ++d)
                {
                    for (size_t y = 0; y < dstHeight; ++y)
                    {
                        for (size_t x = 0; x < dstWidth; ++x)
                        {
                            dst[d * dstWidth * dstHeight + dstWidth * y + x] = src[d * srcWidth * srcHeight + srcWidth * (y + topPad) + x + leftPad];
                        }
                    }
                }
            }
            else
            {
                for (size_t d = 0; d < srcChannels; ++d)
                {
                    for (size_t y = 0; y < dstHeight; ++y)
                    {
                        for (size_t x = 0; x < dstWidth; ++x)
                        {
                            dst[d * dstWidth * dstHeight + dstWidth * y + x] += src[d * srcWidth * srcHeight + srcWidth * (y + topPad) + x + leftPad];
                        }
                    }
                }
            }
        }
    }

    /*
     * paddingWidth, paddingHeight - zero padding added for both sides of the input
     * memory for matrix should be allocated externaly
     */
    template<typename T>
    static void im2col(const T* image,
                       size_t imageWidth,
                       size_t imageHeight,
                       size_t imageChannels,
                       size_t filterWidth,
                       size_t filterHeight,
                       size_t strideWidth,
                       size_t strideHeight,
                       size_t paddingWidth,
                       size_t paddingHeight,
                       T* matrix,
                       bool reversedOrder = false);

    static size_t im2colOutputSize(size_t imageWidth,
                                   size_t imageHeight,
                                   size_t imageChannels,
                                   size_t filterWidth,
                                   size_t filterHeight,
                                   size_t strideWidth,
                                   size_t strideHeight,
                                   size_t paddingWidth,
                                   size_t paddingHeight,
                                   size_t dilationWidth,
                                   size_t dilationHeight);

    /*
     * paddingWidth, paddingHeight - zero padding added for both sides of the input
     * memory for image should be allocated externaly
     */
    template<typename T>
    static void col2im(const T* matrix,
                       size_t imageWidth,
                       size_t imageHeight,
                       size_t imageChannels,
                       size_t filterWidth,
                       size_t filterHeight,
                       size_t strideWidth,
                       size_t strideHeight,
                       size_t paddingWidth,
                       size_t paddingHeight,
                       T* image,
                       bool reversedOrder = false,
                       bool zeroOutput = true);

    /*
     * Rectified Linear Unit
     */
    template<typename T>
    static T ReLU(T x)
    {
        return std::max(static_cast<T>(0), x);
    }
    template<typename T>
    static T ReLU6(T x)
    {
        return std::min(std::max(static_cast<T>(0), x), static_cast<T>(6.0_dt));
    }

    template<typename T>
    static void ReLU(const T& in, T& out)
    {
        std::transform(in.begin(), in.end(), out.begin(), [&](typename T::type val) -> typename T::type { return ReLU(val); });
    }

    template<typename T>
    static void ReLU6(const T& in, T& out)
    {
        std::transform(in.begin(), in.end(), out.begin(), [&](typename T::type val) -> typename T::type { return ReLU6(val); });
    }

    template<typename T>
    static void ReLUBackward(const T& out, const T& delta, T& prevDelta)
    {
#if defined(_OPENMP)
#pragma omp parallel for
#endif
        for (size_t q = 0; q < prevDelta.size(); ++q)
        {
            prevDelta[q] += (out[q] > static_cast<typename T::type>(0)) ? delta[q] : static_cast<typename T::type>(0);
        }
    }

    template<typename T>
    static void ReLU6Backward(const T& out, const T& delta, T& prevDelta)
    {
#if defined(_OPENMP)
#pragma omp parallel for
#endif
        for (size_t q = 0; q < prevDelta.size(); ++q)
        {
            prevDelta[q] += (out[q] > static_cast<typename T::type>(0) && out[q] < static_cast<typename T::type>(6.0f)) ? delta[q] : static_cast<typename T::type>(0);
        }
    }

    /*
     * Gaussian error linear unit
     * @see https://arxiv.org/abs/1606.08415
     */
    static dtype GeLU_Erf(dtype x);
    static dtype GeLU_Tanh(dtype x);

    /*
     * Hard Sigmoid
     */
    template<typename T>
    static T HSigmoid(T x)
    {
        return static_cast<T>(ReLU6(TODTYPE(x) + 3.0_dt) / 6.0_dt);
    }

    /*
     * Hard Swish
     */
    template<typename T>
    static T HSwish(T x)
    {
        return x * HSigmoid(x);
    }
    static dtype sign(dtype x) { return TODTYPE((0.0_dt < x) - (x < 0.0_dt)); }

    template<typename T, typename U>
    static void copyView(const T& view_from, U& view_to, const bool overwrite = false)
    {
        auto retLhs = [](typename T::value_type& lhs, [[maybe_unused]] typename T::value_type& rhs) { return lhs; };

        auto copyViewImpl = [](const T& view_from, U& view_to, auto&& func)
        {
            for (size_t i1 = 0; i1 < view_from.size(0); ++i1)
            {
                for (size_t i2 = 0; i2 < view_from.size(1); ++i2)
                {
                    for (size_t i3 = 0; i3 < view_from.size(2); ++i3)
                    {
                        for (size_t i4 = 0; i4 < view_from.size(3); ++i4)
                        {
                            view_to[i1][i2][i3][i4] = func(view_from[i1][i2][i3][i4], view_to[i1][i2][i3][i4]);
                        }
                    }
                }
            }
        };

        if (overwrite)
        {
            copyViewImpl(view_from, view_to, retLhs);
        }
        else
        {
            copyViewImpl(view_from, view_to, std::plus<typename T::value_type>());
        }
    }

    template<typename T>
    static void unpack4D(const T& src, T& dst, Dimension dir, size_t index, const Name& layerType, const Name& layerName, bool overwrite)
    {
        auto input4d = src.get4DView();
        auto inputDims = yato::dims(src.getDepth(), src.getHeight(), src.getWidth());

        auto outputDims = dst.getShape();

        const typename T::type* startEl = nullptr;
        switch (dir)
        {
            case Dimension::Depth:
                startEl = &input4d[0][index][0][0];
                break;
            case Dimension::Height:
                startEl = &input4d[0][0][index][0];
                break;
            default:
                throw std::runtime_error(layerType + "[" + layerName + "]: unpack4D unknown dim");
        }

        auto srcView = yato::array_view_4d<const typename T::type>(startEl, outputDims, inputDims);
        auto outputView = dst.get4DView();
        Common::copyView(srcView, outputView, overwrite);
    }

    template<typename T>
    static void pack4D(const T& src, T& dst, Dimension dir, size_t index, const Name& layerType, const Name& layerName, bool overwrite)
    {
        auto output4d = dst.get4DView();

        yato::dimensionality<3U, size_t> concatDims(dst.getDepth(), dst.getHeight(), dst.getWidth());

        auto srcView = src.get4DView();
        typename T::type* startEl = nullptr;
        switch (dir)
        {
            case Dimension::Depth:
                startEl = &output4d[0][index][0][0];
                break;
            case Dimension::Height:
                startEl = &output4d[0][0][index][0];
                break;
            default:
                throw std::runtime_error(layerType + "[" + layerName + "]: pack4D unknown dim");
        }

        auto dstView = yato::array_view_4d<typename T::type>(startEl, src.getShape(), concatDims);
        Common::copyView(srcView, dstView, overwrite);
    }

    /*
     * Upper triangle of a rectangular array
     */
    template<typename T>
    static void triu(T* data, size_t nrows, size_t ncols, int diag = 0)
    {
        size_t i = 0;
        int cols = (int)ncols;
        int rows = (int)nrows;
        for (int r = 0; r < rows; ++r)
        {
            for (int c = 0; c < cols; ++c, ++i)
            {
                if (c - r - diag < 0)
                {
                    data[i] = static_cast<T>(0);
                }
            }
        }
    }

    /*
     * Applies a 1D convolution over an input signal composed of several input planes.
     * Supports 2 modes:
     *  1. PyTorch style: Input[N, C, 1, L1] (or [N, 1, C, L1]) -> Output[N, FILTERS, 1, L2] (or [N, 1, FILTERS, L2])
     *  2. TensorFlow style: Input[N, L1, 1, C] (or [N, 1, L1, C]) -> Output[N, L2, 1, FILTERS] (or [N, 1, L2, FILTERS])
     * Output is not zeroed prior to convolution (operator += is used)
     */
    static void conv1d(const dtype* input,
                       dtype* output,
                       const dtype* kernel,
                       const dtype* bias,
                       size_t batchSize,
                       size_t inputSize,
                       size_t inputChannels,
                       size_t outputSize,
                       size_t outputChannels,
                       size_t kernelSize,
                       size_t padding,
                       size_t stride,
                       size_t dilation = 1U,
                       size_t groups = 1U,
                       bool tfStyle = false);

    /*
     * Applies 2D convolution over input tensor, all channels convolved
     * Output is not zeroed prior to convolution (operator += is used)
     */
    template<typename T>
    static void conv2d(const T* input,
                       T* output,
                       const T* kernel,
                       const T* bias,
                       size_t batchSize,
                       size_t inputWidth,
                       size_t inputHeight,
                       size_t inputChannels,
                       size_t outputWidth,
                       size_t outputHeight,
                       size_t outputChannels,
                       size_t kernelWidth,
                       size_t kernelHeight,
                       size_t paddingW,
                       size_t paddingH,
                       size_t strideW,
                       size_t strideH,
                       size_t dilationW = 1U,
                       size_t dilationH = 1U,
                       size_t groups = 1U)
    {
        auto inputs3D = yato::array_view_3d<T>(const_cast<T*>(input), yato::dims(batchSize, inputChannels, inputHeight * inputWidth));
        auto outputs3D = yato::array_view_3d<T>(output, yato::dims(batchSize, outputChannels, outputHeight * outputWidth));
        auto kernelsWeights4D = yato::array_view_4d<T>(const_cast<T*>(kernel), yato::dims(outputChannels, inputChannels / groups, kernelHeight, kernelWidth));

        for (size_t q = 0; q < batchSize; ++q)
        {
            for (size_t d = 0; d < outputChannels; ++d)
            {
                std::fill(outputs3D[q][d].begin(), outputs3D[q][d].end(), static_cast<T>(0.0_dt));
            }

            size_t inputWidthPadded = inputWidth + 2 * paddingW;
            size_t inputHeightPadded = inputHeight + 2 * paddingH;

            std::vector<T> inputPadded(inputChannels * inputHeightPadded * inputWidthPadded);

            Common::addPadding2D(&inputs3D[q][0][0], inputPadded.data(), inputChannels, inputWidth, inputHeight, inputWidthPadded, inputHeightPadded);

            auto inputPadded2D = yato::view(inputPadded).reshape(yato::dims(inputChannels, inputHeightPadded * inputWidthPadded));

            for (size_t group = 0; group < groups; ++group)
            {
                for (size_t kernelIndex = 0; kernelIndex < outputChannels / groups; ++kernelIndex)
                {
                    for (size_t d = 0; d < inputChannels / groups; ++d)
                    {
                        for (size_t oy = 0; oy < outputHeight; ++oy)
                        {
                            for (size_t ox = 0; ox < outputWidth; ++ox)
                            {
                                for (size_t ky = 0; ky < kernelHeight; ++ky)
                                {
                                    for (size_t kx = 0; kx < kernelWidth; ++kx)
                                    {
                                        outputs3D[q][kernelIndex + group * outputChannels / groups][oy * outputWidth + ox] +=
                                            kernelsWeights4D[kernelIndex + group * outputChannels / groups][d][ky][kx] *
                                            inputPadded2D[d + group * inputChannels / groups][oy * inputWidthPadded * strideH + ky * dilationH * inputWidthPadded + ox * strideW + kx * dilationW];
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }

        if (bias)
        {
            for (size_t q = 0; q < batchSize; ++q)
            {
                for (size_t kernelIndex = 0; kernelIndex < outputChannels; ++kernelIndex)
                {
                    for (size_t oy = 0; oy < outputHeight; ++oy)
                    {
                        for (size_t ox = 0; ox < outputWidth; ++ox)
                        {
                            outputs3D[q][kernelIndex][oy * outputWidth + ox] += bias[kernelIndex];
                        }
                    }
                }
            }
        }
    }

    template<typename T = dtype, typename Iterator>
    static void arange(Iterator begin, Iterator end, T start = static_cast<T>(0), T step = static_cast<T>(1))
    {
        auto val = start;
        for (auto p = begin; p != end; ++p)
        {
            *p = static_cast<std::remove_reference_t<decltype(*p)>>(val);
            val += step;
        }
    }

    template<typename T = dtype, typename Iterable>
    static void arange(Iterable& i, T start = static_cast<T>(0), T step = static_cast<T>(1))
    {
        return arange(i.begin(), i.end(), start, step);
    }

    static void replaceAll(std::string& str, const std::string& srcSubstr, const std::string& tgtSubstr)
    {
        size_t start_pos = 0;
        while ((start_pos = str.find(srcSubstr, start_pos)) != std::string::npos)
        {
            str.replace(start_pos, srcSubstr.length(), tgtSubstr);
            start_pos += tgtSubstr.length(); // srcSubstr could be a substring of tgtSubstr
        }
    }

    static bool startsWith(const std::string& str, const std::string& srcSubstr) { return (str.rfind(srcSubstr, 0) == 0); }

    static std::vector<std::string> split(const std::string& string, char delimeter);

    /*
     * @see https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html
     */
    template<typename T>
    static bool shapeIsBroadcastable(const T& from, const T& to)
    {
        const auto n = to.dimensions_num();
        for (size_t i = 0; i < n; ++i)
        {
            if (from[i] != to[i] && from[i] != 1U && to[i] != 1U)
            {
                return false;
            }
        }
        return true;
    }

    static bool endsWith(std::string const& value, std::string const& ending)
    {
        if (ending.size() > value.size())
        {
            return false;
        }
        return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
    }

    static shape getStrides(const shape& tensor_shape);

    static shape offsetToIndexes(size_t offset, const shape& strides);

    static size_t indexesToOffset(const shape& indexes, const shape& strides);
};

template<class T>
bool if_equals(const std::string&& error, const T val1, const T val2)
{
    if (val1 != val2)
    {
        throw(std::runtime_error(error));
    }
    return val1 == val2;
}

} // raul namespace

#endif // COMMON_H

max_stars_repo_name: steelONIONknight/bolt

max_issues_repo_name: steelONIONknight/bolt

max_forks_repo_name: steelONIONknight/bolt


text:
#pragma once

#include "rev/Utilities.h"
#include <algorithm>
#include <cstddef>
#include <glm/glm.hpp>
#include <gsl/gsl_assert>
#include <gsl/span>
#include <limits>
#include <vector>

namespace rev {

template <typename PointType>
struct WeightedControlPoint {
    PointType point;
    typename PointType::value_type weight;
};

template <typename PointType>
class NurbsCurve {
public:
    using ControlPointType = WeightedControlPoint<PointType>;
    using ValueType = typename PointType::value_type;

    NurbsCurve(size_t order, gsl::span<const ValueType> knots,
        gsl::span<const ControlPointType> controlPoints)
        : _order(order)
        , _knots(knots.begin(), knots.end())
        , _controlPoints(controlPoints.begin(), controlPoints.end())
    {
        Expects(order > 1);
        Expects(static_cast<size_t>(controlPoints.size()) >= order);
        Expects(knots.size() == (controlPoints.size() + order));
    }

    ValueType getStart() const { return _knots.front(); }

    ValueType getEnd() const { return _knots.back(); }

    PointType operator[](ValueType position) const
    {
        auto knotsBegin = _knots.begin();
        auto knotsEnd = _knots.end();
        auto spanStart = std::lower_bound(knotsBegin, knotsEnd, position);
        size_t controlPointIndex = std::distance(knotsBegin, spanStart);
        if (controlPointIndex == 0) {
            // Before the first knot. Just return the value of the first control
            // point.
            return _controlPoints.front().point;
        }

        if (spanStart == knotsEnd) {
            // We're past the last knot. Just return the value of the last control
            // point.
            return _controlPoints.back().point;
        }

        size_t controlPointsToProcess = _order;
        size_t controlPointCount = _controlPoints.size();
        if (controlPointIndex > controlPointCount) {
            controlPointsToProcess -= (controlPointIndex - controlPointCount);
        }

        PointType numerator{ 0 };
        ValueType denominator{ 0 };
        do {
            controlPointIndex--;
            auto& controlPoint = _controlPoints[controlPointIndex];
            ValueType weightedBasis
                = basisValue(position, controlPointIndex, _order) * controlPoint.weight;
            numerator += weightedBasis * _controlPoints[controlPointIndex].point;
            denominator += weightedBasis;
            controlPointsToProcess--;
        } while (controlPointIndex && controlPointsToProcess);
        return numerator / denominator;
    }

private:
    ValueType basisValue(ValueType position, size_t controlPointIndex, size_t order) const
    {
        auto& controlPoint = _controlPoints[controlPointIndex];
        ValueType firstKnot = _knots[controlPointIndex];
        ValueType lastKnot = _knots[controlPointIndex + order];
        if (order == 2) {
            // Triangle function
            ValueType middleKnot = _knots[controlPointIndex + 1];
            if (position < middleKnot) {
                return rampUp(position, firstKnot, middleKnot);
            } else {
                return rampDown(position, middleKnot, lastKnot);
            }
        } else {
            size_t lowerOrder = order - 1;
            ValueType secondKnot = _knots[controlPointIndex + 1];
            ValueType secondToLastKnot = _knots[controlPointIndex + lowerOrder];

            ValueType result{ 0 };
            if (position < secondToLastKnot) {
                result += rampUp(position, firstKnot, secondToLastKnot)
                    * basisValue(position, controlPointIndex, lowerOrder);
            }

            if (position > secondKnot) {
                result += rampDown(position, secondKnot, lastKnot)
                    * basisValue(position, controlPointIndex + 1, lowerOrder);
            }

            return result;
        }
    }

    ValueType rampUp(ValueType value, ValueType bottom, ValueType top) const
    {
        if (nearEqual(top, bottom)) {
            return 0.5;
        }
        return (value - bottom) / (top - bottom);
    }

    ValueType rampDown(ValueType value, ValueType bottom, ValueType top) const
    {
        if (nearEqual(top, bottom)) {
            return 0.5;
        }
        return (top - value) / (top - bottom);
    }

    size_t _order;
    std::vector<ValueType> _knots;
    std::vector<ControlPointType> _controlPoints;
};
} // namespace rev

ext: h

max_stars_repo_name: eyebrowsoffire/rev

max_issues_repo_name: eyebrowsoffire/rev

max_forks_repo_name: eyebrowsoffire/rev


text:
//
// Created by karl_ on 2020-5-3.
//

#include <stdlib.h>
#include <math.h>
#include <cblas.h>
#include <core/color-space.h>
#include <imgproc/histogram-equalization.h>

// ref: https://www.math.uci.edu/icamp/courses/math77c/demos/hist_eq.pdf
static void histeq(double const *src, double *dst, size_t area, unsigned range) {
    unsigned *cdf;
    size_t i;
    cdf = calloc(range, sizeof(unsigned));
    for (i = 0; i < area; ++i) {
        cdf[(int) (src[i] * range)] += 1;
    }
    for (i = 1; i < range; ++i) {
        cdf[i] += cdf[i - 1];
    }
    for (i = 0; i < area; ++i) {
        dst[i] = floor((double) (range - 1) * cdf[(int) (src[i] * range)] / area) / range;
    }
    free(cdf);
}

gboolean imgproc_histogram_equalization(CoreImage *src, CoreImage **dst) {
    gdouble *src_data_cast = NULL;
    gdouble *dst_data_cast = NULL;
    gpointer src_data, dst_data;
    CoreSize *size = NULL;
    guint8 channel;
    gsize i, area;
    CoreColorSpace color_space;
    CorePixelType pixel_type;
    gdouble range, inv_range;
    g_return_val_if_fail(src != NULL, FALSE);
    g_return_val_if_fail(dst != NULL, FALSE);

    /* TODO: impl. for multiple channels */
    channel = core_image_get_channel(src);
    g_return_val_if_fail(channel == 1, FALSE);

    src = g_object_ref(src);
    src_data = core_image_get_data(src);
    size = core_image_get_size(src);
    area = core_size_get_area(size);
    color_space = core_image_get_color_space(src);
    pixel_type = core_image_get_pixel_type(src);
    range = core_pixel_get_range(pixel_type);
    inv_range = 1.0 / range;

    /* cast to double and normalize pixel data to [0, 1] */
    if (core_pixel_is_double(pixel_type)) {
        src_data_cast = src_data;
    } else if (core_pixel_is_uint8(pixel_type)) {
        src_data_cast = g_malloc(sizeof(gdouble) * area);
        for (i = 0; i < area; ++i) {
            src_data_cast[i] = ((guint8 *) src_data)[i] * inv_range;
        }
    } else {
        g_return_val_if_fail(FALSE, FALSE);
    }
    dst_data_cast = g_malloc(sizeof(gdouble) * area);

    /* run the algorithm */
    histeq(src_data_cast, dst_data_cast, area, 256);

    /* cast back to uchar */
    if (core_pixel_is_double(pixel_type)) {
        /* exactly do nothing */
    } else if (core_pixel_is_uint8(pixel_type)) {
        /* reuse allocated memory */
        for (i = 0; i < area; ++i) {
            ((guint8 *) dst_data_cast)[i] = dst_data_cast[i] * 255.0;
        }
    } else {
        g_return_val_if_fail(FALSE, FALSE);
    }
    dst_data = dst_data_cast;
    if (core_pixel_is_uint8(pixel_type)) {
        /* release useless memory, can be removed for speed */
        dst_data = g_realloc(dst_data, sizeof(guint8) * area);
    }

    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, color_space, pixel_type, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, color_space, pixel_type, size, FALSE);
    }

    /* free allocated additional memory */
    if (core_pixel_is_uint8(pixel_type)) {
        g_free(src_data_cast);
    }
    g_object_unref(size);
    g_object_unref(src);
    return TRUE;
}

ext: c

max_stars_repo_name: zhengbanbufengbi/EasyPhotoshop

max_issues_repo_name: zhengbanbufengbi/EasyPhotoshop

max_forks_repo_name: zhengbanbufengbi/EasyPhotoshop


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

#ifndef __CCPM_AREATOP_ALLOCATOR_H__
#define __CCPM_AREATOP_ALLOCATOR_H__

#include "atomic_word.h"
#include "list_item.h"
#include <ccpm/interfaces.h>
#include <common/byte_span.h>
#include <gsl/pointers>
#include <array>
#include <cstddef>
#include <ios> // ios_base::fmtflags, ostream
#include <vector>

namespace ccpm
{
	struct level_hints
	{
	private:
		/* Each level has a list for every possible number of contiguous free elements.
		 * of contiguous free elements.
		 * Since "contiguous" free elements do not span a word, that is one list
		 * for every possible contiguous size in a word.
		 * The list at element n locates an area_ctl with has maximal runs of exactly
		 * n+1 free elements,
		 */

		using free_ctls_t = std::array<list_item, alloc_states_per_word>;
		free_ctls_t _free_ctls;
		free_ctls_t::size_type find_free_ctl_ix(unsigned min_run_length) const;
		free_ctls_t::size_type tier_ix_from_run_length(unsigned run_length) const;
	public:
		unsigned _ct_alloc_probe_success;
		unsigned _ct_alloc_probe_failure;
		unsigned _ct_subdivision;

	public:
		static constexpr auto size() { return alloc_states_per_word; }

		/* return smallest tier index which has an area with
		 *   run length >= min_run_length
		 * If no such tier, return sub_states_per_word.
		 */
		auto find_free_ctl(unsigned min_run_length_) const -> const list_item *
		{
			return & _free_ctls[find_free_ctl_ix(min_run_length_)];
		}

		auto find_free_ctl(unsigned min_run_length_) -> list_item *
		{
			return & _free_ctls[find_free_ctl_ix(min_run_length_)];
		}

		const auto *tier_from_run_length(unsigned run_length_) const
		{
			return & _free_ctls[tier_ix_from_run_length(run_length_)];
		}

		auto *tier_from_run_length(unsigned run_length_)
		{
			return & _free_ctls[tier_ix_from_run_length(run_length_)];
		}

		const auto *tier_end() const
		{
			return &_free_ctls[alloc_states_per_word];
		}

		/* Find the longest run length less than mac_run_length.
		 * Used for tracing, when find_free_ctl_ix has failed to
		 * find a long enough run
		 */
		free_ctls_t::size_type find_mac_ctl_ix(unsigned mac_run_length) const;

		level_hints()
			: _free_ctls()
			, _ct_alloc_probe_success(0)
			, _ct_alloc_probe_failure(0)
			, _ct_subdivision(0)
		{}

		using level_ix_t = std::uint8_t;
		void print(
			std::ostream &o_
			, level_ix_t level_
			, std::ios_base::fmtflags // size_format_
		) const;
	};

	struct area_ctl;

	/*
	 * Location of non-persisted items for a single "region" managed by the crash-consistent allocator.
	 * Persisted items are keps in an area_ctl, which is in persistent memory.
	 */
	struct area_top
	{
	private:
		using level_ix_t = std::uint8_t;
		using byte_span = common::byte_span;
		using persist_type = gsl::not_null<ccpm::persister *>;
		area_ctl *_ctl;
		std::size_t _bytes_free;
		bool _all_restored;
		unsigned _trace_level;
		std::ostream *_o;
		using level_hints_vec = std::vector<level_hints>;
		level_hints_vec _level;
		unsigned _ct_allocation;
		byte_span _region; /* for get_region only */

		area_top(
			// persist_type persist
			area_ctl *ctl
			, unsigned trace_level
			, byte_span iov
			, std::ostream &o
		);
		area_top(const area_top &) = delete;
		area_top &operator=(const area_top &) = delete;

		void allocate_strategy_1(
			persist_type persist_
			, void * & ptr_
			, std::size_t bytes
			, std::size_t alignment
			, level_hints_vec::iterator level_
			, unsigned run_length
		);

		bool allocate_recovery_1();
		bool allocate_recovery_2(persist_type persist, level_hints_vec::iterator level);
		bool trace_coarse() const { return 0 < _trace_level; }
		bool trace_fine() const { return 1 < _trace_level; }

	public:
		/* Initial area_ctl */
		explicit area_top(
			persist_type persist
			, byte_span iov
			, unsigned trace_level
			, std::ostream &o
		);
		/* Restored area_ctl */
		explicit area_top(
			persist_type persist
			, byte_span iov
			, const ownership_callback_t &resolver
			, unsigned trace_level
			, std::ostream &o
		);
		~area_top();

		bool includes(const void *addr) const;

		/* Free byte count. Required by users */
		std::size_t bytes_free() const;

		byte_span get_region() const { return _region; }

		void allocate(
			persist_type persist
			, void * & ptr, std::size_t bytes
			, std::size_t alignment
		);

		void deallocate(
			persist_type persist
			, void * & ptr, std::size_t bytes
		);

		void print(
			std::ostream &o
			, std::ios_base::fmtflags size_format
		) const;

		void print_ctls(
			std::ostream *o_
			, std::ios_base::fmtflags format_
		) const;

    level_ix_t height() const { return level_ix_t(_level.size()); }

    void set_root(const byte_span & iov, persist_type persist);
    byte_span get_root() const;

		/*
		 * called by area_ctl to add area_ctl a, at level level_ix, with a longest
		 * free run (consecutive free elements) of free_run, to _level, which is the
		 * non-persistent catalog of area_ctl items.
		 */
		void remove_from_chain(area_ctl *a, level_ix_t level_ix, unsigned longest_run);
		void restore_to_chain(area_ctl *a, level_ix_t level_ix, unsigned run_length);

		bool contains(const void *p) const;

		bool is_in_chain(
			const area_ctl *a
			, level_ix_t level_ix
			, unsigned run_length
		) const;
	};
}

#endif

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
#ifndef __CROSSBOW_BLAS_H_
#define __CROSSBOW_BLAS_H_

#include <cblas.h>

#include <jni.h>

void blas_init (int size, int bufferSize);

void blas_free ();

enum CBLAS_TRANSPOSE getCblasTrans (const char *);

void *getObjectBufferAddress (JNIEnv *, jobject, int);

#endif /* __CROSSBOW_BLAS_H_ */

max_stars_repo_name: lsds/Crossbow

max_issues_repo_name: lsds/Crossbow

max_forks_repo_name: lsds/Crossbow


text:
#ifndef __INS_C_H__
#define __INS_C_H__
/********************************* TRICK HEADER *******************************
PURPOSE:
      (Describe the INS Module On Board, Error equations based on Zipfel, Figure 10.27, space stabilized INS with GPS updates)
LIBRARY DEPENDENCY:
      ((../src/gnc_var.c)
        (../src/Ins_c.c)
        (../../cad/src/global_constants.c)
        (../src/dm_delta_ut.c))
*******************************************************************************/

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include <stdio.h>
#include "cad_utility_c.h"
#include "math_utility_c.h"
#include "time_utility_c.h"

extern const double __DM_sec2r;
extern const double __DM_arcsec2r;
extern const double __WEII3;
extern const double __SMAJOR_AXIS;
extern const double __GM;
extern const double __WEII3;
extern const double ___PI;
extern const double __EPS;
extern const double __DEG;
extern const double __RAD;
/** INS Variables **/
/* Matrix */
extern gsl_matrix *WEII;
extern gsl_matrix *TBIC;
extern gsl_matrix *TDCI;
extern gsl_matrix *TEIC;
extern gsl_matrix *TBDC;
extern gsl_matrix *TBICI;
extern gsl_matrix *TLI;

/* Vector */
extern gsl_vector *EVBI;
extern gsl_vector *EVBID;
extern gsl_vector *ESBI;
extern gsl_vector *ESBID;
extern gsl_vector *RICI;
extern gsl_vector *RICID;
extern gsl_vector *TBIC_Q;
extern gsl_vector *TBIDC_Q;
extern gsl_vector *SBIIC;
extern gsl_vector *VBIIC;
extern gsl_vector *SBEEC;
extern gsl_vector *VBEEC;
extern gsl_vector *WBICI;
extern gsl_vector *EGRAVI;
extern gsl_vector *VBECD;
extern gsl_vector *INS_I_ATT_ERR;
extern gsl_vector *TESTV;
extern gsl_vector *TMP_old;
extern gsl_vector *VBIIC_old;
extern gsl_vector *POS_ERR;
extern gsl_vector *GRAVGI;
extern gsl_vector *TBDCQ;
extern gsl_vector *VBIIC_old_old;
extern gsl_vector *PHI_C;
extern gsl_vector *DELTA_VEL_C;
extern gsl_vector *PHI_LOW_C;
extern gsl_vector *PHI_HIGH_C;

/* Double */
extern double dbic;
extern double dvbec;
extern double alphacx;
extern double betacx;
extern double thtvdcx;
extern double psivdcx;
extern double alppcx;
extern double phipcx;
extern double loncx;
extern double latcx;
extern double altc;
extern double phibdcx;
extern double thtbdcx;
extern double psibdcx;
extern double ins_pos_err;
extern double ins_vel_err;
extern double ins_tilt_err;
extern double ins_pose_err;
extern double ins_vele_err;
extern double ins_phi_err;
extern double ins_tht_err;
extern double ins_psi_err;

/* Unsigned int */
extern unsigned int gpsupdate;
extern unsigned int liftoff;
extern unsigned int ideal;

extern GPS_TIME gpstime;
extern UTC_TIME utctime;
/*******************/
#ifdef __cplusplus
extern "C" {
#endif
    int load_location(double lonx, double latx, double alt);
    int load_angle(double yaw, double roll, double pitch, GPS_TIME gps_time);
    int load_geodetic_velocity(double alpha0x, double beta0x, double dvbe);
    int calculate_INS_derived_TEI(GPS_TIME gps, gsl_matrix *TEIC);
    int AccelHarmonic(const gsl_vector *SBII, double CS[21][21], int n_max, int m_max, const gsl_matrix *TEIC, gsl_vector *acc_out);
    int DCM_2_Euler_angle(const gsl_matrix *TBD, double *phibdc, double *thtbdc, double *psibdc);
    int calculate_INS_derived_phip(gsl_vector *VBECB, double *phipc);
    int calculate_INS_derived_thtvd(gsl_vector *VBECD, double *thtvd);
    int calculate_INS_derived_psivd(gsl_vector *VBECD, double *psivd);
    int calculate_INS_derived_alpp(gsl_vector *VBECB, double *alpp);
    int calculate_INS_derived_beta(gsl_vector *VBECB, double *beta);
    int calculate_INS_derived_alpha(gsl_vector *VBECB, double *alpha);
    int build_VBEB(double _alpha0x, double _beta0x, double _dvbe, gsl_vector *VBEB);
    int INS_update(const double int_step, double *dvbec, unsigned int liftoff, double *alphacx, double *betacx
                    , double *alppcx, double *phipcx, double *loncx, double *latcx, double *altc, double *psivdcx, double *thtvdcx
                    , double *phibdc, double *thtbdc, double *psibdc
                    , gsl_vector *PHI, gsl_vector *DELTA_VEL
                    , gsl_vector *PHI_HIGH, gsl_vector *PHI_LOW, GPS_TIME gps, gsl_matrix *TEIC
                    , gsl_vector *SBIIC, gsl_vector *VBIIC, gsl_vector *VBIIC_old, gsl_vector *GRAVGI, gsl_matrix *TBIC, gsl_vector *SBEEC
                    , gsl_vector *VBEEC, gsl_matrix *WEII, gsl_matrix *TLI, gsl_matrix *TDCI
                    , gsl_matrix *TBICI, gsl_matrix *TBDC, gsl_vector *TBDCQ);
    int INS_init(GPS_TIME gps_time);
    int INS_alloc();
#ifdef __cplusplus
}
#endif
#endif


text:
/*
 * Assignment 4, Question 5 matvec implementation in C using MPI.
 *
 * Copyright (C) 2015, Jonathan Gillett (100437638)
 * All rights reserved.
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include <mpi.h>

// Minimum of two values
#define min(a, b) (((a) < (b)) ? (a) : (b))
// Define process 0 as MASTER
#define MASTER 0
// Message tags
#define REQUEST 1
#define REPLY   2

// Master process, divides up the matrix into rows for computation by slaves
void master(int n_proc, int n_rows, int n_cols) {
    // 0. Allocate space for arrays as needed to compute y = A*x
    double *x = (double *) malloc(n_cols*sizeof(double));
    double *y = (double *) malloc(n_rows*sizeof(double));
    double *z = (double *) malloc(n_rows*sizeof(double));
    double *A = (double *) malloc(n_rows*n_cols*sizeof(double));
    double *sub = (double *) malloc(n_cols*sizeof(double));
    double ans = 0.0;
    int n_sent = 0, row = 0, sender = 0;
    MPI_Status status;

    // 1. Define data A (n_rows x n_cols matrix) and x (n_cols x 1 vector)
    for (int i = 0; i < n_cols; ++i) {
        x[i] = i+1;
    }
    for (int i = 0; i < n_rows; ++i) {
        for (int j = 0; j < n_cols; ++j) {
            A[(i*n_cols)+j] = i+1;
        }
    }

    // 2. Broadcast vector x to all processes
    MPI_Bcast(x, n_cols, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);

    // 3. Send first rows to other processes
    for (int i = 1; i < min(n_rows+1, n_proc); ++i) {
    	memcpy((void *) sub, (void *) &A[n_sent*n_cols], n_cols * sizeof(double));
        MPI_Send(sub, n_cols, MPI_DOUBLE, i, n_sent, MPI_COMM_WORLD);
        n_sent++;
    }

    // Now loop over remaining rows, receiving dot products computed
    for (int i = 0; i < n_rows; ++i) {
        // Receive a computed inner product from one of the processes
        MPI_Recv(&ans, 1, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG,
                 MPI_COMM_WORLD, &status);
        sender = status.MPI_SOURCE;
        row = status.MPI_TAG;
        y[row] = ans;

        // Either send another row back to sender or send final y=Ax
        // along with a tag to signal completion.
        if (n_sent < n_rows) {
        	memcpy((void *) sub, (void *) &A[n_sent*n_cols], n_cols * sizeof(double));
            MPI_Send(sub, n_cols, MPI_DOUBLE, sender, n_sent, MPI_COMM_WORLD);
            n_sent++;
        }
        // Inform sender that there is no more work
        else {
            MPI_Send(y, n_rows, MPI_DOUBLE, sender, n_rows+1, MPI_COMM_WORLD);
        }
    }

    printf("Finished!\n");
    // Calculate the norm to verify ||y-Ax|| is correct
    cblas_dgemv(CblasRowMajor, CblasNoTrans, n_rows, n_cols, 1.0, A, n_cols,
	      		x, 1, 0.0, z, 1);
    for (int i = 0; i < n_rows; ++i) {
    	y[i] -= z[i];
    }
    printf("Check: ||y-Ax|| = %f\n", cblas_dnrm2(n_rows, y, 1));

    free(x);
    free(y);
    free(z);
    free(A);
    free(sub);
}

// Slave process computes inner product and returns result to master process
void slave(int proc_id, int n_rows, int n_cols) {
    // Allocate arrays for receiving and sending rows & dot products
    double *x = (double *) malloc(n_cols*sizeof(double));
    double *my_row = (double *) malloc(n_cols*sizeof(double));
    double ans = 0.0;
    int row = 0;
    MPI_Status status;

    // 1. Receive vector x from master process
    MPI_Bcast(x, n_cols, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);

    // Receive row to calculate from master process
    MPI_Recv(my_row, n_cols, MPI_DOUBLE, MASTER, MPI_ANY_TAG,
             MPI_COMM_WORLD, &status);
    row = status.MPI_TAG;

    while (row < n_rows) {
        // Compute inner product and return result to master process
        ans = cblas_ddot(n_cols, my_row, 1, x, 1);
        MPI_Send(&ans, 1, MPI_DOUBLE, MASTER, row, MPI_COMM_WORLD);
        // Receive new row from master
        MPI_Recv(my_row, n_cols, MPI_DOUBLE, MASTER, MPI_ANY_TAG,
                 MPI_COMM_WORLD, &status);
        row = status.MPI_TAG;
    }
}

int main(int argc, char *argv[]) {
    int    proc_id;       // Process ID number
    int    n_proc;        // Number of processes
    int    n_rows = 1000; // Number of rows
    int    n_cols = 1000; // Number of columns

    // Initialize MPI
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &n_proc);
    MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);

    if (n_proc == 1) {
        puts("Error: requires two or more processes (no work will be done).");
        MPI_Finalize();
        exit(1);
    }

    if (proc_id == MASTER) {
        master(n_proc, n_rows, n_cols);
    }
    else {
        slave(proc_id, n_rows, n_cols);
    }

    MPI_Finalize();
    return 0;
}


text:
#include <stdio.h>
#include <gsl/gsl_combination.h>

void
print_all_combinations ( size_t n, size_t k );

int
main (void)
{
  gsl_combination * c;
  size_t i;

  printf("all subsets of {0,1,2,3} by size (lex. order)\n") ;
  for(i = 0; i <= 4; i++)
    {
      c = gsl_combination_calloc (4, i);
      do
        {
          printf("{");
          gsl_combination_fprintf (stdout, c, " %u");
          printf(" }\n");
        }
      while (gsl_combination_next(c) == GSL_SUCCESS);
      gsl_combination_free(c);
    }
  printf("all subsets of {1,2,3,4} by size (reverse lex. order)\n") ;
  for(i = 0; i <= 4; i++)
    {
      c = gsl_combination_alloc (4, i) ;
      gsl_combination_init_last(c) ;
      do
        {
          printf("{");
          gsl_combination_fprintf (stdout, c, " %u");
          printf(" }\n");
        }
      while (gsl_combination_prev(c) == GSL_SUCCESS);
      gsl_combination_free(c);
    }
  printf("\n");

  print_all_combinations(5, 3);
  print_all_combinations(5, 0);
  print_all_combinations(5, 5);
  print_all_combinations(1, 1);
  print_all_combinations(3, 1);

  return 0;
}

void
print_all_combinations (size_t n, size_t k)
{
  gsl_combination * c = gsl_combination_calloc (n, k);

  printf("combinations %u choose %u\n", n, k);
  do
    {
      gsl_combination_fprintf (stdout, c, " %u");
      printf("\n");
    }
  while (gsl_combination_next(c) == GSL_SUCCESS);
  do
    {
      gsl_combination_fprintf (stdout, c, " %u");
      printf("\n");
    }
  while (gsl_combination_prev(c) == GSL_SUCCESS);
  printf("\n");
  gsl_combination_free(c);
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/*
 * Reinforcement Learning Book / Exercise 4.7
 *
 * $ sudo apt-get install libgsl-dev libgsl2
 * $ gcc car_rental2.c -lgsl -lgslcblas -lm -o car_rental2
 */

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_randist.h>

#define Size 20
#define Iteration 100
#define Rent1_mean 3
#define Rent2_mean 4
#define Return1_mean 3
#define Return2_mean 2

// Global variables
float value[Size + 1][Size + 1];
int policy[Size + 1][Size + 1];

void show_value(){
    printf("[");
    for (int c1 = 0; c1 <= Size; c1++) {
        printf("[");
        for (int c2 = 0; c2 <= Size; c2++) {
            printf("%3.0f", value[c1][c2]);
            if (c2 < Size) {
                printf(", ");
            }
        }
        printf("]");
        if (c1 < Size) {
            printf(",\n ");
        }
    }
    printf("]\n");
}

void show_policy(){
    printf("[");
    for (int c1 = 0; c1 <= Size; c1++) {
        printf("[");
        for (int c2 = 0; c2 <= Size; c2++) {
            printf("%2d", policy[c1][c2]);
            if (c2 < Size) {
                printf(", ");
            }
        }
        printf("]");
        if (c1 < Size) {
            printf(",\n ");
        }
    }
    printf("]\n");
}

float calc_q_val(int c1, int c2, int a){
    float q_val = 0;
    int _c1, _c2;
    float c1_rent_prob, c2_rent_prob, c1_return_prob, c2_return_prob, prob;

    // Moving cars
    c1 -= a;
    c2 += a;
    if (a > 0) {
        q_val -= 2*(a-1);
    } else {
        q_val -= 2*abs(a);
    }
    if (c1 > 10) q_val -= 4;
    if (c2 > 10) q_val -= 4;

    for (int c1_rent = 0; c1_rent <= c1; c1_rent++) {
        for (int c2_rent = 0; c2_rent <= c2; c2_rent++) {
            for (int c1_return = 0; c1_return <= Size-(c1-c1_rent); c1_return++) {
                for (int c2_return = 0; c2_return <= Size-(c2-c2_rent); c2_return++) {
                    _c1 = c1 - c1_rent + c1_return;
                    _c2 = c2 - c2_rent + c2_return;

                    if (c1_rent == c1) {
                        c1_rent_prob = 1;
                        for (int n = 0; n < c1_rent; n++) {
                            c1_rent_prob -= gsl_ran_poisson_pdf(n, Rent1_mean);
                        }
                    } else {
                        c1_rent_prob = gsl_ran_poisson_pdf(c1_rent, Rent1_mean);
                    }

                    if (c2_rent == c2) {
                        c2_rent_prob = 1;
                        for (int n = 0; n < c2_rent; n++) {
                            c2_rent_prob -= gsl_ran_poisson_pdf(n, Rent2_mean);
                        }
                    } else {
                        c2_rent_prob = gsl_ran_poisson_pdf(c2_rent, Rent2_mean);
                    }

                    if (c1_return == Size-(c1-c1_rent)) {
                        c1_return_prob = 1;
                        for (int n = 0; n < c1_return; n++) {
                            c1_return_prob -= gsl_ran_poisson_pdf(n, Return1_mean);
                        }
                    } else {
                        c1_return_prob = gsl_ran_poisson_pdf(c1_return, Return1_mean);
                    }

                    if (c2_return == Size-(c2-c2_rent)) {
                        c2_return_prob = 1;
                        for (int n = 0; n < c2_return; n++) {
                            c2_return_prob -= gsl_ran_poisson_pdf(n, Return2_mean);
                        }
                    } else {
                        c2_return_prob = gsl_ran_poisson_pdf(c2_return, Return2_mean);
                    }

                    prob = c1_rent_prob * c2_rent_prob * c1_return_prob * c2_return_prob;
                    q_val += prob * (10 * (c1_rent + c2_rent) + 0.9 * value[_c1][_c2]);
                }
            }
        }
    }
    return q_val;
}

int min(int a, int b){
    if (a < b) {
        return a;
    } else {
        return b;
    }
}

void run() {
    float q_val;
    float value_update;
    int policy_update;

    for (int c1 = 0; c1 <= Size; c1++) {
        for (int c2 = 0; c2 <= Size; c2++) {
            value[c1][c2] = 0;
            policy[c1][c2] = 0;
        }
    }

    printf("# Iteration 0\n");
    show_policy();
    show_value();

    for (int i = 1; i < Iteration; i++) {
        for (int c1 = 0; c1 <= Size; c1++) {
            for (int c2 = 0; c2 <= Size; c2++) {
                int update = 0;
                for (int a = -min(c2, 5); a <= min(c1, 5); a++) {
                    q_val = calc_q_val(c1, c2, a);
                    if (update == 0 || q_val > value_update) {
                        value_update = q_val;
                        policy_update = a;
                        update = 1;
                    }
                }
                value[c1][c2] = value_update;
                policy[c1][c2] = policy_update;
            }
        }
        printf("# Iteration %d\n", i);
        show_policy();
        show_value();
    }
}

int main() {
    run();
}


text:
#include <glib.h>
#include <gsl/gsl_rng.h>
#include <stdlib.h>
#include "util.h"

hash_t *hash_create () {
  return (char*) g_hash_table_new(g_direct_hash, g_direct_equal);
}

void hash_set (hash_t *table, size_t key, size_t value) {
  g_hash_table_replace((GHashTable*)table, (gpointer)key, (gpointer)value);
}

size_t hash_get (hash_t *table, size_t key, size_t value) {
  size_t orig_key;
  size_t result;
  if (g_hash_table_lookup_extended((GHashTable*)table, (gpointer)key, (gpointer*)&orig_key, (gpointer*)&result)) {
    return result;
  }
  return value;
}

void hash_clear (hash_t *table) {
  return g_hash_table_remove_all((GHashTable*)table);
}

void hash_destroy (hash_t *table) {
  return g_hash_table_destroy((GHashTable*)table);
}

gsl_rng *r = 0;

size_t random_range (size_t lo, size_t hi) {
  if (!r) {
    gsl_rng_env_setup();
    const gsl_rng_type *T = gsl_rng_default;
    r = gsl_rng_alloc (T);
  }
  size_t n = hi - lo;
  return lo + gsl_rng_uniform_int(r, n);//todo this range might not be big enough/type issues may occur here.
}

double random_uniform () {
  if (!r) {
    gsl_rng_env_setup();
    const gsl_rng_type *T = gsl_rng_default;
    r = gsl_rng_alloc (T);
  }
  return gsl_rng_uniform(r);
}

void random_choose (size_t *samples, size_t n, size_t lo, size_t hi) {
  hash_t *table = hash_create();
  for (size_t i = 0; i < n; i++) {
    size_t j = random_range(lo, hi);
    samples[i] = hash_get(table, j, j);
    hash_set(table, j, hash_get(table, lo, lo));
    lo++;
  }
  hash_destroy(table);
}

int size_t_cmp(const void *a, const void *b)
{
    const size_t *ia = (const size_t *)a;
    const size_t *ib = (const size_t *)b;
    return *ia  - *ib;
}

void quicksort(size_t *stuff, size_t n)
{
  if (n < 2) return;

  size_t pivot = stuff[n / 2];

  size_t i, j;
  for (i = 0, j = n - 1; ; i++, j--)
  {
    while (stuff[i] < pivot) i++;
    while (stuff[j] > pivot) j--;

    if (i >= j) break;

    size_t temp = stuff[i];
    stuff[i]     = stuff[j];
    stuff[j]     = temp;
  }

  quicksort(stuff, i);
  quicksort(stuff + i, n - i);
}

void sort (size_t *stuff, size_t n) {
  //quicksort(stuff, n);
  qsort(stuff, n, sizeof(size_t), size_t_cmp);
}

size_t search (const size_t *stuff, size_t lo, size_t hi, size_t key) {
  size_t len = hi - lo;
  size_t half;
  size_t mid;
  while (len > 0) {
    half = len >> 1;
    mid = lo + half;
    if (stuff[mid] >= key) {
      len = half;
    } else {
      lo = mid;
      lo++;
      len = len - half - 1;
    }
  }
  return lo;
}

size_t search_strict (const size_t *stuff, size_t lo, size_t hi, size_t key) {
  size_t len = hi - lo;
  size_t half;
  size_t mid;
  while (len > 0) {
    half = len >> 1;
    mid = lo + half;
    if (stuff[mid] > key) {
      len = half;
    } else {
      lo = mid;
      lo++;
      len = len - half - 1;
    }
  }
  return lo;
}


text:
#pragma once

#include "tino/backends/backends.h"
#include "tino/backends/tensor_wrappers.h"
#include "tino/utils/macros.h"

#ifdef TINO_OPENBLAS_READY
#include <cblas.h>
#endif

namespace tino {
  namespace core {
    class context;
  }
  namespace backends {
    namespace blas {
      enum class layout_t {
#ifdef TINO_OPENBLAS_READY
        RowMajor    = CblasRowMajor,
        ColumnMajor = CblasColMajor
#else
        RowMajor,
        ColumnMajor
#endif
      };
      enum class trans_t {
#ifdef TINO_OPENBLAS_READY
        NoTrans = CblasNoTrans,
        Trans   = CblasTrans
#else
        NoTrans,
        Trans
#endif
      };

      template <typename T>
      struct blasOpts {
        T alpha         = (T)0;
        T beta          = (T)0;
        layout_t layout = layout_t::RowMajor;
        trans_t trans_a = trans_t::NoTrans;
        trans_t trans_b = trans_t::NoTrans;
      };

      template <typename T>
      tensor2d<T>& blas_gemm(core::context& ctx, const blasOpts<T>& blas_opts, tensor2d<T>& A, tensor2d<T>& B, tensor2d<T>& C) {}

      template <>
      tensor2d<float>& blas_gemm(core::context& ctx, const blasOpts<float>& blas_opts, tensor2d<float>& A, tensor2d<float>& B, tensor2d<float>& C) {
#ifdef TINO_OPENBLAS_READY
        const int M   = blas_opts.trans_a == trans_t::NoTrans ? A.template shape<1>() : A.template shape<0>();
        const int N   = blas_opts.trans_b == trans_t::NoTrans ? B.template shape<0>() : B.template shape<1>();
        const int K   = blas_opts.trans_a == trans_t::NoTrans ? A.template shape<0>() : A.template shape<1>();
        const int lda = blas_opts.trans_a == trans_t::NoTrans ? K : M;
        const int ldb = blas_opts.trans_b == trans_t::NoTrans ? N : K;
        const int ldc = N;

        cblas_sgemm(static_cast<CBLAS_LAYOUT>(blas_opts.layout),
                    static_cast<CBLAS_TRANSPOSE>(blas_opts.trans_a),
                    static_cast<CBLAS_TRANSPOSE>(blas_opts.trans_b),
                    M,
                    N,
                    K,
                    blas_opts.alpha,
                    A.data(),
                    lda,
                    B.data(),
                    ldb,
                    blas_opts.beta,
                    C.data(),
                    ldc);
#else
        std::cerr << "invalid backend: OpenBLAS" << std::endl;
        std::exit(1);
#endif

        TINO_MAYBE_UNUSED(ctx);
        TINO_MAYBE_UNUSED(blas_opts);
        TINO_MAYBE_UNUSED(A);
        TINO_MAYBE_UNUSED(B);
        TINO_MAYBE_UNUSED(C);

        return C;
      }

    } // namespace blas
  }   // namespace backends
} // namespace tino

max_stars_repo_name: Senseirn/tino

max_issues_repo_name: Senseirn/tino

max_forks_repo_name: Senseirn/tino


text:
/* multifit_nlinear/fdf.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 * Copyright (C) 2015, 2016 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

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

gsl_multifit_nlinear_workspace *
gsl_multifit_nlinear_alloc (const gsl_multifit_nlinear_type * T,
                            const gsl_multifit_nlinear_parameters * params,
                            const size_t n, const size_t p)
{
  gsl_multifit_nlinear_workspace * w;

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

  w = calloc (1, sizeof (gsl_multifit_nlinear_workspace));
  if (w == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for multifit workspace",
                     GSL_ENOMEM, 0);
    }

  w->x = gsl_vector_calloc (p);
  if (w->x == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
    }

  w->f = gsl_vector_calloc (n);
  if (w->f == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0);
    }

  w->dx = gsl_vector_calloc (p);
  if (w->dx == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
    }

  w->g = gsl_vector_alloc (p);
  if (w->g == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for g", GSL_ENOMEM, 0);
    }

  w->J = gsl_matrix_alloc(n, p);
  if (w->J == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for Jacobian", GSL_ENOMEM, 0);
    }

  w->sqrt_wts_work = gsl_vector_calloc (n);
  if (w->sqrt_wts_work == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for weights", GSL_ENOMEM, 0);
    }

  w->state = (T->alloc)(params, n, p);
  if (w->state == 0)
    {
      gsl_multifit_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for multifit state", GSL_ENOMEM, 0);
    }

  w->type = T;
  w->fdf = NULL;
  w->niter = 0;
  w->params = *params;

  return w;
}

void
gsl_multifit_nlinear_free (gsl_multifit_nlinear_workspace * w)
{
  RETURN_IF_NULL (w);

  if (w->state)
    (w->type->free) (w->state);

  if (w->dx)
    gsl_vector_free (w->dx);

  if (w->x)
    gsl_vector_free (w->x);

  if (w->f)
    gsl_vector_free (w->f);

  if (w->sqrt_wts_work)
    gsl_vector_free (w->sqrt_wts_work);

  if (w->g)
    gsl_vector_free (w->g);

  if (w->J)
    gsl_matrix_free (w->J);

  free (w);
}

gsl_multifit_nlinear_parameters
gsl_multifit_nlinear_default_parameters(void)
{
  gsl_multifit_nlinear_parameters params;

  params.trs = gsl_multifit_nlinear_trs_lm;
  params.scale = gsl_multifit_nlinear_scale_more;
  params.solver = gsl_multifit_nlinear_solver_qr;
  params.fdtype = GSL_MULTIFIT_NLINEAR_FWDIFF;
  params.factor_up = 3.0;
  params.factor_down = 2.0;
  params.avmax = 0.75;
  params.h_df = GSL_SQRT_DBL_EPSILON;
  params.h_fvv = 0.02;

  return params;
}

int
gsl_multifit_nlinear_init (const gsl_vector * x,
                           gsl_multifit_nlinear_fdf * fdf,
                           gsl_multifit_nlinear_workspace * w)
{
  return gsl_multifit_nlinear_winit(x, NULL, fdf, w);
}

int
gsl_multifit_nlinear_winit (const gsl_vector * x,
                            const gsl_vector * wts,
                            gsl_multifit_nlinear_fdf * fdf,
                            gsl_multifit_nlinear_workspace * w)
{
  const size_t n = w->f->size;

  if (n != fdf->n)
    {
      GSL_ERROR ("function size does not match workspace", GSL_EBADLEN);
    }
  else if (w->x->size != x->size)
    {
      GSL_ERROR ("vector length does not match workspace", GSL_EBADLEN);
    }
  else if (wts != NULL && n != wts->size)
    {
      GSL_ERROR ("weight vector length does not match workspace", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* initialize counters for function and Jacobian evaluations */
      fdf->nevalf = 0;
      fdf->nevaldf = 0;
      fdf->nevalfvv = 0;

      w->fdf = fdf;
      gsl_vector_memcpy(w->x, x);
      w->niter = 0;

      if (wts)
        {
          w->sqrt_wts = w->sqrt_wts_work;

          for (i = 0; i < n; ++i)
            {
              double wi = gsl_vector_get(wts, i);
              gsl_vector_set(w->sqrt_wts, i, sqrt(wi));
            }
        }
      else
        {
          w->sqrt_wts = NULL;
        }

      return (w->type->init) (w->state, w->sqrt_wts, w->fdf,
                              w->x, w->f, w->J, w->g);
    }
}

int
gsl_multifit_nlinear_iterate (gsl_multifit_nlinear_workspace * w)
{
  int status =
    (w->type->iterate) (w->state, w->sqrt_wts, w->fdf,
                        w->x, w->f, w->J, w->g, w->dx);

  w->niter++;

  return status;
}

double
gsl_multifit_nlinear_avratio (const gsl_multifit_nlinear_workspace * w)
{
  return (w->type->avratio) (w->state);
}

/*
gsl_multifit_nlinear_driver()
  Iterate the nonlinear least squares solver until completion

Inputs: maxiter  - maximum iterations to allow
        xtol     - tolerance in step x
        gtol     - tolerance in gradient
        ftol     - tolerance in ||f||
        callback - callback function to call each iteration
        callback_params - parameters to pass to callback function
        info     - (output) info flag on why iteration terminated
                   1 = stopped due to small step size ||dx|
                   2 = stopped due to small gradient
                   3 = stopped due to small change in f
                   GSL_ETOLX = ||dx|| has converged to within machine
                               precision (and xtol is too small)
                   GSL_ETOLG = ||g||_inf is smaller than machine
                               precision (gtol is too small)
                   GSL_ETOLF = change in ||f|| is smaller than machine
                               precision (ftol is too small)
        w        - workspace

Return:
GSL_SUCCESS if converged
GSL_MAXITER if maxiter exceeded without converging
GSL_ENOPROG if no accepted step found on first iteration
*/

int
gsl_multifit_nlinear_driver (const size_t maxiter,
                             const double xtol,
                             const double gtol,
                             const double ftol,
                             void (*callback)(const size_t iter, void *params,
                                              const gsl_multifit_nlinear_workspace *w),
                             void *callback_params,
                             int *info,
                             gsl_multifit_nlinear_workspace * w)
{
  int status;
  size_t iter = 0;

  /* call user callback function prior to any iterations
   * with initial system state */
  if (callback)
    callback(iter, callback_params, w);

  do
    {
      status = gsl_multifit_nlinear_iterate (w);

      /*
       * If the solver reports no progress on the first iteration,
       * then it didn't find a single step to reduce the
       * cost function and more iterations won't help so return.
       *
       * If we get a no progress flag on subsequent iterations,
       * it means we did find a good step in a previous iteration,
       * so continue iterating since the solver has now reset
       * mu to its initial value.
       */
      if (status == GSL_ENOPROG && iter == 0)
        {
          *info = status;
          return GSL_EMAXITER;
        }

      ++iter;

      if (callback)
        callback(iter, callback_params, w);

      /* test for convergence */
      status = gsl_multifit_nlinear_test(xtol, gtol, ftol, info, w);
    }
  while (status == GSL_CONTINUE && iter < maxiter);

  /*
   * the following error codes mean that the solution has converged
   * to within machine precision, so record the error code in info
   * and return success
   */
  if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
    {
      *info = status;
      status = GSL_SUCCESS;
    }

  /* check if max iterations reached */
  if (iter >= maxiter && status != GSL_SUCCESS)
    status = GSL_EMAXITER;

  return status;
} /* gsl_multifit_nlinear_driver() */

gsl_matrix *
gsl_multifit_nlinear_jac (const gsl_multifit_nlinear_workspace * w)
{
  return w->J;
}

const char *
gsl_multifit_nlinear_name (const gsl_multifit_nlinear_workspace * w)
{
  return w->type->name;
}

gsl_vector *
gsl_multifit_nlinear_position (const gsl_multifit_nlinear_workspace * w)
{
  return w->x;
}

gsl_vector *
gsl_multifit_nlinear_residual (const gsl_multifit_nlinear_workspace * w)
{
  return w->f;
}

size_t
gsl_multifit_nlinear_niter (const gsl_multifit_nlinear_workspace * w)
{
  return w->niter;
}

int
gsl_multifit_nlinear_rcond (double *rcond, const gsl_multifit_nlinear_workspace * w)
{
  int status = (w->type->rcond) (rcond, w->state);
  return status;
}

const char *
gsl_multifit_nlinear_trs_name (const gsl_multifit_nlinear_workspace * w)
{
  return w->params.trs->name;
}

/*
gsl_multifit_nlinear_eval_f()
  Compute residual vector y with user callback function, and apply
weighting transform if given:

y~ = sqrt(W) y

Inputs: fdf  - callback function
        x    - model parameters
        swts - weight matrix sqrt(W) = sqrt(diag(w1,w2,...,wn))
               set to NULL for unweighted fit
        y    - (output) (weighted) residual vector
               y_i = sqrt(w_i) f_i where f_i is unweighted residual
*/

int
gsl_multifit_nlinear_eval_f(gsl_multifit_nlinear_fdf *fdf,
                            const gsl_vector *x,
                            const gsl_vector *swts,
                            gsl_vector *y)
{
  int s = ((*((fdf)->f)) (x, fdf->params, y));

  ++(fdf->nevalf);

  /* y <- sqrt(W) y */
  if (swts)
    gsl_vector_mul(y, swts);

  return s;
}

/*
gsl_multifit_nlinear_eval_df()
  Compute Jacobian matrix J with user callback function, and apply
weighting transform if given:

J~ = sqrt(W) J

Inputs: x      - model parameters
        f      - residual vector f(x)
        swts   - weight matrix W = diag(w1,w2,...,wn)
                 set to NULL for unweighted fit
        h      - finite difference step size
        fdtype - finite difference method
        fdf    - callback function
        df     - (output) (weighted) Jacobian matrix
                 df = sqrt(W) df where df is unweighted Jacobian
        work   - workspace for finite difference, size n
*/

int
gsl_multifit_nlinear_eval_df(const gsl_vector *x,
                             const gsl_vector *f,
                             const gsl_vector *swts,
                             const double h,
                             const gsl_multifit_nlinear_fdtype fdtype,
                             gsl_multifit_nlinear_fdf *fdf,
                             gsl_matrix *df,
                             gsl_vector *work)
{
  int status;

  if (fdf->df)
    {
      /* call user-supplied function */
      status = ((*((fdf)->df)) (x, fdf->params, df));
      ++(fdf->nevaldf);

      /* J <- sqrt(W) J */
      if (swts)
        {
          const size_t n = swts->size;
          size_t i;

          for (i = 0; i < n; ++i)
            {
              double swi = gsl_vector_get(swts, i);
              gsl_vector_view v = gsl_matrix_row(df, i);

              gsl_vector_scale(&v.vector, swi);
            }
        }
    }
  else
    {
      /* use finite difference Jacobian approximation */
      status = gsl_multifit_nlinear_df(h, fdtype, x, swts, fdf, f, df, work);
    }

  return status;
}

/*
gsl_multifit_nlinear_eval_fvv()
  Compute second direction derivative vector yvv with user
callback function, and apply weighting transform if given:

yvv~ = sqrt(W) yvv

Inputs: h    - step size for finite difference, if needed
        x    - model parameters, size p
        v    - unscaled geodesic velocity vector, size p
        f    - residual vector f(x), size n
        J    - Jacobian matrix J(x), n-by-p
        swts - weight matrix sqrt(W) = sqrt(diag(w1,w2,...,wn))
               set to NULL for unweighted fit
        fdf  - callback function
        yvv  - (output) (weighted) second directional derivative vector
               yvv_i = sqrt(w_i) fvv_i where f_i is unweighted
        work - workspace, size p
*/

int
gsl_multifit_nlinear_eval_fvv(const double h,
                              const gsl_vector *x,
                              const gsl_vector *v,
                              const gsl_vector *f,
                              const gsl_matrix *J,
                              const gsl_vector *swts,
                              gsl_multifit_nlinear_fdf *fdf,
                              gsl_vector *yvv, gsl_vector *work)
{
  int status;

  if (fdf->fvv != NULL)
    {
      /* call user-supplied function */
      status = ((*((fdf)->fvv)) (x, v, fdf->params, yvv));
      ++(fdf->nevalfvv);

      /* yvv <- sqrt(W) yvv */
      if (swts)
        gsl_vector_mul(yvv, swts);
    }
  else
    {
      /* use finite difference approximation */
      status = gsl_multifit_nlinear_fdfvv(h, x, v, f, J,
                                          swts, fdf, yvv, work);
    }

  return status;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl/span>
#include <array>
#include "halley/core/api/audio_api.h"

namespace Halley
{
	class AudioSource {
	public:
		virtual ~AudioSource() {}

		virtual uint8_t getNumberOfChannels() const = 0;
		virtual size_t getSamplesLeft() const = 0;
		virtual bool isReady() const { return true; }
		virtual bool getAudioData(size_t numSamples, AudioMultiChannelSamples dst) = 0;
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <php.h>
#include <cblas.h>
#include <lapacke.h>

#include "php_ext.h"
#include "kernel/operators.h"

void tensor_matmul(zval * return_value, zval * a, zval * b)
{
    unsigned int i, j;
    Bucket * row;
    zval rowC, c;

    zend_zephir_globals_def * zephir_globals = ZEPHIR_VGLOBAL;

    openblas_set_num_threads(zephir_globals->num_threads);

    zend_array * aHat = Z_ARR_P(a);
    zend_array * bHat = Z_ARR_P(b);

    Bucket * ba = aHat->arData;
    Bucket * bb = bHat->arData;

    unsigned int m = zend_array_count(aHat);
    unsigned int p = zend_array_count(bHat);
    unsigned int n = zend_array_count(Z_ARR(bb[0].val));

    double * va = emalloc(m * p * sizeof(double));
    double * vb = emalloc(n * p * sizeof(double));
    double * vc = emalloc(m * n * sizeof(double));

    for (i = 0; i < m; ++i) {
        row = Z_ARR(ba[i].val)->arData;

        for (j = 0; j < p; ++j) {
            va[i * p + j] = zephir_get_doubleval(&row[j].val);
        }
    }

    for (i = 0; i < p; ++i) {
        row = Z_ARR(bb[i].val)->arData;

        for (j = 0; j < n; ++j) {
            vb[i * n + j] = zephir_get_doubleval(&row[j].val);
        }
    }

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, p, 1.0, va, p, vb, n, 0.0, vc, n);

    array_init_size(&c, m);

    for (i = 0; i < m; ++i) {
        array_init_size(&rowC, n);

        for (j = 0; j < n; ++j) {
            add_next_index_double(&rowC, vc[i * n + j]);
        }

        add_next_index_zval(&c, &rowC);
    }

    RETVAL_ARR(Z_ARR(c));

    efree(va);
    efree(vb);
    efree(vc);
}

void tensor_dot(zval * return_value, zval * a, zval * b)
{
    unsigned int i;

    zend_array * aHat = Z_ARR_P(a);
    zend_array * bHat = Z_ARR_P(b);

    Bucket * ba = aHat->arData;
    Bucket * bb = bHat->arData;

    unsigned int n = zend_array_count(aHat);

    double sigma = 0.0;

    for (i = 0; i < n; ++i) {
        sigma += zephir_get_doubleval(&ba[i].val) * zephir_get_doubleval(&bb[i].val);
    }

    RETVAL_DOUBLE(sigma);
}

void tensor_inverse(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval rowB, b;

    zend_zephir_globals_def * zephir_globals = ZEPHIR_VGLOBAL;

    openblas_set_num_threads(zephir_globals->num_threads);

    zend_array * aHat = Z_ARR_P(a);

    Bucket * ba = aHat->arData;

    unsigned int n = zend_array_count(aHat);

    double * va = emalloc(n * n * sizeof(double));

    for (i = 0; i < n; ++i) {
        row = Z_ARR(ba[i].val)->arData;

        for (j = 0; j < n; ++j) {
            va[i * n + j] = zephir_get_doubleval(&row[j].val);
        }
    }

    lapack_int status;

    int pivots[n];

    status = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, va, n, pivots);

    if (status != 0) {
        RETURN_NULL();
    }

    status = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, va, n, pivots);

    if (status != 0) {
        RETURN_NULL();
    }

    array_init_size(&b, n);

    for (i = 0; i < n; ++i) {
        array_init_size(&rowB, n);

        for (j = 0; j < n; ++j) {
            add_next_index_double(&rowB, va[i * n + j]);
        }

        add_next_index_zval(&b, &rowB);
    }

    RETVAL_ARR(Z_ARR(b));

    efree(va);
}

max_stars_repo_name: remicollet/Tensor

max_issues_repo_name: remicollet/Tensor

max_forks_repo_name: remicollet/Tensor


text:
#pragma once
/*
 * (C) Copyright 2020-2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \defgroup ioda_cxx_layout Groups and Data Layout
 * \brief Public API for ioda::Group, ioda::ObsGroup, and data layout policies.
 * \ingroup ioda_cxx_api
 *
 * @{
 * \file Group.h
 * \brief Interfaces for ioda::Group and related classes.
 */

#include <gsl/gsl-lite.hpp>
#include <map>
#include <memory>
#include <string>
#include <vector>

#include "ioda/Attributes/Has_Attributes.h"
#include "ioda/Types/Has_Types.h"
#include "ioda/Variables/FillPolicy.h"
#include "ioda/Variables/Has_Variables.h"
#include "ioda/defs.h"

namespace ioda {
class Group;

namespace Engines {
struct Capabilities;
}  // namespace Engines

namespace detail {
class Group_Base;
class Group_Backend;

/// \brief Hidden base class to prevent constructor confusion.
/// \ingroup ioda_cxx_layout
/// \note enable_shared_from_this is a hint to pybind11.
class IODA_DL Group_Base {
  std::shared_ptr<Group_Backend> backend_;

protected:
  Group_Base(std::shared_ptr<Group_Backend>);

public:
  virtual ~Group_Base();

  std::shared_ptr<Group_Backend> getBackend() const { return backend_; }

  /// Get capabilities of the Engine backing this Group
  virtual ::ioda::Engines::Capabilities getCapabilities() const;

  /// \brief Get the fill value policy used for Variables within this Group
  /// \details The backend has to be consulted for this operation. Storage of this policy is
  ///   backend-dependent.
  virtual FillValuePolicy getFillValuePolicy() const;

  /// \brief List all one-level child groups in this group
  /// \details This function exists to provide the same calling semantics as
  ///   vars.list() and atts.list(). It is useful for human exploration of the contents
  ///   of a Group.
  /// \returns A vector of strings listing the child groups. The strings are unordered.
  /// \see listObjects if you need to enumerate both Groups and Variables, or
  ///   if you want to do a recursive search. This function is more sophisticated, and
  ///   depending on the backend it may be faster than recursive calls to
  ///   list() and vars.list().
  std::vector<std::string> list() const;
  /// Same as list(). Uniform semantics with atts() and vars().
  inline std::vector<std::string> groups() const { return list(); }

  /// \brief List all objects (groups + variables) within this group
  /// \param recurse indicates whether the search should be one-level or recursive. If
  ///   multiple possible paths exist for an object, only one is actually returned.
  /// \param filter allows you to search for only a certain type of object, such as a
  ///   Group or Variable.
  /// \returns a map of ObjectTypes. Each mapping contains a vector of strings indicating
  ///   object names. If a filter is provided, then this map only contains the filtered
  ///   object type. Otherwise, the map always has vectors for each element of ObjectType
  ///   (Group, Variable, etc.), although these vectors may be empty if no object of a
  ///   certain type is found.
  /// \todo This function should list all *distinct* objects. In the future, this will mean
  ///   that 1) hard-linked duplicate objects will only be listed once, 2) soft
  ///   links pointing to the same object will only be traversed once, and
  ///   3) multiple external links to the same object will only be traversed once.
  ///   If you would want to list *indistinct* objects, then there will be a
  ///   listLinks function.
  /// \todo Once links are implemented, add an option to auto-resolve
  ///   soft and external links.
  virtual std::map<ObjectType, std::vector<std::string>> listObjects(ObjectType filter
                                                                     = ObjectType::Ignored,
                                                                     bool recurse = false) const;

  template <ObjectType objectClass>
  std::vector<std::string> listObjects(bool recurse = false) const {
    return listObjects(objectClass, recurse)[objectClass];
  }

  /// Does a group exist at the specified path?
  /// \param name is the group name.
  /// \returns true if a group does exist.
  /// \return false if a group does not exist.
  virtual bool exists(const std::string& name) const;

  /// \brief Create a group
  /// \param name is the group name.
  /// \returns an invalid handle on failure. (i.e. return.isGroup() == false).
  /// \returns a scoped handle to the group on success.
  virtual Group create(const std::string& name);

  /// \brief Open a group
  /// \returns an invalid handle if an error occurred. (i.e. return.isGroup() == false).
  /// \returns a scoped handle to the group upon success.
  /// \note It is possible to have multiple handles opened for the group
  /// simultaneously.
  /// \param name is the name of the child group to open.
  virtual Group open(const std::string& name) const;

  /// Use this to access the metadata for the group / ObsSpace.
  Has_Attributes atts;

  /// Use this to access named data types.
  Has_Types types;

  /// Use this to access variables
  Has_Variables vars;
};

class IODA_DL Group_Backend : public Group_Base {
public:
  virtual ~Group_Backend();
  /// Default fill value policy is NETCDF4. Overridable on a per-backend basis.
  FillValuePolicy getFillValuePolicy() const override;

  /// @}

protected:
  Group_Backend();
};

}  // namespace detail

/** \brief Groups are a new implementation of ObsSpaces.
 * \ingroup ioda_cxx_layout
 *
 * \see \ref Groups for a examples.
 *
 * A group can be thought of as a folder that contains Variables and Metadata.
 * A group can also contain child groups, allowing for our ObsSpaces to exist in a
 * nested tree-like structure, which removes the need for having an ObsSpaceContainer.
 *
 * Groups are implemented in several backends, such as the in-memory IODA store, the HDF5 disk
 * backend, the HDF5 in-memory backend, the ATLAS backend, et cetera. The root Group is mounted
 * using one of these backends (probably as a File object, which is a special type of Group),
 * and additional backends may be mounted into the tree structure.
 *
 * \throws std::logic_error if the backend handle points to something other than a group.
 * \throws std::logic_error if the backend handle is invalid.
 * \throws std::logic_error if any error has occurred.
 * \author Ryan Honeyager (honeyage@ucar.edu)
 **/
class IODA_DL Group : public detail::Group_Base {
public:
  Group();
  Group(std::shared_ptr<detail::Group_Backend>);
  virtual ~Group();
};

}  // namespace ioda

/// @}

max_stars_repo_name: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


text:
#pragma once

#include <gsl\gsl>
#include <winrt\base.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "DynamicLODMaterial.h"
#include "MatrixHelper.h"

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

		int MaxTessellationFactor() const;
		void SetMaxTessellationFactor(int maxTessellationFactor);

		DirectX::XMFLOAT2 TessellationDistances() const;
		void SetTessellationDistances(DirectX::XMFLOAT2 tessellationDistances);
		void SetTessellationDistances(float minTessellationDistance, float maxTessellationDistance);

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

	private:
		DynamicLODMaterial mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::size_t mIndexCount{ 0 };
		bool mUpdateMaterial{ true };
	};
}


text:
/* filter/test.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_filter.h>

/* compare two vectors */
static void
compare_vectors(const double tol, const gsl_vector * v, const gsl_vector * expected,
                const char * desc)
{
  const size_t n = v->size;
  size_t i;

  for (i = 0; i < n; ++i)
    {
      double vi = gsl_vector_get(v, i);
      double ui = gsl_vector_get(expected, i);

      gsl_test_rel(vi, ui, tol, "%s i=%zu", desc, i);
    }
}

/* generate random vector with elements in [-1,1] */
static void
random_vector(gsl_vector * v, gsl_rng * r)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      double vi = 2.0 * gsl_rng_uniform(r) - 1.0; /* in [-1,1] */
      gsl_vector_set(v, i, vi);
    }
}

#include "test_impulse.c"
#include "test_gaussian.c"
#include "test_median.c"
#include "test_rmedian.c"

int
main()
{
  gsl_rng * r = gsl_rng_alloc(gsl_rng_default);

  test_gaussian(r);

  test_impulse(r);
  test_median(r);
  test_rmedian(r);

  gsl_rng_free(r);

  exit (gsl_test_summary());
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/** @file */
#ifndef __CCL_NEUTRINOS_H_INCLUDED__
#define __CCL_NEUTRINOS_H_INCLUDED__

#include <gsl/gsl_spline.h>
#include <gsl/gsl_const_mksa.h>

// maximum number of species
#define CCL_MAX_NU_SPECIES 3
// limits for the precomputed spline of phase
// space diagram in MNU/T
#define CCL_NU_MNUT_MIN 1e-4
#define CCL_NU_MNUT_MAX 500
// and number of points
#define CCL_NU_MNUT_N 1000

// The combination of constants required in Omeganuh2
#define NU_CONST ( \
  8. * pow(M_PI,5) *pow((ccl_constants.KBOLTZ/ ccl_constants.HPLANCK),3)* \
  ccl_constants.KBOLTZ/(15. *pow( ccl_constants.CLIGHT,3))* \
  (8. * M_PI * ccl_constants.GNEWT) / \
  (3. * 100.*100.*1000.*1000. /ccl_constants.MPC_TO_METER /ccl_constants.MPC_TO_METER  * ccl_constants.CLIGHT * ccl_constants.CLIGHT))

CCL_BEGIN_DECLS

typedef enum ccl_neutrino_mass_splits{
  ccl_nu_normal=0,
  ccl_nu_inverted=1,
  ccl_nu_equal=2,
  ccl_nu_sum=3,
  ccl_nu_single=4
} ccl_neutrino_mass_splits;

/**
 * Returns density of one neutrino species at a scale factor a.
 * Users are encouraged to access this quantity via the function ccl_omega_x.
 * @param a Scale factor
 * @param Neff The effective number of species with neutrino mass mnu.
 * @param mnu Pointer to array containing neutrino mass (can be 0).
 * @param T_CMB Temperature of the CMB
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * For specific cases see documentation for ccl_error.c
 * @return OmNuh2 Fractional energy density of neutrions with mass mnu, multiplied by h squared.
 */
double ccl_Omeganuh2(double a, int N_nu_mass, double* mnu, double T_CMB, int * status);

/**
 * Returns mass of one neutrino species at a scale factor a.
 * @param a Scale factor
 * @param Neff The effective number of species with neutrino mass mnu.
 * @param OmNuh2 Fractional energy density of neutrions with mass mnu, multiplied by h squared. (can be 0).
 * @param T_CMB Temperature of the CMB
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * For specific cases see documentation for ccl_error.c
 * @return Mnu Neutrino mass [eV].
 */
double* ccl_nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int * status);

CCL_END_DECLS
#endif

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL