text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */

#pragma once

#include <gsl/gsl-lite.hpp>
#include <nlohmann/json_fwd.hpp>
#include <cstdint>
#include <string>

/**
 * The Event class represents the information needed for a single
 * audit event entry.
 */
class Event {
public:
    Event() = delete;

    /**
     * Construct and initialize a new Event structure based off the
     * provided JSON. See ../README.md for information about the
     * layout of the JSON element.
     *
     * @param entry
     * @throws std::runtime_error for errors accessing the expected
     *                            elements
     */
    explicit Event(const nlohmann::json& json);

    /// The identifier for this entry
    uint32_t id;
    /// The name of the entry
    std::string name;
    /// The full description of the entry
    std::string description;
    /// Set to true if this entry should be handled synchronously
    bool sync;
    /// Set to true if this entry is enabled (or should be dropped)
    bool enabled;
    /// Set to true if the user may enable filtering for the enry
    bool filtering_permitted;
    /// The textual representation of the JSON describing mandatory
    /// fields in the event (NOTE: this is currently not enforced
    /// by the audit daemon)
    std::string mandatory_fields;
    /// The textual representation of the JSON describing the optional
    /// fields in the event (NOTE: this is currently not enforced
    /// by the audit daemon)
    std::string optional_fields;
};


text:
#include <jni.h>
#include <assert.h>
#include <lapacke.h>

/*extern void dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m,
                   int *n, double *alpha, double *A, int *lda, double *B,
                   int *ldb);*/

JNIEXPORT jint Java_JAMAJni_QRDecomposition_dgeqrf (JNIEnv *env, jclass klass, jint matrix_layout,
                                                    jint m, jint n, jdoubleArray a, jint lda,
                                                    jdoubleArray tau){

    double *aElems, *tauElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    tauElems = (*env)-> GetDoubleArrayElements (env, tau, NULL);

    assert(aElems && tauElems);

    info = LAPACKE_dgeqrf ((int) matrix_layout, (lapack_int) m, (lapack_int) n, aElems, (lapack_int) lda, tauElems);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, tau, tauElems, 0);

    return info;

}

JNIEXPORT jint Java_JAMAJni_QRDecomposition_dorgqr (JNIEnv *env, jclass klass, jint matrix_layout,
                                                    jint m, jint n, jint k, jdoubleArray a, jint lda,
                                                    jdoubleArray tau){

    double *aElems, *tauElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    tauElems = (*env)-> GetDoubleArrayElements (env, tau, NULL);

    assert(aElems && tauElems);

    info = LAPACKE_dorgqr((int) matrix_layout, (lapack_int) m, (lapack_int) n, (lapack_int) k, aElems, (lapack_int) lda, tauElems);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, tau, tauElems, 0);

    return info;

}

JNIEXPORT jint Java_JAMAJni_QRDecomposition_dgeqp3  (JNIEnv *env, jclass klass, jint matrix_layout,
                                                     jint m, jint n, jdoubleArray a, jint lda,
                                                     jintArray jpvt, jdoubleArray tau){

    double *aElems, *tauElems;
    lapack_int *jpvtElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    tauElems = (*env)-> GetDoubleArrayElements (env, tau, NULL);
    jpvtElems = (*env)-> GetIntArrayElements (env, jpvt, NULL);

    assert(aElems && tauElems && jpvtElems);

    info = LAPACKE_dgeqp3 ((int) matrix_layout, (lapack_int) m, (lapack_int) n, aElems, (lapack_int) lda, jpvtElems, tauElems);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, tau, tauElems, 0);
    (*env)-> ReleaseIntArrayElements (env, jpvt, jpvtElems, 0);

    return info;
}

JNIEXPORT jint Java_JAMAJni_QRDecomposition_dormqr
(JNIEnv *env, jclass klass, jint matrix_layout, jchar side, jchar trans,
 jint m, jint n, jint k, jdoubleArray a, jint lda, jdoubleArray tau,
 jdoubleArray c, jint ldc){

    double *aElems, *tauElems, *cElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    tauElems = (*env)-> GetDoubleArrayElements (env, tau, NULL);
    cElems = (*env)-> GetDoubleArrayElements (env, c, NULL);

    assert(aElems && tauElems && cElems);

    info = LAPACKE_dormqr ((int) matrix_layout, (char) side, (char) trans,
                           (lapack_int) m, (lapack_int) n, (lapack_int) k,
                           aElems, (lapack_int) lda, tauElems, cElems,
                           (lapack_int) ldc);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, tau, tauElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, c, cElems, JNI_ABORT);

    return info;

}

ext: c


text:
#pragma once

#include <cassert>
#include <stdint.h>

#include <gsl/gsl-lite.hpp> // gsl::span
#include "HypCommands.h" // RECORD_LENGTH

// Utility constexpr function converting a enum class to its underlying type N/U
#include <type_traits>
template <typename E,
          typename = std::enable_if_t<std::is_enum<E>::value>>
constexpr auto ut(E e) noexcept
{
    return static_cast<std::underlying_type_t<E>>(e);
}

class SamplePoint {
public:
    explicit SamplePoint(const gsl::span<uint8_t>& data);

    enum ValueIndex {
        eSeconds = 0,
        eVolts   ,
        eAmps    ,
        emAh_Out , // recalculated
        emAh_In  , // -"- (is rarely present)
        eRPM     , // []
        eAltitude, //  RDU
        eTemp1   , // []
        eTemp2   , // []
        eTemp3   , // []
        eThrottle, // [RDU]
        eTempAmb , //  RDU

        ePowerIn,
        ePowerOut, // [PRM]
        eEfficiency,//[PRM]
        eThrust,   // [RPM]
        eNUM_VALUES
    };

    inline double& operator [](size_t idx) {
        if (idx < eNUM_VALUES) {
            return m_values[idx];
        } else {
            assert(!"Invalid index");
            static double ERROR = 0.0;
            return ERROR;
        }
    }

    inline const double& operator [](size_t idx) const {
        if (idx < eNUM_VALUES) {
            return m_values[idx];
        } else {
            assert(!"Invalid index");
            static const double ERROR = 0.0;
            return ERROR;
        }
    }

//static:
    static const char* SeriesName(size_t idx) {
        if (idx < eNUM_VALUES) { return sSeriesNames[idx]; }
        else                   { return "ERROR"; }
    }

    static const std::array<size_t, eNUM_VALUES> sSeriesPrecision;

private:
    std::array<double, eNUM_VALUES> m_values;
//static
    static const std::array<const char*, eNUM_VALUES> sSeriesNames;
};

ext: h


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

#include "viaio/Vlib.h"
#include "viaio/VImage.h"
#include "viaio/mu.h"
#include "viaio/option.h"

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sort_vector.h>

/* two-pass formula, correct for round-off error */
void VNormalize(gsl_vector *vec)
{
  double ave,var,sd,nx,s,u,tiny=1.0e-6;
  size_t j,nt = vec->size;

  nx = (double)nt;
  ave = 0;
  for (j=0; j<nt; j++) ave += vec->data[j];
  ave /= nx;

  var = u = 0;
  for (j=0; j<nt; j++) {
    s = vec->data[j]-ave;
    u   += s;
    var += s*s;
  }
  var=(var-u*u/nx)/(nx-1);
  sd = sqrt(var);
  if (sd < tiny) {
    for (j=0; j<nt; j++) vec->data[j] = 0;
  }
  else {
    for (j=0; j<nt; j++) {
      u = vec->data[j];
      vec->data[j] = (u-ave)/sd;
    }
  }
}

/* convert to ranks for spearman correlation */
void GetRank(gsl_vector *vec,gsl_permutation *perm,gsl_permutation *rank)
{
  size_t i;
  size_t n = vec->size;
  gsl_sort_vector_index (perm, vec);
  gsl_permutation_inverse (rank, perm);
  for (i=0; i<n; i++) vec->data[i] = (float)rank->data[i];
}

size_t NumVoxels(VImage src)
{
  int b,r,c;
  size_t n=0;
  for (b=0; b<VImageNBands(src); b++) {
    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	if (VGetPixel(src,b,r,c) > 0.01) n++;
      }
    }
  }
  return n;
}

gsl_matrix_float *DataMatrix(VImage *src,int nslices,VImage roi,int metric)
{
  int b,r,c;
  float u=0;
  size_t nt = VImageNBands(src[0]);
  size_t nvox = NumVoxels(roi);
  gsl_matrix_float *X = gsl_matrix_float_calloc(nvox,nt);
  if (X==NULL) VError(" error allocating data matrix, size: %lu x %lu\n",nvox,nt);

  gsl_vector *vec = gsl_vector_calloc(nt);
  gsl_permutation *perm = NULL;
  gsl_permutation *rank = NULL;

  if (metric == 1) {  /* only needed for spearman correlation */
    perm = gsl_permutation_alloc(nt);
    rank = gsl_permutation_alloc(nt);
  }

  size_t i=0,j=0;
  for (b=0; b<VImageNBands(roi); b++) {
    for (r=0; r<VImageNRows(roi); r++) {
      for (c=0; c<VImageNColumns(roi); c++) {

	if (VGetPixel(roi,b,r,c) < 0.01) continue;

	for (j=0; j<nt; j++) {
	  u = (float) VGetPixel(src[b],j,r,c);
	  vec->data[j] = u;
	}

	if (metric == 0) {  /* linear correlation */
	  VNormalize(vec);
	}
	else if (metric == 1) {  /* convert to ranks for Spearman corr */
	  GetRank(vec,perm,rank);
	}
	else {
	  VError(" unknown metric %d",metric);
	}

	for (j=0; j<nt; j++) {
	  gsl_matrix_float_set(X,i,j,(float)vec->data[j]);
	}

	i++;
      }
    }
  }

  gsl_vector_free(vec);
  if (metric==1) {
    gsl_permutation_free (perm);
    gsl_permutation_free (rank);
  }
  return X;
}

/* new voxel map */
VImage VoxelMap(VImage roi)
{
  size_t nvox = NumVoxels(roi);
  VImage map = VCreateImage(1,4,nvox,VShortRepn);
  if (map == NULL) VError(" error allocating addr map");
  VFillImage(map,VAllBands,0);
  int nslices = VImageNBands(roi);
  int nrows = VImageNRows(roi);
  int ncols = VImageNColumns(roi);
  VSetAttr(VImageAttrList(map),"nvoxels",NULL,VLongRepn,(VLong)nvox);
  VSetAttr(VImageAttrList(map),"nslices",NULL,VLongRepn,(VLong)nslices);
  VSetAttr(VImageAttrList(map),"nrows",NULL,VLongRepn,(VLong)nrows);
  VSetAttr(VImageAttrList(map),"ncols",NULL,VLongRepn,(VLong)ncols);
  VPixel(map,0,3,0,VShort) = nslices;
  VPixel(map,0,3,1,VShort) = nrows;
  VPixel(map,0,3,2,VShort) = ncols;

  int b,r,c;
  size_t i = 0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {

	float u = VGetPixel(roi,b,r,c);
	if (u < 0.1) continue;
	VPixel(map,0,0,i,VShort) = b;
	VPixel(map,0,1,i,VShort) = r;
	VPixel(map,0,2,i,VShort) = c;
	i++;
      }
    }
  }
  return map;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef DART_META_H
#define DART_META_H

/*----- System Includes -----*/

#include <memory>
#include <gsl/gsl>
#include <type_traits>

/*----- Local Includes -----*/

#include "shim.h"

/*----- Type Declarations -----*/

namespace dart {

  namespace meta {

    template <class...>
    using void_t = void;

    template <class T>
    struct identity {
      using type = T;
    };

    // Removes one level of "pointer to const", "pointer to volatile"
    // "reference to const", "reference to volatile", etc,
    // to make working with strings easier.
    // Like std::decay except it doesn't decay C arrays,
    // but it otherwise much stronger.
    template <class T>
    struct canonical_type {
      using type = T;
    };
    template <class T>
    struct canonical_type<T const> {
      using type = T;
    };
    template <class T>
    struct canonical_type<T volatile> {
      using type = T;
    };
    template <class T>
    struct canonical_type<T const volatile> {
      using type = T;
    };
    template <class T>
    struct canonical_type<T&> {
      using type = T&;
    };
    template <class T>
    struct canonical_type<T const&> {
      using type = T&;
    };
    template <class T>
    struct canonical_type<T volatile&> {
      using type = T&;
    };
    template <class T>
    struct canonical_type<T const volatile&> {
      using type = T&;
    };
    template <class T>
    struct canonical_type<T*> {
      using type = T*;
    };
    template <class T>
    struct canonical_type<T const*> {
      using type = T*;
    };
    template <class T>
    struct canonical_type<T volatile*> {
      using type = T*;
    };
    template <class T>
    struct canonical_type<T const volatile*> {
      using type = T*;
    };
    template <class T, size_t len>
    struct canonical_type<T (&)[len]> {
      using type = T (&)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T const (&)[len]> {
      using type = T (&)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T volatile (&)[len]> {
      using type = T (&)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T const volatile (&)[len]> {
      using type = T (&)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T (*)[len]> {
      using type = T (*)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T const (*)[len]> {
      using type = T (*)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T volatile (*)[len]> {
      using type = T (*)[len];
    };
    template <class T, size_t len>
    struct canonical_type<T const volatile (*)[len]> {
      using type = T (*)[len];
    };
    template <class T>
    using canonical_type_t = typename canonical_type<T>::type;

    template <class T>
    struct is_ptr_to_const : std::false_type {};
    template <class T>
    struct is_ptr_to_const<T const*> : std::true_type {};

    // Create type traits to check for each type of comparison.
    template <class Lhs, class Rhs, class = void>
    struct are_comparable : std::false_type {};
    template <class Lhs, class Rhs>
    struct are_comparable<
      Lhs,
      Rhs,
      void_t<decltype(std::declval<Lhs>() == std::declval<Rhs>())>
    > : std::true_type {};

    template <class T, class = void>
    struct is_dereferenceable : std::false_type {};
    template <class T>
    struct is_dereferenceable<
      T,
      void_t<decltype(*std::declval<T>())>
    > : std::true_type {};

    template <class...>
    struct conjunction : std::true_type {};
    template <class B1>
    struct conjunction<B1> : B1 {};
    template <class B1, class... Bn>
    struct conjunction<B1, Bn...> :
      std::conditional_t<
        bool(B1::value),
        conjunction<Bn...>,
        B1
      >
    {};

    template <class...>
    struct disjunction : std::false_type {};
    template <class B1>
    struct disjunction<B1> : B1 {};
    template <class B1, class... Bn>
    struct disjunction<B1, Bn...> :
      std::conditional_t<
        bool(B1::value),
        B1,
        disjunction<Bn...>
      >
    {};

    template <class B>
    struct negation : std::integral_constant<bool, !bool(B::value)> {};

    template <class T, class... Ts>
    struct contained : disjunction<std::is_same<T, Ts>...> {};
    template <class T, class... Ts>
    struct not_contained : conjunction<negation<std::is_same<T, Ts>>...> {};

    template <class... Types>
    struct types {
      template <class T, class =
        std::enable_if_t<
          contained<T, Types...>::value
        >
      >
      types(T) {}
    };
    template <class... Types>
    struct not_types {
      template <class T, class =
        std::enable_if_t<
          not_contained<T, Types...>::value
        >
      >
      not_types(T) {}
    };

    template <class...>
    struct first_type {
      using type = void;
    };
    template <class T, class... Ts>
    struct first_type<T, Ts...> {
      using type = T;
    };
    template <class... Ts>
    using first_type_t = typename first_type<Ts...>::type;

    namespace detail {

      template <class Default, class AlwaysVoid, template <class...> class Op, class... Args>
      struct detector {
        using value_t = std::false_type;
        using type = Default;
      };

      template <class Default, template <class...> class Op, class... Args>
      struct detector<Default, void_t<Op<Args...>>, Op, Args...> {
        // Note that std::void_t is a C++17 feature
        using value_t = std::true_type;
        using type = Op<Args...>;
      };

    }

    struct nonesuch {
      nonesuch(nonesuch const&) = delete;
      nonesuch(nonesuch&&) = delete;
      ~nonesuch() = delete;
    };

    template <template <class...> class Op, class... Args>
    using is_detected = typename detail::detector<nonesuch, void, Op, Args...>::value_t;
    template <template <class...> class Op, class... Args>
    using detected_t = typename detail::detector<nonesuch, void, Op, Args...>::type;
    template <class Default, template <class...> class Op, class... Args>
    using detected_or = detail::detector<Default, void, Op, Args...>;

    template <class T>
    using strv_t = decltype(std::declval<T>().strv());
    template <class T>
    using integer_t = decltype(std::declval<T>().integer());
    template <class T>
    using decimal_t = decltype(std::declval<T>().decimal());
    template <class T>
    using boolean_t = decltype(std::declval<T>().boolean());
    template <class T>
    using get_type_t = decltype(std::declval<T>().get_type());

    template <class T>
    struct is_dartlike :
      meta::conjunction<
        is_detected<strv_t, T>,
        is_detected<integer_t, T>,
        is_detected<decimal_t, T>,
        is_detected<boolean_t, T>,
        is_detected<get_type_t, T>
      >
    {};

    template <class T, template <class> class Template>
    struct is_specialization_of : std::false_type {};
    template <class... Args, template <class> class Template>
    struct is_specialization_of<Template<Args...>, Template> : std::true_type {};

    template <class T, template <template <class> class> class Compound>
    struct is_higher_specialization_of : std::false_type {};
    template <template <class> class Inner, template <template <class> class> class Outer>
    struct is_higher_specialization_of<Outer<Inner>, Outer> : std::true_type {};

    template <class T>
    struct is_span : std::false_type {};
#ifdef DART_USING_GSL_LITE
    template <class T>
    struct is_span<gsl::span<T>> : std::true_type {};
#else
    template <class T, std::size_t extent>
    struct is_span<gsl::span<T, extent>> : std::true_type {};
#endif

    template <class T>
    struct is_std_smart_ptr : std::false_type {};
    template <class T, class Del>
    struct is_std_smart_ptr<std::unique_ptr<T, Del>> : std::true_type {};
    template <class T>
    struct is_std_smart_ptr<std::shared_ptr<T>> : std::true_type {};
    template <class T>
    struct is_std_smart_ptr<std::weak_ptr<T>> : std::true_type {};

    template <size_t pos>
    struct priority_tag : priority_tag<pos - 1> {};
    template <>
    struct priority_tag<0> {};

    namespace detail {
      template <class T>
      struct is_builtin_string_impl : std::false_type {};
      template <>
      struct is_builtin_string_impl<char*> : std::true_type {};
      template <>
      struct is_builtin_string_impl<wchar_t*> : std::true_type {};
      template <>
      struct is_builtin_string_impl<char16_t*> : std::true_type {};
      template <>
      struct is_builtin_string_impl<char32_t*> : std::true_type {};
      template <size_t len>
      struct is_builtin_string_impl<char (&)[len]> : std::true_type {};
      template <size_t len>
      struct is_builtin_string_impl<wchar_t (&)[len]> : std::true_type {};
      template <size_t len>
      struct is_builtin_string_impl<char16_t (&)[len]> : std::true_type {};
      template <size_t len>
      struct is_builtin_string_impl<char32_t (&)[len]> : std::true_type {};
    }

    template <class T>
    struct is_builtin_string : detail::is_builtin_string_impl<canonical_type_t<T>> {};

    template <class T>
    struct is_std_string : std::false_type {};
    template <class CharT, class Traits, class Allocator>
    struct is_std_string<std::basic_string<CharT, Traits, Allocator>> : std::true_type {};

    template <class T>
    struct is_std_string_view : std::false_type {};
    template <class CharT, class Traits>
    struct is_std_string_view<shim::basic_string_view<CharT, Traits>> : std::true_type {};

    template <class T>
    struct is_string :
      meta::disjunction<
        is_builtin_string<T>,
        is_std_string<T>,
        is_std_string_view<T>
      >
    {};

    template <class T>
    struct spannable_value_type {
      template <class U>
      static typename U::value_type detect(priority_tag<1>);
      template <class U>
      static typename U::element_type detect(priority_tag<0>);
      template <class U>
      static nonesuch detect(...);
      using type = std::remove_cv_t<decltype(detect<T>(priority_tag<1> {}))>;
    };
    template <class T>
    using spannable_value_type_t = typename spannable_value_type<T>::type;
  }

}

#undef DART_COMPARE_HELPER

#endif


text:
#ifndef _LINEAR_ALGEBRA_HELPERS_H_
#define _LINEAR_ALGEBRA_HELPERS_H_
#include <stddef.h>

#ifdef HAVE_LAPACK
#include <lapacke.h>
#include <cblas.h>
#endif

/* Define constants */
#define LAH_PI 3.14159265358979323846

/*--- Data types ------------------------------------------------------------*/

#ifdef HAVE_LAPACK
typedef size_t lah_index;
#else
typedef size_t lah_index;
#endif

/* Choose type of data */
#ifdef USE_SINGLE  /* Single or double precission */
    typedef float lah_value;
#else
    typedef double lah_value;
#endif

/* Automatically choose BLAS/Lapack routines based on data type */
/* If not OpenBLAS (cblas, lapacke) is used,
 * here is the place to account for another calling convention */
#ifdef HAVE_LAPACK
    #ifdef USE_SINGLE
        #define GEMM cblas_sgemm
        #define POTRF LAPACKE_spotrf
        #define POTRS LAPACKE_spotrs
        #define GETRF LAPACKE_sgetrf
        #define GETRS LAPACKE_sgetrs
        #define GEQRF LAPACKE_sgeqrf
        #define SYEVR LAPACKE_ssyevr
        #define GESDD LAPACKE_sgesdd
        #define LAMCH LAPACKE_slamch
    #else
        #define GEMM cblas_dgemm
        #define POTRF LAPACKE_dpotrf /* Chol Decomposition */
        #define POTRS LAPACKE_dpotrs /* Solving Triangular System */
        #define GETRF LAPACKE_dgetrf /* LU Decomposition */
        #define GETRS LAPACKE_dgetrs /* Solving Triangular System */
        #define GEQRF LAPACKE_dgeqrf /* QR factorization */
        #define SYEVR LAPACKE_dsyevr /* SVD for matSqrt (symmetric matrix) */
        #define GESDD LAPACKE_dgesdd /* General SVD with divide and conquer */
        #define LAMCH LAPACKE_dlamch /* Helper to obtain machine precision */
    #endif /* end Data_Types*/
#endif /* end HAVE_LAPACK*/

/* Check if called from C++ */
#ifdef __cplusplus
	extern "C" {
#endif

/*--- Definitions for Linear Algebra subroutines*/
/* Specifies whether Row- or ColMajor */
typedef enum
{
    lahRowMajor=101,
    lahColMajor=102
} lah_layout;

#define LAH_LAYOUT lahRowMajor
/*#define LAH_LAYOUT(A) (((A)->matType & lahTypeColMajor) ? lahColMajor : lahRowMajor) */
#ifdef HAVE_LAPACK
    #if LAH_LAYOUT == lahColMajor
        #define LAH_CBLAS_LAYOUT CblasColMajor
        #define LAH_LAPACK_LAYOUT LAPACK_COL_MAJOR
        #define LAH_LEADING_DIM(A) ((A)->nR)
    #else
        #define LAH_CBLAS_LAYOUT CblasRowMajor
        #define LAH_LAPACK_LAYOUT LAPACK_ROW_MAJOR
        #define LAH_LEADING_DIM(A) ((A)->nC)
    #endif
#endif

/* Option for LA operations: Transpose Input before Computation or not? */
typedef enum
{
    lahTrans = 'T',
    lahNorm = 'N'
} lah_MatOp;

/* Error Types definition*/
typedef enum
{
    lahReturnOk = 0,          /* Success */
    lahReturnParameterError,  /* Wrong Parameters */
    lahReturnMathError,       /* General Math Error, either error while
                                 * computing or prerequisit not met */
    lahReturnFPError,          /* Error when evaluating function pointer */
    lahReturnExternError          /* Error when evaluating function pointer */
} lah_Return;

typedef enum
{
    /* Properties, which are known in matAlloc */
    lahTypeColMajor = 0x01,     /* 0000 0001 */
    lahTypeNoData = 0x02,       /* 0000 0010 */
    lahTypeVector = 0x04,       /* 0000 0100 */
    lahTypeSquare = 0x08,       /* 0000 1000 */

    /* Properties, which are not known in matAlloc */
    lahTypeDiagonal = 0x10,     /* 0001 0000 */
    lahTypeSymmetric = 0x20,    /* 0010 0000 */
    lahTypePositive = 0x40,     /* 0100 0000 */
    lahTypeOrthogonal = 0x80    /* 1000 0000 */
} lah_MatType;

/* lah_mat: Matrix format */
typedef struct
{
    lah_index nC;          /* Number of Columns */
    lah_index nR;          /* Number of Rows */
    lah_index incRow;      /* increment of rows */
    lah_index incCol;      /* increment of columns */
    lah_value *data;       /* pointer to data */
    lah_MatType matType;   /* flags for several properties of Matrix */
} lah_mat;

/*---------------------------------------------------------------------------*/

/*--- Common Linear Algebra helpers ----------------------------------------*/
/* If provided, these functions are using a CBLAS/LAPACKE library           */
/* else they use a simple handwritten defintion                             */
/*--------------------------------------------------------------------------*/

/* Computes the general Matrix mulitplikation of form:
 * C = alpha*C + beta*op(A)*op(B)
 * similar to the BLAS routine ?gemm.
 * op(A) = A for transA=koopNorm or op(A) = At for transA = koopTrans */
lah_Return lah_matMul(lah_MatOp transA, lah_MatOp transB,
                        lah_value alpha, lah_value beta,
                        lah_mat *C, lah_mat const *A, lah_mat const *B);

/* Adds two matrices C = alpha*A + beta*B */
lah_Return lah_matAdd(lah_mat *C, lah_value alpha, lah_mat const *A,
                      lah_value beta, lah_mat const *B);

/* This will return a lah_mat struct At with inverted order
 * (RowMajor <-> ColMajor). If At is not needed anymore, call free(At)
 * NOT lah_matFree(At), as this will also delete A->data                   */
lah_mat *lah_matTrans(lah_mat const *A);

/* Computes the Cholesky factorization of the form P = L * Lt
 * where P is a symmetric positive-definite matrix, L is a lower triangular
 * matrix and Lt denotes its transpose.                                     */
lah_Return lah_chol(lah_mat *P, lah_index setZero);

/* Computes a rank-1 Update to a Cholesky factorization.
 * Let L is the Cholesky factor of a spd matrix P = L * Lt, then
 * the update is the Cholesky factor of the matrix P_new = P + x * x' */
lah_Return lah_cholUpdate(lah_mat *L, lah_mat *x, lah_value alpha);

/* Computes the Forward Substitution of the linear equation system L * X = B
 * where L is a lower triangular matrix,
 * x and b are matrices with the same number of columns.                    */
lah_Return lah_forwardSub(lah_mat *B, lah_mat const *L);

/* Computes an Update to a matrix C of form D = A*C*A'+B  */
/* workspace needed for special case C = A*C*A'+B         */
lah_Return lah_matUpdate(lah_mat *D, lah_mat const *C,
                           lah_mat const *A, lah_mat const *B,
                           lah_mat *W);
/*--------------------------------------------------------------------------*/

/*--- Function only defined for HAVE_LAPACK --------------------------------*/

#ifdef HAVE_LAPACK

/* Computes the LU factorization of the form A = L * U,
 * as an alternative to Cholesky for general matrices*/
lah_Return lah_LU(lah_mat *A, lah_index setZero, lapack_int *ipiv);

/* Solves the linear equation system L * X = B
 * arising from the LU factorization */
/*lah_Return lah_solveLU(lah_mat *X, const lah_mat *L,
                       const lah_mat *B, lah_index *ipiv);
*/
lah_Return lah_solveLU(lah_mat *B, const lah_mat *L, lapack_int *ipiv);

/* Computes the QR factorization of the form A = L * U,
 * as an alternative to Cholesky for general matrices*/
lah_Return lah_QR(lah_mat *A, lah_value *tau);

/* Computes the eigen values and vectors of a symmetric matrix P*/
lah_Return lah_eigenValue(lah_mat *P, lah_value *S, lah_mat *Z);

/* Computes the Sqrt of a positive Matrix by computing the eigen values */
lah_Return lah_matSqrt(lah_mat *P, lah_mat *Psqrt, lah_mat *Work);

/*General Singular Value Decomposition */
lah_Return lah_SVD(lah_mat *A, lah_mat *U, lah_value *S, lah_mat *Vt);
#endif

/*--------------------------------------------------------------------------*/

/*--- Helper functions -----------------------------------------------------*/

/*--- Utility functions ---*/
/* allocation for the Matrix struct lah_mat */
lah_mat* lah_matAlloc(lah_index nC, lah_index nR, lah_index AllocData);

/* Free all data used bei the Matrix struct lah_mat */
lah_mat* lah_matFree(lah_mat *mat);

/* load a matrix from a file */
lah_mat* lah_matLoad(const char *fname);

/* Print a lah_mat Matrix struct (used for debugging) */
lah_Return lah_matPrint(lah_mat const *A, lah_index brief);

/*--- Noise Generators ---*/
/* Returns a random number from a Gaussian distribution */
lah_value lah_gaussNoise();

/* Returns a random number from a Lorentz (Cauchy) distribution */
lah_value lah_lorentzNoise();

/*--- Property checks ---*/
/* Checks if matrix is symmetric.
 * It compares with $(precision) significant digits */
lah_index lah_isSymmetric(lah_mat *A, lah_index precision);

/* Checks if matrix is positive definite. It tries to compute a Cholesky
 * factorization, which only works for positive definite matrices*/
lah_index lah_isPositive(lah_mat *A);

/*--- Construct different Matrix types ---*/
/* Construct a diagonal matrix from n parameters */
lah_mat *lah_matConstructDiag(lah_index nC, const lah_value *parameter_vec);

/* Construct a symmetric Matrix from n(n+1)/2 parameters */
lah_mat *lah_matConstructSy(lah_index nC, const lah_value *parameter_vec);

/* Construct a orthogonal Matrix from n(n+1)/2 parameters */
lah_mat *lah_matConstructOrtho(lah_index nC, lah_value *parameter_vec);

/* Construct a positive Matrix from n(n+3)/2 parameters */
lah_mat *lah_matConstructPo(lah_index nC, lah_value *parameter_vec);

/* Construct a positive Matrix the simple way from a nxn Maxtrix */
lah_mat *lah_matConstructPo_simple(lah_mat *A);

/*--- Macro definitions ----------------------------------------------------*/

/* Macro to access element of Matrix A[i,j] */
#define LAH_ENTRY(A, i, j) ((A)->data[(i) * (A)->incRow + (j) * (A)->incCol])
#define LAH_MAX(A, B) (((A) > (B)) ? (A) : (B))
#define LAH_MIN(A, B) (((A) < (B)) ? (A) : (B))
#define LAH_SIGN(A) ((A) < 0) ? (-1) : (1))

/*--- Property MACROS---*/
#define LAH_ISTYPE(A, Type) ((A)->matType & (Type))        /* AND */
#define LAH_SETTYPE(A, Type) ((A)->matType |= (Type))      /* OR */
#define LAH_UNSETTYPE(A, Type) ((A)->matType &= (~Type)) /* NAND */
#define LAH_TOGGLETYPE(A, Type) ((A)->matType ^= (Type))   /* XOR */

/*--- LAPACK_ONLY HEADERS ---*/
#ifdef HAVE_LAPACK

#endif

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

#endif


text:
#include <gsl/gsl_spmatrix.h>
#include <munit.h>
#include <osqp/cs.h>
#include <osqp/osqp.h>
#include <unistd.h>

#include <qdm.h>

#include "test.h"

static
MunitResult
test_qdm_vector_seq(
    const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  double min = atof(munit_parameters_get(params, "qdm_vector_seq_min"));
  double max = atof(munit_parameters_get(params, "qdm_vector_seq_max"));
  double delta = atof(munit_parameters_get(params, "qdm_vector_seq_delta"));

  if (min >= max) {
    return MUNIT_SKIP;
  }

  gsl_vector *seq = qdm_vector_seq(min, max, delta);

  munit_assert_size(((max - min) / delta) + 1, ==, seq->size);

  for (size_t i = 0; i < seq->size; i++) {
    double expecting = min + delta * i;
    double found = gsl_vector_get(seq, i);

    munit_logf(MUNIT_LOG_DEBUG, "%.17g", found);
    munit_assert_double_equal(expecting, found, 3);
  }

  return MUNIT_OK;
}

static
MunitResult
test_qdm_vector_set_seq(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  gsl_vector *seq = gsl_vector_alloc(10);

  qdm_vector_set_seq(seq, 0, 9);

  for (size_t i = 0; i < seq->size; i++) {
    double expecting = (double)i;
    double found = gsl_vector_get(seq, i);

    munit_logf(MUNIT_LOG_DEBUG, "e %.17g f %.17g", expecting, found);
    munit_assert_double_equal(expecting, found, 3);
  }

  qdm_vector_set_seq(seq, 9, 0);

  for (size_t i = 0; i < seq->size; i--) {
    double expecting = (double)(9 - i);
    double found = gsl_vector_get(seq, i);

    munit_logf(MUNIT_LOG_DEBUG, "e %.17g f %.17g", expecting, found);
    munit_assert_double_equal(expecting, found, 3);
  }

  qdm_vector_set_seq(seq, 0, 1);
  double data[10] = {
    0.0,
    0.1,
    0.2,
    0.3,
    0.4,
    0.5,
    0.6,
    0.7,
    0.8,
    0.9,
  };

  for (size_t i = 0; i < seq->size; i--) {
    double expecting = data[i];
    double found = gsl_vector_get(seq, i);

    munit_logf(MUNIT_LOG_DEBUG, "e %.17g f %.17g", expecting, found);
    munit_assert_double_equal(expecting, found, 3);
  }

  return MUNIT_OK;
}

static
MunitResult
test_qdm_matrix_det_tmm(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  int status = 0;

  /* Zero */
  double data_zero[4][3] = {
      {1,1,1},
      {2,2,2},
      {3,3,3},
      {4,4,4},
  };

  gsl_matrix_view m = gsl_matrix_view_array((double *)data_zero, 4, 3);

  double det_zero = 0;
  status = qdm_matrix_det_tmm(&m.matrix, &det_zero);
  munit_assert_int(status, ==, 0);
  munit_assert_double_equal(0, det_zero, 3);

  /* Non-zero */
  double data_nonzero[4][3] = {
      {1,2,3},
      {4,5,2},
      {1,1,3},
      {1,1,4},
  };

  m = gsl_matrix_view_array((double *)data_nonzero, 4, 3);

  double det_nonzero = 0;
  status = qdm_matrix_det_tmm(&m.matrix, &det_nonzero);
  munit_assert_int(status, ==, 0);
  munit_assert_double_equal(271, det_nonzero, 3);

  return MUNIT_OK;
}

static
MunitResult
test_qdm_vector_quantile(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  /* Basic */
  double data[10] = {
    1,
    2,
    3,
    4,
    5,
    6,
    7,
    8,
    9,
    10,
  };

  gsl_vector_view v = gsl_vector_view_array((double *)data, 10);
  gsl_vector *probs = qdm_vector_seq(0, 1, 0.25);
  gsl_vector *quantiles = qdm_vector_quantile(&v.vector, probs);

  munit_log(MUNIT_LOG_DEBUG, "basic");
  munit_log_vector(MUNIT_LOG_DEBUG, &v.vector);
  munit_log_vector(MUNIT_LOG_DEBUG, probs);
  munit_log_vector(MUNIT_LOG_DEBUG, quantiles);
  munit_assert_int(quantiles->size, ==, 5);

  double expected[5] = {
    1,
    3.25,
    5.5,
    7.75,
    10,
  };

  for (size_t i = 0; i < 5; i++) {
    munit_assert_double_equal(expected[i], gsl_vector_get(quantiles, i), 3);
  }

  /* Random */
  double rnorm_data[10] = {
    -1.7792761,
    -0.9436401,
    -0.8845596,
    -0.4213213,
    -0.1831428,
    -0.1497455,
    0.1402677,
    0.4509039,
    0.7415657,
    0.8864738,
  };

  v = gsl_vector_view_array((double *)rnorm_data, 10);
  quantiles = qdm_vector_quantile(&v.vector, probs);

  munit_log(MUNIT_LOG_DEBUG, "random");
  munit_log_vector(MUNIT_LOG_DEBUG, &v.vector);
  munit_log_vector(MUNIT_LOG_DEBUG, probs);
  munit_log_vector(MUNIT_LOG_DEBUG, quantiles);
  munit_assert_int(quantiles->size, ==, 5);

  double rnorm_expected[5] = {
    -1.7792761,
    -0.7687500,
    -0.1664441,
    0.3732449,
    0.8864738,
  };

  for (size_t i = 0; i < 5; i++) {
    munit_assert_double_equal(rnorm_expected[i], gsl_vector_get(quantiles, i), 3);
  }

  return MUNIT_OK;
}

static
MunitResult
test_qdm_vector_sum(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  gsl_vector_view v;
  double sum = 0.0;

  /* Basic */
  double basic_data[3] = {
    1,
    2,
    3,
  };

  v = gsl_vector_view_array((double *)basic_data, 3);
  sum = qdm_vector_sum(&v.vector);
  munit_assert_double(6, ==, sum);

  /* Interleaved alternating big sign */
  double iabs_data[4] = {
    1,
    1e100,
    1,
    -1e100,
  };

  v = gsl_vector_view_array((double *)iabs_data, 4);
  sum = qdm_vector_sum(&v.vector);
  munit_assert_double(2, ==, sum);

  /* Long sequence */
  gsl_vector *ls0_data = qdm_vector_seq(100, 1000, 1);
  sum = qdm_vector_sum(ls0_data);
  munit_assert_double(495550, ==, sum);
  gsl_vector_free(ls0_data);

  gsl_vector *ls1_data = qdm_vector_seq(10, 100, 0.1);
  sum = qdm_vector_sum(ls1_data);
  munit_assert_double_equal(49555, sum, 10);
  gsl_vector_free(ls1_data);

  gsl_vector *ls2_data = qdm_vector_seq(1, 10, 0.01);
  sum = qdm_vector_sum(ls2_data);
  munit_assert_double_equal(4955.5, sum, 10);
  gsl_vector_free(ls2_data);

  gsl_vector *ls3_data = qdm_vector_seq(0.000000001, 0.00000001, 0.00000000001);
  sum = qdm_vector_sum(ls3_data);
  munit_assert_double_equal(4.9555e-06, sum, 10);
  gsl_vector_free(ls3_data);

  return MUNIT_OK;
}

static
MunitResult
test_qdm_matrix_select_upper_triangle(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  double grid_data[] = {
      1,2,3,
      4,5,6,
      7,8,9,
  };

  gsl_matrix_view grid = gsl_matrix_view_array(grid_data, 3, 3);

  double expected_data[] = {
      1,2,3,
      0,5,6,
      0,0,9,
  };

  gsl_matrix_view expected = gsl_matrix_view_array(expected_data, 3, 3);

  qdm_matrix_select_upper_triangle(&grid.matrix);

  for (size_t i = 0; i < 3; i++) {
    for (size_t j = 0; j < 3; j++) {
      double e = gsl_matrix_get(&expected.matrix, i, j);
      double r = gsl_matrix_get(&grid.matrix, i, j);

      munit_logf(MUNIT_LOG_DEBUG, "(%zu, %zu)", i, j);
      munit_assert_double(e, ==, r);
    }
  }

  return MUNIT_OK;
}

static
MunitResult
test_qdm_matrix_to_csc_matrix(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  int status = 0;

  int n = 2;

  int P_nnz = 3;
  double P_x[3] = {4.0, 1.0, 2.0};
  int P_i[3] = {0, 0, 1};
  int P_p[3] = {0, 1, 3};

  double input_data[] = {
    4, 1,
    1, 2,
  };

  gsl_matrix_view input = gsl_matrix_view_array(input_data, 2, 2);
  qdm_matrix_select_upper_triangle(&input.matrix);

  munit_log(MUNIT_LOG_DEBUG, "input");
  qdm_matrix_csv_fwrite(stderr, &input.matrix);

  csc *result = NULL;
  status = qdm_matrix_to_csc_matrix(&result, &input.matrix);
  munit_assert_int(status, ==, 0);

  munit_assert_int(P_nnz, ==, result->nzmax);
  munit_assert_int(n, ==, result->n);

  for (size_t i = 0; i < (size_t)(result->n + 1); i++) {
    munit_logf(MUNIT_LOG_DEBUG, "result.p[%zu] = %i", i, result->p[i]);
    munit_assert_int(result->p[i], ==, P_p[i]);
  }

  for (size_t i = 0; i < (size_t)result->nzmax; i++) {
    munit_logf(MUNIT_LOG_DEBUG, "result.i[%zu] = %i", i, result->i[i]);
    munit_assert_int(result->i[i], ==, P_i[i]);
  }

  for (size_t i = 0; i < (size_t)result->nzmax; i++) {
    munit_logf(MUNIT_LOG_DEBUG, "result.x[%zu] = %.17g", i, result->x[i]);
    munit_assert_double(result->x[i], ==, P_x[i]);
  }

  return MUNIT_OK;
}

static
MunitResult
test_qdm_vector_hd5_write(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  int status = 0;

  double v_data[] = {
    1, 2, 3, 4, 5, 6, 7, 8, 9,
  };

  gsl_vector_view v = gsl_vector_view_array(v_data, 9);

  int tf_fd = 0;
  char tf_name[] = "/tmp/qdm-test-hd5-vector-XXXXXX";
  tf_fd = mkstemp(tf_name);
  close(tf_fd);

  hid_t tf_file = qdm_data_create_file(tf_name);
  munit_assert_int(tf_file, >=, 0);
  unlink(tf_name);

  hid_t tf_group = qdm_data_create_group(tf_file, "data");
  munit_assert_int(tf_group, >=, 0);

  status = qdm_vector_hd5_write(tf_group, "v", &v.vector);
  munit_assert_int(status, ==, 0);

  status = H5Gclose(tf_group);
  munit_assert_int(status, ==, 0);

  status = H5Fclose(tf_file);
  munit_assert_int(status, ==, 0);

  return MUNIT_OK;
}

static
MunitResult
test_qdm_matrix_hd5_write(
    MUNIT_UNUSED const MunitParameter params[],
    MUNIT_UNUSED void* fixture
)
{
  int status = 0;

  double m_data[] = {
    1, 2, 3,
    4, 5, 6,
    7, 8, 9,
  };

  gsl_matrix_view m = gsl_matrix_view_array(m_data, 3, 3);

  int tf_fd = 0;
  char tf_name[] = "/tmp/qdm-test-hd5-matrix-XXXXXX";
  tf_fd = mkstemp(tf_name);
  close(tf_fd);

  hid_t tf_file = qdm_data_create_file(tf_name);
  munit_assert_int(tf_file, >=, 0);
  unlink(tf_name);

  hid_t tf_group = qdm_data_create_group(tf_file, "data");
  munit_assert_int(tf_group, >=, 0);

  status = qdm_matrix_hd5_write(tf_group, "m", &m.matrix);
  munit_assert_int(status, ==, 0);

  status = H5Gclose(tf_group);
  munit_assert_int(status, ==, 0);

  status = H5Fclose(tf_file);
  munit_assert_int(status, ==, 0);

  return MUNIT_OK;
}

static char *test_param_qdm_vector_seq_min[] = {
  "1",
  "10",
  "0.1",
  NULL,
};

static char *test_param_qdm_vector_seq_max[] = {
  "10",
  "20",
  "1.1",
  NULL,
};

static char *test_param_qdm_vector_seq_delta[] = {
  "1",
  "0.1",
  "2",
  NULL,
};

static MunitParameterEnum test_params_qdm_vector_seq[] = {
  { "qdm_vector_seq_min"   , test_param_qdm_vector_seq_min   }  ,
  { "qdm_vector_seq_max"   , test_param_qdm_vector_seq_max   }  ,
  { "qdm_vector_seq_delta" , test_param_qdm_vector_seq_delta }  ,
  { NULL                   , NULL                            }  ,
};

MunitTest tests_gsl[] = {
  {
    "/qdm_vector_seq"          , // name
    test_qdm_vector_seq        , // test
    NULL                       , // setup
    NULL                       , // tear_down
    MUNIT_TEST_OPTION_NONE     , // options
    test_params_qdm_vector_seq , // parameters
  },
  {
    "/qdm_vector_set_seq",   // name
    test_qdm_vector_set_seq, // test
    NULL,                    // setup
    NULL,                    // tear_down
    MUNIT_TEST_OPTION_NONE,  // options
    NULL,                    // parameters
  },
  {
    "/qdm_matrix_det_tmm"   , // name
    test_qdm_matrix_det_tmm , // test
    NULL                    , // setup
    NULL                    , // tear_down
    MUNIT_TEST_OPTION_NONE  , // options
    NULL                    , // parameters
  },
  {
    "/qdm_vector_quantile"   , // name
    test_qdm_vector_quantile , // test
    NULL                     , // setup
    NULL                     , // tear_down
    MUNIT_TEST_OPTION_NONE   , // options
    NULL                     , // parameters
  },
  {
    "/qdm_vector_sum"      , // name
    test_qdm_vector_sum    , // test
    NULL                   , // setup
    NULL                   , // tear_down
    MUNIT_TEST_OPTION_NONE , // options
    NULL                   , // parameters
  },
  {
    "/qdm_matrix_select_upper_triangle"   , // name
    test_qdm_matrix_select_upper_triangle , // test
    NULL                                  , // setup
    NULL                                  , // tear_down
    MUNIT_TEST_OPTION_NONE                , // options
    NULL                                  , // parameters
  },
  {
    "/qdm_matrix_to_csc_matrix"   , // name
    test_qdm_matrix_to_csc_matrix , // test
    NULL                          , // setup
    NULL                          , // tear_down
    MUNIT_TEST_OPTION_NONE        , // options
    NULL                          , // parameters
  },
  {
    "/qdm_vector_hd5_write"   , // name
    test_qdm_vector_hd5_write , // test
    NULL                      , // setup
    NULL                      , // tear_down
    MUNIT_TEST_OPTION_NONE    , // options
    NULL                      , // parameters
  },
  {
    "/qdm_matrix_hd5_write"   , // name
    test_qdm_matrix_hd5_write , // test
    NULL                      , // setup
    NULL                      , // tear_down
    MUNIT_TEST_OPTION_NONE    , // options
    NULL                      , // parameters
  },
  { NULL, NULL, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
};

static const MunitSuite test_suite_gsl = {
  "/gsl"                  , // name
  tests_gsl               , // tests
  NULL                    , // suites
  1                       , // iterations
  MUNIT_SUITE_OPTION_NONE , // options
};

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
static char help[] =
   "Newton's method for a two-variable system.  Implements analytical Jacobian.\n"
   "Adds struct to hold parameter.\n";

#include <petsc.h>

typedef struct {
  double  b;
} AppCtx;

extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);
extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void*);

//STARTMAIN
int main(int argc,char **argv) {
  SNES   snes;         // nonlinear solver context
  Vec    x,r;          // solution, residual vectors
  Mat    J;
  AppCtx user;
  PetscErrorCode ierr;

  PetscInitialize(&argc,&argv,NULL,help);
  user.b = 2.0;

  ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,2); CHKERRQ(ierr);
  ierr = VecSetFromOptions(x); CHKERRQ(ierr);
  ierr = VecDuplicate(x,&r); CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&J); CHKERRQ(ierr);
  ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2); CHKERRQ(ierr);
  ierr = MatSetFromOptions(J); CHKERRQ(ierr);
  ierr = MatSetUp(J); CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
  ierr = SNESSetFunction(snes,r,FormFunction,&user);CHKERRQ(ierr);
  ierr = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  ierr = VecSet(x,1.0);CHKERRQ(ierr);
  ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  VecDestroy(&x);  VecDestroy(&r);  SNESDestroy(&snes);  MatDestroy(&J);
  return PetscFinalize();
}
//ENDMAIN

//STARTFORM
PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void *ctx) {
    PetscErrorCode ierr;
    AppCtx       *user = (AppCtx*)ctx;
    const double b = user->b, *ax;
    double       *aF;

    ierr = VecGetArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF);CHKERRQ(ierr);
    aF[0] = (1.0 / b) * PetscExpReal(b * ax[0]) - ax[1];
    aF[1] = ax[0] * ax[0] + ax[1] * ax[1] - 1.0;
    ierr = VecRestoreArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecRestoreArray(F,&aF);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode FormJacobian(SNES snes, Vec x, Mat J, Mat P, void *ctx) {
    PetscErrorCode ierr;
    AppCtx       *user = (AppCtx*)ctx;
    const double b = user->b, *ax;
    double       v[4];
    int          row[2] = {0,1}, col[2] = {0,1};

    ierr = VecGetArrayRead(x,&ax); CHKERRQ(ierr);
    v[0] = PetscExpReal(b * ax[0]);  v[1] = -1.0;
    v[2] = 2.0 * ax[0];              v[3] = 2.0 * ax[1];
    ierr = VecRestoreArrayRead(x,&ax); CHKERRQ(ierr);
    ierr = MatSetValues(P,2,row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDFORM


text:
//
// Created by L.Jonathan Feldstein
//

#ifndef CIMPLE_SYSTEM_H
#define CIMPLE_SYSTEM_H

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
#include "cimple_polytope_library.h"

// EXAMPLE ALLOCATION FOR STRUCT
// Try to allocate structure.
// struct struct_name *return_**** = malloc(sizeof(struct));
// Try to allocate data, free structure if fail.
// return_****->**** = malloc(***);
// Free structure if fail, return NULL
// if (return_****->**** == NULL){
//    return NULL;
//}

/**
 * Functions that help perform several mathematical operations
 * on the way to calculating the input.
 */
typedef struct auxiliary_matrices{

    gsl_matrix * A_K;
    gsl_matrix * A_N;
    gsl_matrix * Ct;

    /*
    B_diag = B
    E_diag = E
    for i in range(N-1):
    B_diag = _block_diag2(B_diag, B)
    E_diag = _block_diag2(E_diag, E)
     */
    gsl_matrix * B_diag;
    gsl_matrix * E_diag;

    gsl_matrix * L_default;
    gsl_matrix * E_default;

    //LU = ([Uset->size1*N, n+N*m])

    gsl_matrix * LU;

    //MU = tile(Uset->G, N)

    gsl_vector * MU;

    //GU = ([Uset->size1*N, p*N])

    gsl_matrix * GU;

    //K_hat = tile(K, N)

    gsl_vector * K_hat;

}auxiliary_matrices;

/**
 * Cost function to be minimized:
 *
 * |Rx|_{ord} + |Qu|_{ord} + r'x + mid_weight * |xc - x(N)|_{ord}
 *
 *      R: state cost matrix for:
 *      x = [x(1)' x(2)' .. x(N)']'
 *      If empty, zero matrix is used.
 *      dim[N*n x N*n] (n = x.size)
 *
 *      r: cost vector for state trajectory:
 *      x = [x(1)' x(2)' .. x(N)']'
 *      r: size (N*xdim x 1)
 *
 *      Q: input cost matrix for control input::
 *      u = [u(0)' u(1)' .. u(N-1)']'
 *      If empty, identity matrix is used.
 *      dim[N*m x N*m] m = u.size
 *
 *      distance_error_weight: cost weight for |x(N)-xc|_{ord}
 *
 *      ord: norm used for cost function: ord in {1, 2, INFINITY}
*/
typedef struct cost_function{

    gsl_matrix *R;
    gsl_matrix *Q;
    gsl_vector *r;
    double distance_error_weight;

}cost_function;

/**
 * Abstraction of the system
 * contains all the polytopes
 *
 * number_of_region = total number of possible states that can be reached
 * regions: array of region, each containing one or multiple polytopes
 *
 * closed_loop: if 'true' only first input of next N calculated inputs is applied.
 *              All others are discarded.
 * conservative: if 'true' x(0)...x(N-1) are in starting polytope, x(N) is in final polytope
 *               if false x(1)...x(N-1) can be anywhere
 *
 * ord: norm that is used for minimizing cost function in {1, 2, INFINITY}
 *
 * time_horizon: number of next steps that are taken into account in calculation of the path
 */
typedef struct discrete_dynamics{

    int abstract_states_count;
    int number_of_original_regions;
    abstract_state **abstract_states_set;
    abstract_state **original_regions;
    int closed_loop;
    int conservative;
    int ord;
    size_t time_horizon;

}discrete_dynamics;

/**
 * Dynamics of the system
 *
 * x(t+1) = Ax(t) + Bu(t) + E
 *
 * A system dynamics
 * B input dynamics
 * E disturbance
 * U_set, W_set constraints on states and inputs
 * aux_matrices: auxiliary matrices to fasten calculation of next input
 */
typedef struct system_dynamics{

    gsl_matrix *A;
    gsl_matrix *B;
    gsl_matrix *E;
    gsl_vector *K;
    polytope *W_set;
    polytope *U_set;
    auxiliary_matrices *aux_matrices;

}system_dynamics;

/**
 * State plant is in at the beginning of the calculation
 *
 * current_abs_state starting region
 */
typedef struct current_state{

    int current_abs_state;
    gsl_vector *x;

} current_state;

/**
 * @brief "Constructor" Dynamically allocates the memory all auxiliary matrices need
 * @param n s_dyn.A.size2
 * @param p s_dyn.E.size2
 * @param m s_dyn.B.size2
 * @param u_set_size Uset.size1
 * @param N time horizon
 * @return
 */
struct auxiliary_matrices *aux_matrices_alloc(size_t n,
                                              size_t p,
                                              size_t m,
                                              size_t u_set_size,
                                              size_t N);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the auxiliary matrices
 * @param aux_matrices
 */
void aux_matrices_free(auxiliary_matrices *aux_matrices);

/**
 * @brief "Constructor" Dynamically allocates the memory the complete system dynamics need
 * @param n
 * @param m
 * @param p
 * @param w_set_size
 * @param u_set_size
 * @param N
 * @return
 */
struct system_dynamics *system_dynamics_alloc (size_t n,
                                               size_t m,
                                               size_t p,
                                               size_t w_set_size,
                                               size_t u_set_size,
                                               size_t N);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the system dynamics
 * @param system_dynamics
 */
void system_dynamics_free(system_dynamics * system_dynamics);

/**
 * @brief "Constructor" Dynamically allocates the memory for the state of the plant
 * @param n
 * @param abstract_state
 * @return
 */
struct current_state *state_alloc(size_t n, int abstract_state);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the state of the plant
 * @param state
 */
void state_free(current_state *state);

/**
 * @brief "Constructor" Dynamically allocates the memory for cost function matrices
 * @param n
 * @param m
 * @param N
 * @param distance_error_weight
 * @return
 */
struct cost_function *cost_function_alloc(size_t n,size_t m, size_t N, double distance_error_weight);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the cost function matrices
 * @param cost_function
 */
void cost_function_free(cost_function *cost_function);

/**
 * @brief "Constructor" Dynamically allocates the memory for the discrete abstraction of the system
 * @param polytopes_in_region
 * @param polytope_sizes
 * @param hull_sizes
 * @param orig_polytopes_in_region
 * @param orig_polytope_sizes
 * @param orig_hull_sizes
 * @param n
 * @param abstract_states_count
 * @param number_of_original_regions
 * @param closed_loop
 * @param conservative
 * @param ord
 * @param time_horizon
 * @return
 */
struct discrete_dynamics *discrete_dynamics_alloc(int *polytopes_in_region,
                                                  size_t *polytope_sizes,
                                                  size_t *hull_sizes,
                                                  int *orig_polytopes_in_region,
                                                  size_t *orig_polytope_sizes,
                                                  size_t *orig_hull_sizes,
                                                  int *transitions_in_sizes,
                                                  int *transitions_out_sizes,
                                                  size_t n,
                                                  int abstract_states_count,
                                                  int number_of_original_regions,
                                                  int closed_loop,
                                                  int conservative,
                                                  int ord,
                                                  size_t time_horizon);
/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the discrete dynamics
 * @param d_dyn
 */
void discrete_dynamics_free(discrete_dynamics *d_dyn);

/**
 *
 */
typedef struct control_computation_arguments{

    gsl_matrix *u;
    current_state * now;
    discrete_dynamics *d_dyn;
    system_dynamics *s_dyn;
    int target_abs_state;
    cost_function * f_cost;
    size_t current_time_horizon;
    polytope **polytope_list_backup;

}control_computation_arguments;

typedef struct total_safemode_computation_arguments{

    current_state *now;
    gsl_matrix *u;
    polytope *current;
    polytope *safe;
    system_dynamics *s_dyn;
    size_t time_horizon;
    cost_function *f_cost;
    polytope **polytope_list_backup;

}total_safemode_computation_arguments;

typedef struct next_safemode_computation_arguments{

    current_state *now;
    gsl_matrix *u;
    system_dynamics *s_dyn;
    size_t time_horizon;
    cost_function *f_cost;
    polytope **polytope_list_backup;

}next_safemode_computation_arguments;

/**
 * "Constructor" Dynamically allocates the space for the get_input thread
 */
struct control_computation_arguments *cc_arguments_alloc(current_state *now,
                                                         gsl_matrix* u,
                                                         system_dynamics *s_dyn,
                                                         discrete_dynamics *d_dyn,
                                                         cost_function *f_cost,
                                                         size_t current_time_horizon,
                                                         int target_abs_state,
                                                         polytope **polytope_list);

/**
 * "Constructor" Dynamically allocates the space for the arguments of the safemode computation thread
 */
struct total_safemode_computation_arguments *sm_arguments_alloc(current_state *now,
                                                                gsl_matrix * u,
                                                                polytope *current,
                                                                polytope *safe,
                                                                system_dynamics * s_dyn,
                                                                size_t time_horizon,
                                                                cost_function *f_cost,
                                                                polytope **polytope_list_safemode);

/**
 * "Constructor" Dynamically allocates the space for the arguments of the next step towards safemode computation thread
 */
struct next_safemode_computation_arguments *next_sm_arguments_alloc(current_state *now,
                                                                    gsl_matrix * u,
                                                                    system_dynamics * s_dyn,
                                                                    size_t time_horizon,
                                                                    cost_function *f_cost,
                                                                    polytope **polytope_list_safemode);
#endif //CIMPLE_SYSTEM_H

max_stars_repo_name: shaesaert/TuLiPXML

max_issues_repo_name: shaesaert/TuLiPXML

max_forks_repo_name: shaesaert/TuLiPXML


text:
#include <stdio.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_rng.h>

#define NUM 5

int main() {
  int i;
  double data[NUM];
  double mean, variance, largest, smallest;

  gsl_rng* r = gsl_rng_alloc(gsl_rng_default);
  for (i = 0; i < NUM; ++i) data[i] = gsl_rng_uniform(r);

  mean = gsl_stats_mean(data, 1, NUM);
  variance = gsl_stats_variance(data, 1, NUM);
  largest = gsl_stats_max(data, 1, NUM);
  smallest = gsl_stats_min(data, 1, NUM);

  printf("The data set is");
  for (i = 0; i < NUM; ++i) printf(" %6.3f", data[i]);
  printf("\n");

  printf("Mean:     %6.3f\n", mean);
  printf("Variance: %6.3f\n", variance);
  printf("Max value:%6.3f\n", largest);
  printf("Min value:%6.3f\n", smallest);

  return 0;
}

max_stars_repo_name: julnamoo/practice-linux

max_issues_repo_name: julnamoo/practice-linux

max_forks_repo_name: julnamoo/practice-linux


text:
/*
We are using the OpenCV implementation and adding gridding to the
feature detection.

Modified : computeKeyPoints
Added methods : keepBestKeyPoints and subimageToImageCoordinates
*/

/*********************************************************************
* Software License Agreement (BSD License)
*
*  Copyright (c) 2009, Willow Garage, Inc.
*  All rights reserved.
*
*  Redistribution and use in source and binary forms, with or without
*  modification, are permitted provided that the following conditions
*  are met:
*
*   * Redistributions of source code must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*   * Redistributions in binary form must reproduce the above
*     copyright notice, this list of conditions and the following
*     disclaimer in the documentation and/or other materials provided
*     with the distribution.
*   * Neither the name of the Willow Garage nor the names of its
*     contributors may be used to endorse or promote products derived
*     from this software without specific prior written permission.
*
*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
*  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
*  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
*  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
*  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
*  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
*  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
*  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
*  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
*  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
*  POSSIBILITY OF SUCH DAMAGE.
*********************************************************************/

/** Authors: Ethan Rublee, Vincent Rabaud, Gary Bradski */

#pragma once

#include "Utils\thread_memory.h"
#include "ImageData.h"

#include <iterator>

#include <opencv2\core\types.hpp>

#include <gsl\gsl>

namespace UnitTests
{
    class OpenCVModifiedUnitTest;
    class OrbFeatureDetectorUnitTest;
}

class OrbDetector
{
    using ImageData = mage::ImageData<mage::ImageAllocator>;
public:
    OrbDetector(
        unsigned int gaussianKernelSize,
        unsigned int nfeatures,
        float scaleFactor,
        unsigned int nlevels,
        unsigned int patchSize,
        unsigned int fastThreshold,
        bool useOrientation,
        float featureFactorANMS,
        float featureStrengthANMS,
        int strongResponseANMS,
        float minRobustFactor,
        float maxRobustFactor,
        int numCellsX,
        int numCellsY);

    // Compute the ORB_Impl features and descriptors on an image
    void DetectAndCompute(
        mage::thread_memory memory,
        ImageData& imageData,
        const cv::Mat& image);
private:

    void ComputeKeyPoints(const cv::Mat& imagePyramid,
        gsl::span<const cv::Rect> layerInfo,
        gsl::span<const float> layerScale,
        mage::thread_memory memory,
        ImageData& result);

    static void ResizeAndComputeOrbDescriptorsPrerotated(
        const cv::Mat& image,
        gsl::span<const cv::KeyPoint> keypoints,
        gsl::span<mage::ORBDescriptor> descriptors,
        const std::vector<std::vector<signed char>>& scaledPatterns,
        size_t dsize,
        float upScale);    // amount to upscale the source image

    static void ComputeOrbDescriptorsPrerotated(
        const cv::Mat& imagePyramid,
        gsl::span<const cv::Rect> layerInfo,
        gsl::span<const float> layerScale,
        gsl::span<const cv::KeyPoint> keypoints,
        gsl::span<mage::ORBDescriptor> descriptors,
        const signed char* _patternRotated,
        size_t dsize);

    static void ComputeOrbDescriptors(
        const cv::Mat& imagePyramid,
        gsl::span<const cv::Rect> layerInfo,
        gsl::span<const float> layerScale,
        gsl::span<const cv::KeyPoint> keypoints,
        gsl::span<mage::ORBDescriptor> descriptors,
        gsl::span<const cv::Point> _pattern,
        size_t dsize);

    static void ICAngles(
        const cv::Mat& img,
        gsl::span<const cv::Rect> layerinfo,
        gsl::span<const int> u_max,
        int half_k,
        gsl::span<cv::KeyPoint> pts);

    static void HarrisResponses(
        const cv::Mat& img,
        gsl::span<const cv::Rect> layerinfo,
        gsl::span<cv::KeyPoint> pts,
        int blockSize,
        float harris_k);

    static void MakeRandomPattern(int patchSize, cv::Point* pattern, int npoints);

    enum class FeatureType
    {
        TYPE_5_8 = 0, TYPE_7_12 = 1, TYPE_9_16 = 2
    };

    void AdaptiveNonMaximalSuppresion(mage::temp::vector<cv::KeyPoint>& keypoints,
        unsigned int numToKeep, int threshold, mage::thread_memory memory);
    void FAST(
        const cv::Mat& img,
        mage::temp::vector<cv::KeyPoint>& keypoints,
        int threshold,
        bool nonmax_suppression,
        uchar* threshold_tab,
        mage::thread_memory& memory,
        FeatureType type = FeatureType::TYPE_9_16);

    const unsigned int m_gaussianKernelSize;
    const unsigned int m_nfeatures;
    const float m_scaleFactor;
    const unsigned int m_nlevels;
    const unsigned int m_patchSize;
    const unsigned int m_fastThreshold;
    const bool m_useOrientation;
    const float m_featureFactorANMS;
    const float m_featureStrengthANMS;
    const int m_strongResponseANMS;
    const float m_minRobustFactorANMS;
    const float m_maxRobustFactorANMS;
    const int m_numCellsX;
    const int m_numCellsY;

    friend class ::UnitTests::OpenCVModifiedUnitTest;
    friend class ::UnitTests::OrbFeatureDetectorUnitTest;
};

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
/*
 * The MIT License is a permissive free software license, which permits reuse
 * within both open source and proprietary software. The software is licensed
 * as is, and no warranty is given as to fitness for purpose or absence of
 * infringement of third parties' rights, such as patents. Generally use for
 * research activities is allowed regardless of any third party patents, but
 * commercial use may be subject to a separate license.
 *
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2015 Tuomo Raitio
 *
 * 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.
 *
 *
 *
 *
 * <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
 *              GlottHMM Speech Parameter Extractor
 * <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
 *
 * This program reads a speech file and extracts speech
 * parameters using glottal inverse filtering.
 *
 * This program has been written in Aalto University,
 * Department of Signal Processign and Acoustics, Espoo, Finland
 *
 * Author: Tuomo Raitio
 * Acknowledgements: Antti Suni, Paavo Alku, Martti Vainio
 *
 * File Analysis.c
 * Version: 1.1
 *
 */

/***********************************************/
/*                 INCLUDE                     */
/***********************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sndfile.h> 				/* Read and write wav */
#include <gsl/gsl_vector.h>		/* GSL, Vector */
#include <gsl/gsl_matrix.h>		/* GSL, Matrix */
#include <gsl/gsl_fft_real.h>		/* GSL, FFT */
#include <libconfig.h>				/* Configuration file */
#include "AnalysisFunctions.h"

/*******************************************************************/
/*                          MAIN                                   */
/*******************************************************************/

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

	/* Check command line format */
	if(Check_command_line(argc) == EXIT_FAILURE)
		return EXIT_FAILURE;

	/* Read default configuration file and assign parameters */
	struct config_t *conf_def = Read_config(argv[2]);
	if(conf_def == NULL)
		return EXIT_FAILURE;
	PARAM params;
	if(Assign_config_parameters(conf_def,&params,DEF_CONF) == EXIT_FAILURE) {
		config_destroy(conf_def);
		free(conf_def);
		return EXIT_FAILURE;
	}

	/* Read user configuration file and assign parameters */
	if(argc == 4) {
		struct config_t *conf_usr = Read_config(argv[3]);
		if(conf_usr == NULL)
			return EXIT_FAILURE;
		if(Assign_config_parameters(conf_usr,&params,USR_CONF) == EXIT_FAILURE) {
			config_destroy(conf_usr);
			free(conf_usr);
			return EXIT_FAILURE;
		}
	}

	/* Check the validity of the parameters */
	if(Check_parameter_validity(&params) == EXIT_FAILURE)
		return EXIT_FAILURE;

	/* Read soundfile, define sampling frequency, invert signal if requested */
	gsl_vector *signal = Read_soundfile(argv[1], &params);
	if(signal == NULL)
		return EXIT_FAILURE;

	/* Read external F0 file if requested */
	gsl_vector *fundf;
	if(Read_external_f0(&fundf,&params) == EXIT_FAILURE) {
		gsl_vector_free(signal);
		return EXIT_FAILURE;
	}

	/* High-pass filtering */
	if(HighPassFilter(signal,&params) == EXIT_FAILURE) {
		gsl_vector_free(signal);
		return EXIT_FAILURE;
	}

	/*******************************************************************/
	/*                    ALLOCATE MEMORY                              */
	/*******************************************************************/

	/* Allocate memory for variables */
	int i,j,index;
	double f0 = 0;
	gsl_vector *frame,*frame0,*glottal,*gain,*uvgain,*f0_frame,*glottal_f0,*f0_frame0,*glottsig,*glottsig_f0;
	gsl_vector *source_signal,*h1h2,*naq,*ph1h2,*pnaq;
	gsl_matrix *fundf_candidates, *LSF, *LSF2, *bp_gain, *spectral_tilt, *HNR, *waveform, *harmonics,
			*fftmatrix_vt, *fftmatrix_src, *fftmatrix_uv;
	Allocate_variables(&frame,&frame0,&glottal,&gain,&uvgain,&f0_frame,&glottal_f0,&f0_frame0,&glottsig,
			&glottsig_f0,&source_signal,&fundf_candidates,&LSF,&LSF2,&bp_gain,&spectral_tilt,&HNR,
			&waveform,&harmonics,&h1h2,&naq,&fftmatrix_vt,&fftmatrix_src,&fftmatrix_uv,&params);

	/*******************************************************************/
	/*                   ALLOCATE PULSE LIBRARY                        */
	/*******************************************************************/

	/* Allocate memory for pulse library variables */
	gsl_vector *pulse_inds,*pulse_lengths,*pulse_pos,*pgain;
	gsl_matrix *gpulses,*gpulses_rs,*plsf,*ptilt,*pharm,*phnr,*pwaveform;
	Allocate_pulselib_variables(&gpulses,&gpulses_rs,&pulse_inds,&pulse_pos,&pulse_lengths,
			&plsf,&ptilt,&pharm,&phnr,&pwaveform,&pgain,&ph1h2,&pnaq,&params);

	/********************************************************************/
	/*                    EXTRACT PARAMETERS                            */
	/********************************************************************/

	/* Analysis: Start reading signal vector */
	for(index=0; index<params.n_frames; index++) {

		/* Print progress */
		Print_progress(index,params.n_frames,argv[1],params.FS);

		/* Get samples to frames */
		Get_samples_to_frames(signal,frame,frame0,f0_frame,f0_frame0,params.shift,index);

		/* Gain extraction */
		Gain(frame, gain, index, USE_WINDOWING);
		UnvoicedGain(frame, uvgain, index, USE_WINDOWING, params.unvoiced_frame_length);
		BandPassGain(frame, bp_gain, &params, index);

		/* Preliminary inverse filtering for GCI detection (e.g. in GCI-weighted SWLP) */
		if(params.lp_method == LP_METHOD_ID_WLP && params.lp_weighting == LP_WEIGHTING_ID_GCI) {
			InverseFiltering_long(f0_frame, f0_frame0, glottsig_f0, fundf, glottsig, index, &params);
			FundF(glottsig_f0, frame, fundf_candidates, fundf, bp_gain, index, &params);
			InverseFiltering_long(frame, frame0, glottsig, fundf, glottsig, index, &params);
		}

		/* Inverse filtering */
		InverseFiltering(frame, frame0, glottal, LSF, LSF2, spectral_tilt, fundf, glottsig, fftmatrix_vt,
				fftmatrix_src, fftmatrix_uv, index, &params);
		InverseFiltering_long(f0_frame, f0_frame0, glottal_f0, fundf, glottsig_f0, index, &params);

		/* Fundamental frequency estimation. Define a f0 value for the frame (even if unvoiced) */
		FundF(glottal_f0, frame, fundf_candidates, fundf, bp_gain, index, &params);
		f0 = Define_current_f0(fundf,f0,index);

		/* Estimate Harmonic-to-Noise Ratio (HNR) and magnitudes of the first N harmonics */
		Harmonic_analysis(glottal_f0, harmonics, HNR, h1h2, f0, params.FS, index);

		/* Extract pulses, pitch-synchronous spectrum estimation */
		Extract_pulses(f0_frame,glottal_f0,fundf,naq,gpulses,gpulses_rs,pulse_pos,pulse_inds,pulse_lengths,
				plsf,ptilt,pharm,phnr,pgain,ph1h2,pnaq,LSF,spectral_tilt,harmonics,h1h2,HNR,gain,pwaveform,waveform,index,&params);

		/* Construct source signal */
		Construct_source(source_signal, glottal, glottal_f0, index, &params);
	}

	/********************************************************************/
	/*                POSTPROCESSING OF PARAMETERS                      */
	/********************************************************************/

	/* Add missing frames to the beginning and the end of parameter vectors/matrices.
	 * This is due to the analysis scheme, where the analysis starts
	 * and ends with whole frames instead of zero-padding */
	Add_missing_frames(&fundf,&fundf_candidates,&gain,&uvgain,&HNR,&spectral_tilt,&LSF,&LSF2,&harmonics,&waveform,&naq,&h1h2,
			&fftmatrix_vt,&fftmatrix_src,&fftmatrix_uv,&params);

	/* Postprocess F0 */
	F0_postprocess(fundf,fundf_candidates,&params);

	/* Median filtering */
	MedFilt5_matrix(HNR);
	MedFilt5_matrix(harmonics);
	MedFilt5(h1h2);
	MedFilt5(naq);

	/* Replace unvoiced fetures */
	for(i=0; i<fundf->size; i++) {
		if(gsl_vector_get(fundf, i) == 0) {
			if(params.sep_vuv_spectrum == 0) {
				for(j=0; j<params.lpc_order_vt; j++)
					gsl_matrix_set(LSF, i, j, gsl_matrix_get(LSF2, i, j));
			}
			for(j=0; j<fftmatrix_vt->size2; j++)
				gsl_matrix_set(fftmatrix_vt, i, j, gsl_matrix_get(fftmatrix_uv, i, j));
			gsl_vector_set(gain,i,gsl_vector_get(uvgain,i));
		}
	}

	/* Noise reduction */
	Noise_reduction(gain,&params);

	/* Fill possible gaps in waveform and NAQ */
	Fill_waveform_gaps(fundf,waveform);
	Fill_naq_gaps(fundf,naq);

	/* Select new values to pulse parameters according to refined parameters */
	Select_new_refined_values(fundf,gain,LSF,spectral_tilt,HNR,harmonics,pgain,plsf,ptilt,
			phnr,pharm,pulse_pos,pulse_lengths,waveform,h1h2,ph1h2,naq,pnaq,&params);

	/* Select only unique pulses */
	Select_unique_pulses(&gpulses,&gpulses_rs,&pulse_lengths,&pulse_pos,&pulse_inds,&plsf,
			&ptilt,&phnr,&pharm,&pwaveform,&pgain,&ph1h2,&pnaq,&params);

	/* Select only one pulse per frame */
	Select_one_pulse_per_frame(&gpulses,&gpulses_rs,&pulse_lengths,&pulse_pos,&pulse_inds,&plsf,
				&ptilt,&phnr,&pharm,&pwaveform,&pgain,&ph1h2,&pnaq,&params);

	/********************************************************************/
	/*                      FORMANT ENHANCEMENT                         */
	/********************************************************************/

	/* Formant enhancement by LSFs*/
	if(params.formant_enh_method == FORMANT_ENH_ID_LSF)
		LSF_Postfilter(LSF, &params);

	/* Formant Enhancement by re-estimating LPC and modifying the autocorrelation */
	if(params.formant_enh_method == FORMANT_ENH_ID_LPC)
		LPC_Postfilter(LSF, &params);

	/* Differentiate LSFs if requested */
	if(params.differential_lsf == 1) {
		Differentiate_LSFs(&LSF);
		Differentiate_LSFs(&LSF2);
		Differentiate_LSFs(&spectral_tilt);
	}

	/* Convert F0 to logarithmic scale */
	Convert_F0_to_log(fundf,&params);

	/********************************************************************/
	/*                    WRITE PARAMETERS TO FILE                      */
	/********************************************************************/

	/* Open files for writing parameters, free memory */
	Write_parameters_to_file(argv[1],LSF,LSF2,spectral_tilt,HNR,harmonics,waveform,fundf,gain,h1h2,naq,source_signal,
			fftmatrix_vt,fftmatrix_src,&params);
	Free_variables(frame,frame0,signal,glottal,glottsig,glottsig_f0,uvgain,f0_frame,f0_frame0,glottal_f0,bp_gain,
			fundf_candidates,fftmatrix_uv);

	/**************************************************************/
	/*            WRITE  PULSE LIBRARY TO FILE                    */
	/**************************************************************/

	/* Select unique pulses and write pulse library data to file, free memory */
	Write_pulselibrary_to_file(argv[1],gpulses,gpulses_rs,pulse_lengths,pulse_pos,pulse_inds,plsf,
			ptilt,phnr,pharm,pwaveform,pgain,ph1h2,pnaq,&params);

	/* Finish */
	printf("\nFinished analysis.\n\n");
	return EXIT_SUCCESS;
}

/***********/
/*   EOF   */
/***********/

max_stars_repo_name: mjansche/GlottHMM

max_issues_repo_name: mjansche/GlottHMM

max_forks_repo_name: mjansche/GlottHMM


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "sph_data_types.h"
#include "sph_linked_list.h"
#include "sph_compute.h"

int compute_density_3d_naive(int N,double h,
                             double* restrict x, double* restrict y,
                             double* restrict z,double* restrict nu,
                             double* restrict Fx){
  #pragma omp parallel for
  for(int64_t ii=0;ii<N;ii+=1){
    Fx[ii] = 0;
    for(int64_t jj=0;jj<N;jj+=1){
      double dist = 0.;

      dist += (x[ii]-x[jj])*(x[ii]-x[jj]);
      dist += (y[ii]-y[jj])*(y[ii]-y[jj]);
      dist += (z[ii]-z[jj])*(z[ii]-z[jj]);

      dist = sqrt(dist);

      Fx[ii] += nu[jj]*w_bspline_3d(dist,h);
    }
  }

  return 0;
}

int compute_density_3d_ref(int N,double h,
                           double* restrict x, double* restrict y,
                           double* restrict z, double* restrict nu,
                           double* restrict Fx){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);
  #pragma omp parallel for
  for(int64_t ii=0;ii<N;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];
    double rhoii = 0.0;

    #pragma omp simd reduction(+:rhoii) aligned(x,y,z,nu)
    for(int64_t jj=0;jj<N;jj+=1){
      double q = 0.;

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      //q = sqrt(q);//*inv_h;
      q = sqrt(q)*inv_h;

      rhoii += nu[jj]*w_bspline_3d_simd(q);
    }
    Fx[ii] = kernel_constant*rhoii;
  }

  return 0;
}

int compute_density_3d_tiled(int N,double h,
                             double* restrict x, double* restrict y,
                             double* restrict z, double* restrict nu,
                             double* restrict Fx){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);
  const int64_t STRIP = 1000;
  const int64_t N_tiles = N/STRIP;
  const int64_t N_prime = N - N%STRIP;

  #pragma omp parallel for num_threads(24)
  for(int64_t ii=0;ii<N;ii+=1)
    Fx[ii] = 0.;

  #pragma omp parallel for
  for(int64_t i=0;i<N_prime;i+=STRIP){
    for(int64_t j=0;j<N_prime;j+=STRIP){
      for(int64_t ii=i;ii<i+STRIP;ii+=1){
        double xii = x[ii];
        double yii = y[ii];
        double zii = z[ii];
        double rhoii = 0.0;

        #pragma omp simd reduction(+:rhoii) aligned(x,y,z,nu)
        for(int64_t jj=j;jj<j+STRIP;jj+=1){
          double q = 0.;

          double xij = xii-x[jj];
          double yij = yii-y[jj];
          double zij = zii-z[jj];

          q += xij*xij;
          q += yij*yij;
          q += zij*zij;

          q = sqrt(q)*inv_h;

          rhoii += nu[jj]*w_bspline_3d_simd(q);
        }
        Fx[ii] += kernel_constant*rhoii;
      }
    }
  }

  #pragma omp parallel for
  for(int64_t j=0;j<N_prime;j+=STRIP){
    for(int64_t ii=N_prime;ii<N;ii+=1){
      double xii = x[ii];
      double yii = y[ii];
      double zii = z[ii];
      double rhoii = 0.0;

      #pragma omp simd reduction(+:rhoii) aligned(x,y,z,nu)
      for(int64_t jj=j;jj<j+STRIP;jj+=1){
        double q = 0.;

        double xij = xii-x[jj];
        double yij = yii-y[jj];
        double zij = zii-z[jj];

        q += xij*xij;
        q += yij*yij;
        q += zij*zij;

        q = sqrt(q)*inv_h;

        rhoii += nu[jj]*w_bspline_3d_simd(q);
      }
      Fx[ii] += kernel_constant*rhoii;
    }
  }

  return 0;
}

void swap_int64_t2(int64_t* a, int64_t* b) {
  int64_t t0 = *a;
  int64_t t1 = *(a+1);

  *a     = *b;
  *(a+1) = *(b+1);
  *b     = t0;
  *(b+1) = t1;
}

int64_t partition_int64_t2(int64_t arr[], int64_t low, int64_t high)
{
  int64_t pivot = arr[2*high]; // pivot
  int64_t i     = (low - 1);   // Index of smaller element and indicates the right position of pivot found so far

  for (int64_t j = low; j <= high - 1; j++){
    // If current element is smaller than the pivot
    if (arr[2*j] < pivot){
      i++; // increment index of smaller element
      swap_int64_t2(&arr[2*i], &arr[2*j]);
    }
  }

  swap_int64_t2(&arr[2*(i + 1)], &arr[2*high]);
  return (i + 1);
}

/*
void insertionSort_int64_t2(int64_t *arr, int64_t n){
    // https://www.geeksforgeeks.org/insertion-sort/

    int64_t i=0, key0=0, key1=0, j=0;
    for (i = 1; i < n; i++){
        key0 = arr[2*i+0];
        key1 = arr[2*i+1];
        j = i - 1;

        // Move elements of arr[0..i-1], that are
        // greater than key, to one position ahead
        // of their current position
        while (j >= 0 && arr[2*j] > key0){
            arr[2*(j+1)+0] = arr[2*j+0];
            arr[2*(j+1)+1] = arr[2*j+1];

            j = j-1;
        }

        arr[2*(j+1)+0] = key0;
        arr[2*(j+1)+1] = key1;
    }
}*/

void insertionSort_int64_t2(int64_t *arr, int64_t low, int64_t high)
{
    int64_t i, key0, key1, j;
    for (i = low+1; i < high; i++){
        key0 = arr[2*i+0];
        key1 = arr[2*i+1];

        j = i - 1;

        /* Move elements of arr[0..i-1], that are
        greater than key, to one position ahead
        of their current position */
        while (j >= 0 && arr[2*j] > key0){
            arr[2*(j+1)+0] = arr[2*j+0];
            arr[2*(j+1)+1] = arr[2*j+1];

            j = j-1;
        }
        arr[2*(j+1)+0] = key0;
        arr[2*(j+1)+1] = key1;
    }
}

#define MIN_DEPTH 120

void quickSort_int64_t2(int64_t *arr, int64_t low, int64_t high){
  // https://www.geeksforgeeks.org/quick-sort/

  if (low < high){
    // pi is partitioning index, arr[p] is now
    //   at right place
    int64_t pi = partition_int64_t2(arr, low, high);

    // Separately sort elements before
    // partition and after partition
    //quickSort_int64_t2(arr, low, pi-1);
    //quickSort_int64_t2(arr, pi + 1, high);

    if(pi-low-1 >= MIN_DEPTH){
      #pragma omp task shared(arr)
      quickSort_int64_t2(arr, low, pi-1);
    }
    else{
      //quickSort_int64_t2(arr, low, pi-1);
      //#pragma omp task shared(arr)
      insertionSort_int64_t2(arr,low,pi-1);
    }

    if(high - (pi+1) >= MIN_DEPTH){
      #pragma omp task shared(arr)
      quickSort_int64_t2(arr, pi+1, high);
    }
    else{
      //quickSort_int64_t2(arr, low, pi-1);
      //#pragma omp task shared(arr)
      insertionSort_int64_t2(arr,pi+1,high);
    }
  }
}

int main(){

  int err,dbg=0;
  int64_t N = 1000000;
  double h=0.05;
  linkedListBox *box;
  SPHparticle *lsph;

  err = SPHparticle_SoA_malloc(N,&lsph);
  if(err)
    printf("error in SPHparticle_SoA_malloc\n");

  err = gen_unif_rdn_pos( N,123123123,lsph);
  if(err)
    printf("error in gen_unif_rdn_pos\n");

  box = (linkedListBox*)malloc(1*sizeof(linkedListBox));
  if(box==NULL)
    printf("error alocating the linkedListBox\n");

  box->Xmin = -1.0; box->Ymin = -1.0; box->Zmin = -1.0;
  box->Xmax =  2.0; box->Ymax =  2.0; box->Zmax =  2.0;
  box->Nx = (int)( (box->Xmax-box->Xmin)/(2*h) );
  box->Ny = (int)( (box->Ymax-box->Ymin)/(2*h) );
  box->Nz = (int)( (box->Zmax-box->Zmin)/(2*h) );
  box->N  = (box->Nx)*(box->Ny)*(box->Nz);
  double min_val = fmin((box->Xmax-box->Xmin)/box->Nx,fmin((box->Ymax-box->Ymin)/box->Ny,(box->Zmax-box->Zmin)/box->Nz));
  box->width = (int)( 0.5 + 2*h/min_val );
  box->hbegin = kh_init(0);
  box->hend = kh_init(1);

  double t0,t1,t2,t3,t4,t5,t6,t7;
  t0 = omp_get_wtime();

  if(dbg)
    printf("hello - 3\n");
  err = compute_hash_MC3D(N,lsph,box);

  t1 = omp_get_wtime();

  if(dbg)
    printf("hello - 4\n");

  //qsort(lsph->hash,N,2*sizeof(int64_t),compare_int64_t);
  //quickSort_int64_t2(lsph->hash,0,N-1);

  #pragma omp parallel num_threads(1)
  {
    #pragma omp single
    quickSort_int64_t2(lsph->hash,0,N-1);
  }

  /*
  #pragma omp parallel num_threads(1)
  {
    #pragma omp single
    quicksort_omp(lsph->hash,0,N);
  }*/

  t2 = omp_get_wtime();

  if(dbg)
    printf("hello - 5\n");
  void *swap_arr = malloc(N*sizeof(double));
  err = reorder_lsph_SoA(N,lsph,swap_arr);
  if(err)
    printf("error in reorder_lsph_SoA\n");

  t3 = omp_get_wtime();

  if(dbg)
    printf("hello - 6\n");
  err = setup_interval_hashtables(N,lsph,box);
  if(err)
    printf("error in setup_interval_hashtables\n");

  t4 = omp_get_wtime();

  if(dbg)
    printf("hello - 7\n");

  //err = compute_density_3d(N,h,lsph,box);
  //err = compute_density_3d_load_ballanced(N,h,lsph,box);
  err = compute_density_3d_symmetrical_load_ballance(N,h,lsph,box);
  //err = compute_density_3d_symmetrical_lb_branching(N,h,lsph,box);
  //err = compute_density_3d_symmetrical_lb(N,h,lsph,box);
  if(err)
    printf("error in setup_interval_hashtables\n");

  t5 = omp_get_wtime();

  if(dbg)
    printf("hello - 7\n");

  //err = compute_density_3d(N,h,lsph,box);
  //err = compute_density_3d_innerOmp(N,h,lsph,box);
  //err = compute_density_3d_loopswapped(N,h,lsph,box);
  //
  if(err)
    printf("error in setup_interval_hashtables\n");

  t6 = omp_get_wtime();

  if(dbg)
    printf("hello - 8\n");

  //err = compute_density_3d_tiled(N,h,lsph->x,lsph->y,lsph->z,lsph->nu,lsph->Fx);
  //err = compute_density_3d_ref(N,h,lsph->x,lsph->y,lsph->z,lsph->nu,lsph->Fx);
  //err = compute_density_3d_naive(N,h,lsph->x,lsph->y,lsph->z,lsph->nu,lsph->Fx);
  err = compute_density_3d_load_ballanced(N,h,lsph,box);
  if(err)
    printf("error in compute_density_3d_ref\n");

  t7 = omp_get_wtime();

  printf("compute_hash_MC3D calculation time : %lf : %lf\n",t1-t0,100*(t1-t0)/(t5-t0));
  printf("qsort calculation time : %lf : %lf\n",t2-t1,100*(t2-t1)/(t5-t0));
  printf("reorder_lsph_SoA calculation time : %lf : %lf\n",t3-t2,100*(t3-t2)/(t5-t0));
  printf("setup_interval_hashtables calculation time : %lf : %lf\n",t4-t3,100*(t4-t3)/(t5-t0));
  printf("compute_density_3d base calculation time : %lf : %lf\n",t5-t4,100*(t5-t4)/(t5-t0));
  printf("compute_density_3d load balanced calculation time : %lf : %lf\n",t6-t5,0.);
  printf("Total Linked-List compute_density_3d calculation time : %lf : %lf\n",t5-t0,100*(t5-t0)/(t5-t0));
  printf("Reference tiled compute_density_3d calculation time : %lf\n",t7-t6);

  if(dbg)
    printf("hello - 9\n");

  FILE *fp = fopen("data/sph_density_compute_ref.csv","w");
  if(fp!=NULL){
    for(int64_t i=0;i<N;i+=1)
      fprintf(fp,"%ld %.12lf %.12lf %.6lg\n",i,
                                              lsph->rho[i],
                                              lsph->Fx[i],
                                              fabs(lsph->rho[i]-lsph->Fx[i]));
    fclose(fp);
  }

  if(dbg)
    printf("hello - 10\n");
  SPHparticleSOA_safe_free(N,&lsph);
  safe_free_box(box);
  free(swap_arr);

  return 0;
}

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
/* cdf/binomial.c
 *
 * Copyright (C) 2004 Jason H. Stover.
 * Copyright (C) 2004 Giulio Bottazzi
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf_gamma.h>

#include "error.h"

/* Computes the cumulative distribution function for a binomial
   random variable. For a binomial random variable X with n trials
   and success probability p,

           Pr( X <= k ) = Pr( Y >= p )

   where Y is a beta random variable with parameters k+1 and n-k.

   The binomial distribution has the form,

   prob(k) =  n!/(k!(n-k)!) *  p^k (1-p)^(n-k) for k = 0, 1, ..., n

   The cumulated distributions can be expressed in terms of normalized
   incomplete beta functions (see Abramowitz & Stegun eq. 26.5.26 and
   eq. 6.6.3).

   Reference:

   W. Feller, "An Introduction to Probability and Its
   Applications," volume 1. Wiley, 1968. Exercise 45, page 173,
   chapter 6.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>

double
gsl_cdf_binomial_P (const unsigned int k, const double p, const unsigned int n)
{
  double P;
  double a;
  double b;

  if (p > 1.0 || p < 0.0)
    {
      CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
    }

  if (k >= n)
    {
      P = 1.0;
    }
  else
    {
      a = (double) k + 1.0;
      b = (double) n - k;
      P = gsl_cdf_beta_Q (p, a, b);
    }

  return P;
}

double
gsl_cdf_binomial_Q (const unsigned int k, const double p, const unsigned int n)
{
  double Q;
  double a;
  double b;

  if (p > 1.0 || p < 0.0)
    {
      CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
    }

  if (k >= n)
    {
      Q = 0.0;
    }
  else
    {
      a = (double) k + 1.0;
      b = (double) n - k;
      Q = gsl_cdf_beta_P (p, a, b);
    }

  return Q;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE