text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_INV_H
#define SCILIB_LINALG_INV_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <lapacke.h>
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <experimental/mdspan>
#include <exception>
#include <cassert>
#include <type_traits>

namespace Sci {
namespace Linalg {

// Matrix inversion.
template <class T_a, class T_res, class Layout>
    requires(std::is_same_v<std::remove_cv_t<T_a>, double>)
inline void inv(Sci::Matrix_view<T_a, Layout> a,
                Sci::Matrix_view<T_res, Layout> res)
{
    namespace stdex = std::experimental;

    static_assert(a.is_contiguous());
    static_assert(res.is_contiguous());

    assert(a.extent(0) == a.extent(1));

    if (det(a) == 0.0) {
        throw std::runtime_error("inv: matrix not invertible");
    }
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(0));

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
    }

    Sci::copy(a, res);

    Sci::Vector<BLAS_INT, Layout> ipiv(n);
    Sci::Linalg::lu(res, ipiv.view()); // perform LU factorization

    BLAS_INT info =
        LAPACKE_dgetri(matrix_layout, n, res.data(), lda, ipiv.data());
    if (info != 0) {
        throw std::runtime_error("dgetri: matrix inversion failed");
    }
}

template <class Layout, class Allocator>
inline Sci::Matrix<double, Layout, Allocator>
inv(const Sci::Matrix<double, Layout, Allocator>& a)
{
    Sci::Matrix<double, Layout, Allocator> res(a.extent(0), a.extent(1));
    inv(a.view(), res.view());
    return res;
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_INV_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


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

#include <php.h>
#include <cblas.h>
#include <lapacke.h>
#include "kernel/operators.h"

/**
 * Matrix-matrix multiplication i.e. linear transformation of matrices A and B.
 *
 * @param return_value
 * @param a
 * @param b
 */
void tensor_matmul(zval * return_value, zval * a, zval * b)
{
    unsigned int i, j;
    Bucket * row;
    zval rowC, c;

    zend_array * aa = Z_ARR_P(a);
    zend_array * ab = Z_ARR_P(b);

    Bucket * ba = aa->arData;
    Bucket * bb = ab->arData;

    unsigned int m = zend_array_count(aa);
    unsigned int p = zend_array_count(ab);
    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);
}

/**
 * Dot product between vectors A and B.
 *
 * @param return_value
 * @param a
 * @param b
 */
void tensor_dot(zval * return_value, zval * a, zval * b)
{
    unsigned int i;

    zend_array * aa = Z_ARR_P(a);
    zend_array * ab = Z_ARR_P(b);

    Bucket * ba = aa->arData;
    Bucket * bb = ab->arData;

    unsigned int n = zend_array_count(aa);

    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);
}

/**
 * Return the multiplicative inverse of a square matrix A.
 *
 * @param return_value
 * @param a
 */
void tensor_inverse(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval rowB, b;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int n = zend_array_count(aa);

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

    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;

    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);
    efree(pivots);
}

/**
 * Return the (Moore-Penrose) pseudoinverse of a general matrix A.
 *
 * @param return_value
 * @param a
 */
void tensor_pseudoinverse(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval b, rowB;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int m = zend_array_count(aa);
    unsigned int n = zend_array_count(Z_ARR(ba[0].val));
    unsigned int k = MIN(m, n);

    double * va = emalloc(m * n * sizeof(double));
    double * vu = emalloc(m * m * sizeof(double));
    double * vs = emalloc(k * sizeof(double));
    double * vvt = emalloc(n * n * sizeof(double));
    double * vb = emalloc(n * m * sizeof(double));

    for (i = 0; i < m; ++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 = LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'A', m, n, va, n, vs, vu, m, vvt, n);

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

    for (i = 0; i < k; ++i) {
        cblas_dscal(m, 1.0 / vs[i], &vu[i], m);
    }

    cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, n, m, m, 1.0, vvt, n, vu, m, 0.0, vb, m);

    array_init_size(&b, n);

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

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

        add_next_index_zval(&b, &rowB);
    }

    RETVAL_ARR(Z_ARR(b));

    efree(va);
    efree(vu);
    efree(vs);
    efree(vvt);
    efree(vb);
}

/**
 * Return the row echelon form of matrix A.
 *
 * @param return_value
 * @param a
 */
void tensor_ref(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval rowB, b;
    zval tuple;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int m = zend_array_count(aa);
    unsigned int n = zend_array_count(Z_ARR(ba[0].val));

    double * va = emalloc(m * n * sizeof(double));
    int * pivots = emalloc(MIN(m, n) * sizeof(int));

    for (i = 0; i < m; ++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 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, m, n, va, n, pivots);

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

    array_init_size(&b, m);

    long swaps = 0;

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

        for (j = 0; j < i; ++j) {
            add_next_index_double(&rowB, 0.0);
        }

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

        add_next_index_zval(&b, &rowB);

        if (i + 1 != pivots[i]) {
            ++swaps;
        }
    }

    array_init_size(&tuple, 2);

    add_next_index_zval(&tuple, &b);
    add_next_index_long(&tuple, swaps);

    RETVAL_ARR(Z_ARR(tuple));

    efree(va);
    efree(pivots);
}

/**
 * Compute the Cholesky decomposition of matrix A and return the lower triangular matrix.
 *
 * @param return_value
 * @param a
 */
void tensor_cholesky(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval rowB, b;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int n = zend_array_count(aa);

    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 = LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', n, va, n);

    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 <= i; ++j) {
            add_next_index_double(&rowB, va[i * n + j]);
        }

        for (j = i + 1; j < n; ++j) {
            add_next_index_double(&rowB, 0.0);
        }

        add_next_index_zval(&b, &rowB);
    }

    RETVAL_ARR(Z_ARR(b));

    efree(va);
}

/**
 * Compute the LU factorization of matrix A and return a tuple with lower, upper, and permutation matrices.
 *
 * @param return_value
 * @param a
 */
void tensor_lu(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval rowL, l, rowU, u, rowP, p;
    zval tuple;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int n = zend_array_count(aa);

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

    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 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, va, n, pivots);

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

    array_init_size(&l, n);
    array_init_size(&u, n);
    array_init_size(&p, n);

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

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

        add_next_index_double(&rowL, 1.0);

        for (j = i + 1; j < n; ++j) {
            add_next_index_double(&rowL, 0.0);
        }

        add_next_index_zval(&l, &rowL);
    }

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

        for (j = 0; j < i; ++j) {
            add_next_index_double(&rowU, 0.0);
        }

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

        add_next_index_zval(&u, &rowU);
    }

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

        for (j = 0; j < n; ++j) {
            if (j == pivots[i] - 1) {
                add_next_index_long(&rowP, 1);
            } else {
                add_next_index_long(&rowP, 0);
            }
        }

        add_next_index_zval(&p, &rowP);
    }

    array_init_size(&tuple, 3);

    add_next_index_zval(&tuple, &l);
    add_next_index_zval(&tuple, &u);
    add_next_index_zval(&tuple, &p);

    RETVAL_ARR(Z_ARR(tuple));

    efree(va);
    efree(pivots);
}

/**
 * Compute the eigendecomposition of a general matrix A and return the eigenvalues and eigenvectors in a tuple.
 *
 * @param return_value
 * @param a
 */
void tensor_eig(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval eigenvalues;
    zval eigenvectors;
    zval eigenvector;
    zval tuple;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int n = zend_array_count(aa);

    double * va = emalloc(n * n * sizeof(double));
    double * wr = emalloc(n * sizeof(double));
    double * wi = emalloc(n * sizeof(double));
    double * vr = 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 = LAPACKE_dgeev(LAPACK_ROW_MAJOR, 'N', 'V', n, va, n, wr, wi, NULL, n, vr, n);

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

    array_init_size(&eigenvalues, n);
    array_init_size(&eigenvectors, n);

    for (i = 0; i < n; ++i) {
        add_next_index_double(&eigenvalues, wr[i]);

        array_init_size(&eigenvector, n);

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

        add_next_index_zval(&eigenvectors, &eigenvector);
    }

    array_init_size(&tuple, 2);

    add_next_index_zval(&tuple, &eigenvalues);
    add_next_index_zval(&tuple, &eigenvectors);

    RETVAL_ARR(Z_ARR(tuple));

    efree(va);
    efree(wr);
    efree(wi);
    efree(vr);
}

/**
 * Compute the eigendecomposition of a symmetric matrix A and return the eigenvalues and eigenvectors in a tuple.
 *
 * @param return_value
 * @param a
 */
void tensor_eig_symmetric(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval eigenvalues;
    zval eigenvectors;
    zval eigenvector;
    zval tuple;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int n = zend_array_count(aa);

    double * va = emalloc(n * n * sizeof(double));
    double * wr = emalloc(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 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', n, va, n, wr);

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

    array_init_size(&eigenvalues, n);
    array_init_size(&eigenvectors, n);

    for (i = 0; i < n; ++i) {
        add_next_index_double(&eigenvalues, wr[i]);

        array_init_size(&eigenvector, n);

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

        add_next_index_zval(&eigenvectors, &eigenvector);
    }

    array_init_size(&tuple, 2);

    add_next_index_zval(&tuple, &eigenvalues);
    add_next_index_zval(&tuple, &eigenvectors);

    RETVAL_ARR(Z_ARR(tuple));

    efree(va);
    efree(wr);
}

/**
 * Compute the singular value decomposition of a matrix A and return the singular values, and unitary matrices U and VT in a tuple.
 *
 * @param return_value
 * @param a
 */
void tensor_svd(zval * return_value, zval * a)
{
    unsigned int i, j;
    Bucket * row;
    zval u, rowU;
    zval s;
    zval vt, rowVt;
    zval tuple;

    zend_array * aa = Z_ARR_P(a);

    Bucket * ba = aa->arData;

    unsigned int m = zend_array_count(aa);
    unsigned int n = zend_array_count(Z_ARR(ba[0].val));
    unsigned int k = MIN(m, n);

    double * va = emalloc(m * n * sizeof(double));
    double * vu = emalloc(m * m * sizeof(double));
    double * vs = emalloc(k * sizeof(double));
    double * vvt = emalloc(n * n * sizeof(double));

    for (i = 0; i < m; ++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 = LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'A', m, n, va, n, vs, vu, m, vvt, n);

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

    array_init_size(&u, m);
    array_init_size(&s, k);
    array_init_size(&vt, n);

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

        for (j = 0; j < m; ++j) {
            add_next_index_double(&rowU, vu[i * m + j]);
        }

        add_next_index_zval(&u, &rowU);
    }

    for (i = 0; i < k; ++i) {
        add_next_index_double(&s, vs[i]);
    }

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

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

        add_next_index_zval(&vt, &rowVt);
    }

    array_init_size(&tuple, 3);

    add_next_index_zval(&tuple, &u);
    add_next_index_zval(&tuple, &s);
    add_next_index_zval(&tuple, &vt);

    RETVAL_ARR(Z_ARR(tuple));

    efree(va);
    efree(vu);
    efree(vs);
    efree(vvt);
}

ext: c


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef thdpool_3a298e4f_792c_44ff_bc62_c639d1eb697c_h
#define thdpool_3a298e4f_792c_44ff_bc62_c639d1eb697c_h

#include <thread>
#include <mutex>
#include <condition_variable>
#include <future>
#include <functional>
#include <gslib/std.h>

__gslib_begin__

using std::thread;
using std::function;

class thread_pool
{
public:
    thread_pool(int threads)
    {
        for(int i = 0; i < threads; i ++) {
            _workers.emplace_back(
                [this]() {
                    for(;;) {
                        function<void()> task;
                        {
                            std::unique_lock<std::mutex> lock(_mtx_worker);
                            _cv_worker.wait(lock, [this] { return _stop || !_tasks.empty(); });
                            if(_stop && _tasks.empty())
                                return;
                            task = std::move(_tasks.front());
                            _tasks.pop();
                        }
                        task();
                        {
                            std::unique_lock<std::mutex> lock(_mtx_main);
                            -- _undone;
                        }
                        _cv_main.notify_one();
                    }
                }
            );
        }
    }
    ~thread_pool()
    {
        {
            std::unique_lock<std::mutex> lock(_mtx_worker);
            _stop = true;
        }
        _cv_worker.notify_all();
        for(thread& worker : _workers)
            worker.join();
    }

protected:
    vector<thread>          _workers;
    queue<function<void()>> _tasks;
    std::mutex              _mtx_worker;
    std::mutex              _mtx_main;
    std::condition_variable _cv_worker;
    std::condition_variable _cv_main;
    volatile bool           _stop = false;
    volatile int            _undone = 0;

public:
    template<class _func, class... _args>
    auto add(_func&& f, _args&&... args) ->
        std::future<typename std::result_of<_func(_args...)>::type>
    {
        using return_type = typename std::result_of<_func(_args...)>::type;
        auto task = std::make_shared<std::packaged_task<return_type()>>(
            std::bind(std::forward<_func>(f), std::forward<_args>(args)...)
            );
        std::future<return_type> res = task->get_future();
        {
            std::unique_lock<std::mutex> lock(_mtx_worker);
            if(_stop)
                throw std::runtime_error("enqueue on stopped thread pool.");
            _tasks.emplace([task]() { (*task)(); });
        }
        {
            std::unique_lock<std::mutex> lock(_mtx_main);
            ++ _undone;
        }
        _cv_worker.notify_one();
        return res;
    }
    void join()
    {
        std::unique_lock<std::mutex> lock(_mtx_main);
        _cv_main.wait(lock, [this] { return !_undone; });
    }
};

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/* multifit_nlinear/test.c
 *
 * Copyright (C) 2007, 2013, 2015, 2016 Brian Gough, 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.
 */

/* These tests are based on the NIST Statistical Reference Datasets
   See http://www.nist.gov/itl/div898/strd/index.html for more
   information. */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_ieee_utils.h>

#include "test_fdf.c"

static void
test_proc(const gsl_multifit_nlinear_trs *trs,
          const gsl_multifit_nlinear_scale *scale,
          const gsl_multifit_nlinear_solver *solver,
          const int fdtype)
{
  gsl_multifit_nlinear_parameters fdf_params =
    gsl_multifit_nlinear_default_parameters();

  fdf_params.trs = trs;
  fdf_params.scale = scale;
  fdf_params.solver = solver;
  fdf_params.fdtype = fdtype;

  test_fdf_main(&fdf_params);
}

int
main (void)
{
  const gsl_multifit_nlinear_trs **nlinear_trs[6];
  const gsl_multifit_nlinear_solver **nlinear_solvers[5];
  const gsl_multifit_nlinear_scale **nlinear_scales[3];
  const gsl_multifit_nlinear_trs **trs;
  const gsl_multifit_nlinear_solver **solver;
  const gsl_multifit_nlinear_scale **scale;
  size_t i = 0;

  gsl_ieee_env_setup();

  /* initialize arrays */

  nlinear_trs[0] = &gsl_multifit_nlinear_trs_lm;
  nlinear_trs[1] = &gsl_multifit_nlinear_trs_lmaccel;
  nlinear_trs[2] = &gsl_multifit_nlinear_trs_dogleg;
  nlinear_trs[3] = &gsl_multifit_nlinear_trs_ddogleg;
  nlinear_trs[4] = &gsl_multifit_nlinear_trs_subspace2D;
  nlinear_trs[5] = NULL;

  nlinear_solvers[0] = &gsl_multifit_nlinear_solver_cholesky;
  nlinear_solvers[1] = &gsl_multifit_nlinear_solver_mcholesky;
  nlinear_solvers[2] = &gsl_multifit_nlinear_solver_qr;
  nlinear_solvers[3] = &gsl_multifit_nlinear_solver_svd;
  nlinear_solvers[4] = NULL;

  /* skip Marquardt scaling since it won't pass */
  nlinear_scales[0] = &gsl_multifit_nlinear_scale_levenberg;
  nlinear_scales[1] = &gsl_multifit_nlinear_scale_more;
  nlinear_scales[2] = NULL;

  /* run testsuite over all parameter combinations */

  for (trs = nlinear_trs[i]; trs != NULL; trs = nlinear_trs[++i])
    {
      size_t j = 0;

      fprintf(stderr, "trs = %s\n", (*trs)->name);

      for (solver = nlinear_solvers[j]; solver != NULL; solver = nlinear_solvers[++j])
        {
          size_t k = 0;

          /* don't use Cholesky solver with dogleg methods */
          if (i > 1 && *solver == gsl_multifit_nlinear_solver_cholesky)
            continue;

          fprintf(stderr, "solver = %s\n", (*solver)->name);
          for (scale = nlinear_scales[k]; scale != NULL; scale = nlinear_scales[++k])
            {
              test_proc(*trs, *scale, *solver, GSL_MULTIFIT_NLINEAR_FWDIFF);
              test_proc(*trs, *scale, *solver, GSL_MULTIFIT_NLINEAR_CTRDIFF);
            }
        }
    }

  exit (gsl_test_summary ());
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once

#include <utility>
#include <array>
#include <vector>
#include <unordered_map>
#include <DirectXMathExtend.h>
#include <gsl.h>
#include <VertexTraits.h>
#include <minmax>
#include "bvh.h"
#include <hlslm\hlsl.hpp>
#include <iostream>

namespace Geometrics
{
	template <typename IndexType>
	struct Triangle
	{
		static constexpr size_t VertexCount = 3;
		typedef std::array<IndexType, VertexCount> container_type;
		container_type v;

		typedef IndexType value_type;
		typedef size_t size_type;
		typedef ptrdiff_t difference_type;
		typedef IndexType *pointer;
		typedef const IndexType *const_pointer;
		typedef IndexType& reference;
		typedef const IndexType& const_reference;

		typedef typename container_type::iterator iterator;
		typedef typename container_type::const_iterator const_iterator;

		typedef typename container_type::reverse_iterator reverse_iterator;
		typedef typename container_type::const_reverse_iterator const_reverse_iterator;

		Triangle() {}
		Triangle(const IndexType &v0, const IndexType &v1, const IndexType &v2)
		{
			v[0] = v0; v[1] = v1; v[2] = v2;
		}
		Triangle(const IndexType *v_)
		{
			v[0] = v_[0]; v[1] = v_[1]; v[2] = v_[2];
		}
		template <class S> explicit Triangle(const S &x)
		{
			v[0] = x[0];  v[1] = x[1];  v[2] = x[2];
		}
		IndexType &operator[] (int i) { return v[i]; }
		const IndexType &operator[] (int i) const { return v[i]; }
		operator const IndexType * () const { return &(v[0]); }
		operator IndexType * () { return &(v[0]); }

		void rewind() {
			std::swap(v[0], v[2]);
		}

		int indexof(IndexType v_) const
		{
			return (v[0] == v_) ? 0 :
				(v[1] == v_) ? 1 :
				(v[2] == v_) ? 2 : -1;
		}

		// contrainer access
		auto begin() { return v.begin(); }
		auto end() { return v.end(); }
		auto begin() const { return v.begin(); }
		auto end() const { return v.end(); }

		auto rbegin() { return v.rbegin(); }
		auto rend() { return v.rend(); }
		auto rbegin() const { return v.rbegin(); }
		auto rend() const { return v.rend(); }

		size_t size() const { return VertexCount; }
	};

	struct submesh_data
	{
		uint32_t vertex_offset; // start vertices index of this geometry
		uint32_t vertex_count;
		uint32_t index_offset; // start index position of this geometry
		uint32_t index_count;
	};

	// generate per vertex normal for given triangle mesh
	template <typename VertexType, typename IndexType>
	bool generate_normal(gsl::span<VertexType> vertices, gsl::span<const Triangle<IndexType>> facets)
	{
		using namespace DirectX::VertexTraits;
		using namespace DirectX;

		if (!has_normal<VertexType>::value)
			return false;

		//static_assert(has_normal<VertexType>::value, "The vertex type dose not contains normal field");

		std::vector<XMVECTOR, XMAllocator> normals(vertices.size());

		// set zero
		std::memset(normals.data(), 0, normals.size() * sizeof(XMVECTOR));

		XMVECTOR v0, v1, v2, n;
		for (const auto& face : facets)
		{
			v0 = get_position(vertices[face[0]]);
			v1 = get_position(vertices[face[1]]);
			v2 = get_position(vertices[face[2]]);
			v1 -= v0;
			v2 -= v0;

			v1 = XMVector3Normalize(v1);
			v2 = XMVector3Normalize(v2);
			n = XMVector3Cross(v1, v2); // weighted normal

			normals[face[0]] += n;
			normals[face[1]] += n;
			normals[face[2]] += n;
		}

		for (size_t i = 0; i < vertices.size(); i++)
		{
			n = XMVector3Normalize(normals[i]);
			set_normal(vertices[i], n);
		}

		return true;
	}

	// generate per vertex tangent for given triangle mesh
	template <typename VertexType, typename IndexType>
	bool generate_tangent(gsl::span<VertexType> vertices, gsl::span<const Triangle<IndexType>> facets)
	{
		using namespace DirectX::VertexTraits;
		using namespace DirectX;

		if (!has_tangent<VertexType>::value)
			return false;
		//static_assert(has_tangent<VertexType>::value, "The vertex type dose not contains tangent field");

		std::vector<XMVECTOR, XMAllocator> tan1(vertices.size() * 2);
		XMVECTOR* tan2 = &tan1[vertices.size()];
		std::memset(tan1.data(), 0, tan1.size() * sizeof(XMVECTOR));

		for (const auto& face : facets)
		{
			{
				XMVECTOR v0, v1, v2, w0, w1, w2;
				v0 = get_position(vertices[face[0]]);
				w0 = get_uv(vertices[face[0]]);

				v1 = get_position(vertices[face[1]]);
				w1 = get_uv(vertices[face[1]]);

				v2 = get_position(vertices[face[2]]);
				w2 = get_uv(vertices[face[2]]);
			}
			XMFLOAT4A v_1, v_2, w_1, w_2;

			XMStoreA(v_1, v1 - v0);
			XMStoreA(v_2, v2 - v0);
			XMStoreA(w_1, w1 - w0);
			XMStoreA(w_2, w2 - w0);

			float x1 = v_1.x; float x2 = v_2.x;
			float y1 = v_1.y; float y2 = v_2.y;
			float z1 = v_1.z; float z2 = v_2.z;

			float s1 = w_1.x; float s2 = w_2.x;
			float t1 = w_1.y; float t2 = w_2.y;

			float r = 1.0F / (s1 * t2 - s2 * t1);
			XMVECTOR sdir = XMVectorSet(
				(t2 * x1 - t1 * x2) * r,
				(t2 * y1 - t1 * y2) * r,
				(t2 * z1 - t1 * z2) * r, .0f);
			XMVECTOR tdir = XMVectorSet(
				(s1 * x2 - s2 * x1) * r,
				(s1 * y2 - s2 * y1) * r,
				(s1 * z2 - s2 * z1) * r, .0f);

			tan1[face[0]] += sdir;
			tan1[face[1]] += sdir;
			tan1[face[2]] += sdir;

			tan2[face[0]] += tdir;
			tan2[face[1]] += tdir;
			tan2[face[2]] += tdir;
		}

		for (size_t i = 0; i < vertices.size(); i++)
		{
			XMVECTOR n = get_normal(vertices[i], n);
			XMVECTOR t = tan1[i];

			// Gram-Schmidt orthogonalize
			XMVECTOR nt = XMVector3Normalize(t - n * XMVector3Dot(n, t));

			XMVECTOR w = XMVectorLess(XMVector3Dot(XMVector3Cross(n, t), tan2[i]), XMVectorZero());
			w = XMVectorSelect(g_XMNegativeOne.v, g_XMOne.v, w);

			nt = _DXMEXT XMVectorSelect<0, 0, 0, 1>(nt, w);
			set_tangent(vertices[i], nt);
		}

		return true;
	}

	namespace Detail
	{
		XM_HAS_MEMBER(vertices, has_vertices);
		XM_HAS_MEMBER(indices, has_indices);
	}

	template <typename MeshType>
	struct concept_mesh_type : public std::integral_constant<bool, Detail::has_vertices<MeshType>::value && Detail::has_indices<MeshType>::value> {};

	template <typename _MeshType>
	class SubMesh : public submesh_data
	{
	public:
		typedef _MeshType MeshType; // Parent MeshTypw
		typedef typename MeshType::VertexType	VertexType;
		typedef typename MeshType::IndexType	IndexType;
		typedef typename MeshType::FaceType		FaceType;

		static constexpr size_t PolygonVertex = FaceType::VertexCount;

		SubMesh(MeshType& mesh, const submesh_data& metric)
			: parent(mesh), submesh_data(metric)
		{
			vertices = gsl::span<VertexType>(&mesh.vertices[metric.vertex_offset], metric.vertex_count);
			indices = gsl::span<IndexType>(&mesh.indices[metric.index_offset], metric.index_count);
		}

		MeshType&						parent;
		gsl::span<VertexType>			vertices;
		gsl::span<IndexType>			indices;

		gsl::span<Triangle<IndexType>>	facets()
		{
			return gsl::span<const FaceType>(
				reinterpret_cast<const FaceType*>(indices.data()),
				indices.size() / FaceType::VertexCount);
		}

		bool empty() const {
			return vertices.empty() || indices.empty();
		}

	};

	template <typename _VertexType, typename _IndexType = uint16_t, typename _FaceType = Triangle<_IndexType>>
	class PolygonSoup
	{
	public:
		typedef _VertexType VertexType;
		typedef _IndexType	IndexType;
		typedef _FaceType	FaceType;

		static constexpr size_t PolygonVertex = FaceType::VertexCount;

		PolygonSoup() {}

		bool empty() const {
			return vertices.empty() || indices.empty();
		}

		void clear()
		{
			vertices.clear();
			indices.clear();
		}

		inline const FaceType& facet(int idx) const
		{
			return reinterpret_cast<const FaceType&>(indices[idx * FaceType::VertexCount]);
		}

		inline FaceType& facet(int idx)
		{
			return reinterpret_cast<FaceType&>(indices[idx * FaceType::VertexCount]);
		}

		void add_facet(const FaceType& new_facet)
		{
			indices.insert(indices.end(), new_facet.begin(), new_facet.end());
		}

		template <class... _TIndecies>
		std::enable_if_t<sizeof...(_TIndecies) == FaceType::VertexCount> add_facet(_TIndecies... _indices)
		{
			indices.insert(indices.end(), { static_cast<IndexType>(_indices)... });
		}

		gsl::span<FaceType> facets() {
			return gsl::span<FaceType>(
				reinterpret_cast<FaceType*>(indices.data()),
				indices.size() / FaceType::VertexCount);
		}

		gsl::span<const FaceType> facets() const {
			return gsl::span<const FaceType>(
				reinterpret_cast<const FaceType*>(indices.data()),
				indices.size() / FaceType::VertexCount);
		}

		const VertexType& vertex(int facet, int vidx) const
		{
			return this->vertices[this->indices[facet * FaceType::VertexCount + vidx]];
		}

		VertexType& vertex(int facet, int vidx)
		{
			return this->vertices[this->indices[facet * FaceType::VertexCount + vidx]];
		}

		// flip all the polygons' winding and vertices' normal
		void flip()
		{
			using namespace DirectX::VertexTraits;

			if (has_normal<VertexType>::value)
			{
				for (auto& v : vertices)
					set_normal(v, -get_normal(v));
			}

			static constexpr int vc = FaceType::VertexCount;
			for (int i = 0; i < indices.size() / vc; i++)
				std::reverse(indices.begin() + i * vc, indices.begin() + (i + 1)* vc);
		}

		void generate_normal()
		{
			static_assert(PolygonVertex == 3, "This mesh is not trigle mesh, please trianglize it first");
			generate_normal<VertexType, IndexType>(this->vertices, this->facets());
		}

		void generate_tangent()
		{
			static_assert(PolygonVertex == 3, "This mesh is not trigle mesh, please trianglize it first");
			generate_tangent<VertexType, IndexType>(this->vertices, this->facets());
		}

		// applies an uniform transform to all vertices in the mesh
		void XM_CALLCONV transform(DirectX::FXMMATRIX M)
		{
			using namespace DirectX::VertexTraits;
			for (auto& v : vertices)
			{
				XMVECTOR p = get_position(v);
				p = _DXMEXT XMVector3TransformCoord(p, M);
				set_position(v, p);

				if (has_normal<VertexType>::value)
				{
					p = get_normal(v);
					p = _DXMEXT XMVector3TransformNormal(p, M);
					set_normal(v, p);

					if (has_tangent<VertexType>::value)
					{
						p = get_tangent(v);
						p = _DXMEXT XMVector3TransformNormal(p, M);
						set_tangent(v, p);
					}
				}
			}
		}

		std::vector<VertexType, DirectX::AlignedAllocator<VertexType>> vertices;
		std::vector<IndexType>	indices;
	};

	template <class _VertexType, class _IndexType>
	int XM_CALLCONV intersect(const PolygonSoup<_VertexType, _IndexType, Triangle<_IndexType>> &Mesh, DirectX::FXMVECTOR Origin, DirectX::FXMVECTOR Direction, std::vector<float>* distances = nullptr)
	{
		using namespace DirectX;
		using namespace DirectX::VertexTraits;
		int count = 0;
		XMVECTOR vDir = XMVector3Normalize(Direction);
		for (const auto& tri : Mesh.facets())
		{
			float distance;

			XMVECTOR v0 = get_position(&Mesh.vertices[tri[0]]);
			XMVECTOR v1 = get_position(&Mesh.vertices[tri[1]]);
			XMVECTOR v2 = get_position(&Mesh.vertices[tri[2]]);

			bool hr = DirectX::TriangleTests::Intersects(Origin, vDir, v0, v1, v2, distance);
			if (hr) {
				++count;
				if (distances) {
					distances->push_back(distances);
				}
			}
		}
		return count;
	}

	template <typename _TVertex>
	struct MeshRayIntersectedVertex : public _TVertex
	{
		using Vector3 = DirectX::Vector3;

		float	distance;
		Vector3 barycentric;
		int		facet; // facet id

		inline bool operator < (const MeshRayIntersectedVertex& rhs)
		{
			return this->distance < rhs.distance;
		}

		inline bool is_valiad() const {return  facet >= 0; }
		inline operator bool() const { return is_valiad(); }
	};

	template <size_t _Dim>
	struct PolygonContainmentInfo
	{
		static constexpr size_t Dim = _Dim;
		DirectX::ContainmentType containment; // overall polygon containment type
		DirectX::ContainmentType vertex_containments[Dim];
		operator DirectX::ContainmentType() const { return containment; }
	};

	struct XM_ALIGNATTR aabb_t {
		::hlsl::xmvector3f min;
		::hlsl::xmvector3f max;
		mutable int	 intersected;

		template <size_t dim>
		static bool less(const aabb_t& lhs, const aabb_t & rhs)
		{
			using namespace hlsl;
			return hlsl::less(lhs.min, rhs.min).get<dim>();
		}

		aabb_t operator+(const aabb_t& rhs) const
		{
			const aabb_t& lhs = *this;
			using namespace ::hlsl;
			aabb_t ret{ ::hlsl::min(lhs.min,rhs.min),::hlsl::max(lhs.max,rhs.max) };
			return ret;
		}

		// Convert min-max representation of this box to center-extent representation
		DirectX::BoundingBox get_dxbox() const
		{
			const float epsilon = 1e-3;
			using namespace ::hlsl;
			auto& aabb = *this;
			DirectX::BoundingBox box;
			xmvector3f extent = (aabb.max - aabb.min) * 0.5f;
			extent = ::hlsl::max(extent, xmfloat(epsilon));
			box.Center << (aabb.max + aabb.min) * 0.5f;
			box.Extents << extent;
			return box;
		}

		operator DirectX::BoundingBox() const { return get_dxbox(); }
	};

	// A struct stores fixed size polygon (position only)
	template <size_t _Size>
	struct PlainPolygon {
		hlsl::xmvector3f v[_Size];
	};

	/// <summary>
	/// Basic triangle mesh, each index represent an edge, which is the edge oppsite to the vertex in it's owner triangle
	/// </summary>
	template <typename _VertexType, typename _IndexType = uint16_t>
	class TriangleMesh : public PolygonSoup<_VertexType, _IndexType, Triangle<_IndexType>>
	{
	public:
		static const size_t VertexCount = FaceType::VertexCount;
		using triangle_handle = IndexType;
		using PlainTriangle = PlainPolygon<3>;
		using ContainmentType = DirectX::ContainmentType;
		using BaseType = PolygonSoup<_VertexType, _IndexType, Triangle<_IndexType>>;
		using BaseType::PolygonVertex;
		using typename BaseType::VertexType;
		using typename BaseType::IndexType;
		using IntersectionVertex = MeshRayIntersectedVertex<_VertexType>;

	public:
		// edge's reverse edge
		// stores the adjacent edges of a edge in a triangle
		std::vector<IndexType>	revedges;
		DirectX::BoundingBox	aabb;
		DirectX::BoundingOrientedBox obb;

		PlainTriangle XM_CALLCONV get_triangle_position(IndexType fid) const
		{
			using namespace DirectX::VertexTraits;
			using namespace hlsl;
			auto f = this->indices.begin() + fid * 3;
			PlainTriangle tri = {
			xmvector3f(get_position(this->vertices[f[0]]))
			,xmvector3f(get_position(this->vertices[f[1]]))
			,xmvector3f(get_position(this->vertices[f[2]])) };
			return tri;
		}

#define _EXPANDTRI(tri) tri.v[0].v, tri.v[1].v, tri.v[2].v

		aabb_t XM_CALLCONV get_triangle_aabb(triangle_handle t) const
		{
			using namespace hlsl;
			PlainTriangle tri = get_triangle_position(t);
			xmvector3f vmax = hlsl::max(tri.v[0], tri.v[1]);
			vmax = hlsl::max(vmax, tri.v[2]);
			xmvector3f vmin = hlsl::min(tri.v[0], tri.v[1]);
			vmin = hlsl::min(vmin, tri.v[2]);
			aabb_t ret;
			ret.min = vmin;
			ret.max = vmax;
			ret.intersected = true;
			return ret;
		}

		using triangle_bvh_t = kd_bvh<triangle_handle, float, 3, aabb_t>;
		triangle_bvh_t				triangle_bvh;

		using vertex_bvh_t = kd_bvh<::hlsl::xmvector3f, float, 3, aabb_t>;
		vertex_bvh_t				vertex_bvh;
		// An Lookup table to find the 'Ienditification' of a vertex
		// Make sense when multiple vertex shares the same spatial location but have different index
		std::vector<IndexType>		vertex_id;

		TriangleMesh() : triangle_bvh([this](const triangle_handle& t) {
			return this->get_triangle_aabb(t);
		}) {}

		// vertex's first adjacant edge
		// std::vector<IndexType> vedges;

		union EdgeType
		{
			_IndexType v[2];
			struct {
				_IndexType v0, v1;
			};
		};

		inline EdgeType edge(int facet, int edge) const
		{
			EdgeType e;
			auto& tri = this->facet(facet);
			switch (edge)
			{
			case 0:
				e.v0 = tri[1];
				e.v1 = tri[2];
				break;
			case 1:
				e.v0 = tri[2];
				e.v1 = tri[0];
				break;
			case 2:
				e.v0 = tri[0];
				e.v1 = tri[1];
				break;
			default:
				e.v0 = -1;
				e.v1 = -1;
				break;
			}
			return e;
		}

		inline EdgeType edge(int eid) const
		{
			return edge(eid / VertexCount, eid % VertexCount);
		}

		inline int adjacentFacet(int facet, int edge) const
		{
			return revedges[facet * VertexCount + edge] / VertexCount;
		}

		inline int adjacentFacet(int eid) const
		{
			return revedges[eid] / VertexCount;
		}

		void build_reverse_edges()
		{
			// intialize all adjacant edges to -1
			revedges.resize(this->indices.size(), -1);
			_IndexType vsize = this->vertices.size();
			int esize = this->indices.size();
			int fsize = this->indices.size() / VertexCount;

			using hash_type = uint64_t;
			// make sure int for hash is enough
			assert(vsize*vsize < std::numeric_limits<hash_type>::max());
			std::unordered_map<hash_type, _IndexType> edges(esize * 2);
			// max items inside this table should be less than (esize / 2)

			for (_IndexType fid = 0; fid < fsize; fid++)
			{
				for (_IndexType i = 0; i < VertexCount; i++)
				{
					_IndexType eid = i + fid * VertexCount;
					auto e = this->edge(fid, i);
					e.v0 = vertex_id[e.v0];
					e.v1 = vertex_id[e.v1];

					hash_type ehash = e.v0 * vsize + e.v1;
					hash_type revehash = e.v0 + e.v1 * vsize;

					auto revItr = edges.find(revehash);
					if (revItr == edges.end())
						edges[ehash] = eid;
					else // find reversed edge, remove from edges map
					{
						auto revEid = revItr->second;
						revedges[eid] = revEid;
						revedges[revEid] = eid;
						auto reve = this->edge(revEid);
						reve.v0 = vertex_id[reve.v0];
						reve.v1 = vertex_id[reve.v1];
						assert(reve.v0 == e.v1 && reve.v1 == e.v0);
						edges.erase(revItr);
					}
				}
			}
		}

		// build the adjacent map so we can access all the 1 rings
		void build()
		{
			int vsize = this->vertices.size();
			int esize = this->indices.size();
			int fsize = this->indices.size() / VertexCount;

			// rebuild the triangle-bvh
			triangle_bvh.resize(fsize);
			std::iota(triangle_bvh.begin(), triangle_bvh.end(), 0);
			triangle_bvh.rebuild();

			// build the vertex-bvh
			const aabb_t& bb = this->triangle_bvh.get_volumn(this->triangle_bvh.root());
			::hlsl::xmvector3f voxel_ext(1e-4f);

			vertex_bvh.set_volumn_getter([voxel_ext](const ::hlsl::xmvector3f v) {
				return aabb_t{ v - voxel_ext , v + voxel_ext, 0 };
			});
			vertex_bvh.resize(vsize);
			for (int i = 0; i < vsize; i++)
				vertex_bvh[i] = ::hlsl::xmvector3f(::DirectX::VertexTraits::get_position(this->vertices[i]));
			vertex_bvh.rebuild();

			vertex_id.resize(vsize);
			std::iota(vertex_id.begin(), vertex_id.end(), 0);
			std::unordered_map<uint64_t, IndexType> vids(vsize);

			// Use 20-bits to hash a float component, plus 1 bit signbit
			// Generate a total 63-bits
			auto hash_position = [](const ::hlsl::xmvector3f vp) -> uint64_t{
				auto vf = vp * 3e4f; // 30 * 3e4 = 9e5 < 1e6 = 2^20
				auto vi = vf.cast<uint32_t>();
				uint32_t signbits32 = ::hlsl::detail::move_mask(vi);
				uint64_t signbits = signbits32;
				vi &= ::hlsl::xmuint((1ULL << 20ULL) - 1ULL);
				uint32_t vis32[3];
				vis32 << vi;
				uint64_t vis[3] = { vis32[0], vis32[1], vis32[2]};
				uint64_t vhash = ((signbits&7) << 60ULL) | (vis[0] << 40ULL) | (vis[1] << 20ULL) | (vis[2]);
				return vhash;
			};

			for (int i = 0; i < vsize; i++)
			{
				auto vp = *vertex_bvh.get_object(i);
				std::pair<uint64_t, IndexType> p = { hash_position(vp), i };
				auto rib = vids.insert(p);
				IndexType vid = rib.first->second;
				if (!rib.second && vid != i)
				{
					// query-merge-set merge algorithm
					int root = vid;
					// find the ultimate parent (root)
					while (vertex_id[root] != root) root = vertex_id[root];

					while (vid != root)
					{
						int cvid = vid;
						vid = vertex_id[vid];
						vertex_id[cvid] = root;
					}
					vertex_id[i] = root;
					//if (vid != i)
					//	std::cout << "Glued vertex pair (" << i << ',' << vid << ')' << std::endl;
				}
			}

			build_reverse_edges();

			//for (int fid = 0; fid < this->indices.size() / 3; ++fid)
			//{
			//	auto tri = this->facet(fid);
			//	using namespace DirectX::VertexTraits;
			//	DirectX::XMVECTOR v0 = get_position(this->vertices[tri[0]]);
			//	DirectX::XMVECTOR v1 = get_position(this->vertices[tri[1]]);
			//	DirectX::XMVECTOR v2 = get_position(this->vertices[tri[2]]);
			//	auto&& vol = triangle_bvh.get_volumn(fid);
			//	bool cond = true;
			//	assert(DirectX::XMVector3LessOrEqual(v0, vol.max.v));
			//	assert(DirectX::XMVector3LessOrEqual(v1, vol.max.v));
			//	assert(DirectX::XMVector3LessOrEqual(v2, vol.max.v));
			//	assert(DirectX::XMVector3LessOrEqual(vol.min.v, v0));
			//	assert(DirectX::XMVector3LessOrEqual(vol.min.v, v1));
			//	assert(DirectX::XMVector3LessOrEqual(vol.min.v, v2));
			//}
		}

		// generate a persoude vertex from the interpolation of the trianlge
		VertexType persudo_vertex(int fid, DirectX::FXMVECTOR baycentric) const
		{
			const auto& tri = this->facet(fid);

			const VertexType& v0 = this->vertices[tri[0]];
			const VertexType& v1 = this->vertices[tri[1]];
			const VertexType& v2 = this->vertices[tri[2]];

			using namespace DirectX::VertexTraits;
			return interpolate_vertex(baycentric, v0, v1, v2);
		}

		protected:
			IntersectionVertex XM_CALLCONV get_ray_intersected_vertex(DirectX::FXMVECTOR vOri, DirectX::FXMVECTOR vDir, IndexType fid) const
			{
				using DirectX::XMVECTOR;
				IntersectionVertex info;
				float distance = .0f;
				info.facet = -1;
				PlainTriangle tri = this->get_triangle_position(fid);
				bool hr = DirectX::TriangleTests::Intersects(vOri, vDir, _EXPANDTRI(tri), distance);
				if (!hr) return info; // No intersection, return with info.facet == -1
				// Fill the persudo-vertex
				XMVECTOR pos = distance * vDir + vOri;
				XMVECTOR baryc = DirectX::TriangleTests::BarycentricCoordinate(pos, _EXPANDTRI(tri));
				info.barycentric = baryc;
				info.facet = fid;
				info.distance = distance;
				(VertexType&)info = this->persudo_vertex(fid, baryc);
				return info;
			}

			template <typename _Ty>
			inline static void push_back_output_itr(nullptr_t itr, _Ty&& element)
			{
			}

			template <typename _OutItr, typename _Ty>
			inline static std::enable_if_t<!std::is_same<_OutItr,nullptr_t>::value>
				push_back_output_itr(_OutItr& itr, _Ty&& element)
			{
				static_assert(std::is_same<std::decay_t<decltype(*itr)>,std::decay_t<_Ty>>::value,"Output Iterator's Type is not expected");
				*(itr++) = element;
			}

		public:

		IntersectionVertex XM_CALLCONV first_intersect(DirectX::FXMVECTOR Origin, DirectX::FXMVECTOR Direction) const
		{
			using namespace DirectX;
			assert(this->triangle_bvh.valiad() && "This triangle mesh is modified but not rebuild");

			XMVECTOR vOri = Origin;
			XMVECTOR vDir = XMVector3Normalize(Direction);
			auto ray_aabb_inter = [vOri, vDir](const aabb_t& aabb) -> float
			{
				float f;
				auto box = aabb.get_dxbox();
				bool succ = box.Intersects(vOri, vDir, f);
				aabb.intersected = succ;
				return succ ? f : -1.0f;
			};
			auto ray_triangle_inter = [this, vOri, vDir](const int fid) -> float
			{
				PlainTriangle tri = this->get_triangle_position(fid);
				float f;
				bool succ = DirectX::TriangleTests::Intersects(vOri, vDir, _EXPANDTRI(tri), f);
				return succ ? f : -1.0f;
			};

			//std::cout << counter << std::endl;
			auto pfid = find_first_of(this->triangle_bvh, ray_triangle_inter, ray_aabb_inter);
			IntersectionVertex info;
			info.facet = -1;
			if (pfid)
				info = get_ray_intersected_vertex(vOri, vDir, *pfid);
			return info;
		}

		// Ray intersection test with advanced infomation
		template <typename _TOutItr>
		int XM_CALLCONV all_intersections(DirectX::FXMVECTOR Origin, DirectX::FXMVECTOR Direction, _TOutItr outItr) const
		{
			using namespace DirectX;
			assert(this->triangle_bvh.valiad() && "This triangle mesh is modified but not rebuild");

			size_t count = 0;
			XMVECTOR vOri = Origin;
			XMVECTOR vDir = XMVector3Normalize(Direction);

			struct XM_ALIGNATTR ray_aabb_intersector {
				DirectX::XMVECTOR origin;
				DirectX::XMVECTOR direction;
				bool XM_CALLCONV operator ()(const aabb_t& aabb) const
				{
					float f;
					auto box = aabb.get_dxbox();
					bool succ = box.Intersects(origin, direction, f);
					aabb.intersected = succ;
					return succ;
				}
			} intersector{ vOri , vDir };

			for (auto triangles = find_all_of(triangle_bvh, intersector);
				triangles;
				++triangles)
			{
				int fid = *triangles;
				IntersectionVertex info = get_ray_intersected_vertex(vOri, vDir, fid);
				if (info.is_valiad()) {
					++count;
					push_back_output_itr(outItr, std::move(info));
				}
			}
			return count;
		}

		template <typename _TPred, typename _TContainer>
		void find_adjacant_facets_of(IndexType facet_id, _TPred&& pred, _TContainer& result_map) const
		{
			auto containment = pred(facet_id);
			result_map.insert(std::make_pair(facet_id, containment));
			if (!containment) return;
			for (int e = 0; e < PolygonVertex; e++)
			{
				int eid = facet_id * PolygonVertex + e;
				auto rev_edge = this->revedges[eid];
				auto adj_face = rev_edge / PolygonVertex;
				// if not visited yet
				if (rev_edge != -1 && result_map.find(adj_face) == result_map.end())
					this->find_adjacant_facets_of(adj_face, pred, result_map);
			}
		}

	};

	namespace Internal
	{
#ifndef min
		template <typename T>
		inline T clamp(T value, T minV, T maxV)
		{
			return std::max(std::min(value, maxV), minV);
		}
#else
		template <typename T>
		inline T clamp(T value, T minV, T maxV)
		{
			return max(min(value, maxV), minV);
		}
#endif
	}

	/// <summary>
	/// Compute the closest projection point from P0 to triangle(V0,V1,V2).
	/// </summary>
	/// <param name="P0">The p0.</param>
	/// <param name="V0">The v0.</param>
	/// <param name="V1">The v1.</param>
	/// <param name="V2">The v2.</param>
	/// <returns></returns>
	inline DirectX::XMVECTOR XM_CALLCONV Projection(DirectX::FXMVECTOR P0, DirectX::FXMVECTOR V0, DirectX::FXMVECTOR V1, DirectX::GXMVECTOR V2)
	{
		using namespace DirectX;
		using namespace Geometrics::Internal;
		XMVECTOR edge0 = V1 - V0;
		XMVECTOR edge1 = V2 - V0;
		XMVECTOR p0 = V0 - P0;

		//XMVectorBaryCentric();

		float a = XMVectorGetX(_DXMEXT XMVector3LengthSq(edge0));
		float b = XMVectorGetX(_DXMEXT XMVector3Dot(edge0, edge1));
		float c = XMVectorGetX(_DXMEXT XMVector3LengthSq(edge1));
		float d = XMVectorGetX(_DXMEXT XMVector3Dot(edge0, p0));
		float e = XMVectorGetX(_DXMEXT XMVector3Dot(edge1, p0));

		float det = a*c - b*b;
		float s = b*e - c*d;
		float t = b*d - a*e;

		if (s + t < det)
		{
			if (s < 0.f)
			{
				if (t < 0.f)
				{
					if (d < 0.f)
					{
						s = clamp(-d / a, 0.f, 1.f);
						t = 0.f;
					}
					else
					{
						s = 0.f;
						t = clamp(-e / c, 0.f, 1.f);
					}
				}
				else
				{
					s = 0.f;
					t = clamp(-e / c, 0.f, 1.f);
				}
			}
			else if (t < 0.f)
			{
				s = clamp(-d / a, 0.f, 1.f);
				t = 0.f;
			}
			else
			{
				float invDet = 1.f / det;
				s *= invDet;
				t *= invDet;
			}
		}
		else
		{
			if (s < 0.f)
			{
				float tmp0 = b + d;
				float tmp1 = c + e;
				if (tmp1 > tmp0)
				{
					float numer = tmp1 - tmp0;
					float denom = a - 2 * b + c;
					s = clamp(numer / denom, 0.f, 1.f);
					t = 1 - s;
				}
				else
				{
					t = clamp(-e / c, 0.f, 1.f);
					s = 0.f;
				}
			}
			else if (t < 0.f)
			{
				if (a + d > b + e)
				{
					float numer = c + e - b - d;
					float denom = a - 2 * b + c;
					s = clamp(numer / denom, 0.f, 1.f);
					t = 1 - s;
				}
				else
				{
					s = clamp(-e / c, 0.f, 1.f);
					t = 0.f;
				}
			}
			else
			{
				float numer = c + e - b - d;
				float denom = a - 2 * b + c;
				s = clamp(numer / denom, 0.f, 1.f);
				t = 1.f - s;
			}
		}

		return V0 + s * edge0 + t * edge1;
	}

	inline float XM_CALLCONV Distance(DirectX::FXMVECTOR P0, DirectX::FXMVECTOR V0, DirectX::FXMVECTOR V1, DirectX::GXMVECTOR V2)
	{
		using namespace DirectX;
		XMVECTOR vProj = Projection(P0, V0, V1, V2);
		vProj -= P0;
		return XMVectorGetX(_DXMEXT XMVector3Length(vProj));
	}

	template <typename _VertexType, typename _IndexType>
	inline float XM_CALLCONV Distance(const TriangleMesh<_VertexType, _IndexType> &Mesh, DirectX::FXMVECTOR Point)
	{
		using namespace DirectX;
		float minDis = numeric_limits<float>::max();
		for (const auto& tri : Mesh.facets())
		{
			float dis;

			XMVECTOR v0 = XMLoadFloat3(&Mesh.vertices[tri[0]].position);
			XMVECTOR v1 = XMLoadFloat3(&Mesh.vertices[tri[1]].position);
			XMVECTOR v2 = XMLoadFloat3(&Mesh.vertices[tri[2]].position);

			dis = Distance(Point, v0, v1, v2);
			if (dis < minDis) minDis = dis;
			//minDis = std::min(dis,minDis);
		}

		//XMVECTOR vDis = XMVectorReplicate(numeric_limits<float>::max());
		//for (const auto& vertex : Mesh.vertices)
		//{
		//	XMVECTOR vPos = XMLoadFloat3(&vertex.position);
		//	vPos -= Point;
		//	XMVECTOR dis = XMVector3Length(vPos);
		//	vDis = XMVectorMin(vDis,dis);
		//}

		//float minDisV = XMVectorGetX(vDis);
		//assert(minDisV >= minDis);
		return minDis;
		//return XMVectorGetX(minDis);
	}

	template <typename _VertexType, typename _IndexType>
	bool XM_CALLCONV Inside(const TriangleMesh<_VertexType, _IndexType> &Mesh, DirectX::FXMVECTOR Point)
	{
		XMFLOAT3A Direction;
		Direction.x = (float)std::rand() / (RAND_MAX + 1);
		Direction.y = (float)std::rand() / (RAND_MAX + 1);
		Direction.z = (float)std::rand() / (RAND_MAX + 1);
		XMVECTOR vDir = XMLoadFloat3A(&Direction);
		auto count = intersect(Mesh, Point, vDir, nullptr);
		return count & 1; //count % 2
	}

}

max_stars_repo_name: ArcEarth/SrInspection

max_issues_repo_name: ArcEarth/SrInspection

max_forks_repo_name: ArcEarth/SrInspection


text:
#ifndef LDA_TOPIC_WORD_PROB_MULTI_SELECT_H
#define LDA_TOPIC_WORD_PROB_MULTI_SELECT_H

#include "Lambda.h"
#include "LambdaCreationFunctions.h"
#include "MultiSelectionComp.h"
#include "PDBVector.h"
#include "LDATopicWordProb.h"
#include "TopicAssignment.h"
#include "LDATopicWordProb.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <random>
#include <gsl/gsl_vector.h>

/* This class generates the entries for the topic probability for each word */
using namespace pdb;
class LDATopicWordProbMultiSelection
    : public MultiSelectionComp<LDATopicWordProb, TopicAssignment> {

private:
    Vector<double> prior;
    Handle<Vector<char>> myMem;
    int numTopics;

public:
    ENABLE_DEEP_COPY

    LDATopicWordProbMultiSelection() {}

    LDATopicWordProbMultiSelection(Vector<double>& fromPrior, unsigned numTopics) {

        this->prior = fromPrior;
        this->numTopics = numTopics;

	/* Set up the random number generator */
        gsl_rng* src = gsl_rng_alloc(gsl_rng_mt19937);
        std::random_device rd;
        std::mt19937 gen(rd());
        gsl_rng_set(src, gen());

        int spaceNeeded = sizeof(gsl_rng) + src->type->size;
        myMem = makeObject<Vector<char>>(spaceNeeded, spaceNeeded);

        memcpy(myMem->c_ptr(), src, sizeof(gsl_rng));
        memcpy(myMem->c_ptr() + sizeof(gsl_rng), src->state, src->type->size);

        gsl_rng_free(src);
    }

    Lambda<bool> getSelection(Handle<TopicAssignment> checkMe) override {
        return makeLambda(checkMe, [](Handle<TopicAssignment>& checkMe) { return true; });
    }

    Lambda<Vector<Handle<LDATopicWordProb>>> getProjection(
        Handle<TopicAssignment> checkMe) override {
        return makeLambda(checkMe, [&](Handle<TopicAssignment>& checkMe) {

            gsl_rng* rng = getRng();
            Vector<unsigned>& wordCount = checkMe->getVector();

            unsigned wordNum = prior.size();
            double* mySamples = new double[wordNum];
            double* totalProb = new double[wordNum];
            double* myPrior = prior.c_ptr();
            for (int i = 0; i < wordNum; i++) {
                totalProb[i] = myPrior[i] + wordCount[i];
            }

	    /* Sample the word probability */
            gsl_ran_dirichlet(rng, wordNum, totalProb, mySamples);

	    /* Create the topic probability for each word */
            Vector<Handle<LDATopicWordProb>> result(wordNum);
            for (int i = 0; i < wordNum; i++) {
                Handle<LDATopicWordProb> myTWP =
                    makeObject<LDATopicWordProb>(numTopics, i, checkMe->getKey(), mySamples[i]);
                result.push_back(myTWP);
            }

            delete[] mySamples;
            delete[] totalProb;

            return result;

        });
    }

    gsl_rng* getRng() {
        gsl_rng* dst = (gsl_rng*)myMem->c_ptr();
        dst->state = (void*)(myMem->c_ptr() + sizeof(gsl_rng));
        dst->type = gsl_rng_mt19937;
        return dst;
    }
};

#endif


text:
/* movstat/test_median.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 <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_statistics.h>

/* compute filtered data by explicitely constructing window, sorting it and finding median */
int
slow_movmedian(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y,
               const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      double yi = gsl_stats_median(window, 1, wsize);

      gsl_vector_set(y, i, yi);
    }

  free(window);

  return GSL_SUCCESS;
}

/* test root sequence */
static void
test_median_root(const double tol, const size_t n, const size_t K, const gsl_movstat_end_t etype)
{
  gsl_movstat_workspace *w = gsl_movstat_alloc(K);
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  char buf[2048];
  size_t i;

  /* test a root sequence (square input): x = [zero one zero] */
  gsl_vector_set_all(x, 0.0);

  for (i = n / 3; i <= n / 2; ++i)
    gsl_vector_set(x, i, 1.0);

  /* compute y = median(x) and test y = x */
  gsl_movstat_median(etype, x, y, w);

  sprintf(buf, "n=%zu K=%zu endtype=%u SMF root sequence", n, K, etype);
  compare_vectors(tol, y, x, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_movstat_free(w);
}

static double
func_median(const size_t n, double x[], void * params)
{
  (void) params;
  return gsl_stats_median(x, 1, n);
}

static void
test_median_proc(const double tol, const size_t n, const size_t H, const size_t J,
                 const gsl_movstat_end_t etype, gsl_rng *rng_p)
{
  gsl_movstat_workspace *w;
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *z = gsl_vector_alloc(n);
  gsl_movstat_function F;
  char buf[2048];

  F.function = func_median;
  F.params = NULL;

  if (H == J)
    w = gsl_movstat_alloc(2*H + 1);
  else
    w = gsl_movstat_alloc2(H, J);

  /* test moving median with random input */
  random_vector(x, rng_p);

  /* y = median(x) with slow brute force algorithm */
  slow_movmedian(etype, x, y, H, J);

  /* z = median(x) */
  gsl_movstat_median(etype, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u median random", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = median(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_movstat_median(etype, z, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u median random in-place", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = median(x) with user-defined function */
  gsl_movstat_apply(etype, &F, x, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u median user", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
  gsl_movstat_free(w);
}

static void
test_median(gsl_rng * rng_p)
{
  test_median_root(GSL_DBL_EPSILON, 1000, 3, GSL_MOVSTAT_END_PADZERO);
  test_median_root(GSL_DBL_EPSILON, 200, 15, GSL_MOVSTAT_END_PADVALUE);
  test_median_root(GSL_DBL_EPSILON, 100, 5, GSL_MOVSTAT_END_TRUNCATE);

  test_median_proc(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 1, 1, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 100, 150, 150, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 8, 8, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 0, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 5, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 15, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 10, 15, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 100, 150, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 150, 100, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 100, 100, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_median_proc(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 1, 1, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 100, 150, 150, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 8, 8, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 0, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 5, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 15, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 10, 15, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 100, 150, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 150, 100, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 100, 100, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_median_proc(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 1, 1, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 100, 150, 150, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 8, 8, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 0, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 5, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 15, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 10, 15, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 100, 150, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 150, 100, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 100, 100, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}


text:
/*
 * Copyright (c) 2011-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 *
 * @precisions normal z -> s d c
 *
 */
#ifndef _DPLASMA_CORES_H_
#define _DPLASMA_CORES_H_

#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include "dplasma_complex.h"
#include "parsec/data_dist/matrix/matrix.h"
#include <core_blas.h>

#define PLASMA_BLKADDR(A, type, m, n)  (type *)plasma_getaddr(A, m, n)
#define PLASMA_BLKLDD(A, k) ( ( (k) + (A).i/(A).mb) < (A).lm1 ? (A).mb : (A).lm%(A).mb )

#endif /* _DPLASMA_CORES_H_ */

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_randist.h>

int
main (void)
{
  double x;
  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for (x = 0.1; x < 2; x+= 0.1)
    {
      double y0 = exp (x);
      double sigma = 0.1 * y0;
      double dy = gsl_ran_gaussian (r, sigma);

      printf ("%g %g %g\n", x, y0 + dy, sigma);
    }

  gsl_rng_free(r);

  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 * @copyright (c) 2017 King Abdullah University of Science and Technology (KAUST).
 *                     All rights reserved.
 **/
/**
 * @file time_zpotrf_tile.c
 *
 * This file shows how to generate tile low-rank (TLR) matrix and factorize it using Cholesky factorization.
 *
 * HiCMA is a software package provided by King Abdullah University of Science and Technology (KAUST)
 *
 * @version 0.1.1
 * @author Kadir Akbudak
 * @date 2018-11-08
 **/

/*
 * @copyright (c) 2009-2014 The University of Tennessee and The University
 *                          of Tennessee Research Foundation.
 *                          All rights reserved.
 * @copyright (c) 2012-2016 Inria. All rights reserved.
 * @copyright (c) 2012-2014 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria, Univ. Bordeaux. All rights reserved.
 */

/**
 * The meaning of the descriptors:
 * - AUV: U and V, side by side
 * - AD : U*V
 * - A  : the original, non-approximated problem
 * - Ark: rank of U and V, each tile of the matrix is a single integer in fact a double.
 *
 **/

#include "morse.h"
#include "timing.h"
#include "hicma_constants.h"
#include "hicma_struct.h"
#include "hicma_z.h"
#include <stdio.h>
#include <time.h>
#include <sys/time.h>
//#include <mpi.h> //MPI_Wtime()

#include "starpu.h"
#ifdef MKL
#include <mkl.h>
//#pragma message("MKL is used")
#else
#include <cblas.h>
#ifdef LAPACKE_UTILS
#include <lapacke_utils.h>
#endif
#include <lapacke.h>
//#pragma message("MKL is NOT used")
#endif

#include "starsh-spatial.h"

#include <assert.h>
#include "hicma_z.h"
#include "auxcompute_z.h"
#include "auxdescutil.h"
#include "hicma.h"
#include <math.h>
#include <time.h>

#undef  CBLAS_SADDR
#define CBLAS_SADDR(_val) (_val)

// zgytlr uses starsh in MPI mode.
STARSH_blrf *mpiF;

int print_progress = 1;   // Print progress about the execution
char datebuf[128];
time_t timer;
struct tm* tm_info;
#define PROGRESS(str) \
    if(print_progress){ \
        int myrank = MORSE_My_Mpi_Rank();\
        time(&timer); \
        tm_info = localtime(&timer); \
        strftime(datebuf, 26, "%Y-%m-%d %H:%M:%S",tm_info); \
        fprintf(stderr, "%d:%s\t%d\t%s\t%s\n", myrank, datebuf, __LINE__, __func__, str);\
        fflush(stderr);\
    }
//#undef PROGRESS
//#define PROGRESS(str)

int store_only_diagonal_tiles = 0;
int global_check = 0;
int global_always_fixed_rank = 0;
int global_fixed_rank = 0;
int global_omit_computation = 1;
int num_mpi_ranks;
int run_potrf = 1;
int diag_nrows = 0;
int main_print_index = 0;
int print_index = 0;
int print_index_end = 0;
int main_print_mat = 0;
int print_mat = 0;
int use_scratch = 1; // Use scratch memory provided by starpu
int calc_rank_stat = 1;

void fwrite_array(int m, int n, int ld, double* arr, char* file){
        FILE* fp = fopen(file, "w");
        if(fp == NULL){
            fprintf(stderr, "File %s cannot be opened to write\n", file);
            exit(1);
        }
        int i, j;
        fprintf(fp, "%d %d\n", m, n);
        for(i = 0; i < m; i++){
            for(j = 0; j < n; j++){
                fprintf(fp, "%d\t", (int)arr[ld*j+i] );
            }
            fprintf(fp, "\n" );
        }
        fclose(fp);
}

double timediff(struct timeval begin, struct timeval end){
    double elapsed = (end.tv_sec - begin.tv_sec) +
        ((end.tv_usec - begin.tv_usec)/1000000.0);
    return elapsed;
}

    int
RunTest(int *iparam, double *dparam, morse_time_t *t_, char* rankfile)
{
    // print progress info only on ROOT process
    if(MORSE_My_Mpi_Rank() != 0)
        print_progress = 0;
    PROGRESS("RunTest started");

    // this paramater enables storing only diagonal tiles in a tall and skinny matrix
    store_only_diagonal_tiles = 1;
    //chameleon/runtime/starpu/control/runtime_descriptor.c
    //MORSE_user_tag_size(31,26);
    //MORSE_user_tag_size(31,29);
    MORSE_user_tag_size(31,27);// When I added tile_to_lapack for descArk, I got not enough number of desc error

    // get parameters coming from command line
    PASTE_CODE_IPARAM_LOCALS( iparam );

    // set global variable so that p.. files can fill dense matrix
    global_check = check;
    // calculate total number of mpi processes (it is not used for now)
    num_mpi_ranks = P*Q;
    print_index = iparam[IPARAM_HICMA_PRINTINDEX];
    print_index_end = iparam[IPARAM_HICMA_PRINTINDEXEND];
    print_mat   = iparam[IPARAM_HICMA_PRINTMAT];
    int64_t _nb = iparam[IPARAM_NB];
    LDA = chameleon_max(M, iparam[IPARAM_LDA]);
    int hicma_maxrank  = iparam[IPARAM_HICMA_MAXRANK];
    global_always_fixed_rank = iparam[IPARAM_HICMA_ALWAYS_FIXED_RANK];

    int saveNB = NB;
    NB = MB;
    size_t ncols_AD;
    int saveP = P;
    int saveQ = Q;
    if (store_only_diagonal_tiles == 1) {
        ncols_AD = MB;
    } else {
        ncols_AD = M;
    }
    int saveN = N;
    N = ncols_AD;
    PASTE_CODE_ALLOCATE_MATRIX_TILE( descAD, 1, double, MorseRealDouble, LDA, M, N );
    N = saveN;
    P = saveP;
    Q = saveQ;
    PROGRESS("descAD is allocated");

    size_t ncols_Dense;
    size_t ld_Dense;
    int saveMB = MB;
    if(check == 0) {
        ncols_Dense = MT;
        MB = NB = 1;
        ld_Dense = MT;
    } else {
        ncols_Dense = M;
        ld_Dense = M;
    }
    /*descDense is full matrix if numerical accuracy will be checked.
     * Otherwise it is MB-by-MB matrix with 1-by-1 tiles */
    PASTE_CODE_ALLOCATE_MATRIX_TILE( descDense, 1, double, MorseRealDouble, ld_Dense, ncols_Dense, ncols_Dense );
    if(check == 0) {
        MB = saveMB;
    } else {
    }
    PROGRESS("descDense is allocated");
    NB = saveNB;

    int MTMB = MT * MB; // roundup number of rows/columns for AUV
    int nrows_AUV = MTMB;
    int ld_AUV = MTMB;
    // allocate descUV
    saveN = N;
    N = N * 2;
    saveNB = NB;
    NB = NB * 2;
    //printf("N:%d NB:%d\n", N, NB);
    PASTE_CODE_ALLOCATE_MATRIX_TILE( descAUV, 1, double, MorseRealDouble, ld_AUV, nrows_AUV, N );
    N = saveN;
    NB = saveNB;
    PROGRESS("descAUV is allocated");

    /* tile dimension of rank descriptor must be 1 */
    /* when LD for rk matrices is 1, program exits*/
    int bigMB = MB;
    int bigNB = NB;
    MB = NB = 1;
    PASTE_CODE_ALLOCATE_MATRIX_TILE( descArk, 1, double, MorseRealDouble, MT, MT, NT);
    PROGRESS("descA's are allocated");
    MB = bigMB;
    NB = bigNB;

    int diag_dense = 1;
    int fixedrank = iparam[IPARAM_RK]; //genargs->k
    double fixedacc = pow(10, -1.0*iparam[IPARAM_ACC]);

    char sym;
    if (run_potrf)
        sym = 'S';
    else
        sym = 'N';
    int probtype = iparam[IPARAM_HICMA_STARSH_PROB];
    int maxrank  = iparam[IPARAM_HICMA_STARSH_MAXRANK];
    //double ddecay = pow(10, -1.0*iparam[IPARAM_HICMA_STARSH_DECAY]);
    double ddecay = dparam[IPARAM_HICMA_STARSH_DECAY];

    int  initial_maxrank, final_maxrank;
    double initial_avgrank, final_avgrank;
    HICMA_problem_t hicma_problem;

    hicma_problem.ndim = 2;

    //BEGIN: rndtiled
    if(iparam[IPARAM_HICMA_STARSH_PROB] == HICMA_STARSH_PROB_RND) {
        hicma_problem.noise    = 1.0; //value added to diagonal
    }
    //END:   rndtiled

    //BEGIN: geostat
    //double theta[3] = {1, 0.1, 0.5}; //initially
    double theta[3] = {
        1.0, //sigma
        0.01, //beta
        10.0 //nu Aleks used 10.0 in his paper
    };
    if(iparam[IPARAM_HICMA_STARSH_PROB] == HICMA_STARSH_PROB_GEOSTAT) {
        hicma_problem.theta = theta;
        hicma_problem.noise = 0.0;
        hicma_problem.noise = 1.e-2;
        hicma_problem.kernel_type = STARSH_SPATIAL_MATERN2_SIMD;
    }
    //END: geostat

    //BEGIN: ss
    if(iparam[IPARAM_HICMA_STARSH_PROB] == HICMA_STARSH_PROB_SS) {
        //sigma=1.0 default value line 193 of stars-h/src/applications/spatial.c
        // Correlation length
        hicma_problem.beta  = 0.1;
        //If fixed rank is required set beta=1 and a sample case will be like this nb=25 maxrank=10 m=2500 So ranks will decrease.

        // Smoothing parameter for Matern kernel
        hicma_problem.nu    = 0.5;
        // Shift added to diagonal elements
        hicma_problem.noise = 1.e-4; //not enough for matrices larger than 600K
        hicma_problem.noise = 5.e-4; //works for 640K but does not work for 10M
        hicma_problem.noise = 1.e-2; //
    }
    //END: ss

    //BEGIN: edsin
    if(iparam[IPARAM_HICMA_STARSH_PROB] == HICMA_STARSH_PROB_EDSIN) {
        // Wave number, >= 0
        hicma_problem.wave_k = dparam[IPARAM_HICMA_STARSH_WAVE_K];
        hicma_problem.diag = M;
    //printf("%s %d: %g\n", __FILE__, __LINE__, hicma_problem.wave_k);
    }
    //END: edsin

    PROGRESS("generating coordinates started");
    struct timeval tvalBefore, tvalAfter;  // removed comma
    gettimeofday (&tvalBefore, NULL);
    HICMA_zgenerate_problem(probtype, sym, ddecay, M, MB, MT, NT, &hicma_problem);
    gettimeofday (&tvalAfter, NULL);
    if(MORSE_My_Mpi_Rank()==0){
        printf("Tproblem:%g\n",
                (tvalAfter.tv_sec - tvalBefore.tv_sec)
                 +(tvalAfter.tv_usec - tvalBefore.tv_usec)/1000000.0
              );
        fflush(stderr);
        fflush(stdout);
    }
    PROGRESS("generating coordinates ended");
    mpiF = hicma_problem.starsh_format; // This is assignment will be hidden from user in release

	// DO NOT enforce compression of diagonal tiles
	int compress_diag = 0;
    PROGRESS("nompi zgytlr starting");
    //descDense original problem
    gettimeofday (&tvalBefore, NULL);
    HICMA_zgytlr_Tile(MorseLower, descAUV, descAD, descArk, 0, maxrank, fixedacc, compress_diag, descDense);
    gettimeofday (&tvalAfter, NULL);
    if(MORSE_My_Mpi_Rank()==0){
        printf("Tcompress:%g\n",
                (tvalAfter.tv_sec - tvalBefore.tv_sec)
                 +(tvalAfter.tv_usec - tvalBefore.tv_usec)/1000000.0
                );
        fflush(stderr);
        fflush(stdout);
    }
    PROGRESS("nompi zgytlr finished");
    fflush(stderr);
    fflush(stdout);
    /*return 0; //TODO*/

    if(calc_rank_stat == 1) {
        PASTE_TILE_TO_LAPACK( descArk, Ark_initial, 1, double, MT, NT );
        if(MORSE_My_Mpi_Rank()==0){

            sprintf(rankfile, "%s-1", rankfile);
            fwrite_array(descArk->m, descArk->n, descArk->m, Ark_initial, rankfile);

            HICMA_stat_t hicma_statrk_initial;
            zget_stat(MorseLower, Ark_initial, MT, NT, MT,  &hicma_statrk_initial);
            printf("initial_ranks:");
            zprint_stat(hicma_statrk_initial);
            fflush(stderr);
            fflush(stdout);
        }
    }

    if (global_always_fixed_rank == 1) {
        fprintf(stderr, "%s %d Fixed rank: %d\n", __FILE__, __LINE__, global_fixed_rank);
    }

    if(0 && num_mpi_ranks == 1 && initial_maxrank > N){ //FIXME Enable for distributed mem
        fprintf(stderr,"%s %d %d\t|N:%d is less than actual maxrank:%d\n", __FILE__, __LINE__, MORSE_My_Mpi_Rank(), N, initial_maxrank);
        exit(1);
    }
    int set_diag = 0;

    /* Save A for check */
    PROGRESS("pasting original dense descAD into Adense and Adense2 started");
    // Adense: original dense problem.
    PASTE_TILE_TO_LAPACK( descDense, Adense, check, double, LDA, M );
    double one = 1.0, zero = 0.0, minusone = -1.0, diagVal = M;
    double* swork = NULL;
    //double* cp_L_Adense =  calloc(LDA*M, sizeof(double));
    if(check){
        swork  = calloc(2*M, sizeof(double));
        {size_t i, j;
            double* orgAdense = calloc(LDA*M, sizeof(double));
            for(j = 0; j < M; j++){
                for(i = 0; i < M; i++){
                    orgAdense[j*LDA+i] = Adense[j*LDA+i];
                }
            }
            int info = LAPACKE_dpotrf_work(
                    LAPACK_COL_MAJOR,
                    'L',
                    M, orgAdense, LDA);
            if(0 && info != 0){ //FIXME
                fprintf(stderr, "%s\t|%d\t|Error in LAPACK potrf. info:%d, This errors means "
                        "that the matrix generated is not positive definite\n", __FILE__, __LINE__, info);
            }
            for(j = 0; j < M; j++){
                for(i = 0; i < j; i++){
                    orgAdense[j*LDA+i] = zero;
                }
            }
            /*for(j = 0; j < M; j++) {                */
                /*for(i = 0; i < M; i++){*/
                    /*cp_L_Adense[j*LDA+i] = orgAdense[j*LDA+i];*/
                /*}*/
            /*}*/
            if(main_print_mat ){printf("L of Adense\n");printmat(orgAdense,M,M,LDA,MB, MB);}
            double normOrgAdense = 0.0;
            /*HICMA_znormest(M, M, orgAdense, &normOrgAdense, swork);*/
            /*printf("norm_L_OrgAdense:%e\n",normOrgAdense);*/
            free(orgAdense);
        }
    }
    PASTE_TILE_TO_LAPACK( descDense, Adense2, check, double, LDA, M );
    PROGRESS("pasting original dense descAD into Adense and Adense2 finished");
    PROGRESS("potrf started");
    START_TIMING();
    HICMA_zpotrf_Tile(MorseLower, descAUV, descAD, descArk, fixedrank, maxrank, fixedacc );
    STOP_TIMING();
    fflush(stderr);
    fflush(stdout);
    PROGRESS("potrf finished");
    if(check){
        HICMA_zuncompress(MorseLower, descAUV, descDense, descArk);
        HICMA_zdiag_vec2mat(descAD, descDense);
        PASTE_CODE_FREE_MATRIX( descAD ); //@KADIRLBL001
        descAD = descDense; // descAD was only diagonals.
        // After this line, descAD is dense matrix containing approximate L
        // So no need to adapt below code for descAD containg only diagonals.
    }
    if(calc_rank_stat == 1) {
        PASTE_TILE_TO_LAPACK( descArk, Ark_final, 1, double, MT, NT );
        if(MORSE_My_Mpi_Rank()==0){
            sprintf(rankfile, "%s-2", rankfile);
            fwrite_array(descArk->m, descArk->n, descArk->m, Ark_final, rankfile);
            HICMA_stat_t hicma_statrk_final;
            zget_stat(MorseLower, Ark_final, MT, NT, MT,  &hicma_statrk_final);
            printf("final_ranks:");
            zprint_stat(hicma_statrk_final);
            fflush(stderr);
            fflush(stdout);
        }
    }

    int check_dense = 0;
    int check_app = 1;
    if(check == 0){
        check_dense = check_app = 0;
    }
    if(check_app ) {
        PROGRESS("checking accuracy");
        if( MORSE_My_Mpi_Rank()==0){
#ifndef COMPLEX
            if(main_print_mat){printf("Adense2\n");printmat(Adense2,M,M,LDA,MB, MB);}
            double normA;
            {size_t i, j;
                for(j = 0; j < M; j++){
                    for(i = 0; i < j; i++){
                        Adense2[j*LDA+i] = zero;
                    }
                }
            }
            PROGRESS("normaA started");
            HICMA_znormest(M, M, Adense2, &normA, swork);
            // Ahicma: result of TLR potrf
            PASTE_TILE_TO_LAPACK( descAD, Ahicma, check, double, LDA, M );
            /*if(0){size_t i,j;*/
                /*for(j = 0; j < M; j++) {                */
                    /*for(i = 0; i < M; i++){*/
                        /*Ahicma[j*LDA+i] = cp_L_Adense[j*LDA+i];*/
                    /*}*/
                /*}*/
            /*}*/
            double normAhicma = 0.0;
            {size_t i, j;
                for(j = 0; j < M; j++){
                    for(i = 0; i < j; i++){
                        Ahicma[j*LDA+i] = zero;
                    }
                }
                double* orgAhicma = calloc(LDA*M, sizeof(double));
                for(j = 0; j < M; j++){
                    for(i = 0; i < M; i++){
                        orgAhicma[j*LDA+i] = Ahicma[j*LDA+i];
                    }
                }
                HICMA_znormest(M, M, orgAhicma, &normAhicma, swork);
                free(orgAhicma);
            }
            if(set_diag){size_t j; for(j = 0; j < M; j++){ Ahicma[j*LDA+j] = diagVal; } }
            if(main_print_mat){printf("Ahicma\n");printmat(Ahicma,M,M,LDA, MB, MB);}
            //LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', M, Ahicma, LDA);
            // AhicmaT: transpose of Ahicma
            PROGRESS("copy descAd into AhicmaT started");
            PASTE_TILE_TO_LAPACK( descAD,  AhicmaT, check, double, LDA, M );

            {size_t i, j;
                for(j = 0; j < M; j++){
                    for(i = 0; i < j; i++){
                        Adense[j*LDA+i] = zero;
                    }
                }
            }

            if(main_print_mat){printf("Ahicma-upperzero\n");printmat(Ahicma,M,M,LDA, MB, MB);}
            PROGRESS("Transpose A started");
            LAPACKE_dge_trans(LAPACK_COL_MAJOR, M, M, Ahicma, LDA, AhicmaT, LDA);
            if(main_print_mat){printf("AhicmaT\n");printmat(AhicmaT,M,M,LDA, MB, MB);}
            PROGRESS("TRMM started");
            cblas_dtrmm (CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, M, M, one, Ahicma, LDA, AhicmaT, LDA);
            if(main_print_mat){printf("Ahicma*AhicmaT\n");printmat(AhicmaT,M,M,LDA, MB, MB);}
            //double tmpnorm;normest(M, M, AhicmaT, &tmpnorm, swork);printf("tmpnorm:%e\n",tmpnorm);
            {size_t i, j;
                for(j = 0; j < M; j++){
                    for(i = 0; i < j; i++){
                        AhicmaT[j*LDA+i] = zero;
                    }
                }
            }

            size_t nelm = M * M;
            if(main_print_mat)printf("nelm:%zu M:%d N:%d\n", nelm, M, N);
            PROGRESS("DAXPY started");
            cblas_daxpy(nelm, minusone, AhicmaT, 1, Adense, 1);
            if(main_print_mat){printf("Adense-(Ahicma*AhicmaT)\n");printmat(Adense,M,M,LDA, MB, MB);}

            double normDenseAppDiff;
            PROGRESS("Norm of difference started");
            HICMA_znormest(M, M, Adense, &normDenseAppDiff, swork);
            double accuracyDenseAppDiff = normDenseAppDiff/normA;
            //printf("normA:%.2e normDenseAppdiff:%.2e Accuracy: %.2e\n", normA, normDenseAppDiff,  accuracyDenseAppDiff);
            dparam[IPARAM_RES] = normDenseAppDiff;
            dparam[IPARAM_ANORM] = normA;
            dparam[IPARAM_XNORM] = normA;
            dparam[IPARAM_BNORM] = normAhicma;
#endif
        } else {
            PASTE_TILE_TO_LAPACK( descAD, Ahicma, check, double, LDA, M );
            PASTE_TILE_TO_LAPACK( descAD,  AhicmaT, check, double, LDA, M );
        }
        PROGRESS("checking accuracy is finished");
    }

    PASTE_CODE_FREE_MATRIX( descAUV );
    PROGRESS("descAUV is freed");
    if(check == 0) { // If there is no check, then descAD and descDense are different. Refer to @KADIRLBL001
        PASTE_CODE_FREE_MATRIX( descAD );
        PROGRESS("descAD is freed");
    }
    PASTE_CODE_FREE_MATRIX( descArk );
    PROGRESS("descArk is freed");
    PASTE_CODE_FREE_MATRIX( descDense );
    PROGRESS("descDense is freed");
    PROGRESS("freed descs");
    return 0;
}

max_stars_repo_name: Quansight/hicma

max_issues_repo_name: Quansight/hicma


text:
#ifndef __GSL_VERSION_H__
#define __GSL_VERSION_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif
__BEGIN_DECLS

#define GSL_VERSION "2.4"
#define GSL_MAJOR_VERSION 2
#define GSL_MINOR_VERSION 4

GSL_VAR const char * gsl_version;

__END_DECLS

#endif /* __GSL_VERSION_H__ */

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:
#include <stdio.h>
#include <cblas.h>

int main(void)
{
	float a[] = { 1, 2, 3, 4, 5, 6 };
	float b[] = { 0, 0, 0, 0, 0, 0 };
	const int n = sizeof a / sizeof a[0];
	int i;

	cblas_scopy(n/2, a, 2, b, 2);	/* copy every other element */
	for (i = 0; i < n; ++i)
		printf("%f ", b[i]);
	putchar('\n');

	return 0;
}

max_stars_repo_name: ecashin/mmm-blog

max_issues_repo_name: ecashin/mmm-blog

max_forks_repo_name: ecashin/mmm-blog


text:
/**
 * @license    BSD 3-Clause
 * @copyright  Pawel Okas
 * @version    $Id$
 * @brief
 *
 * @authors    Pawel Okas
 * created on: 30-03-2019
 *
 * @copyright Copyright (c) 2019, Pawel Okas
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 *     1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *     2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the
 *        documentation and/or other materials provided with the distribution.
 *     3. Neither the name of the copyright holder 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 HOLDER
 * 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.
 */

#ifndef _AT24MAC_H_
#define _AT24MAC_H_

#include <charconv>
#include <gsl/span>
#include <optional>
#include <string_view>
#include "I2CDevice/I2CDevice.h"
#include "microhal.h"

/**
 * \addtogroup Devices
 * @{
 * @class AT24MAC
 * @}
 */
class AT24MAC {
 public:
    using I2C = microhal::I2C;
    enum class Error {
        AcknowledgeFailure = static_cast<int>(I2C::Error::AcknowledgeFailure),
        ArbitrationLost = static_cast<int>(I2C::Error::ArbitrationLost),
        Bus = static_cast<int>(I2C::Error::Bus),
        None = static_cast<int>(I2C::Error::None),
        Overrun = static_cast<int>(I2C::Error::Overrun),
        Timeout = static_cast<int>(I2C::Error::Timeout),
        Unknown = static_cast<int>(I2C::Error::Unknown),
        DataOverflow,
        Addres
    };

    struct SerialNumber {
        uint8_t serial[128 / 8];

        bool operator!=(SerialNumber b) { return std::equal(std::begin(serial), std::begin(serial) + sizeof(serial), std::begin(b.serial)); }
        std::string toString() const {
            std::string str;
            str.reserve(18);
            for (uint_fast8_t i = 0; i < sizeof(serial); i++) {
                std::array<char, 2> buff;
                auto [p, ec] = std::to_chars(buff.data(), buff.data() + buff.size(), serial[i], 16);
                str.append(std::string_view(buff.data(), p - str.data()));
                if (i + 1 != sizeof(serial)) {
                    str.append(":");
                }
            }
            return str;
        }
    };

    static constexpr const size_t pageSize = 16;
    static constexpr const size_t memorySizeInBytes = 256;

 private:
    using Endianness = microhal::Endianness;
    using Access = microhal::Access;
    using span = gsl::span<uint8_t>;  // Todo change to std::span when it will be available in gcc

    // create alias to microhal::Address, we just want to type less
    template <typename T, T i>
    using Address = microhal::Address<T, i>;

    struct Register {
        static constexpr auto SerialNumberReg =
            microhal::makeRegister<SerialNumber, Access::ReadOnly, Endianness::Little>(Address<uint8_t, 0b1000'0000>{});
        static constexpr auto EUIAddress = microhal::makeRegister<uint64_t, Access::ReadOnly, Endianness::Big>(Address<uint8_t, 0b1001'1000>{});
    };

    static constexpr uint8_t getMACi2cAddrFromMemoryI2cAddr(uint8_t memoryAddress) { return (memoryAddress & 0x0F) | 0xB0; }

 public:
    AT24MAC(I2C &i2c, uint8_t address) : memory(i2c, address), mac(i2c, getMACi2cAddrFromMemoryI2cAddr(address)) {}

    static std::string_view toString(Error error);

    Error readEUI(uint64_t &eui) { return static_cast<Error>(mac.readRegister(Register::EUIAddress, eui)); }
    Error readSerialNumber(SerialNumber &serial) { return static_cast<Error>(mac.readRegister(Register::SerialNumberReg, serial)); }

    // Memory access functions
    Error readByte(uint8_t address, uint8_t &data) { return static_cast<Error>(memory.read(address, data)); }
    Error read(uint8_t address, span data) { return static_cast<Error>(memory.read(address, data)); }

    Error writeByte(uint8_t address, uint8_t data) { return static_cast<Error>(memory.write(address, data)); }
    Error writePage(uint8_t pageAddress, span data) {
        if (data.size_bytes() > 16) return Error::DataOverflow;
        if ((pageAddress % 16) != 0) return Error::Addres;
        return static_cast<Error>(memory.write(pageAddress, data));
    }
    Error write(uint8_t address, span data) {
        if (data.size_bytes() > memorySizeInBytes) return Error::DataOverflow;

        uint_fast8_t bytesInFirstPage = pageSize - (address % pageSize);
        span firstPage(data.data(), bytesInFirstPage);
        writePage(address, firstPage);

        address += bytesInFirstPage;
        uint_fast8_t bytesToWrite = (data.size_bytes() - bytesInFirstPage);
        uint8_t *dataPtr = data.data() + bytesInFirstPage;
        while (bytesToWrite / pageSize) {
            span page(dataPtr, pageSize);
            writePage(address, page);
            address += pageSize;
            dataPtr += pageSize;
            bytesToWrite -= pageSize;
        }

        uint_fast8_t lastPageSize = bytesToWrite % pageSize;
        if (lastPageSize) {
            span page(dataPtr, lastPageSize);
            writePage(address, page);
        }
        return Error::None;
    }

    void writeWait() { std::this_thread::sleep_for(std::chrono::milliseconds{5}); }

 private:
    microhal::I2CDevice memory;
    microhal::I2CDevice mac;
};

#endif /* _AT24MAC_H_ */

max_stars_repo_name: microHAL/microhal-drivers

max_issues_repo_name: microHAL/microhal-drivers

max_forks_repo_name: microHAL/microhal-drivers


text:
/*
 * PIController.h
 *
 *  Created on: Mar 31, 2013
 *      Author: vgomez, Sep Thijssen
 */

#ifndef PICONTROLLER_H_
#define PICONTROLLER_H_

#define VERB_LEVEL 1

#include "global.h"
#include <boost/property_tree/ptree.hpp>

#include <iostream>
#include <fstream>
#include <algorithm>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <Eigen/Dense>
using namespace Eigen;

class PIController {
public:
	static double dt;		// time increment (required in Model)
	static int units;		// number of UAVs (required in PIController and Model)
	static double R;	    // control (required in Model)
	static double nu;		// noise level
	static double lambda;  	// PI scalar ratio = nu*R
	static int dtperstep;	// number of dt per control step
	static int H;			// number of steps in horizon
	static double dS;		// duration of a control step = dt*dtperstep
	static double stdv;		// local standard deviation = sqrt(nu/dS)
	static int dimUAVx;
	static int dimUAVu;
	static int seed;

protected:
	gsl_rng *r;					// Random seed.
	int dimX;					// dimension state
	int dimU;					// dimension control
	int N;						// number of rollouts applied each step
	double stdv_dt;				// new time scale for noise
	vvec u_exp;					// [H][dimU] current best, used for exploring
	std::string outfile_matlab;	// name of the output .m file
	std::string outfile_costs;	// name of a log file with costs (comma separated values format)
	boost::property_tree::ptree pt;

	// needed for feedback
	vector<MatrixXd> A,lhs,dlhs,rhs,drhs;
	int dimF;

public :

	// nested class Model used by the PIController to simulate dynamics
	class Model {
	protected:
		vec state;			// [dimX] running state
		PIController *ctrl; // pointer to PIController object (outer class)

	public:

		Model() {};
		virtual ~Model() {};

		// apply control uuu on the current state
		virtual void step(const vec&) = 0;

		// replaces the current state
		void setState(const vec& X0) { state = X0; }

		// returns the current state
		vec getState() const { return state; }

		// replaces the current state
		void setController(PIController *_ctrl) { ctrl = _ctrl; }

		// sets the concrete model (implemented by subclass)
		virtual void setProperties() = 0;

		// Output is the immediate state reward
		// writes in of the components of the reward to be logged
		virtual double immediateStateReward(std::ofstream &of) = 0;

		// Output is the immediate state reward
		// sets different attributes of the subclass to be logged if necessary
		virtual double immediateStateReward() = 0;

		// Output is the end state reward. Input is a state.
		double endStateReward() { return immediateStateReward(); };

	};

private:

protected:
	Model *model;			// The dynamical model the PIController uses

	public:
	PIController() {};
	virtual ~PIController() { unsetModel(); };

	void Init(const boost::property_tree::ptree &pt);

	void setSeed(const double &s) { seed = s; };

	// sets the concrete model (implemented by subclass)
	void setModel();

	// unsets the concrete model (implemented by subclass)
	void unsetModel();

	// Returns state based reward of a rollout
	double runningStateReward(const vec&, const vvec&);

	// Returns the cost of a control sequence.
	double runningControlCost(const vvec&) const;

	vec predictState(const vec& X, const vec& A);

	// Output is the immediate control cost. Input is a control action.
	double immediateControlCost(const vec&) const;

//	double immediateStateReward(const vec& X)	{	model->setState(X);	return model->immediateStateReward(); }

	void printTime() const;

	void plotCurrent(const vec&, const vec&) const;

	void plotSetup() const;

	boost::property_tree::ptree getProperties() { return pt; }

	// The PI-control algorithm, input is a state, output is control (VelHeight)
	vec computeControl(const vec &);

	// PI-Controller with feedback, input is a state, output is control (VelHeight)
	// property A is updated
	vec computeControlFeedback(const vec &);

};

#endif /* PICONTROLLER_H_ */


text:
#pragma once
#include "Ds4Color.h"
#include <gsl/span>

/**
 * \brief PoD for setting \c Ds4Device parameters like motor speed, light color, etc.
 */
struct Ds4Output
{
	uint8_t  rightMotor    = 0;
	uint8_t  leftMotor     = 0;
	Ds4Color lightColor    = {};
	uint8_t  flashOnDur    = 0;
	uint8_t  flashOffDur   = 0;
	uint8_t  volumeLeft    = 50;
	uint8_t  volumeRight   = 50;
	uint8_t  volumeMic     = 50;
	uint8_t  volumeSpeaker = 50;

	/**
	 * \brief Updates a buffer with the stored device parameters.
	 * \param buffer Buffer to update.
	 * \return \c true if changes have been made to \a buffer.
	 */
	bool update(const gsl::span<uint8_t>& buffer) const;
};