text:
#include "std_includes.h"

#ifndef MATRIX_H
#define MATRIX_H

#ifdef CRANIUM_USE_MKL
#define CRANIUM_USE_BLAS
#include "mkl.h"
#else
#ifdef CRANIUM_USE_CBLAS
#define CRANIUM_USE_BLAS
#include <cblas.h>
#endif
#endif

// represents user-supplied training data
typedef struct DataSet_ {
    size_t rows;
    size_t cols;
    float** data;
} DataSet;

// represents a matrix of data in row-major order
typedef struct Matrix_ {
    size_t rows;
    size_t cols;
    float* data;
} Matrix;

#ifndef RAND_MOD
#define RAND_MOD 2147483647
#endif

// create dataset given user data
static DataSet* createDataSet(size_t rows, size_t cols, float** data);

// uses memory of the original data to split dataset into batches
static DataSet** createBatches(DataSet* allData, int numBatches);

// split a dataset into row matrices
static Matrix** splitRows(DataSet* dataset);

// destroy dataset
static void destroyDataSet(DataSet* dataset);

// convert dataset to matrix
static Matrix* dataSetToMatrix(DataSet* dataset);

// creates a matrix given data
static Matrix* createMatrix(size_t rows, size_t cols, float* data);

// creates a matrix zeroed out
static Matrix* createMatrixZeroes(size_t rows, size_t cols);

// get an element of a matrix
static float getMatrix(Matrix* mat, size_t row, size_t col);

// set an element of a matrix
static void setMatrix(Matrix* mat, size_t row, size_t col, float val);

// sets the values in $to equal to values in $from
static void copyValuesInto(Matrix* from, Matrix* to);

// prints the entries of a matrix
static void printMatrix(Matrix* input);

// sets each entry in matrix to 0
static void zeroMatrix(Matrix* orig);

// returns transpose of matrix
static Matrix* transpose(Matrix* orig);

// transposes matrix and places data into $origT
static void transposeInto(Matrix* orig, Matrix* origT);

// adds two matrices and returns result
static Matrix* add(Matrix* A, Matrix* b);

// adds $from to $to and places result in $to
static void addTo(Matrix* from, Matrix* to);

// adds $B, a row vector, to each row of $A
static Matrix* addToEachRow(Matrix* A, Matrix* B);

// multiplies every element of $orig by $C
static void scalarMultiply(Matrix* orig, float c);

// multiplies $A and $B (ordering: AB) and returns product matrix
static Matrix* multiply(Matrix* A, Matrix* B);

// multiplies $A and $B (ordering: AB) and places values into $into
static void multiplyInto(Matrix* A, Matrix* B, Matrix* into);

// element-wise multiplcation
static Matrix* hadamard(Matrix* A, Matrix* B);

// places values of hadamard product of $A and $B into $into
static void hadamardInto(Matrix* A, Matrix* B, Matrix* into);

// returns a shallow copy of input matrix
static Matrix* copy(Matrix* orig);

// returns 1 if matrices are equal, 0 otherwise
static int equals(Matrix* A, Matrix* B);

// frees a matrix and its data
static void destroyMatrix(Matrix* matrix);

/*
    Begin functions.
*/

static DataSet* createDataSet(size_t rows, size_t cols, float** data){
    DataSet* dataset = (DataSet*)malloc(sizeof(DataSet));
    dataset->rows = rows;
    dataset->cols = cols;
    dataset->data = data;
    return dataset;
}

DataSet** createBatches(DataSet* allData, int numBatches){
    DataSet** batches = (DataSet**)malloc(sizeof(DataSet*) * numBatches);
    int remainder = allData->rows % numBatches;
    int i;
    int curRow = 0;
    for (i = 0; i < numBatches; i++){
        size_t batchSize = allData->rows / numBatches;
        if (remainder-- > 0){
            batchSize++;
        }
        batches[i] = createDataSet(batchSize, allData->cols, allData->data + curRow);
        curRow += batchSize;
    }
    return batches;
}

static Matrix** splitRows(DataSet* dataset){
    Matrix** rows = (Matrix**)malloc(sizeof(Matrix*) * dataset->rows);
    int i;
    for (i = 0; i < dataset->rows; i++){
        rows[i] = createMatrix(1, dataset->cols, dataset->data[i]);
    }
    return rows;
}

void shuffleTogether(DataSet* A, DataSet* B, std::minstd_rand &gen){
    assert(A->rows == B->rows);
    int i;
    for (i = 0; i < A->rows - 1; i++){
        size_t j = i + gen() / (RAND_MOD / (A->rows - i) + 1);
        float* tmpA = A->data[j];
        A->data[j] = A->data[i];
        A->data[i] = tmpA;
        float* tmpB = B->data[j];
        B->data[j] = B->data[i];
        B->data[i] = tmpB;
    }
}

static void destroyDataSet(DataSet* dataset){
    int i;
    for (i = 0; i < dataset->rows; i++){
        free(dataset->data[i]);
    }
    free(dataset->data);
    free(dataset);
}

static Matrix* dataSetToMatrix(DataSet* dataset){
    Matrix* convert = (Matrix*)malloc(sizeof(Matrix));
    convert->rows = dataset->rows;
    convert->cols = dataset->cols;
    convert->data = (float*)malloc(sizeof(float) * dataset->rows * dataset->cols);
    int i, j;
    for (i = 0; i < dataset->rows; i++){
        for (j = 0; j < dataset->cols; j++){
            setMatrix(convert, i, j, dataset->data[i][j]);
        }
    }
    return convert;
}

Matrix* createMatrix(size_t rows, size_t cols, float* data){
    assert(rows > 0 && cols > 0);
    Matrix* matrix = (Matrix*)malloc(sizeof(Matrix));
    matrix->rows = rows;
    matrix->cols = cols;
    matrix->data = data;
    return matrix;
}

Matrix* createMatrixZeroes(size_t rows, size_t cols){
    assert(rows > 0 && cols > 0);
    Matrix* matrix = (Matrix*)malloc(sizeof(Matrix));
    matrix->rows = rows;
    matrix->cols = cols;
    float* data = (float*)calloc(rows * cols, sizeof(float));
    matrix->data = data;
    return matrix;
}

inline float getMatrix(Matrix* mat, size_t row, size_t col){
    return mat->data[row * mat->cols + col];
}

inline void setMatrix(Matrix* mat, size_t row, size_t col, float val){
    mat->data[row * mat->cols + col] = val;
}

void copyValuesInto(Matrix* from, Matrix* to){
    assert(from->rows == to->rows && from->cols == to->cols);
    memcpy(to->data, from->data, sizeof(float) * to->rows * to->cols);
}

void printMatrix(Matrix* input){
    int i, j;
    for (i = 0; i < input->rows; i++){
        printf("\n");
        for (j = 0; j < input->cols; j++){
            printf("%.2f ", getMatrix(input, i, j));
        }
    }
    printf("\n");
}

void zeroMatrix(Matrix* orig){
    memset(orig->data, 0, orig->rows * orig->cols * sizeof(float));
}

Matrix* transpose(Matrix* orig){
    float* data = (float*)malloc(sizeof(float) * orig->rows * orig->cols);
    Matrix* transpose = createMatrix(orig->cols, orig->rows, data);
    int i, j;
    for (i = 0; i < orig->rows; i++){
        for (j = 0; j < orig->cols; j++){
            setMatrix(transpose, i, j, getMatrix(orig, i, j));
        }
    }
    return transpose;
}

void transposeInto(Matrix* orig, Matrix* origT){
    assert(orig->rows == origT->cols && orig->cols == origT->rows);
    int i, j;
    for (i = 0; i < orig->rows; i++){
        for (j = 0; j < orig->cols; j++){
            setMatrix(origT, j, i, getMatrix(orig, i, j));
        }
    }
}

Matrix* add(Matrix* A, Matrix* B){
    assert(A->rows == B->rows && A->cols == B->cols);
    float* data = (float*)malloc(sizeof(float) * A->rows * B->rows);
    Matrix* result = createMatrix(A->rows, A->cols, data);
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < A->cols; j++){
            setMatrix(result, i, j, getMatrix(B, i, j) + getMatrix(A, i, j));
        }
    }
    return result;
}

void matrix_mul_add(Matrix * to, float a, Matrix * x, float b, bool square = false) {
	assert(to->rows == x->rows && to->cols == x->cols);
	int i, j;
	for (i = 0; i < to->rows; i++) {
		if (!square) {
			for (j = 0; j < to->cols; j++)
				setMatrix(to, i, j, getMatrix(to, i, j) * a + getMatrix(x, i, j) * b);
		}
		else for (j = 0; j < to->cols; j++) {
			float val = getMatrix(x, i, j);
			setMatrix(to, i, j, getMatrix(to, i, j) * a + val * val * b);
		}
	}
}

void matrix_replace(Matrix * to, Matrix * x, float b) {
	assert(to->rows == x->rows && to->cols == x->cols);
	int i, j;
	for (i = 0; i < to->rows; i++) {
		for (j = 0; j < to->cols; j++)
			setMatrix(to, i, j, getMatrix(x, i, j) * b);
	}
}

void adam_update(Matrix * to, Matrix * m, Matrix * v, float a) {
	assert(to->rows == m->rows && to->cols == m->cols);
	assert(to->rows == v->rows && to->cols == v->cols);
	int i, j;
	for (i = 0; i < to->rows; i++) {
		for (j = 0; j < to->cols; j++)
			setMatrix(to, i, j, (getMatrix(m, i, j) * a) / (sqrt(getMatrix(v, i, j)) + 0.000000001));
	}
}

void addTo(Matrix* from, Matrix* to){
    assert(from->rows == to->rows && from->cols == to->cols);
    int i, j;
    for (i = 0; i < from->rows; i++){
        for (j = 0; j < from->cols; j++){
            setMatrix(to, i, j, getMatrix(from, i, j) + getMatrix(to, i, j));
        }
    }
}

// add B to each row of A
Matrix* addToEachRow(Matrix* A, Matrix* B){
    assert(A->cols == B->cols && B->rows == 1);
    float* data = (float*)malloc(sizeof(float) * A->rows * A->cols);
    Matrix* result = createMatrix(A->rows, A->cols, data);
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < A->cols; j++){
            setMatrix(result, i, j, getMatrix(A, i, j) + getMatrix(B, 0, j));
        }
    }
    return result;
}

void scalarMultiply(Matrix* orig, float c){
    int i, j;
    for (i = 0; i < orig->rows; i++){
        for (j = 0; j < orig->cols; j++){
            setMatrix(orig, i, j, getMatrix(orig, i, j) * c);
        }
    }
}

Matrix* multiply(Matrix* A, Matrix* B){
    assert(A->cols == B->rows);
    float* data = (float*)malloc(sizeof(float) * A->rows * B->cols);
    Matrix* result = createMatrix(A->rows, B->cols, data);
#ifdef CRANIUM_USE_BLAS
    zeroMatrix(result);
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, A->rows, B->cols
    , A->cols, 1, A->data, A->cols, B->data, B->cols, 1, result->data, result->cols);
    return result;
#endif
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < B->cols; j++){
            float sum = 0;
            int k;
            for (k = 0; k < B->rows; k++){
                sum += getMatrix(A, i, k) * getMatrix(B, k, j);
            }
            setMatrix(result, i, j, sum);
        }
    }
    return result;
}

void multiplyInto(Matrix* A, Matrix* B, Matrix* into){
#ifdef CRANIUM_USE_BLAS
    zeroMatrix(into);
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, A->rows, B->cols
    , A->cols, 1, A->data, A->cols, B->data, B->cols, 1, into->data, into->cols);
    return;
#endif
    assert(A->cols == B->rows);
    assert(A->rows == into->rows && B->cols == into->cols);
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < B->cols; j++){
            float sum = 0;
            int k;
            for (k = 0; k < B->rows; k++){
                sum += getMatrix(A, i, k) * getMatrix(B, k, j);
            }
            setMatrix(into, i, j, sum);
        }
    }
}

Matrix* hadamard(Matrix* A, Matrix* B){
    assert(A->rows == B->rows && A->cols == B->cols);
    float* data = (float*)malloc(sizeof(float) * A->rows * A->cols);
    Matrix* result = createMatrix(A->rows, A->cols, data);
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < A->cols; j++){
            setMatrix(result, i, j, getMatrix(A, i, j) * getMatrix(B, i, j));
        }
    }
    return result;
}

void hadamardInto(Matrix* A, Matrix* B, Matrix* into){
    assert(A->rows == B->rows && A->cols == B->cols);
    assert(A->rows == into->rows && A->cols == into->cols);
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < A->cols; j++){
            setMatrix(into, i, j, getMatrix(A, i, j) * getMatrix(B, i, j));
        }
    }
}

Matrix* copy(Matrix* orig){
    float* data = (float*)malloc(sizeof(float) * orig->rows * orig->cols);
    memcpy(data, orig->data, sizeof(float) * orig->cols * orig->rows);
    return createMatrix(orig->rows, orig->cols, data);
}

int equals(Matrix* A, Matrix* B){
    if (A->rows != B->rows){
        return 0;
    }
    if (A->cols != B->cols){
        return 0;
    }
    int i, j;
    for (i = 0; i < A->rows; i++){
        for (j = 0; j < A->cols; j++){
            if (getMatrix(A, i, j) != getMatrix(B, i, j)){
                return 0;
            }
        }
    }
    return 1;
}

void destroyMatrix(Matrix* matrix){
    free(matrix->data);
    free(matrix);
}

#endif

max_stars_repo_name: ewail/DiaNN

max_issues_repo_name: ewail/DiaNN

max_forks_repo_name: ewail/DiaNN


text:
#pragma once

#include "BasicTypes.h"
#include "CharInfo.h"
#include "Distance.h"
#include "FontTypes.h"
#include "types/ND0.h"
#include "types/ND1_Node.h"
#include "types/ND1_Box.h"
#include "types/ND2_Char.h"
#include "types/ND2_Box.h"
#include "types/ND2_Rule.h"
#include "types/ND2_Glue.h"
#include "types/ND2_Kern.h"
#include "types/ND2_Penalty.h"
#include <gsl.h>
#include <list>

ext: h


text:
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 Microsoft Corporation. All rights reserved.
//
// This code is licensed under the MIT License (MIT).
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
//
///////////////////////////////////////////////////////////////////////////////

#ifndef BITVEC_H
#define BITVEC_H

#include <gsl/gsl_assert> // for Expects
#include <gsl/gsl_util>   // for narrow_cast, narrow

using gsl::narrow_cast;

#include <cstddef>  // for ptrdiff_t, size_t, nullptr_t
#include <iterator> // for reverse_iterator, distance, random_access_...
#include <memory>

#ifdef _MSC_VER
#pragma warning(push)

// turn off some warnings that are noisy about our Expects statements
#pragma warning(disable : 4127) // conditional expression is constant
#pragma warning(disable : 4702) // unreachable code

// Turn MSVC /analyze rules that generate too much noise. TODO: fix in the tool.
#pragma warning(                                                               \
    disable : 26495) // uninitalized member when constructor calls constructor
#pragma warning(                                                               \
    disable : 26446) // parser bug does not allow attributes on some templates

#endif // _MSC_VER

namespace gb {

template <std::ptrdiff_t Extent> class bitvec {
  public:
	// constants and types
	using element_type = unsafe_bit_t;
	using value_type = std::remove_cv_t<element_type>;
	using index_type = std::ptrdiff_t;
	using pointer = shared_ptr_unsynchronized<element_type>;

	using iterator = details::span_iterator<bitvec<Extent>, false>;
	using const_iterator = details::span_iterator<bitvec<Extent>, true>;
	using reverse_iterator = std::reverse_iterator<iterator>;
	using const_reverse_iterator = std::reverse_iterator<const_iterator>;

	using size_type = index_type;

	// [bitvec.cons], bitvec constructors, copy, assignment, and destructor
	template <bool Dependent = false,
	          // "Dependent" is needed to make "std::enable_if_t<Dependent ||
	          // Extent <= 0>" SFINAE, since "std::enable_if_t<Extent <= 0>" is
	          // ill-formed when Extent is greater than 0.
	          class = std::enable_if_t<(Dependent || Extent <= 0)>>
	bitvec() noexcept : storage_(nullptr, details::extent_type<0>()) {};

	bitvec(pointer ptr, index_type count) : storage_(ptr, count) {}

	bitvec(const bitvec &other) noexcept = default;

	template <std::ptrdiff_t OtherExtent,
	          class = std::enable_if_t<OtherExtent == Extent ||
	                                   Extent == dynamic_extent>>
	bitvec(const bitvec<OtherExtent> &other)
	    : storage_(other.data(),
	               details::extent_type<OtherExtent>(other.size())) {}

	~bitvec() noexcept = default;
	bitvec &operator=(const bitvec &other) noexcept = default;
	inline void operator=(bool value) {
		write(*this, value);
	};

	template <std::ptrdiff_t OtherExtent>
	void operator&=(const bitvec<OtherExtent> &other) noexcept {
		static_assert(Extent == OtherExtent || Extent == dynamic_extent || OtherExtent == dynamic_extent,
			"Mismatching operand sizes");
		_and(*this, *this, other);
	};
	template <std::ptrdiff_t OtherExtent>
	void operator|=(const bitvec<OtherExtent> &other) noexcept {
		static_assert(Extent == OtherExtent || Extent == dynamic_extent || OtherExtent == dynamic_extent,
			"Mismatching operand sizes");
		_or(*this, *this, other);
	};
	template <std::ptrdiff_t OtherExtent>
	void operator^=(const bitvec<OtherExtent> &other) noexcept {
		static_assert(Extent == OtherExtent || Extent == dynamic_extent || OtherExtent == dynamic_extent,
			"Mismatching operand sizes");
		_xor(*this, *this, other);
	};

	// [bitvec.sub], bitvec subviews
	template <std::ptrdiff_t Count> bitvec<Count> first() const {
		Expects(Count >= 0 && Count <= size());
		return {data(), Count};
	}

	template <std::ptrdiff_t Count = 1>
	GSL_SUPPRESS(bounds .1) // NO-FORMAT: attribute
	bitvec<Count> last() const {
		Expects(Count >= 0 && size() - Count >= 0);
		return {data_plus(size() - Count), Count};
	}

	template <std::ptrdiff_t Offset, std::ptrdiff_t Count = dynamic_extent>
	GSL_SUPPRESS(bounds .1) // NO-FORMAT: attribute
	auto subspan() const ->
	    typename details::calculate_subspan_type<Extent, Offset, Count>::type {
		Expects((Offset >= 0 && size() - Offset >= 0) &&
		        (Count == dynamic_extent ||
		         (Count >= 0 && Offset + Count <= size())));

		return {data_plus(Offset),
		        Count == dynamic_extent ? size() - Offset : Count};
	}

	template <std::ptrdiff_t Count>
	bitvec<Count> subspan(std::ptrdiff_t offset) const {
		Expects((offset >= 0 && size() - offset >= 0) &&
		        (Count == dynamic_extent ||
		         (Count >= 0 && offset + Count <= size())));

		return {data_plus(offset),
		        Count == dynamic_extent ? size() - offset : Count};
	}

	bitvec<dynamic_extent> first(index_type count) const {
		Expects(count >= 0 && count <= size());
		return {data(), count};
	}

	bitvec<dynamic_extent> last(index_type count) const {
		return make_subspan(size() - count, dynamic_extent,
		                    subspan_selector<Extent>{});
	}

	bitvec<dynamic_extent> subspan(index_type offset,
	                             index_type count = dynamic_extent) const {
		return make_subspan(offset, count, subspan_selector<Extent>{});
	}

	// [bitvec.obs], bitvec observers
	index_type size() const noexcept { return storage_.size(); }
	index_type size_bytes() const noexcept {
		return size() * narrow_cast<index_type>(sizeof(element_type));
	}
	bool empty() const noexcept { return size() == 0; }

	bitvec<1> operator[](index_type idx) const { return {data_plus(idx), 1}; }

	pointer at(index_type idx) const {
		Expects(CheckRange(idx, idx < storage_.size()));
		return data()[idx];
	}
	pointer operator()(index_type idx) const { return this->operator[](idx); }
	pointer data() const noexcept { return storage_.data(); }
	element_type *cptr() const noexcept { return storage_.data().get(); }

	// [bitvec.iter], bitvec iterator support
	iterator begin() const noexcept { return {this, 0}; }
	iterator end() const noexcept { return {this, size()}; }

	const_iterator cbegin() const noexcept { return {this, 0}; }
	const_iterator cend() const noexcept { return {this, size()}; }

	reverse_iterator rbegin() const noexcept { return reverse_iterator{end()}; }
	reverse_iterator rend() const noexcept { return reverse_iterator{begin()}; }

	const_reverse_iterator crbegin() const noexcept {
		return const_reverse_iterator{cend()};
	}
	const_reverse_iterator crend() const noexcept {
		return const_reverse_iterator{cbegin()};
	}

#ifdef _MSC_VER
	// Tell MSVC how to unwrap spans in range-based-for
	pointer _Unchecked_begin() const noexcept { return data(); }
	pointer _Unchecked_end() const noexcept {
		GSL_SUPPRESS(bounds .1) // NO-FORMAT: attribute
		return data() + size();
	}
#endif // _MSC_VER

  private:
	static constexpr bool CheckRange(index_type idx, index_type size) {
		// Optimization:
		//
		// idx >= 0 && idx < size
		// =>
		// static_cast<size_t>(idx) < static_cast<size_t>(size)
		//
		// because size >=0 by bitvec construction, and negative idx will
		// wrap around to a value always greater than size when casted.

		// check if we have enough space to wrap around
		if (sizeof(index_type) <= sizeof(size_t)) {
			return narrow_cast<size_t>(idx) < narrow_cast<size_t>(size);
		} else {
			return idx >= 0 && idx < size;
		}
	}

	// Needed to remove unnecessary null check in subspans
	struct KnownNotNull {
		pointer p;
	};

	// this implementation detail class lets us take advantage of the
	// empty base class optimization to pay for only storage of a single
	// pointer in the case of fixed-size spans
	template <class ExtentType> class storage_type : public ExtentType {
	  public:
		// KnownNotNull parameter is needed to remove unnecessary null check
		// in subspans and constructors from arrays
		template <class OtherExtentType>
		storage_type(KnownNotNull data, OtherExtentType ext)
		    : ExtentType(ext), data_(data.p) {
			Expects(ExtentType::size() >= 0);
		}

		template <class OtherExtentType>
		storage_type(pointer data, OtherExtentType ext)
		    : ExtentType(ext), data_(data) {
			Expects(ExtentType::size() >= 0);
			Expects(data || ExtentType::size() == 0);
		}

		pointer data() const noexcept { return data_; }

	  private:
		pointer data_;
	};

	storage_type<details::extent_type<Extent>> storage_;

	// The rest is needed to remove unnecessary null check
	// in subspans and constructors from arrays
	bitvec(KnownNotNull ptr, index_type count) : storage_(ptr, count) {}

	template <std::ptrdiff_t CallerExtent> class subspan_selector {};

	template <std::ptrdiff_t CallerExtent>
	bitvec<dynamic_extent> make_subspan(index_type offset, index_type count,
	                                  subspan_selector<CallerExtent>) const {
		const bitvec<dynamic_extent> tmp(*this);
		return tmp.subspan(offset, count);
	}

	GSL_SUPPRESS(bounds .1) // NO-FORMAT: attribute
	bitvec<dynamic_extent> make_subspan(index_type offset, index_type count,
	                                  subspan_selector<dynamic_extent>) const {
		Expects(offset >= 0 && size() - offset >= 0);

		pointer p = data_plus(offset);

		if (count == dynamic_extent) {
			return {KnownNotNull{p}, size() - offset};
		}

		Expects(count >= 0 && size() - offset >= count);
		return {KnownNotNull{p}, count};
	}

  private:
	// Implements `data() + offset` taking care of ownership.
	pointer data_plus(index_type offset) const {
		Expects(offset >= 0 && size() - offset >= 0);

		// Uses same ref counter as data()
		return pointer(data(), data().get() + offset);
	}
};

} // namespace gb

#ifdef _MSC_VER
#pragma warning(pop)
#endif // _MSC_VER

#endif // BITVEC_H

max_stars_repo_name: CapacitorSet/FHE-tools

max_issues_repo_name: CapacitorSet/Glovebox

max_forks_repo_name: CapacitorSet/Glovebox


text:
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <string.h>
#include <float.h>
#include <time.h>
#include <errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

#include "function.h"
#include "norm.h"

// Exception system configuration
#define ERR_T int
#define ERR_V ERRNO
#define OK 0
#include "exceptions.h"

// Input
#define MAX_BUF 256
#define MATCH(line, str) !strncmp(line, str, strlen(str))

// User options
#define DEFAULT_TOLERANCE 1.0e-12
#define DEFAULT_REL_TOL 1
#define DEFAULT_MAX_ZERO_DIST 1.0e-12
#define DEFAULT_USER_NORM 0
#define DEFAULT_MAX_ITER 25
#define DEFAULT_MAX_DIV_ITER 10
#define DEFAULT_JX_REUSE 5

struct options {
	double tolerance;
	char rel_tol;
	double max_zero_dist;
	char user_norm;
	unsigned int max_iter;
	unsigned int max_div_iter;
	unsigned int jx_reuse;
};

// Result
struct additional_data {
	double max_error;
	double* fx;
	unsigned int iter_count;
	double delta_t;
};

// Function declarations

/*
 * The result will be stored in x and x0 will contain the difference between x
 * and the result of the previous iteration.
 */
ERR_T findroot(int dim, double* x0, double* x, struct options* options,
		struct additional_data* data);

ERR_T norm(int dim, double* x, double* n, struct options* options);

ERR_T input_data(char* path, int* dim, double** x0, struct options* options);

void output_result(int dim, double* result, struct additional_data* data);

void output_vector(int dim, double* x);

void handler(const char* reason, const char* file, int line, int gsl_errno);

// Program

int main(int argc, char** argv) {
	int dim;
	char* path;
	double* x = NULL;
	double* x0 = NULL;
	struct options options;
	struct additional_data additional_data;

	// Save the conf file path
	TRY(1, argc == 2)
	path = argv[1];

	// Initialize options with default data
	options.tolerance = DEFAULT_TOLERANCE;
	options.rel_tol = DEFAULT_REL_TOL;
	options.max_zero_dist = DEFAULT_MAX_ZERO_DIST;
	options.user_norm = DEFAULT_USER_NORM;
	options.max_iter = DEFAULT_MAX_ITER;
	options.max_div_iter = DEFAULT_MAX_DIV_ITER;
	options.jx_reuse = DEFAULT_JX_REUSE;

	// Get conf file data
	TRY(2, input_data(path, &dim, &x0, &options) == OK)

	// Find root
	TRY(3, (x = (double*) malloc(dim * sizeof(double))) != NULL)
	TRY(4, findroot(dim, x0, x, &options, &additional_data) == OK)

	output_result(dim, x, &additional_data);

	RETURN(0)

	EXCEPT(
		case 1:
			printf("Usage: findroot path\n");
			break;
		case 2:
			fprintf(stderr,
					"[x] Ezin izan da konfigurazio fitxategia ondo "
					"irakurri.\n");
			break;
		case 3:
			fprintf(stderr,
					"[x] Ezin izan da memoria nahikoa erreserbatu.\n");
			printf("[?] MEMORIA KOPURUA: %lu\n", dim * sizeof(double));
			break;
		case 4:
			fprintf(stderr, "[x] Ezin izan da emaitza kalkulatu.\n");
			break;
	)

	FINALLY(
		free(x);
		free(x0);
	)
}

ERR_T findroot(int dim, double* x0, double* x, struct options* options,
		struct additional_data* data) {
	int s;
	unsigned int iter_count, iter_div_count, jx_reuse_count;
	double max_err, max_err_prev, zero_dist;
	double* fx = NULL;
	double* jx = NULL;
	clock_t begin, end;

	// INITIALIZATION

	TRY(1, (fx = (double*) malloc(dim * sizeof(double))) != NULL)
	TRY(2, (jx = (double*) malloc(dim * dim * sizeof(double))) != NULL)

	gsl_set_error_handler(&handler);

	gsl_vector_view x_gsl = gsl_vector_view_array(x, dim);
	gsl_vector_view x0_gsl = gsl_vector_view_array(x0, dim);
	gsl_vector_view fx_gsl = gsl_vector_view_array(fx, dim);
	gsl_matrix_view jx_gsl = gsl_matrix_view_array(jx, dim, dim);
	gsl_permutation* p = gsl_permutation_alloc(dim);

	// NEWTON-RAPHSON LOOP
	begin = clock();

	/*
	 * a, b and c are vectors of dim dimensions
	 * FX and JX are square matrixes of dim dimensions
	 *
	 * b = a - c
	 * WHERE JX * c = FX
	 *
	 * a beign the result of the previous iteration (or the initial point)
	 * b beign the result of the current iteration (or the final result)
	 */

	// x0 == a

	memcpy(x, x0, dim * sizeof(double));

	// x == a
	// x0 == a

	iter_count = 0;
	iter_div_count = 0;
	max_err = DBL_MAX;
	max_err_prev = DBL_MAX;
	jx_reuse_count = options->jx_reuse;

	do {
		f(dim, x, fx);

		// Reusage of the Jacobian matrix
		if (jx_reuse_count == options->jx_reuse) {
			jakobiarra(dim, x, jx);
			TRY(3, gsl_linalg_LU_decomp(&jx_gsl.matrix, p, &s) == OK)
			jx_reuse_count = 0;
		} else {
			++jx_reuse_count;
		}

		TRY(4,
			gsl_linalg_LU_solve(
				&jx_gsl.matrix, p, &fx_gsl.vector, &x0_gsl.vector) == OK
		)

		// x == a
		// x0 == c

		TRY(5, gsl_vector_sub(&x_gsl.vector, &x0_gsl.vector) == OK)

		TRY(6, norm(dim, x0, &max_err, options) == OK)

		// Relative error
		if (options->rel_tol) {
			// zero_dist is reused to save memory
			TRY(7, norm(dim, x, &zero_dist, options) == OK)
			max_err /= zero_dist;
		}

		// Zero dist
		TRY(8, norm(dim, fx, &zero_dist, options) == OK)

		// x == b
		// x0 == c

		// UPDATES

		++iter_count;

		// Track divergence iterations
		if (max_err > max_err_prev) {
			++iter_div_count;
		} else {
			iter_div_count = 0;
		}
		max_err_prev = max_err;

	} while ((max_err > options->tolerance ||
			zero_dist > options->max_zero_dist) &&
			iter_count < options->max_iter &&
			iter_div_count < options->max_div_iter);

	end = clock();

	if (iter_count == options->max_iter)
		printf("[!] Iterazio mugara iritsi da.\n");
	if (iter_div_count == options->max_div_iter)
		printf("[!] Iterazio dibergente mugara iritsi da.\n");

	data->max_error = max_err;
	data->fx = fx;
	data->iter_count = iter_count;
	data->delta_t = (double) (end - begin) / (double) CLOCKS_PER_SEC;

	EXCEPT(
		case 1:
			fprintf(stderr,
					"[x] Ezin izan da memoria nahikoa erreserbatu.\n");
			printf("[?] MEMORIA KOPURUA: %lu\n", dim * sizeof(double));
			break;
		case 2:
			fprintf(stderr,
					"[x] Ezin izan da memoria nahikoa erreserbatu.\n");
			printf("[?] MEMORIA KOPURUA: %lun", dim * dim * sizeof(double));
			break;
		case 3:
			fprintf(stderr,
					"[x] Ezin izan da JX * x = FX ekuazio sistema "
					"linealaren LU deskonposaketa egin.\n");
			printf("[?] HURBILPEN PARTZIALA: ");
			output_vector(dim, x);
			printf("[?] ERRORE MAXIMOA: %.*g\n", DBL_DIG, max_err);
			printf("[?] F(X): ");
			output_vector(dim, fx);
			printf("[?] ITERAZIO_KOPURUA: %u\n", iter_count);
			break;
		case 4:
			fprintf(stderr,
					"[x] Ezin izan da JX * x = FX ekuazio sistema lineala "
					"ebatzi LU deskonposaketa erabiliz.\n");
			printf("[?] HURBILPEN PARTZIALA: ");
			output_vector(dim, x);
			printf("[?] ERRORE MAXIMOA: %.*g\n", DBL_DIG, max_err);
			printf("[?] F(X): ");
			output_vector(dim, fx);
			printf("[?] ITERAZIO_KOPURUA: %u\n", iter_count);
			break;
		case 5:
			fprintf(stderr,
					"[x] Bektoreen arteko kenketa egitean errore kritiko "
					"bat egon da.");
			printf("[?] F(X): ");
			output_vector(dim, fx);
			printf("[?] ITERAZIO_KOPURUA: %u\n", iter_count);
			break;
		case 6:
			fprintf(stderr,
					"[x] JX * x = FX ekuazio sistema linealaren emaitzaren "
					"norma kalkulatzean errore kritiko bat egon da.\n");
			printf("[?] F(X): ");
			output_vector(dim, fx);
			printf("[?] ITERAZIO_KOPURUA: %u\n", iter_count);
			break;
		case 7:
			fprintf(stderr,
					"[x] x-ren norma kalkulatzean errore kritiko bat egon "
					"da.\n");
			printf("[?] HURBILPEN PARTZIALA: ");
			output_vector(dim, x);
			printf("[?] F(X): ");
			output_vector(dim, fx);
			printf("[?] ITERAZIO_KOPURUA: %u\n", iter_count);
			break;
		case 8:
			fprintf(stderr,
					"[x] f(x)-ren norma kalkulatzean errore kritiko bat egon "
					"da.\n");
			printf("[?] HURBILPEN PARTZIALA: ");
			output_vector(dim, x);
			printf("[?] ERRORE MAXIMOA: %.*g\n", DBL_DIG, max_err);
			output_vector(dim, fx);
			printf("[?] ITERAZIO_KOPURUA: %u\n", iter_count);
			break;
	)

	FINALLY(
		free(fx);
		free(jx);
		gsl_permutation_free(p);
	)
}

ERR_T norm(int dim, double* x, double* n, struct options* options) {
	if (options->user_norm) {
		norma(dim, x, n);
	} else {
		gsl_vector_view x_gsl = gsl_vector_view_array(x, dim);
		TRY(1, (*n = gsl_blas_dnrm2(&x_gsl.vector)) >= 0.0)
	}

	EXCEPT(
		case 1:
			fprintf(stderr, "[x] Ezin izan da norma kalkulatu, posible da "
					"bektorea handiegia izatea.\n");
			break;
	)

	FINALLY()
}

ERR_T input_data(char* path, int* dim, double** x0, struct options* options) {
	FILE* f;
	char line[MAX_BUF];
	int count = 0;

	TRY(1, (f = fopen(path, "r")) != NULL)

	// Read dimension
	*dim = 0;
	while (fgets(line, MAX_BUF, f) != NULL) {
		if (line[0] != '#' && (line[0] != '\0')) {
			if (MATCH(line, "dimentsioa")) {
				TRY(2, sscanf(line, "dimentsioa %d", dim) == 1)
				TRY(3, *dim > 0)
				break;
			}
		}
	}

	TRY(4, *dim != 0)

	TRY(5, (*x0 = (double*) malloc((*dim) * sizeof(double))) != NULL)

	// Read x0 and optional parameters
	while (fgets(line, MAX_BUF, f) != NULL) {
		if (line[0] != '#' && (line[0] != '\0')) {
			if (MATCH(line, "tolerantzia")) {
				TRY(6,
					sscanf(line, "tolerantzia %lf", &(options->tolerance)) == 1)
				TRY(7, options->tolerance > 0.0)
			} else if (MATCH(line, "erlatiboa")) {
				if (MATCH(line, "erlatiboa bai")) {
					options->rel_tol = 1;
				} else if (MATCH(line, "erlatiboa ez")) {
					options->rel_tol = 0;
				} else {
					TRY(8, 0);
				}
			} else if (MATCH(line, "max_zero_distantzia")) {
				TRY(9,
					sscanf(line, "max_zero_distantzia %lf",
							&(options->max_zero_dist)) == 1)
				TRY(10, options->max_zero_dist > 0.0)
			} else if (MATCH(line, "ordezko_norma")) {
				if (MATCH(line, "ordezko_norma bai")) {
					options->user_norm = 1;
				} else if (MATCH(line, "ordezko_norma ez")) {
					options->user_norm = 0;
				} else {
					TRY(11, 0);
				}
			} else if (MATCH(line, "iterazio_maximoa")) {
				TRY(12,	sscanf(line, "iterazio_maximoa %u",
							&(options->max_iter)) == 1)
				TRY(13, options->max_iter != 0)
			} else if (MATCH(line, "dibergentzia_iterazio_maximoa")) {
				TRY(14, sscanf(line, "dibergentzia_iterazio_maximoa %u",
						&(options->max_div_iter)) == 1)
			} else if (MATCH(line, "jakobiar_berrerabilpena")) {
				TRY(15, sscanf(line, "jakobiar_berrerabilpena %u",
							&(options->jx_reuse)) == 1)
			} else {
				TRY(16, count <= *dim)
				TRY(17, sscanf(line, "%lf", (*x0) + count) == 1)
				++count;
			}
		}
	}

	TRY(18, count == *dim)

	EXCEPT(
		case 1:
			fprintf(stdout,
					"[x] Ezin izan da konfigurazio fitxategia ireki.\n");
			break;
		case 2:
			fprintf(stdout,
					"[x] Sintaxi desegokia dimentsioa emateko lerroan\n");
			printf("[?] LERROA: %s", line);
			break;
		case 3:
			fprintf(stdout,
					"[x] Dimentsioak zero baina handiagoa izan behar du.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 4:
			fprintf(stdout, "[x] Ez da aurkitu dimentsioa.\n");
			break;
		case 5:
			fprintf(stdout,
					"[x] Ezin izan da X0-rentzat memoria erreserbatu\n.");
			printf("[?] MEMORIA KOPURUA: %lu", (*dim) * sizeof(double));
			break;
		case 6:
			fprintf(stdout,
					"[x] Sintaxi desegokia tolerantzia emateko lerroan\n.");
			printf("[?] LERROA: %s", line);
			break;
		case 7:
			fprintf(stdout,
					"[x] Tolerantziak zero baina handiagoa izan behar du.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 8:
			fprintf(stdout,
					"[x] Sintaxi desegokia tolerantzia erlatiboa aukeratzeko "
					"lerroan. Aukera honek 'bai' eta 'ez' balioak hartu "
					"ditzake soilik.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 9:
			fprintf(stdout,
					"[x] Sintaxi okerra zerotik distantzia maximoa emateko "
					"lerroan.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 10:
			fprintf(stdout,
					"[x] Zerotik distantzia maximoak zero baina handiagoa izan "
					"behar du.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 11:
			fprintf(stdout,
					"[x] Sintaxi desegokia ordezko norma aukeratzeko "
					"lerroan. Aukera honek 'bai' eta 'ez' balioak hartu "
					"ditzake soilik.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 12:
			fprintf(stdout,
					"[x] Sintaxi okerra iterazio kopuru maximoa emateko "
					"lerroan. Kopuruak ezin du negatiboa izan.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 13:
			fprintf(stdout,
					"[x] Iterazio kopuru maximoak ezin du zero izan.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 14:
			fprintf(stdout,
					"[x] Sintaxi okerra iterazio dibergente kopuru maximoa "
					"emateko lerroan. Kopuruak ezin du negatiboa izan.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 15:
			fprintf(stdout,
					"[x] Sintaxi okerra jakobiarraren berrerabilpen ratioa "
					"emateko lerroan. Ratioak ezin du negatiboa izan.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 16:
			fprintf(stdout,
					"[x] X0-ren elementu kopurura dimentsioa baina handiagoa "
					"da.\n");
			printf("[?] DIMENTSIOA: %i\n", *dim);
			break;
		case 17:
			fprintf(stdout,
					"[x] Sintaxi desegokia.\n");
			printf("[?] LERROA: %s", line);
			break;
		case 18:
			fprintf(stdout,
					"[x] Dimentsioa eta X0-ren tamaina ez datoz bat.\n");
			printf("[?] DIMENTSIOA: %i\n", *dim);
			printf("[?] TAMAINA: %i\n", count);
			break;
		case 19:
			fprintf(stdout,
					"[x] Konfigurazio fitxategia ixtean errore kritiko bat "
					"egon da.\n");
			break;
	)

	FINALLY(
		if (f != NULL) {
			TRY(19, fclose(f) == OK);

			// On fail, so that it doesn't try again and again
			f = NULL;
		}
	)
}

void output_result(int dim, double* result, struct additional_data* data) {
	printf("\n");

	printf("Erroa: ");
	output_vector(dim, result);

	printf("Errore maximoa: %.*g\n", DBL_DIG, data->max_error);

	printf("Erroaren irudia: ");
	output_vector(dim, data->fx);

	printf("Iterazio kopurua: %u\n", data->iter_count);

	printf("Denbora: %.*g seg\n", DBL_DIG, data->delta_t);

	printf("\n");
}

void output_vector(int dim, double* x) {
	int i;

	printf("(");
	printf("%.*g", DBL_DIG, x[0]);
	for (i = 1; i < dim; ++i)
		printf(", %.*g", DBL_DIG, x[i]);
	printf(")\n");
}

void handler(const char* reason, const char* file, int line, int gsl_errno) {
	fprintf(stderr, "[x] ");

	switch(gsl_errno) {
		case GSL_EDOM:
			fprintf(stderr, "GSL DOMAIN ERROR: ");
			break;
		case GSL_ERANGE:
			fprintf(stderr, "GSL RANGE ERROR: ");
			break;
		case GSL_ENOMEM:
			fprintf(stderr, "GSL NO MEMORY AVAILABLE: ");
			break;
		case GSL_EINVAL:
			fprintf(stderr, "GSL INVALID ARGUMENT: ");
			break;
		default:
			fprintf(stderr, "GSL ERROR: ");
			break;
	}

	fprintf(stderr, "%s\n", reason);
}

ext: c

max_stars_repo_name: oersted/find-root

max_issues_repo_name: oersted/find-root

max_forks_repo_name: oersted/find-root


text:
/// @file
///
/// @brief test of gsl pointer wrappers
///
/// @copyright (c) 2007 CSIRO
/// Australia Telescope National Facility (ATNF)
/// Commonwealth Scientific and Industrial Research Organisation (CSIRO)
/// PO Box 76, Epping NSW 1710, Australia
/// atnf-enquiries@csiro.au
///
/// This file is part of the ASKAP software distribution.
///
/// The ASKAP software distribution is free software: you can redistribute it
/// and/or modify it under the terms of the GNU General Public License as
/// published by the Free Software Foundation; either version 2 of the License,
/// or (at your option) any later version.
///
/// This program is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with this program; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
///
/// @author Max Voronkov <maxim.voronkov@csiro.au>

// unit test include
#include <cppunit/extensions/HelperMacros.h>

// own includes
#include "utils/SharedGSLTypes.h"

// gsl includes
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

namespace askap {

namespace utility {

/// @brief helper specialisation to allow testing of the destructor call
/// @details It is not used anywhere outside this test suit
template<> struct CustomGSLDeleter<bool> {
    /// @brief this method just sets the passed object to true
    /// @param[in] obj pointer to a boolean variable
    void operator()(bool *obj) const { *obj = true; }
};

class SharedGSLTypesTest : public CppUnit::TestFixture
{
   CPPUNIT_TEST_SUITE(SharedGSLTypesTest);
   CPPUNIT_TEST(testVector);
   CPPUNIT_TEST(testMatrix);
   CPPUNIT_TEST_EXCEPTION(testNullPointer,AskapError);
   CPPUNIT_TEST(testDestruction);
   CPPUNIT_TEST_SUITE_END();
public:
   void testVector() {
      const size_t nElements = 10;
      SharedGSLVector vec = createGSLVector(nElements);
      for (size_t el = 0; el < nElements; ++el) {
           gsl_vector_set(vec.get(), el, double(el));
      }
      // check the content (can't really check that the destructor is
      // called)
      for (size_t el = 0; el < nElements; ++el) {
          CPPUNIT_ASSERT_DOUBLES_EQUAL(double(el), gsl_vector_get(vec.get(), el), 1e-6);
      }
      // destructor should be called on leaving this method
   }

   void testMatrix() {
      const size_t nRow = 10;
      const size_t nCol = 12;
      SharedGSLMatrix matr = createGSLMatrix(nRow, nCol);
      for (size_t row = 0; row < nRow; ++row) {
           for (size_t col = 0; col < nCol; ++col) {
                gsl_matrix_set(matr.get(), row, col, double(row*col));
           }
      }

      // check the content (can't really check that the destructor is
      // called)
      for (size_t row = 0; row < nRow; ++row) {
           for (size_t col = 0; col < nCol; ++col) {
                CPPUNIT_ASSERT_DOUBLES_EQUAL(double(row*col), gsl_matrix_get(matr.get(), row,col), 1e-6);
           }
      }
      // destructor should be called on leaving this method
   }

   void testNullPointer() {
      gsl_vector *nullVec = NULL;
      // the following should throw an exception
      createGSLObject(nullVec);
   }

   void testDestruction() {
      bool destructorCalledBuf = false;
      {
        const boost::shared_ptr<bool> destructorCalledPtr = createGSLObject(&destructorCalledBuf);
        CPPUNIT_ASSERT(destructorCalledPtr);
        CPPUNIT_ASSERT_EQUAL(false, *destructorCalledPtr);
        CPPUNIT_ASSERT_EQUAL(false, destructorCalledBuf);
      }
      CPPUNIT_ASSERT_EQUAL(true, destructorCalledBuf);
   }
};

} // namespace utility

} // namespace askap

max_stars_repo_name: rtobar/askapsoft

max_issues_repo_name: ATNF/askapsoft

max_forks_repo_name: ATNF/askapsoft


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

/*
 * This generator is taken from
 *
 * Donald E. Knuth
 * The Art of Computer Programming
 * Volume 2
 * Third Edition
 * Addison-Wesley
 * Page 108
 *
 * This implementation copyright (C) 2001 Brian Gough, Carlo Perassi
 * and (C) 2003 Heiko Bauke.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#define AAA 40692
#define MMM 2147483399UL
#define QQQ 52774
#define RRR 3791

static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);

typedef struct
{
  unsigned long int x;
}
ran_state_t;

static inline unsigned long int
ran_get (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  long int y = state->x;
  long int r = RRR * (y / QQQ);

  y = AAA * (y % QQQ) - r;
  if (y < 0)
    y += MMM;

  state->x = y;

  return state->x;
}

static double
ran_get_double (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 2147483399.0;
}

static void
ran_set (void *vstate, unsigned long int s)
{
  ran_state_t *state = (ran_state_t *) vstate;

  if ((s%MMM) == 0)
    s = 1;                      /* default seed is 1 */

  state->x = s % MMM;

  return;
}

static const gsl_rng_type ran_type = {
  "lecuyer21",                  /* name */
  MMM-1,                        /* RAND_MAX */
  1,                            /* RAND_MIN */
  sizeof (ran_state_t),
  &ran_set,
  &ran_get,
  &ran_get_double
};

const gsl_rng_type *gsl_rng_lecuyer21 = &ran_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef libceed_solids_examples_utils_h
#define libceed_solids_examples_utils_h

#include <ceed.h>
#include <petsc.h>

// Translate PetscMemType to CeedMemType
static inline CeedMemType MemTypeP2C(PetscMemType mem_type) {
  return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST;
}

#endif // libceed_solids_examples_utils_h

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
/* Copyright (C) 2021 Atsushi Togo */
/* All rights reserved. */

/* This file is part of phonopy. */

/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions */
/* are met: */

/* * Redistributions of source code must retain the above copyright */
/*   notice, this list of conditions and the following disclaimer. */

/* * Redistributions in binary form must reproduce the above copyright */
/*   notice, this list of conditions and the following disclaimer in */
/*   the documentation and/or other materials provided with the */
/*   distribution. */

/* * Neither the name of the phonopy project 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 __phono3py_H__
#define __phono3py_H__

#ifdef MKL_LAPACKE
#include <mkl.h>
#else
#include <lapacke.h>
#endif
#include "phonoc_array.h"

long ph3py_get_interaction(Darray *fc3_normal_squared,
                           const char *g_zero,
                           const Darray *frequencies,
                           const lapack_complex_double *eigenvectors,
                           const long (*triplets)[3],
                           const long num_triplets,
                           const long (*bz_grid_addresses)[3],
                           const long D_diag[3],
                           const long Q[3][3],
                           const double *fc3,
                           const long is_compact_fc3,
                           const double (*svecs)[3],
                           const long multi_dims[2],
                           const long (*multi)[2],
                           const double *masses,
                           const long *p2s_map,
                           const long *s2p_map,
                           const long *band_indices,
                           const long symmetrize_fc3_q,
                           const double cutoff_frequency);
long ph3py_get_pp_collision(double *imag_self_energy,
                            const long relative_grid_address[24][4][3], /* thm */
                            const double *frequencies,
                            const lapack_complex_double *eigenvectors,
                            const long (*triplets)[3],
                            const long num_triplets,
                            const long *triplet_weights,
                            const long (*bz_grid_addresses)[3], /* thm */
                            const long *bz_map,                 /* thm */
                            const long bz_grid_type,
                            const long D_diag[3],
                            const long Q[3][3],
                            const double *fc3,
                            const long is_compact_fc3,
                            const double (*svecs)[3],
                            const long multi_dims[2],
                            const long (*multi)[2],
                            const double *masses,
                            const long *p2s_map,
                            const long *s2p_map,
                            const Larray *band_indices,
                            const Darray *temperatures,
                            const long is_NU,
                            const long symmetrize_fc3_q,
                            const double cutoff_frequency);
long ph3py_get_pp_collision_with_sigma(
    double *imag_self_energy,
    const double sigma,
    const double sigma_cutoff,
    const double *frequencies,
    const lapack_complex_double *eigenvectors,
    const long (*triplets)[3],
    const long num_triplets,
    const long *triplet_weights,
    const long (*bz_grid_addresses)[3],
    const long D_diag[3],
    const long Q[3][3],
    const double *fc3,
    const long is_compact_fc3,
    const double (*svecs)[3],
    const long multi_dims[2],
    const long (*multi)[2],
    const double *masses,
    const long *p2s_map,
    const long *s2p_map,
    const Larray *band_indices,
    const Darray *temperatures,
    const long is_NU,
    const long symmetrize_fc3_q,
    const double cutoff_frequency);
void ph3py_get_imag_self_energy_at_bands_with_g(
    double *imag_self_energy,
    const Darray *fc3_normal_squared,
    const double *frequencies,
    const long (*triplets)[3],
    const long *triplet_weights,
    const double *g,
    const char *g_zero,
    const double temperature,
    const double cutoff_frequency,
    const long num_frequency_points,
    const long frequency_point_index);
void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
    double *detailed_imag_self_energy,
    double *imag_self_energy_N,
    double *imag_self_energy_U,
    const Darray *fc3_normal_squared,
    const double *frequencies,
    const long (*triplets)[3],
    const long *triplet_weights,
    const long (*bz_grid_addresses)[3],
    const double *g,
    const char *g_zero,
    const double temperature,
    const double cutoff_frequency);
void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
                                         const Darray *fc3_normal_squared,
                                         const long *band_indices,
                                         const double *frequencies,
                                         const long (*triplets)[3],
                                         const long *triplet_weights,
                                         const double epsilon,
                                         const double temperature,
                                         const double unit_conversion_factor,
                                         const double cutoff_frequency);
void ph3py_get_real_self_energy_at_frequency_point(
    double *real_self_energy,
    const double frequency_point,
    const Darray *fc3_normal_squared,
    const long *band_indices,
    const double *frequencies,
    const long (*triplets)[3],
    const long *triplet_weights,
    const double epsilon,
    const double temperature,
    const double unit_conversion_factor,
    const double cutoff_frequency);
void ph3py_get_collision_matrix(double *collision_matrix,
                                const Darray *fc3_normal_squared,
                                const double *frequencies,
                                const long (*triplets)[3],
                                const long *triplets_map,
                                const long *map_q,
                                const long *rotated_grid_points,
                                const double *rotations_cartesian,
                                const double *g,
                                const long num_ir_gp,
                                const long num_gp,
                                const long num_rot,
                                const double temperature,
                                const double unit_conversion_factor,
                                const double cutoff_frequency);
void ph3py_get_reducible_collision_matrix(double *collision_matrix,
                                          const Darray *fc3_normal_squared,
                                          const double *frequencies,
                                          const long (*triplets)[3],
                                          const long *triplets_map,
                                          const long *map_q,
                                          const double *g,
                                          const long num_gp,
                                          const double temperature,
                                          const double unit_conversion_factor,
                                          const double cutoff_frequency);
void ph3py_get_isotope_scattering_strength(
    double *gamma,
    const long grid_point,
    const double *mass_variances,
    const double *frequencies,
    const lapack_complex_double *eigenvectors,
    const long num_grid_points,
    const long *band_indices,
    const long num_band,
    const long num_band0,
    const double sigma,
    const double cutoff_frequency);
void ph3py_get_thm_isotope_scattering_strength(
    double *gamma,
    const long grid_point,
    const long *ir_grid_points,
    const long *weights,
    const double *mass_variances,
    const double *frequencies,
    const lapack_complex_double *eigenvectors,
    const long num_ir_grid_points,
    const long *band_indices,
    const long num_band,
    const long num_band0,
    const double *integration_weights,
    const double cutoff_frequency);
void ph3py_distribute_fc3(double *fc3,
                          const long target,
                          const long source,
                          const long *atom_mapping,
                          const long num_atom,
                          const double *rot_cart);
void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
                            const double (*delta_fc2s)[3][3],
                            const double *inv_U,
                            const double (*site_sym_cart)[3][3],
                            const long *rot_map_syms,
                            const long num_atom,
                            const long num_site_sym,
                            const long num_disp);
void ph3py_get_permutation_symmetry_fc3(double *fc3, const long num_atom);
void ph3py_get_permutation_symmetry_compact_fc3(double *fc3,
                                                const long p2s[],
                                                const long s2pp[],
                                                const long nsym_list[],
                                                const long perms[],
                                                const long n_satom,
                                                const long n_patom);
void ph3py_transpose_compact_fc3(double *fc3,
                                 const long p2s[],
                                 const long s2pp[],
                                 const long nsym_list[],
                                 const long perms[],
                                 const long n_satom,
                                 const long n_patom,
                                 const long t_type);
long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                             long *map_q,
                                             const long grid_point,
                                             const long mesh[3],
                                             const long is_time_reversal,
                                             const long num_rot,
                                             const long (*rec_rotations)[3][3],
                                             const long swappable);
long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
                                const long grid_point,
                                const long (*bz_grid_addresses)[3],
                                const long *bz_map,
                                const long *map_triplets,
                                const long num_map_triplets,
                                const long D_diag[3],
                                const long Q[3][3],
                                const long bz_grid_type);
long ph3py_get_integration_weight(double *iw,
                                  char *iw_zero,
                                  const double *frequency_points,
                                  const long num_band0,
                                  const long relative_grid_address[24][4][3],
                                  const long mesh[3],
                                  const long (*triplets)[3],
                                  const long num_triplets,
                                  const long (*bz_grid_addresses)[3],
                                  const long *bz_map,
                                  const long bz_grid_type,
                                  const double *frequencies1,
                                  const long num_band1,
                                  const double *frequencies2,
                                  const long num_band2,
                                  const long tp_type,
                                  const long openmp_per_triplets,
                                  const long openmp_per_bands);
void ph3py_get_integration_weight_with_sigma(double *iw,
                                             char *iw_zero,
                                             const double sigma,
                                             const double sigma_cutoff,
                                             const double *frequency_points,
                                             const long num_band0,
                                             const long (*triplets)[3],
                                             const long num_triplets,
                                             const double *frequencies,
                                             const long num_band,
                                             const long tp_type);
long ph3py_get_grid_index_from_address(const long address[3],
                                       const long mesh[3]);
void ph3py_get_gr_grid_addresses(long gr_grid_addresses[][3],
                                 const long D_diag[3]);
long ph3py_get_reciprocal_rotations(long rec_rotations[48][3][3],
                                    const long (*rotations)[3][3],
                                    const long num_rot,
                                    const long is_time_reversal);
long ph3py_transform_rotations(long (*transformed_rots)[3][3],
                               const long (*rotations)[3][3],
                               const long num_rot,
                               const long D_diag[3],
                               const long Q[3][3]);
long ph3py_get_snf3x3(long D_diag[3],
                      long P[3][3],
                      long Q[3][3],
                      const long A[3][3]);
long ph3py_transform_rotations(long (*transformed_rots)[3][3],
                               const long (*rotations)[3][3],
                               const long num_rot,
                               const long D_diag[3],
                               const long Q[3][3]);
long ph3py_get_ir_grid_map(long *ir_grid_map,
                           const long D_diag[3],
                           const long PS[3],
                           const long (*grg_rotations)[3][3],
                           const long num_rot);
long ph3py_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
                                 long *bz_map,
                                 long *bzg2grg,
                                 const long D_diag[3],
                                 const long Q[3][3],
                                 const long PS[3],
                                 const double rec_lattice[3][3],
                                 const long type);
long ph3py_rotate_bz_grid_index(const long bz_grid_index,
                                const long rotation[3][3],
                                const long (*bz_grid_addresses)[3],
                                const long *bz_map,
                                const long D_diag[3],
                                const long PS[3],
                                const long bz_grid_type);
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
                                       const long num_column,
                                       const long num_temp,
                                       const long num_sigma);
void ph3py_expand_collision_matrix(double *collision_matrix,
                                   const long *rot_grid_points,
                                   const long *ir_grid_points,
                                   const long num_ir_gp,
                                   const long num_grid_points,
                                   const long num_rot,
                                   const long num_sigma,
                                   const long num_temp,
                                   const long num_band);
long ph3py_get_neighboring_gird_points(long *relative_grid_points,
                                       const long *grid_points,
                                       const long (*relative_grid_address)[3],
                                       const long mesh[3],
                                       const long (*bz_grid_addresses)[3],
                                       const long *bz_map,
                                       const long bz_grid_type,
                                       const long num_grid_points,
                                       const long num_relative_grid_address);
long ph3py_get_thm_integration_weights_at_grid_points(
    double *iw,
    const double *frequency_points,
    const long num_band0,
    const long num_band,
    const long num_gp,
    const long (*relative_grid_address)[4][3],
    const long D_diag[3],
    const long *grid_points,
    const long (*bz_grid_addresses)[3],
    const long *bz_map,
    const long bz_grid_type,
    const double *frequencies,
    const long *gp2irgp_map,
    const char function);

#endif


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

#pragma once

#include "engine_error.h"
#include "protocol_binary.h"
#include "rbac.h"
#include "types.h"

#include <mcbp/protocol/opcode.h>
#include <nlohmann/json_fwd.hpp>
#include <gsl/gsl>
#include <string>

/**
 * Commands to operate on a specific cookie.
 */
struct ServerCookieIface {
    virtual ~ServerCookieIface() = default;

    /**
     * Store engine-specific session data on the given cookie.
     *
     * The engine interface allows for a single item to be
     * attached to the connection that it can use to track
     * connection-specific data throughout duration of the
     * connection.
     *
     * @param cookie The cookie provided by the frontend
     * @param engine_data pointer to opaque data
     */
    virtual void store_engine_specific(gsl::not_null<const void*> cookie,
                                       void* engine_data) = 0;

    /**
     * Retrieve engine-specific session data for the given cookie.
     *
     * @param cookie The cookie provided by the frontend
     *
     * @return the data provied by store_engine_specific or NULL
     *         if none was provided
     */
    virtual void* get_engine_specific(gsl::not_null<const void*> cookie) = 0;

    /**
     * Check if datatype is supported by the connection.
     *
     * @param cookie The cookie provided by the frontend
     * @param datatype The datatype to test
     *
     * @return true if connection supports the datatype or else false.
     */
    virtual bool is_datatype_supported(gsl::not_null<const void*> cookie,
                                       protocol_binary_datatype_t datatype) = 0;

    /**
     * Check if mutation extras is supported by the connection.
     *
     * @param cookie The cookie provided by the frontend
     *
     * @return true if supported or else false.
     */
    virtual bool is_mutation_extras_supported(
            gsl::not_null<const void*> cookie) = 0;

    /**
     * Check if collections are supported by the connection.
     *
     * @param cookie The cookie provided by the frontend
     *
     * @return true if supported or else false.
     */
    virtual bool is_collections_supported(
            gsl::not_null<const void*> cookie) = 0;

    /**
     * Retrieve the opcode of the connection, if
     * ewouldblock flag is set. Please note that the ewouldblock
     * flag for a connection is cleared before calling into
     * the engine interface, so this method only works in the
     * notify hooks.
     *
     * @param cookie The cookie provided by the frontend
     *
     * @return the opcode from the binary_header saved in the
     * connection.
     */
    virtual cb::mcbp::ClientOpcode get_opcode_if_ewouldblock_set(
            gsl::not_null<const void*> cookie) = 0;

    /**
     * Validate given ns_server's session cas token against
     * saved token in memached, and if so incrment the session
     * counter.
     *
     * @param cas The cas token from the request
     *
     * @return true if session cas matches the one saved in
     * memcached
     */
    virtual bool validate_session_cas(uint64_t cas) = 0;

    /**
     * Decrement session_cas's counter everytime a control
     * command completes execution.
     */
    virtual void decrement_session_ctr() = 0;

    /**
     * Let a connection know that IO has completed.
     * @param cookie cookie representing the connection
     * @param status the status for the io operation
     */
    virtual void notify_io_complete(gsl::not_null<const void*> cookie,
                                    ENGINE_ERROR_CODE status) = 0;

    /**
     * Notify the core that we're holding on to this cookie for
     * future use. (The core guarantees it will not invalidate the
     * memory until the cookie is invalidated by calling release())
     */
    virtual ENGINE_ERROR_CODE reserve(gsl::not_null<const void*> cookie) = 0;

    /**
     * Notify the core that we're releasing the reference to the
     * The engine is not allowed to use the cookie (the core may invalidate
     * the memory)
     */
    virtual ENGINE_ERROR_CODE release(gsl::not_null<const void*> cookie) = 0;

    /**
     * Set the priority for this connection
     */
    virtual void set_priority(gsl::not_null<const void*> cookie,
                              CONN_PRIORITY priority) = 0;

    /**
     * Get the priority for this connection
     */
    virtual CONN_PRIORITY get_priority(gsl::not_null<const void*> cookie) = 0;

    /**
     * Get the bucket the connection is bound to
     *
     * @cookie The connection object
     * @return the bucket identifier for a cookie
     */
    virtual bucket_id_t get_bucket_id(gsl::not_null<const void*> cookie) = 0;

    /**
     * Get connection id
     *
     * @param cookie the cookie sent to the engine for an operation
     * @return a unique identifier for a connection
     */
    virtual uint64_t get_connection_id(gsl::not_null<const void*> cookie) = 0;

    /**
     * Check if the cookie have the specified privilege in it's
     * active set.
     *
     * @todo We should probably add the key we want to access as part
     *       of the API. We're going to need that when we're adding
     *       support for collections. For now let's assume that it
     *       won't be a big problem to fix that later on.
     * @param cookie the cookie sent to the engine for an operation
     * @param privilege the privilege to check for
     * @return true if the cookie have the privilege in its active set,
     *         false otherwise
     */
    virtual cb::rbac::PrivilegeAccess check_privilege(
            gsl::not_null<const void*> cookie,
            cb::rbac::Privilege privilege) = 0;

    /**
     * Method to map an engine error code to the appropriate mcbp response
     * code (the client may not support all error codes so we may have
     * to remap some).
     *
     * @param cookie the client cookie (to look up the client connection)
     * @param code the engine error code to get the mcbp response code.
     * @return the mcbp response status to use
     * @throws std::engine_error if the error code results in being
     *                           ENGINE_DISCONNECT after remapping
     *         std::logic_error if the error code doesn't make sense
     *         std::invalid_argument if the code doesn't exist
     */
    virtual cb::mcbp::Status engine_error2mcbp(
            gsl::not_null<const void*> cookie, ENGINE_ERROR_CODE code) = 0;

    /**
     * Get the log information to be used for a log entry.
     *
     * The typical log entry from the core is:
     *
     *  `id> message` - Data read from ta client
     *  `id: message` - Status messages for this client
     *  `id< message` - Data sent back to the client
     *
     * If the caller wants to dump more information about the connection
     * (like socket name, peer name, user name) the pair returns this
     * info as the second field. The info may be invalidated by the core
     * at any time (but not while the engine is operating in a single call
     * from the core) so it should _not_ be cached.
     */
    virtual std::pair<uint32_t, std::string> get_log_info(
            gsl::not_null<const void*> cookie) = 0;

    virtual std::string get_authenticated_user(
            gsl::not_null<const void*> cookie) = 0;

    virtual in_port_t get_connected_port(gsl::not_null<const void*> cookie) = 0;

    /**
     * Set the error context string to be sent in response. This should not
     * contain security sensitive information. If sensitive information needs to
     * be preserved, log it with a UUID and send the UUID.
     *
     * Note this has no affect for the following response codes.
     *   cb::mcbp::Status::Success
     *   cb::mcbp::Status::SubdocSuccessDeleted
     *   cb::mcbp::Status::SubdocMultiPathFailure
     *   cb::mcbp::Status::Rollback
     *   cb::mcbp::Status::NotMyVbucket
     *
     * @param cookie the client cookie (to look up client connection)
     * @param message the message string to be set as the error context
     */
    virtual void set_error_context(gsl::not_null<void*> cookie,
                                   cb::const_char_buffer message) = 0;

    /**
     * Set a JSON object to be included in an error response (along side
     * anything set by set_error_context).
     *
     * The json object cannot include "error" as a top-level key
     *
     * Note this has no affect for the following response codes.
     *   cb::mcbp::Status::Success
     *   cb::mcbp::Status::SubdocSuccessDeleted
     *   cb::mcbp::Status::SubdocMultiPathFailure
     *   cb::mcbp::Status::Rollback
     *   cb::mcbp::Status::NotMyVbucket
     *
     * @param cookie the client cookie (to look up client connection)
     * @param json extra json object to include in a error response.
     */
    virtual void set_error_json_extras(gsl::not_null<void*> cookie,
                                       const nlohmann::json& json) = 0;
};


text:
/* multimin/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_ieee_utils.h>

#include "test_funcs.h"

int
test_fdf(const char * desc, gsl_multimin_function_fdf *f, initpt_function initpt, const gsl_multimin_fdfminimizer_type *T);

int
main (void)
{
  const gsl_multimin_fdfminimizer_type *fdfminimizers[5];
  const gsl_multimin_fdfminimizer_type ** T;

  gsl_ieee_env_setup ();

  fdfminimizers[0] = gsl_multimin_fdfminimizer_steepest_descent;
  fdfminimizers[1] = gsl_multimin_fdfminimizer_conjugate_pr;
  fdfminimizers[2] = gsl_multimin_fdfminimizer_conjugate_fr;
  fdfminimizers[3] = gsl_multimin_fdfminimizer_vector_bfgs;
  fdfminimizers[4] = 0;

  T = fdfminimizers;

  while (*T != 0)
    {
      test_fdf("Roth", &roth, roth_initpt,*T);
      test_fdf("Wood", &wood, wood_initpt,*T);
      test_fdf("Rosenbrock", &rosenbrock, rosenbrock_initpt,*T);
      T++;
    }

  T = fdfminimizers;

  while (*T != 0)
    {
      test_fdf("NRoth", &Nroth, roth_initpt,*T);
      test_fdf("NWood", &Nwood, wood_initpt,*T);
      test_fdf("NRosenbrock", &Nrosenbrock, rosenbrock_initpt,*T);
      T++;
    }

  exit (gsl_test_summary());
}

int
test_fdf(const char * desc,
         gsl_multimin_function_fdf *f,
	 initpt_function initpt,
	 const gsl_multimin_fdfminimizer_type *T)
{
  int status;
  size_t iter = 0;
  double step_size;

  gsl_vector *x = gsl_vector_alloc (f->n);

  gsl_multimin_fdfminimizer *s;

  (*initpt) (x);

  step_size = 0.1 * gsl_blas_dnrm2 (x);

  s = gsl_multimin_fdfminimizer_alloc(T, f->n);

  gsl_multimin_fdfminimizer_set (s, f, x, step_size, 0.1);

#ifdef DEBUG
  printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
  printf("g "); gsl_vector_fprintf (stdout, s->gradient, "%g");
#endif

  do
    {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate(s);

#ifdef DEBUG
      printf("%i: \n",iter);
      printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
      printf("g "); gsl_vector_fprintf (stdout, s->gradient, "%g");
      printf("f(x) %g\n",s->f);
      printf("dx %g\n",gsl_blas_dnrm2(s->dx));
      printf("\n");
#endif

      status = gsl_multimin_test_gradient(s->gradient,1e-3);
    }
  while (iter < 5000 && status == GSL_CONTINUE);

  status |= (fabs(s->f) > 1e-5);

  gsl_test(status, "%s, on %s: %i iterations, f(x)=%g",
	   gsl_multimin_fdfminimizer_name(s),desc, iter, s->f);

  gsl_multimin_fdfminimizer_free(s);
  gsl_vector_free(x);

  return status;
}


text:
#pragma once

#include <gsl/gsl>
#include "halley/utils/utils.h"
#include <vector>

namespace Halley {
	class String;
	class Path;

	class FileSystem
	{
	public:
		static bool exists(const Path& p);
		static bool createDir(const Path& p);
		static bool createParentDir(const Path& p);

		static int64_t getLastWriteTime(const Path& p);
		static bool isFile(const Path& p);
		static bool isDirectory(const Path& p);

		static void copyFile(const Path& src, const Path& dst);
		static bool remove(const Path& path);

		static void writeFile(const Path& path, gsl::span<const gsl::byte> data);
		static void writeFile(const Path& path, const Bytes& data);
		static Bytes readFile(const Path& path);

		static std::vector<Path> enumerateDirectory(const Path& path);

		static Path getRelative(const Path& path, const Path& parentPath);
		static Path getAbsolute(const Path& path);

		static size_t fileSize(const Path& path);

		static Path getTemporaryPath();

		static int runCommand(const String& command);
	};
}

max_stars_repo_name: Healthire/halley

max_issues_repo_name: Healthire/halley

max_forks_repo_name: Healthire/halley