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

#ifndef SALESMAN_EXAMPLE_PATH_MERGER_H
#define SALESMAN_EXAMPLE_PATH_MERGER_H

#include <algorithm>
#include <functional>
#include <iterator>
#include <stack>
#include <utility>
#include <vector>
#include <gsl/gsl_util>
#include <gsl/span>
#include <repeat.h>
#include "disjoint_set_data_structure.h"
#include "permutation.h"

template<class UniformRandomBitGenerator>
class path_merger {
public:
	explicit path_merger(UniformRandomBitGenerator& g) noexcept;
	permutation operator()(const permutation& lhs, const permutation& rhs);
private:
	using edge_type = std::pair<unsigned, unsigned>;
	using edge_vector = std::vector<edge_type>;
	edge_vector to_edges(const permutation& perm) const;
	permutation to_permutation(const edge_vector& edges) const;
	UniformRandomBitGenerator& rand;
};

template<class UniformRandomBitGenerator>
path_merger<UniformRandomBitGenerator>::path_merger(UniformRandomBitGenerator& g) noexcept
	: rand(g) {}

template<class UniformRandomBitGenerator>
permutation path_merger<UniformRandomBitGenerator>::operator()(const permutation& lhs, const permutation& rhs) {
	Expects(lhs.size() == rhs.size());
	Expects(lhs.size() > 0);
	const std::size_t size = lhs.size();
	edge_vector lhs_edges = to_edges(lhs);
	edge_vector rhs_edges = to_edges(rhs);
	edge_vector result_edges;
	result_edges.reserve(size);
	std::set_intersection(lhs_edges.begin(), lhs_edges.end(), rhs_edges.begin(), rhs_edges.end(), std::back_inserter(result_edges));
	disjoint_set_data_structure components(size);
	std::vector<unsigned> missing_edges(size, 2);
	for (const auto& [left, right] : result_edges) {
		components.merge(left, right);
		gsl::at(missing_edges, left)--;
		gsl::at(missing_edges, right)--;
	}
	std::shuffle(lhs_edges.begin(), lhs_edges.end(), rand);
	std::shuffle(rhs_edges.begin(), rhs_edges.end(), rand);
	gsl::span<edge_type> spans[2] {lhs_edges, rhs_edges};
	while (!std::all_of(std::begin(spans), std::end(spans), std::mem_fn(&gsl::span<edge_type>::empty))) {
		for (auto&& span : spans) {
			auto it = std::find_if(span.begin(), span.end(), [&](const edge_type& edge) {
				const auto& [lhs, rhs] = edge;
				return gsl::at(missing_edges, lhs) && gsl::at(missing_edges, rhs) && components.find(lhs) != components.find(rhs);
			});
			if (it != span.end()) {
				const auto& [left, right] = *it;
				result_edges.emplace_back(left, right);
				components.merge(left, right);
				gsl::at(missing_edges, left)--;
				gsl::at(missing_edges, right)--;
			}
			span = span.subspan(it - span.begin());
		}
	}
	std::stack<unsigned, std::vector<unsigned>> s;
	for (unsigned i = 0; i < size; i++) {
		repeat(gsl::at(missing_edges, i), [&] {
			if (!s.empty() && components.find(s.top()) != components.find(i)) {
				const unsigned& left = s.top();
				const unsigned& right = i;
				result_edges.emplace_back(left, right);
				components.merge(left, right);
				gsl::at(missing_edges, left)--;
				gsl::at(missing_edges, right)--;
				s.pop();
			} else {
				s.push(i);
			}
		});
	}
	if (!s.empty()) {
		const unsigned left = s.top();
		s.pop();
		const unsigned& right = s.top();
		result_edges.emplace_back(left, right);
	}
	Ensures(result_edges.size() == lhs.size());
	return to_permutation(result_edges);
}

template<class UniformRandomBitGenerator>
auto path_merger<UniformRandomBitGenerator>::to_edges(const permutation& perm) const -> edge_vector {
	Expects(perm.size() > 0);
	edge_vector result {std::minmax({perm.front(), perm.back()})};
	result.reserve(perm.size());
	std::transform(std::next(perm.begin()), perm.end(), perm.begin(), std::inserter(result, result.end()), [](unsigned lhs, unsigned rhs) {
		return std::minmax({lhs, rhs});
	});
	std::sort(result.begin(), result.end());
	Ensures(result.size() == perm.size());
	return result;
}

template<class UniformRandomBitGenerator>
permutation path_merger<UniformRandomBitGenerator>::to_permutation(const edge_vector& edges) const {
	Expects(edges.size() > 0);
	const std::size_t size = edges.size();
	std::vector<std::vector<unsigned>> graph(size);
	for (auto&& node : graph) {
		node.reserve(2);
	}
	for (const auto& [lhs, rhs] : edges) {
		gsl::at(graph, lhs).push_back(rhs);
		gsl::at(graph, rhs).push_back(lhs);
	}
	permutation result(size);
	unsigned prev = 0;
	unsigned cur = graph.front().front();
	std::for_each(std::next(result.begin()), result.end(), [&](unsigned& element) {
		element = cur;
		const auto& node = gsl::at(graph, cur);
		unsigned next = node.front() == prev ? node.back() : node.front();
		prev = cur;
		cur = next;
	});
	Ensures(result.size() == edges.size());
	return result;
}

#endif

max_stars_repo_name: SirEmentaler/Eugenics-Wars

max_issues_repo_name: SirEmentaler/Eugenics-Wars

max_forks_repo_name: SirEmentaler/Eugenics-Wars


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 utility_77d83a8d_bf54_4fbe_ab3c_b06d0be9e40e_h
#define utility_77d83a8d_bf54_4fbe_ab3c_b06d0be9e40e_h

#include <gslib/std.h>
#include <gslib/type.h>

/*
 * This file include a series of intersection operations.
 * Here are all the details:
 *                      point       ratio
 * linear - linear:     available   available
 * linear - quad:       none        available
 * linear - cubic:      none        available
 * quad - quad:         none        available
 * quad - cubic:        available   none
 * cubic - cubic:       available   none
 * To convert from ratio to point was quite simple, this was why I suggest ratio first.
 * other cases limited by the algorithm where the point was the first result, I hope you
 * notice this and use certain methods to get your result, so I didn't offer an altered
 * way, but you can simply achieve that by calling the reparameterize functions.
 * Caution that the precision of intersection may cause loss of roots,
 * set tolerance around 0.2 - 0.3 might be good.
 */

__gslib_begin__

gs_export extern void linear_interpolate(vec2 c[], const vec2& p1, const vec2& p2, int step);
gs_export extern void quadratic_interpolate(vec2 c[], const vec2& p1, const vec2& p2, const vec2& p3, int step);
gs_export extern void cubic_interpolate(vec2 c[], const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4, int step);
gs_export extern void hermite_interpolate(vec2 c[], const vec2& p1, const vec2& p2, const vec2& t1, const vec2& t2, int step);
gs_export extern float linear_reparameterize(const vec2& p1, const vec2& p2, const vec2& p);
gs_export extern float quad_reparameterize(const vec3 para[2], const vec2& p);
gs_export extern float cubic_reparameterize(const vec4 para[2], const vec2& p);
gs_export extern float best_quad_reparameterize(const vec3 para[2], const vec2& p);
gs_export extern float best_cubic_reparameterize(const vec4 para[2], const vec2& p);
gs_export extern bool is_concave_angle(const vec2& p1, const vec2& p2, const vec2& p3);
gs_export extern bool is_concave_angle(const vec2& p1, const vec2& p2, const vec2& p3, bool cw);
gs_export extern bool is_approx_line(const vec2& p1, const vec2& p2, const vec2& p3, float tolerance);
gs_export extern void get_linear_coefficient(vec3& coef, const vec2& p, const vec2& d);
gs_export extern void get_linear_parameter_equation(vec2 para[2], const vec2& a, const vec2& b);
gs_export extern void get_quad_parameter_equation(vec3 para[2], const vec2& a, const vec2& b, const vec2& c);
gs_export extern void get_cubic_parameter_equation(vec4 para[2], const vec2& a, const vec2& b, const vec2& c, const vec2& d);
gs_export extern void get_first_derivate_factor(vec2 df1[2], const vec3 para[2]);
gs_export extern void get_first_derivate_factor(vec3 df1[2], const vec4 para[2]);
gs_export extern void get_second_derivate_factor(vec2 df2[2], const vec4 para[2]);
gs_export extern void get_first_derivate_factor(vec3 df1[2], const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4);
gs_export extern void get_second_derivate_factor(vec2 df2[2], const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4);
gs_export extern void get_first_derivate(vec2& d, const vec3 df1[2], float t);
gs_export extern void get_second_derivate(vec2& d, const vec2 df2[2], float t);
gs_export extern float get_cubic_curvature(const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4, float t);
gs_export extern int solve_univariate_quadratic(float t[2], const vec3& coef);
gs_export extern int solve_univariate_cubic(float t[3], const vec4& coef);
gs_export extern int solve_univariate_quartic(float t[4], const float coef[5]);
gs_export extern int get_cubic_inflection(float t[2], const vec2& a, const vec2& b, const vec2& c, const vec2& d);
gs_export extern int get_cubic_inflection(float t[2], const vec3 ff[2], const vec2 sf[2]);
gs_export extern int get_quad_extrema(float t[], const vec2& ff);
gs_export extern int get_cubic_extrema(float t[2], const vec3& ff);
gs_export extern void get_quad_bound_box(rectf& rc, const vec2& p1, const vec2& p2, const vec2& p3);
gs_export extern void get_cubic_bound_box(rectf& rc, const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4);
gs_export extern int intersection_quad_linear(float t[2], const vec3 quad[2], const vec3& linear);
gs_export extern int intersection_cubic_linear(float t[3], const vec4 cubic[2], const vec3& linear);
gs_export extern int intersection_quad_quad(float ts[4][2], const vec3 quad1[2], const vec3 quad2[2]);
gs_export extern void intersectp_linear_linear(vec2& ip, const vec2& p1, const vec2& p2, const vec2& d1, const vec2& d2);
gs_export extern int intersectp_cubic_quad(vec2 ip[6], const vec2 cp1[4], const vec2 cp2[3], float tolerance);
gs_export extern int intersectp_cubic_cubic(vec2 ip[9], const vec2 cp1[4], const vec2 cp2[4], float tolerance);
gs_export extern bool get_self_intersection(float ts[2], const vec2& a, const vec2& b, const vec2& c, const vec2& d);
gs_export extern void split_quad_bezier(vec2 c[5], const vec2 p[3], float t);
gs_export extern void split_quad_bezier(vec2 c[7], const vec2 p[3], float t1, float t2);
gs_export extern void split_cubic_bezier(vec2 c[7], const vec2 p[4], float t);
gs_export extern void split_cubic_bezier(vec2 c[10], const vec2 p[4], float t1, float t2);
gs_export extern void offset_cubic_bezier(vec2 c[4], const vec2 p[4], float d);
gs_export extern float quad_bezier_length(const vec2 cp[3]);
gs_export extern float cubic_bezier_length(const vec2 cp[4], float tolerance);
gs_export extern int cubic_to_quad_bezier(vector<vec2>& quadctl, const vec2 cp[4], float tolerance);
gs_export extern float quad_control_length(const vec2& a, const vec2& b, const vec2& c);
gs_export extern float cubic_control_length(const vec2& a, const vec2& b, const vec2& c, const vec2& d);
gs_export extern float point_line_distance(const vec2& p, const vec2& p1, const vec2& p2);
gs_export extern int get_interpolate_step(const vec2& a, const vec2& b, const vec2& c, float step_len = -1.f);
gs_export extern int get_interpolate_step(const vec2& a, const vec2& b, const vec2& c, const vec2& d, float step_len = -1.f);
gs_export extern int get_rough_interpolate_step(const vec2& a, const vec2& b, const vec2& c);
gs_export extern int get_rough_interpolate_step(const vec2& a, const vec2& b, const vec2& c, const vec2& d);
gs_export extern float get_cubic_klmcoords(vec3 m[4], const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4,
    const vec2& ac1, const vec2& ac2, const vec2& ac3, const vec2& ac4  /* absolute coord for error control */
    );
gs_export extern float get_include_angle(const vec2& d1, const vec2& d2);
gs_export extern float get_include_angle(const vec2& a, const vec2& b, const vec2& c);
gs_export extern float get_rotate_angle(const vec2& d1, const vec2& d2);
gs_export extern float get_rotate_angle(const vec2& a, const vec2& b, const vec2& c);
gs_export extern void get_reduce_point(vec2& p, const vec2& a, const vec2& b, const vec2& c, const vec2& d);
gs_export extern bool point_in_triangle(const vec2& p, const vec2& p1, const vec2& p2, const vec2& p3);
gs_export extern int point_in_polygon(const vec2& p, const vector<vec2>& poly);   /* 0: outside; 1 : inside; -1 : coincide */
gs_export extern float get_triangle_area(const vec2& p1, const vec2& p2, const vec2& p3);
gs_export extern vec3 get_barycentric_coords(const vec2& p, const vec2& p1, const vec2& p2, const vec2& p3);
gs_export extern bool clip_line_rect(float t[2], const vec2& p1, const vec2& p2, const rect& clipbox);
gs_export extern bool clip_line_rectf(float t[2], const vec2& p1, const vec2& p2, const rectf& clipbox);
gs_export extern void trace_quad_strip(const vec2 cp[], int size);
gs_export extern void trace_cubic_strip(const vec2 cp[], int size);

inline bool is_isosigned(float a, float b) { return (int)((*(uint*)&a ^ *(uint*)&b)) >= 0; }
inline bool is_isosigned(double a, double b) { return (long long)((*(unsigned long long*)&a ^ *(unsigned long long*)&b)) >= 0; }
inline bool fuzzy_zero(float f) { return (f <= 1e-4f && f >= -1e-4f); }
inline bool fuzzy_zero(float f, float tol) { return f <= tol && f >= -tol; }
inline bool fuzzy_between(float a, float b, float c, float tol) { return b > a - tol && b < c + tol; }
inline bool fuzzy_greater_inclusive(float m, float n, float tol) { return m > (n - tol); }
inline bool fuzzy_less_inclusive(float m, float n, float tol) { return m < (n + tol); }
inline bool fuzzy_greater_exclusive(float m, float n, float tol) { return m > (n + tol); }
inline bool fuzzy_less_exclusive(float m, float n, float tol) { return m < (n - tol); }
inline bool is_parallel(const vec2& a, const vec2& b) { return fuzzy_zero(a.ccw(b)); }

inline vec2& eval_linear(vec2& v, const vec2 para[2], float t)
{
    vec2 sv(t, 1.f);
    v.x = para[0].dot(sv);
    v.y = para[1].dot(sv);
    return v;
}

inline vec2& eval_quad(vec2& v, const vec3 para[2], float t)
{
    vec3 sv(t * t, t , 1.f);
    v.x = para[0].dot(sv);
    v.y = para[1].dot(sv);
    return v;
}

inline vec2& eval_cubic(vec2& v, const vec4 para[2], float t)
{
    vec4 sv(t * t * t, t * t, t, 1.f);
    v.x = para[0].dot(sv);
    v.y = para[1].dot(sv);
    return v;
}

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/* linalg/hermtd.c
 *
 * Copyright (C) 2001, 2007, 2009 Brian Gough
 *
 * 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.
 */

/* Factorise a hermitian matrix A into
 *
 * A = U T U'
 *
 * where U is unitary and T is real symmetric tridiagonal.  Only the
 * diagonal and lower triangular part of A is referenced and modified.
 *
 * On exit, T is stored in the diagonal and first subdiagonal of
 * A. Since T is symmetric the upper diagonal is not stored.
 *
 * U is stored as a packed set of Householder transformations in the
 * lower triangular part of the input matrix below the first subdiagonal.
 *
 * The full matrix for U can be obtained as the product
 *
 *       U = U_N ... U_2 U_1
 *
 * where
 *
 *       U_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [0, ..., 0, 1, A(i+2,i), A(i+3,i), ... , A(N,i)]
 *
 * This storage scheme is the same as in LAPACK.  See LAPACK's
 * chetd2.f for details.
 *
 * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex_math.h>

#include <gsl/gsl_linalg.h>

int
gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("hermitian tridiagonal decomposition requires square matrix",
                 GSL_ENOTSQR);
    }
  else if (tau->size + 1 != A->size1)
    {
      GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i;

      const gsl_complex zero = gsl_complex_rect (0.0, 0.0);
      const gsl_complex one = gsl_complex_rect (1.0, 0.0);
      const gsl_complex neg_one = gsl_complex_rect (-1.0, 0.0);

      for (i = 0 ; i < N - 1; i++)
        {
          gsl_vector_complex_view c = gsl_matrix_complex_column (A, i);
          gsl_vector_complex_view v = gsl_vector_complex_subvector (&c.vector, i + 1, N - (i + 1));
          gsl_complex tau_i = gsl_linalg_complex_householder_transform (&v.vector);

          /* Apply the transformation H^T A H to the remaining columns */

          if ((i + 1) < (N - 1)
              && !(GSL_REAL(tau_i) == 0.0 && GSL_IMAG(tau_i) == 0.0))
            {
              gsl_matrix_complex_view m =
                gsl_matrix_complex_submatrix (A, i + 1, i + 1,
                                              N - (i+1), N - (i+1));
              gsl_complex ei = gsl_vector_complex_get(&v.vector, 0);
              gsl_vector_complex_view x = gsl_vector_complex_subvector (tau, i, N-(i+1));
              gsl_vector_complex_set (&v.vector, 0, one);

              /* x = tau * A * v */
              gsl_blas_zhemv (CblasLower, tau_i, &m.matrix, &v.vector, zero, &x.vector);

              /* w = x - (1/2) tau * (x' * v) * v  */
              {
                gsl_complex xv, txv, alpha;
                gsl_blas_zdotc(&x.vector, &v.vector, &xv);
                txv = gsl_complex_mul(tau_i, xv);
                alpha = gsl_complex_mul_real(txv, -0.5);
                gsl_blas_zaxpy(alpha, &v.vector, &x.vector);
              }

              /* apply the transformation A = A - v w' - w v' */
              gsl_blas_zher2(CblasLower, neg_one, &v.vector, &x.vector, &m.matrix);

              gsl_vector_complex_set (&v.vector, 0, ei);
            }

          gsl_vector_complex_set (tau, i, tau_i);
        }

      return GSL_SUCCESS;
    }
}

/*  Form the orthogonal matrix U from the packed QR matrix */

int
gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
                          const gsl_vector_complex * tau,
                          gsl_matrix_complex * U,
                          gsl_vector * diag,
                          gsl_vector * sdiag)
{
  if (A->size1 !=  A->size2)
    {
      GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
    }
  else if (tau->size + 1 != A->size1)
    {
      GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
    }
  else if (U->size1 != A->size1 || U->size2 != A->size1)
    {
      GSL_ERROR ("size of U must match size of A", GSL_EBADLEN);
    }
  else if (diag->size != A->size1)
    {
      GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
    }
  else if (sdiag->size + 1 != A->size1)
    {
      GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;

      size_t i;

      /* Initialize U to the identity */

      gsl_matrix_complex_set_identity (U);

      for (i = N - 1; i-- > 0;)
        {
          gsl_complex ti = gsl_vector_complex_get (tau, i);

          gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (A, i);

          gsl_vector_complex_const_view h =
            gsl_vector_complex_const_subvector (&c.vector, i + 1, N - (i+1));

          gsl_matrix_complex_view m =
            gsl_matrix_complex_submatrix (U, i + 1, i + 1, N-(i+1), N-(i+1));

          gsl_linalg_complex_householder_hm (ti, &h.vector, &m.matrix);
        }

      /* Copy diagonal into diag */

      for (i = 0; i < N; i++)
        {
          gsl_complex Aii = gsl_matrix_complex_get (A, i, i);
          gsl_vector_set (diag, i, GSL_REAL(Aii));
        }

      /* Copy subdiagonal into sdiag */

      for (i = 0; i < N - 1; i++)
        {
          gsl_complex Aji = gsl_matrix_complex_get (A, i+1, i);
          gsl_vector_set (sdiag, i, GSL_REAL(Aji));
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
                            gsl_vector * diag,
                            gsl_vector * sdiag)
{
  if (A->size1 !=  A->size2)
    {
      GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
    }
  else if (diag->size != A->size1)
    {
      GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
    }
  else if (sdiag->size + 1 != A->size1)
    {
      GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;

      size_t i;

      /* Copy diagonal into diag */

      for (i = 0; i < N; i++)
        {
          gsl_complex Aii = gsl_matrix_complex_get (A, i, i);
          gsl_vector_set (diag, i, GSL_REAL(Aii));
        }

      /* Copy subdiagonal into sd */

      for (i = 0; i < N - 1; i++)
        {
          gsl_complex Aji = gsl_matrix_complex_get (A, i+1, i);
          gsl_vector_set (sdiag, i, GSL_REAL(Aji));
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 * @file piecewise_regimes.c
 * @author Carl Boettiger <cboettig@gmail.com>
 * @date 17 May 2010
 *
 * @section LICENSE
 *
 * GPLv3
 *
 * @section DESCRIPTION
 *
 * Code for computing the likelihood of a multitype OU process on a phylogenetic tree
 *
 */

#include "likelihood.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include "optimizers.h"
#include "mvn.h"

/**
* @brief Function which wraps general format used by siman and multimin
*        into the likelihood format of calc_lik
*
* @param v vector of values being optimized by the algorithm
* @param params specifies the tree and observed data
*
* @return minus log likelihood
*/
double optim_func (const gsl_vector *v, void *params)
{
  tree * mytree = (tree *) params;
  int i, n_regimes = mytree->n_regimes;
// Force root to match Theta of regime assigned to it, should also reduce vector size!
  mytree->Xo[0] = mytree->theta[mytree->regimes[0]];
 // mytree->Xo[0] = gsl_vector_get(v, 0);
  for(i = 0; i < n_regimes; i++){
    mytree->alpha[i] = v->data[i];
    mytree->theta[i] = v->data[n_regimes+i];
    mytree->sigma[i] = v->data[2*n_regimes+i];

    if( (mytree->alpha[i] < 0) | (mytree->sigma[i] < 0) ){ return GSL_POSINF; }
    if(  mytree->theta[i] < 0 ){ return GSL_POSINF; }
  }

  double llik = 0;
    calc_lik(mytree->Xo, mytree->alpha, mytree->theta, mytree->sigma,
         mytree->regimes, mytree->ancestor, mytree->branch_length,
             mytree->traits, &(mytree->n_nodes), mytree->lca_matrix, &llik);
    return -llik;
}

/** Fit the model by maximum likelihood, called by R */
void fit_model(double * Xo,
  double * alpha,
  double * theta,
  double * sigma,
  int * regimes,
  int * ancestor,
  double * branch_length,
  double * traits,
  int * n_nodes,
  int * n_regimes,
  double * llik,
  int * use_siman)
{

//  Perhaps don't want to see gsl_errors once we're handling in R?
//  Still want better/smarter behavior of the optimizer
  gsl_set_error_handler_off ();

    /** First, initialize the tree structure to facilitate passing all these elements */
  int i,j;
  tree * mytree = (tree  *) malloc(sizeof(tree));
  mytree->Xo = Xo;
  mytree->alpha = alpha;
  mytree->theta = theta;
  mytree->sigma = sigma;
  mytree->regimes = regimes;
  mytree->ancestor = ancestor;
  mytree->branch_length = branch_length;
  mytree->traits = traits;
  mytree->n_nodes = *n_nodes;
  mytree->n_regimes = *n_regimes;

    /** Second, save time by storing the least common ancestor matrix */
  /* Save time by calculating least common ancestor ahead of time.
   * Though this does the calc for all nodes, only tip ones are used.
   * Current get_lca algorithm is only designed for tips anyway.  */
  double sep; // separation time, not actually used
  mytree->lca_matrix = (int *) malloc(gsl_pow_2(*n_nodes) * sizeof(int) );
  for(i=0; i < *n_nodes; i++){
    for(j=0; j < *n_nodes; j++){
      mytree->lca_matrix[*n_nodes*i+j] = get_lca(i,j, *n_nodes, ancestor, branch_length, &sep);
    }
  }

    /** Then initialize the vector of parameters to be searched */
  gsl_vector *x = gsl_vector_alloc(3 * *n_regimes);
//  gsl_vector_set(x, 0, *Xo);  /* Forces Xo to match theta (?) */
  for(i = 0; i < *n_regimes; i++){
    gsl_vector_set(x, i, alpha[i]);
    gsl_vector_set(x, *n_regimes+i, theta[i]);
    gsl_vector_set(x, 2 * *n_regimes+i, mytree->sigma[i]);
  }

    /** Select and call a likelihood optimization routine
     *  given the tree (with data) and parameters */
  gsl_rng * rng = gsl_rng_alloc(gsl_rng_default);
  if(*use_siman)
    *llik = siman(x, mytree, rng);
  else
    *llik = multimin(x, mytree);

  gsl_rng_free(rng);
  gsl_vector_free(x);
  free(mytree->lca_matrix);
  free(mytree);
}

/* Externally exported function to determine LCA matrix */
/**
* @brief Determine the last common ancestor matrix of the tips,
* used in the likelihood calculation
*
* @param ancestor ancestor of the specified nide
* @param branch_length double branch lengths below each of the nodes
* @param n_nodes: integer number of nodes (including tips) in the tree
* @param lca_matrix must be an integer array of size n_nodes^2
*/
void calc_lca(int * ancestor, double * branch_length,
              int * n_nodes, int * lca_matrix)
{
  int i,j;
  /* Save time by calculating least common ancestor ahead of time.
   * Though this does the calc for all nodes, only tip ones are used.
   * Current get_lca algorithm is only designed for tips anyway.  */
  double sep; // separation time, not actually used
  for(i=0; i < *n_nodes; i++){
    for(j=0; j < *n_nodes; j++){
      lca_matrix[*n_nodes*i+j] = get_lca(i,j, *n_nodes, ancestor,
                                               branch_length, &sep);
    }
  }

}

/** Simulate under the model.  Called by R */
void simulate_model (double * Xo, double * alpha, double * theta,
                     double * sigma, int * regimes, int * ancestor,
                     double * branch_length, double * traits, int * n_nodes,
                     int * n_regimes, double * llik, double * seed)
{

  /* Turn off likelihood calcuation on the simulated data for speed */
  int return_likelihood = 0; // consider making an option at the .C level

  int i,j;
  tree * mytree = (tree  *) malloc(sizeof(tree));
  mytree->Xo = Xo;
  mytree->alpha = alpha;
  mytree->theta = theta;
  mytree->sigma = sigma;
  mytree->regimes = regimes;
  mytree->ancestor = ancestor;
  mytree->branch_length = branch_length;
  mytree->traits = traits;
  mytree->n_nodes = *n_nodes;
  mytree->n_regimes = *n_regimes;

  /* Save time by calculating least common ancestor ahead of time. */
  double sep; // separation time, not actually used
  mytree->lca_matrix = (int *) malloc(gsl_pow_2(*n_nodes) * sizeof(int) );
  for(i=0; i < *n_nodes; i++){
    for(j=0; j < *n_nodes; j++){
      mytree->lca_matrix[*n_nodes*i+j] =
                get_lca(i,j, *n_nodes, ancestor, branch_length, &sep);
    }
  }

    /* Allocate the parameter vector */
  gsl_vector *x = gsl_vector_alloc(1 + 3 * *n_regimes);
  gsl_vector_set(x, 0, *Xo);
  for(i = 0; i < *n_regimes; i++){
    gsl_vector_set(x, 1+i, alpha[i]);
    gsl_vector_set(x, 1+*n_regimes+i, theta[i]);
    gsl_vector_set(x, 1+2 * *n_regimes+i, mytree->sigma[i]);
  }

  gsl_rng * rng = gsl_rng_alloc(gsl_rng_default);
  gsl_rng_set(rng, *seed);

   /* These are the only lines that differ from the fit_models wrapper */
  simulate (rng, mytree);

  if(return_likelihood){
    calc_lik(mytree->Xo, mytree->alpha, mytree->theta, mytree->sigma,
             mytree->regimes, mytree->ancestor, mytree->branch_length,
                   mytree->traits, &(mytree->n_nodes), mytree->lca_matrix, llik);
  }

//  for(i=0; i< mytree->n_regimes; i++) printf("%g %g %g\n", mytree->alpha[i], mytree->theta[i], mytree->sigma[i] );
//  printf("sim llik: %g\n", *llik);

  gsl_rng_free(rng);
  gsl_vector_free(x);
  free(mytree->lca_matrix);
  free(mytree);
}

/** Test code by running in pure C */
int main(void)
{
  /** Define the tree and the trait data (Anoles tree from "ouch" R package) */
  int n_nodes = 45;

  double branch_length[] = { 0.0, 12./38, 20./38,  2./38,  2./38,  4./38,
               8./38,  5./38,  5./38,  5./38,  5./38, 10./38,
               9./38,  4./38,  8./38,  2./38, 20./38,  2./38,
               4./38,  2./38, 1./38,  2./38, 26./38, 4./38,
               2./38, 2./38,  2./38,  2./38, 15./38, 10./38,
               10./38, 10./38, 10./38, 16./38, 12./38,  4./38,
               2./38,  2./38, 10./38,  8./38,  2./38,  1./38,
               1./38,  2./38,  2./38};

/* Almost a star tree */
/*  double branch_length[] = { 0.000, .001, .001,  .001,  .001,
			     0.001, .001, .001,  .001,  .001,
			     0.001, .001, .001,  .001,  .001,
			     0.001, .001, .001,  .001,  .001,
			     0.001, .001, .999,  .999,  .999,
			     0.999, .999, .999,  .999,  .999,
			     0.999, .999, .999,  .999,  .999,
			     0.999, .999, .999,  .999,  .999,
			     0.999, .999, .999,  .999,  .999 };
*/

  int ancestor[] = {-1,  0,  1,  2,  3,  2,  0,  6,  7,  8,  9,
                     8,  7, 12, 13,14,  6, 16, 17, 18, 19, 18, 1,
                     3,  4,  4,  5,  5,  9, 10, 10, 11, 11, 12,
                    13, 14, 15, 15, 16, 17, 19, 20, 20, 21, 21};

  double traits[] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                     0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                     0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                     0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                     0.000000, 0.000000, 2.602690, 2.660260, 2.660260,
                     2.653242, 2.674149, 2.701361, 7.161247, 13.299534,
                     3.328627, 3.353407, 3.360375, 3.049273, 2.906901,
                     2.980619, 2.933857, 3.975530, -3.104587, 3.346389,
                     12.928524, -12.939162, 7.990720, 8.058707, -5.068053};

/*
  int regimes[] =  {1, 1, 2, 2, 2,
                    2, 1, 1, 0, 0,
                    0, 0, 1, 1, 1,
                    1, 1, 1, 1, 1,
                    1, 1, 2, 2, 2,
                    2, 2, 2, 0, 0,
                    0, 0, 0, 1, 1,
                    1, 1, 0, 1, 1,
                    1, 1, 1, 1, 1};
*/

  int regimes[] =   {0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0};

  double Xo = 3;

  double alpha[3] = {2.6, 2.6, 2.6};
  double theta[3] = {3., 3.,3};
  double sigma[3] = {0.2, 0.2, 0.2 };
  int n_regimes = 3;

  double llik = 0;

  /* Initialize the lca matrix */
  int * lca = (int *) calloc(n_nodes*n_nodes, sizeof(int));

  /* test the pure likelihood function first */
  calc_lca(ancestor, branch_length, &n_nodes, lca);

  calc_lik(&Xo, alpha, theta, sigma, regimes, ancestor,
           branch_length, traits, &n_nodes, lca, &llik);
  printf("log likelihood: %g\n", llik);

  /* We'll use the faster Nelder-Mead algorithm for the check */
  int use_siman = 0;
  /* Test the fitting routine */

  fit_model(&Xo, alpha, theta, sigma, regimes, ancestor, branch_length,
            traits, &n_nodes, &n_regimes, &llik, &use_siman);
  printf("Xo = %g\n", Xo);
  printf("alphas: %g %g %g\n", alpha[0], alpha[1], alpha[2]);
  printf("thetas: %g %g %g\n", theta[0], theta[1], theta[2]);
  printf("sigmas: %g %g %g\n", sigma[0], sigma[1], sigma[2]);
  printf("log likelihood: %g\n", llik);

  /* Test the simulation code */
/*
  double seed = 1.0;
  printf("seed: %lo\n", (unsigned long int) seed);
  simulate_model(&Xo, alpha, theta, sigma, regimes, ancestor, branch_length,
                 traits, &n_nodes, &n_regimes, &llik, &seed);
  int i;
  printf("simulated traits:\n");
  for(i=0; i < 45; i++)
    printf("%g ", traits[i]);
  printf("\n");
*/

  free(lca);
  return 0;
}

max_stars_repo_name: cboettig/wrightscape

max_issues_repo_name: cboettig/wrightscape

max_forks_repo_name: cboettig/wrightscape


text:
/* Cholesky Decomposition
 *
 * Copyright (C) 2000  Thomas Walter
 *
 * 03 May 2000: Modified for GSL by Brian Gough
 * 29 Jul 2005: Additions by Gerard Jungman
 * Copyright (C) 2000,2001, 2002, 2003, 2005, 2007 Brian Gough, Gerard Jungman
 *
 * This 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, or (at your option) any
 * later version.
 *
 * This source 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.
 */

/*
 * Cholesky decomposition of a symmetrix positive definite matrix.
 * This is useful to solve the matrix arising in
 *    periodic cubic splines
 *    approximating splines
 *
 * This algorithm does:
 *   A = L * L'
 * with
 *   L  := lower left triangle matrix
 *   L' := the transposed form of L.
 *
 */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

static inline
double
quiet_sqrt (double x)
     /* avoids runtime error, for checking matrix for positive definiteness */
{
  return (x >= 0) ? sqrt(x) : GSL_NAN;
}

int
gsl_linalg_cholesky_decomp (gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t i,j,k;
      int status = 0;

      /* Do the first 2 rows explicitly.  It is simple, and faster.  And
       * one can return if the matrix has only 1 or 2 rows.
       */

      double A_00 = gsl_matrix_get (A, 0, 0);

      double L_00 = quiet_sqrt(A_00);

      if (A_00 <= 0)
        {
          status = GSL_EDOM ;
        }

      gsl_matrix_set (A, 0, 0, L_00);

      if (M > 1)
        {
          double A_10 = gsl_matrix_get (A, 1, 0);
          double A_11 = gsl_matrix_get (A, 1, 1);

          double L_10 = A_10 / L_00;
          double diag = A_11 - L_10 * L_10;
          double L_11 = quiet_sqrt(diag);

          if (diag <= 0)
            {
              status = GSL_EDOM;
            }

          gsl_matrix_set (A, 1, 0, L_10);
          gsl_matrix_set (A, 1, 1, L_11);
        }

      for (k = 2; k < M; k++)
        {
          double A_kk = gsl_matrix_get (A, k, k);

          for (i = 0; i < k; i++)
            {
              double sum = 0;

              double A_ki = gsl_matrix_get (A, k, i);
              double A_ii = gsl_matrix_get (A, i, i);

              gsl_vector_view ci = gsl_matrix_row (A, i);
              gsl_vector_view ck = gsl_matrix_row (A, k);

              if (i > 0) {
                gsl_vector_view di = gsl_vector_subvector(&ci.vector, 0, i);
                gsl_vector_view dk = gsl_vector_subvector(&ck.vector, 0, i);

                gsl_blas_ddot (&di.vector, &dk.vector, &sum);
              }

              A_ki = (A_ki - sum) / A_ii;
              gsl_matrix_set (A, k, i, A_ki);
            }

          {
            gsl_vector_view ck = gsl_matrix_row (A, k);
            gsl_vector_view dk = gsl_vector_subvector (&ck.vector, 0, k);

            double sum = gsl_blas_dnrm2 (&dk.vector);
            double diag = A_kk - sum * sum;

            double L_kk = quiet_sqrt(diag);

            if (diag <= 0)
              {
                status = GSL_EDOM;
              }

            gsl_matrix_set (A, k, k, L_kk);
          }
        }

      /* Now copy the transposed lower triangle to the upper triangle,
       * the diagonal is common.
       */

      for (i = 1; i < M; i++)
        {
          for (j = 0; j < i; j++)
            {
              double A_ij = gsl_matrix_get (A, i, j);
              gsl_matrix_set (A, j, i, A_ij);
            }
        }

      if (status == GSL_EDOM)
        {
          GSL_ERROR ("matrix must be positive definite", GSL_EDOM);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
                           const gsl_vector * b,
                           gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve for c using forward-substitution, L c = b */

      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
                         gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Solve for c using forward-substitution, L c = b */

      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_cholesky_invert()
  Compute the inverse of a symmetric positive definite matrix in
Cholesky form.

Inputs: LLT - matrix in cholesky form on input
              A^{-1} = L^{-t} L^{-1} on output

Return: success or error
*/

int
gsl_linalg_cholesky_invert(gsl_matrix * LLT)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      size_t N = LLT->size1;
      size_t i, j;
      double sum;
      gsl_vector_view v1, v2;

      /* invert the lower triangle of LLT */
      for (i = 0; i < N; ++i)
        {
          double ajj;

          j = N - i - 1;

          gsl_matrix_set(LLT, j, j, 1.0 / gsl_matrix_get(LLT, j, j));
          ajj = -gsl_matrix_get(LLT, j, j);

          if (j < N - 1)
            {
              gsl_matrix_view m;

              m = gsl_matrix_submatrix(LLT, j + 1, j + 1,
                                       N - j - 1, N - j - 1);
              v1 = gsl_matrix_subcolumn(LLT, j, j + 1, N - j - 1);

              gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit,
                             &m.matrix, &v1.vector);

              gsl_blas_dscal(ajj, &v1.vector);
            }
        } /* for (i = 0; i < N; ++i) */

      /*
       * The lower triangle of LLT now contains L^{-1}. Now compute
       * A^{-1} = L^{-t} L^{-1}
       *
       * The (ij) element of A^{-1} is column i of L^{-1} dotted into
       * column j of L^{-1}
       */

      for (i = 0; i < N; ++i)
        {
          for (j = i + 1; j < N; ++j)
            {
              v1 = gsl_matrix_subcolumn(LLT, i, j, N - j);
              v2 = gsl_matrix_subcolumn(LLT, j, j, N - j);

              /* compute Ainv_{ij} = sum_k Linv_{ki} Linv_{kj} */
              gsl_blas_ddot(&v1.vector, &v2.vector, &sum);

              /* store in upper triangle */
              gsl_matrix_set(LLT, i, j, sum);
            }

          /* now compute the diagonal element */
          v1 = gsl_matrix_subcolumn(LLT, i, i, N - i);
          gsl_blas_ddot(&v1.vector, &v1.vector, &sum);
          gsl_matrix_set(LLT, i, i, sum);
        }

      /* copy the transposed upper triangle to the lower triangle */

      for (j = 1; j < N; j++)
        {
          for (i = 0; i < j; i++)
            {
              double A_ij = gsl_matrix_get (LLT, i, j);
              gsl_matrix_set (LLT, j, i, A_ij);
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_cholesky_invert() */

int
gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
{
  const size_t N = A->size1;
  size_t i, j;

  /* initial Cholesky */
  int stat_chol = gsl_linalg_cholesky_decomp(A);

  if(stat_chol == GSL_SUCCESS)
  {
    /* calculate D from diagonal part of initial Cholesky */
    for(i = 0; i < N; ++i)
    {
      const double C_ii = gsl_matrix_get(A, i, i);
      gsl_vector_set(D, i, C_ii*C_ii);
    }

    /* multiply initial Cholesky by 1/sqrt(D) on the right */
    for(i = 0; i < N; ++i)
    {
      for(j = 0; j < N; ++j)
      {
        gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
      }
    }

    /* Because the initial Cholesky contained both L and transpose(L),
       the result of the multiplication is not symmetric anymore;
       but the lower triangle _is_ correct. Therefore we reflect
       it to the upper triangle and declare victory.
       */
    for(i = 0; i < N; ++i)
      for(j = i + 1; j < N; ++j)
        gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
  }

  return stat_chol;
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#ifdef HAVE_PETSC
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#include <mpi.h>
#include <petsc.h>

#include "common_c.h"
#include "FCMangle.h"

#define fillsparsecpetscs FortranCInterface_GLOBAL_(fillsparsecpetscs, FILLSPARSECPETSCS)
#define fillsparsecpetscc FortranCInterface_GLOBAL_(fillsparsecpetscc, FILLSPARSECPETSCC)

#define COLMAJ2D(row,col,numrow) (row-1)+(col-1)*numrow
#define COLMAJ3D(a,b,c,amax,bmax,cmax) (a-1)+amax*((b-1)+bmax*(c-1))
#define ROWMAJ2D_ONE(row,col,numcol) (row-1)*numcol+(col-1)
typedef long long int gcorp_t;

void fillsparsecpetscs(gcorp_t* ieng, double* EGmass, Mat* lhsP)
{
        int npro = propar.npro;
        int nshl = shpdat.nshl;
        double* mb = (double*) malloc(sizeof(double)*nshl*nshl); //block to insert
        int e,i,j,aa; //following along with fillsparse.f
        PetscInt* locat = (PetscInt*) malloc(sizeof(PetscInt)*nshl);
        for(e=0;e<npro;e++)
        {
         for(aa=0;aa<nshl;aa++) locat[aa]=ieng[e+npro*aa]-1;
//         for(aa=0;aa<nshl;aa++) assert(locat[aa]>=0);
         for (i=0; i<nshl; i++)  {   // fill up Ke with respective egmass
           for (j=0; j<nshl; j++)  {
            mb[nshl*i + j] = EGmass[e + npro*(i + nshl*j)];
           }
         }
         //MatSetValuesBlocked(*lhsP, nshl , locat, nshl, locat, mb, ADD_VALUES);
         PetscInt petsc_nshl;
         petsc_nshl = (PetscInt) nshl;
         MatSetValues(*lhsP, petsc_nshl , locat, petsc_nshl, locat, mb, ADD_VALUES);
        }
        free(mb);
	free(locat);
}
void fillsparsecpetscc(gcorp_t* ieng, double* EGmass, Mat* lhsP)
{
        int npro = propar.npro;
        int nshl = shpdat.nshl;
        int nedof = conpar.nedof;
        double* mb = (double*) malloc(sizeof(double)*nedof*nedof); //block to insert
        int e,i,j,aa; //following along with fillsparse.f
        //int* locat = (int*) malloc(sizeof(int)*nshl);
        PetscInt* locat = (PetscInt*) malloc(sizeof(PetscInt)*nshl);
        for(e=0;e<npro;e++)
        {
         for(aa=0;aa<nshl;aa++) locat[aa]=ieng[e+npro*aa]-1;
//         for(aa=0;aa<nshl;aa++) assert(locat[aa]>=0);
         for (i=0; i<nedof; i++)  {   /* fill up Ke with respective egmass */
           for (j=0; j<nedof; j++)  {
            mb[nedof*i + j] = EGmass[e + npro*(i + nedof*j)];
           }
         }
         //MatSetValuesBlocked(*lhsP, nshl , locat, nshl, locat, mb, ADD_VALUES);
         PetscInt petsc_nshl;
         petsc_nshl = (PetscInt) nshl;
         MatSetValuesBlocked(*lhsP, petsc_nshl , locat, petsc_nshl, locat, mb, ADD_VALUES);
        }
        free(mb);
	free(locat);
}
#endif


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once
#include <iosfwd>
#include <vector>
#include <algorithm>
#include <string>
#include <cstring>
#include <gsl/gsl>
#include "onnxruntime_config.h"

namespace onnxruntime {
#ifdef __GNUC__
#pragma GCC diagnostic push
#ifdef HAS_NULL_DEREFERENCE
#pragma GCC diagnostic ignored "-Wnull-dereference"
#endif
#endif
class TensorShape {
  // We use negative numbers for unknown symbolic dimension. Each negative
  // number represents a unique symbolic dimension.
 public:
  TensorShape() = default;

  TensorShape(const TensorShape& other) : TensorShape(other.GetDims()) {}
  TensorShape& operator=(const TensorShape& other);

  TensorShape(TensorShape&& other) { operator=(std::move(other)); }
  TensorShape& operator=(TensorShape&& other);

  TensorShape(gsl::span<const int64_t> dims);
  TensorShape(const std::vector<int64_t>& dims) : TensorShape(gsl::make_span(dims)) {}
  TensorShape(const std::initializer_list<int64_t>& dims) : TensorShape(gsl::make_span(dims)) {}
  TensorShape(const int64_t* dimension_sizes, size_t dimension_count) : TensorShape(gsl::span<const int64_t>(dimension_sizes, dimension_count)) {}
  TensorShape(const std::vector<int64_t>& dims, size_t start, size_t end) : TensorShape(gsl::span<const int64_t>(&dims[start], end - start)) {}

  // Create a TensorShape that points to an existing buffer internally. As no copy is made, 'data' must remain valid for the life of the TensorShape
  static const TensorShape FromExistingBuffer(const std::vector<int64_t>& data) {
    return TensorShape(External{}, gsl::span<int64_t>(const_cast<int64_t*>(data.data()), data.size()));
  }

  /**
     Return the dimension specified by <idx>.
  */
  int64_t operator[](size_t idx) const { return values_[idx]; }
  int64_t& operator[](size_t idx) { return values_[idx]; }

  bool operator==(const TensorShape& other) const noexcept { return GetDims() == other.GetDims(); }
  bool operator!=(const TensorShape& other) const noexcept { return GetDims() != other.GetDims(); }

  size_t NumDimensions() const noexcept {
    return values_.size();
  }

  /**
     Copy dims into an array with given size
  */
  void CopyDims(int64_t* dims, size_t num_dims) const {
    memcpy(dims, values_.begin(), sizeof(int64_t) * std::min(num_dims, NumDimensions()));
  }

  /**
     Copy dims from a specific start dim into an array with given size
     `start_dim` is expected to be in the inclusive range [0, NumDimensions() - 1]
     and this function does no checks to ensure that
  */
  void CopyDims(int64_t* dims, size_t start_dim, size_t num_dims) const {
    memcpy(dims, values_.begin() + start_dim, sizeof(int64_t) * std::min(num_dims, NumDimensions() - start_dim));
  }

  /**
     Return underlying vector representation.
  */
  gsl::span<const int64_t> GetDims() const { return values_; }
  std::vector<int64_t> GetDimsAsVector() const { return std::vector<int64_t>(values_.begin(), values_.end()); }

  /**
   * Return the total number of elements. Returns 1 for an empty (rank 0) TensorShape.
   *
   * May return -1
   */
  int64_t Size() const;

  /**
     Return the total number of elements up to the specified dimension.
     If the dimension interval is empty (dimension == 0), return 1.
     @param dimension Return size up to this dimension. Value must be between 0 and this->NumDimensions(), inclusive.
  */
  int64_t SizeToDimension(size_t dimension) const;

  /**
     Return the total number of elements from the specified dimension to the end of the tensor shape.
     If the dimension interval is empty (dimension == this->NumDimensions()), return 1.
     @param dimension Return size from this dimension to the end. Value must be between 0 and this->NumDimensions(),
                      inclusive.
  */
  int64_t SizeFromDimension(size_t dimension) const;

  /**
     Return a new TensorShape of the dimensions from dimstart to dimend.
  */
  TensorShape Slice(size_t dimstart, size_t dimend) const;

  /**
     Return a new TensorShape of the dimensions from dimstart to end.
  */
  TensorShape Slice(size_t dimstart) const { return Slice(dimstart, values_.size()); }

  /**
     output dimensions nicely formatted
  */
  std::string ToString() const;

  /**
     Calculate size between start and end.
     Assumes start and end are between 0 and this->NumDimensions(), inclusive, and that
     start < end.
  */
  int64_t SizeHelper(size_t start, size_t end) const;

  /**
     empty shape or 1D shape (1) is regarded as scalar tensor
  */
  bool IsScalar() const {
    size_t len = values_.size();
    return len == 0 || (len == 1 && values_[0] == 1);
  }

 private:

  struct External {};
  TensorShape(External, gsl::span<int64_t> buffer) : values_{buffer} {}

  void Allocate(size_t size);

  gsl::span<int64_t> values_;
  int64_t small_buffer_[5];
  std::unique_ptr<int64_t[]> allocated_buffer_;

  friend struct ProviderHostImpl; // So that the shared provider interface can access Allocate
};
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif
// operator<< to nicely output to a stream
std::ostream& operator<<(std::ostream& out, const TensorShape& shape);

}  // namespace onnxruntime

max_forks_repo_name: jthelin/onnxruntime


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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>

#include "design.h"

/*This function computes the p-distance between 2 points in D dimensions:
the exponent p is tunable*/

double distance(double *x,double *y,int D,double p){

	double dist = 0.0;
	int i;

	for(i=0;i<D;i++){
		dist += pow(fabs(x[i]-y[i]),p);
	}

	return pow(dist,1.0/p);

}

/*This is the cost function of the problem: it is a sum of all the reciprocal of
the pairs distances, with some tunable exponent lambda; note that the final
1/lambda exponentiation is not performed here!*/

double cost(gsl_matrix *data,int Npoints,int D,double p,double lambda){

	double sum = 0.0;
	int i,j;

	for(i=0;i<Npoints;i++){
		for(j=i+1;j<Npoints;j++){

			//Add the contribution of pair (i,j) to the cost function
			sum += pow(pow(D,1.0/p)/distance(gsl_matrix_ptr(data,i,0),gsl_matrix_ptr(data,j,0),D,p),lambda);

		}
	}

	return (2.0/(Npoints*(Npoints-1)))*sum;

}

/*This computes the cost function in the particular case in which all the points are
equally spaced on the diagonal of the hypercube*/

double diagonalCost(int Npoints,double lambda){

	double sum = 0.0;
	int i,j;

	for(i=0;i<Npoints;i++){
		for(j=i+1;j<Npoints;j++){

			//Add the contribution of pair (i,j) to the cost function
			sum += pow((Npoints-1)*1.0/(j-i),lambda);

		}
	}

	return (2.0/(Npoints*(Npoints-1)))*sum;

}

/*This function computes the variation of the cost when a pair of coordinates is exchanged;
this function also performs an in-place swap of the coordinates. More specifically:
this function swaps coordinate d of points i1 and i2 and returns the variation of the cost
function due to this swap*/

/**/

double swap(gsl_matrix *data,int Npoints,int D,double p,double lambda,int i1,int i2, int d){

	double costBefore,costAfter,temp;
	int i;

	//initialize to 0
	costBefore = costAfter = 0.0;

	/*compute the contribution of points i1 and i2 to the cost function, before swapping;
	sum over all the particles except i1 and i2*/
	for(i=0;i<Npoints;i++){

		if(i!=i1 && i!=i2){
			costBefore += pow(pow(D,1.0/p)/distance(gsl_matrix_ptr(data,i,0),gsl_matrix_ptr(data,i1,0),D,p),lambda) + pow(pow(D,1.0/p)/distance(gsl_matrix_ptr(data,i,0),gsl_matrix_ptr(data,i2,0),D,p),lambda);
		}

	}

	//perform the coordinate swap
	temp = gsl_matrix_get(data,i1,d);
	gsl_matrix_set(data,i1,d,gsl_matrix_get(data,i2,d));
	gsl_matrix_set(data,i2,d,temp);

	/*compute the contribution of points i1 and i2 to the cost function, after swapping;
	sum over all the particles except i1 and i2*/
	for(i=0;i<Npoints;i++){

		if(i!=i1 && i!=i2){
			costAfter += pow(pow(D,1.0/p)/distance(gsl_matrix_ptr(data,i,0),gsl_matrix_ptr(data,i1,0),D,p),lambda) + pow(pow(D,1.0/p)/distance(gsl_matrix_ptr(data,i,0),gsl_matrix_ptr(data,i2,0),D,p),lambda);
		}

	}

	//return the cost difference
	return (2.0/(Npoints*(Npoints-1)))*(costAfter - costBefore);

}

/*This function swaps the particle pair back: needed if the original swap didn't improve the cost function*/

void swapBack(gsl_matrix *data,int i1,int i2,int d){

	double temp;

	//perform the coordinate swap
	temp = gsl_matrix_get(data,i1,d);
	gsl_matrix_set(data,i1,d,gsl_matrix_get(data,i2,d));
	gsl_matrix_set(data,i2,d,temp);

}

/*Main design sampler*/

double sample(int Npoints,int D,double p,double lambda,int seed,int maxIterations,gsl_matrix *data,double *costValues){

	int i,d,i1,i2,iterCount;

	const gsl_rng_type *T;
	gsl_rng *r;
	gsl_permutation *perm = gsl_permutation_alloc(Npoints);

	double currentCost,deltaCost,deltaPerc;

	//Initialize random number generator
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);

	//Initialize permutation
	gsl_permutation_init(perm);

	//Initialize random number generator with provided seed
	gsl_rng_set(r,seed);

	//Initialize the point coordinates in data with random permutations of (1..Npoints) to enforce latin hypercube structure
	for(d=0;d<D;d++){

		//Shuffle the numbers
		gsl_ran_shuffle(r,perm->data,Npoints,sizeof(size_t));

		//Permute coordinates
		for(i=0;i<Npoints;i++){
			gsl_matrix_set(data,i,d,(double)perm->data[i]/(Npoints-1));
		}

	}

	/*The loop does the following: it swaps a random coordinate of a random pair,
	checks if the cost is lower. If so, it keeps the configuration, otherwise it
	reverses it and tries a new one.*/

	iterCount = 0;
	currentCost = cost(data,Npoints,D,p,lambda);
	deltaPerc = 0.0;

	while(1){

		//Decide which coordinate to swap of which pair

		i1 = gsl_rng_uniform_int(r,Npoints);
		while((i2=gsl_rng_uniform_int(r,Npoints))==i1);
		d = gsl_rng_uniform_int(r,D);

		//Compute the change in the cost function
		deltaCost = swap(data,Npoints,D,p,lambda,i1,i2,d);

		/*Check if gain in cost is positive or negative: if positive, revert the swap, if negative keep it;
		anyway, log the result*/
		if(deltaCost>=0){
			swapBack(data,i1,i2,d);
		} else{
			currentCost += deltaCost;
			deltaPerc = deltaCost/currentCost;
		}

		//Save the current cost to array
		costValues[iterCount] = currentCost;

		//Criterion to break the loop
		if(++iterCount == maxIterations) break;

	}

	//Release resources for random number generator and permutations
	gsl_rng_free(r);
	gsl_permutation_free(perm);

	//Return the relative cost change due to the last iteration
	return deltaPerc;

}

max_stars_repo_name: asabyr/LensTools

max_issues_repo_name: asabyr/LensTools

max_forks_repo_name: asabyr/LensTools


text:
#ifndef QCQP_B_H
#define QCQP_B_H

#include "Matrix.h"
#include "ALM_L_Katyusha.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

/*This class solves problem of the form
min 1/2 x^\top Q_0 x+b_0^\top x
s.t. 1/2 x^\top Q_i x+b_i^\top x -1\leq 0,\enspace \forall i=1,\dots,m
-b\leq x\leq b

using IPALM with loopless Katyusha as inner solver.
*/

template<typename L, typename D>
class QCQP_b: public ALM_L_Katyusha<L, D>
{
private:

  vector<D> Lambda; // maximal eigenvalues of Q0, Q1,...Qm.

  vector<D> norm_c; //||[b0;b1;...;bm] ||_2

  vector<D> norm_1_c; //||[b0;b1;...;bm] ||_1

  D b;  // absolute bound on each variable

  D nb_q_c;  //number of quadratic constraints

  D nb_v; //number of variables

  Matrix<L,D> data_Q; // Q=[Q0;Q1;...;Qm]; b=[b0;b1;...;bm].

protected:

public:

  QCQP_b(const char* Matrix_file,D val_b)
  :ALM_L_Katyusha<L,D>(), data_Q(Matrix_file)
  {
    b= val_b;
    nb_v=data_Q.nfeatures;
    nb_q_c=data_Q.nsamples/data_Q.nfeatures- 1;
    this->mu_g=0;
  }

  inline D get_n(){
    return data_Q.nfeatures;
  }

  inline D get_m(){
    return data_Q.nsamples/data_Q.nfeatures- 1;
  }

 void set_dimension(){
    this->d= nb_v;
    this->m= nb_q_c;
  }

  void set_is_feasibility_constraint(){
    for(L i=0;i<nb_q_c;i++)
    this->is_feasiblity_constraint[i]=1;
  }

  inline D value_of_P(vector<D> & x){
    D res=0;
    for(L i=0;i<this->d;i++)
    {
      if (fabs(x[i])> b){
        cout<< "error: the variable x_"<< i<<"="<<x[i]<<". It should be in [-"<<b<<", "<<b<<" ]"<< endl;
        system("pause");
        return std::numeric_limits<double>::max();
      }
    }
    return res;
  }

  inline void prox_of_P(D alpha,vector<D> & x, vector<D> &  y){
    for(L i=0;i<this->d;i++)
    {
      D x2=x[i];
      if (x2>= b){
        y[i]=b;
      }
      else if(x2<=-b){
        y[i]=-b;
      }
      else{
        y[i]=x2;
      }
    }
  }

  D prox_of_h_star_j(D x,D beta, L j){
    if(x>=0) return x;
    else return 0;
  }

  D value_of_h_star_j(D x,L j){
    if(x>=0) return 0;
    else {
      cout<< "Error! x="<<x<<" should be nonnegative."<<endl;
      system("pause");
      return std::numeric_limits<double>::max();
    }
  }

  void compute_gradient_f0(vector<D> & x, vector<D>& g){
    for(L j=0;j<nb_v;j++){
      g[j]=0;
      for(L k = data_Q.ptr[j]; k < data_Q.ptr[j + 1];k++){
        L kj=data_Q.row_idx[k];
        g[j]+=x[kj]*data_Q.A[k];
      }
      g[j]+=data_Q.b[j];
    }
  }

  D compute_gradient_and_value_f_i(vector<D> & x ,vector<D> & g, L i){
    D fi=0;
    for(L j=0;j<nb_v;j++){
      g[j]=0;
      L j2=i*nb_v+j;
      for(L k = data_Q.ptr[j2]; k < data_Q.ptr[j2 + 1];k++){
        L kj=data_Q.row_idx[k];
        g[j]+=x[kj]*data_Q.A[k];
      }
      D bj=data_Q.b[j2];
      fi+=x[j]*g[j]/2+x[j]*bj;
      g[j]+=bj;
    }
    return fi-1;
  }

  D value_of_f0(vector<D> & x){
    return my_value_of_f_i( x, 0);
  }

  D value_of_f_i(vector<D> & x, L i){
    return my_value_of_f_i( x, i)-1;
  }

  D my_value_of_f_i(vector<D> & x, L i){
    D fi=0;
    for(L j=0;j<nb_v;j++){
      D Qjx=0;
      L j2=i*nb_v+j;
      for(L k = data_Q.ptr[j2]; k < data_Q.ptr[j2 + 1];k++){
        L kj=data_Q.row_idx[k];
        Qjx+=x[kj]*data_Q.A[k];
      }
      fi+=x[j]*Qjx/2+x[j]*data_Q.b[j2];
    }
    return fi;
  }

  D value_of_h_j(D x, L i){
    cout<< "Error! The "<<i<<"th function should be a feasibility constraint"<< endl;
    system("pause");
    return std::numeric_limits<double>::max();
  }

  D distance_to_domain_of_h_j(D x, L i){
    if(x<=0)
    return 0;
    else
    return x;
  }

  D compute_Qi_x(vector<D> & x ,vector<D> & g, L i){
    D normg=0;
    for(L j=0;j<nb_v;j++){
      g[j]=0;
      L j2=i*nb_v+j;
      for(L k = data_Q.ptr[j2]; k < data_Q.ptr[j2 + 1];k++){
        L kj=data_Q.row_idx[k];
        g[j]+=x[kj]*data_Q.A[k];
      }
      normg+=g[j]*g[j];
    }
    return sqrt(normg);
  }

  void compute_lambda(){
    D maxv=0;
    D minv=std::numeric_limits<double>::max();
    Lambda.resize(nb_q_c+1);
    norm_c.resize(nb_q_c+1,0);
    norm_1_c.resize(nb_q_c+1,0);
    std::vector<D> xk(nb_v,1);
    std::vector<D> yk(nb_v);
    D normk;
    D tmp;
    for (L i= 0; i<nb_q_c+1 ; i++){
      for(L kk=0;kk<20;kk++){
        normk= compute_Qi_x(xk , yk ,i);
        for (L j=0;j<nb_v;j++)
        xk[j]=yk[j]/normk;
      }
      D res=0;
      normk= compute_Qi_x(xk, yk ,i);
      for (L j=0;j<nb_v;j++)
      res+=xk[j]*yk[j];
      Lambda[i]= res;
      if (res> maxv){
        maxv= res;
      }
      if (res< minv){
        minv= res;
      }
      D res2= 0;
      D res3= 0;
      for (L j=0; j< nb_v; j++){
        res2+= data_Q.b[i*nb_v+j]*data_Q.b[i*nb_v+j];
        res3+= fabs(data_Q.b[i*nb_v+j]);
      }
      norm_c[i]= sqrt(res2);
      norm_1_c[i]=res3;
    }
    cout<< "max Lambda= "<< maxv<< " min Lambda= "<< minv<< endl;
  }

  inline void set_Li_Lf(){
    this->Li.clear();
    this->Li.resize(this->nsamples,0);
    this->Li[0]=this->nsamples*Lambda[0];
    this->Lf=Lambda[0];
    D maxLi=this->Li[0];
    D minLi=this->Li[0];
    for (L i=0; i< nb_q_c; i++){
      D M_p_j= nb_v*Lambda[i+ 1]*b*b+ b*norm_1_c[i+1]+ 1;
      D M_dp_j= sqrt(nb_v)*Lambda[i+ 1]*b+ norm_c[i+1];
      D L_dp_j= Lambda[i+1];
      D d_s_j= M_p_j+ this->beta_s*this->lambda_s[i];
      D tmp=L_dp_j*d_s_j/this->beta_s+ M_dp_j*M_dp_j/this->beta_s;
      this->Li[i+1]= this->nsamples*tmp;
      this->Lf+=tmp;
      if (this->Li[i+1]> maxLi){
        maxLi= this->Li[i+ 1];
      }
      if (this->Li[i+ 1]< minLi){
        minLi= this->Li[i+ 1];
      }
    }
    this->sumLi=this->Lf*this->nsamples;
    cout<<"  max of Li: "<<maxLi<<"  min of Li: "<<minLi<<" sumLi="<<this->sumLi<<" Lf= "<<this->Lf<<endl;
  }

  void ALM_QCQP_solver(D beta_0, D epsilon_0,  D eta, D rho,vector<D> & x0,vector<D> & y0,L val_tau, L max_nb_outer, L p_N_1, L p_N_2,string filename1, string filename2, D time){
    compute_lambda();
    this->ALM_solve_with_L_Katyusha(beta_0, epsilon_0,  eta, rho,x0,y0, val_tau, max_nb_outer,  p_N_1,  p_N_2,  filename1,  filename2, time);
  }

};

#endif


text:
#ifndef SEARCHPATTERNHEXMATCHER_HHHHH
#define SEARCHPATTERNHEXMATCHER_HHHHH

#include <gsl/string_span>
#include <vector>
#include <string>
#include <algorithm>

#include "SearchPatternBase.h"

namespace MPSig {
    namespace detail {
        class SearchPatternHexMatcher final : public SearchPatternBase {
        private:
            std::vector<std::pair<char, bool>> m_compiledPattern;
        public:
            SearchPatternHexMatcher(gsl::cstring_span<-1> hexPattern);

            virtual SearchPatternBase::ExecFirstResult<gsl_span_cit_t> ExecFirst(const MPSig::SigScannerMemoryData& data, gsl_span_cit_t begin, gsl_span_cit_t end) const override;
            virtual SearchPatternBase::ExecFirstResult<gsl_span_crit_t> ExecFirst(const MPSig::SigScannerMemoryData& data, gsl_span_crit_t begin, gsl_span_crit_t end) const override;

            virtual SearchPatternBase::ExecFirstResult<gsl_span_cit_t> ExecDepend(const MPSig::SigScannerMemoryData& data, gsl_span_cit_t begin, gsl_span_cit_t end) const override;
            virtual SearchPatternBase::ExecFirstResult<gsl_span_crit_t> ExecDepend(const MPSig::SigScannerMemoryData& data, gsl_span_crit_t begin, gsl_span_crit_t end) const override;

            virtual std::unique_ptr<SearchPatternBase> Clone() const override;
        };

    }
}

#endif


text:
/* AUTORIGHTS
Copyright (C) 2007 Princeton University

This file is part of Ferret Toolkit.

Ferret Toolkit is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "LSH.h"

extern int LSH_release (cass_table_t *table);
extern int LSH_dump (cass_table_t *table);

int LSH_init_private (cass_table_t *table, const char *param)

{
	LSH_t *lsh;
	cass_size_t D, M, L, H;
	cass_size_t n;
	const gsl_rng_type *T;
	gsl_rng *r;
	int i, j;
	float w;

	lsh = type_calloc(LSH_t, 1);
	assert(table->parent_cfg != NULL);
	D = lsh->D = table->parent_cfg->vec_dim;
	M = lsh->M = param_get_int(param, "-M", 10);
	assert(M <= sizeof(uint32_t) * 8);
	L = lsh->L = param_get_int(param, "-L", 1);
	H = lsh->H = param_get_int(param, "-H", 1017881);

	lsh->W = type_calloc(float, L);

	w = param_get_float(param, "-w", 1.0);

	for (i = 0; i < L; i++) lsh->W[i] = w;

	if (strstr(param, "-W") && strstr(param, "-recall")) debug("-W and -recall conflicts.\n");

	param_get_float_array(param, "-W", &n, lsh->W);

	lsh->alphas = type_matrix_alloc(float, L * M, D);
	lsh->rnd = type_matrix_alloc(uint32_t, L, M);
	lsh->betas = (float *)malloc(L * M * sizeof(float));

	lsh->tmp = type_matrix_alloc(uint32_t, L, M);
	lsh->tmp2 = type_calloc(uint32_t, L);

	gsl_rng_env_setup();

	T = gsl_rng_default;
	r = gsl_rng_alloc(T);

	for (i = 0; i < L * M; i++)
	{
		for (j = 0; j < D; j++)
		{
			lsh->alphas[i][j] = gsl_ran_gaussian(r, 1.0);
		}
	}

	for (i = 0; i < L * M; i++)
	{
		lsh->betas[i] = gsl_rng_uniform(r) * lsh->W[i / M];
	}
	gsl_rng_free(r);

	srand(0);
	for (i = 0; i < L; i++)
		for (j = 0; j < M; j++)
			lsh->rnd[i][j] = rand();

	lsh->est = NULL;

	if (strstr(param, "-recall"))
	{
		cass_size_t d_step, max_T;
		float d_min, d_max;

		d_step = param_get_int(param, "-d_step", 100);
		max_T = param_get_int(param, "-T_max", 100);
		d_min = param_get_float(param, "-d_min", 0.0);
		d_max = param_get_float(param, "-d_max", 100);
		LSH_recall_init(&lsh->recall, d_step, d_min, d_max, lsh->M, lsh->L, max_T, lsh->W[0]);
	}

	table->__private = lsh;
	return 0;
}

int LSH_free_private (cass_table_t *table)
{
	LSH_t *lsh = table->__private;
	if (table->loaded) LSH_release(table);
	matrix_free(lsh->rnd);
	free(lsh->betas);
	matrix_free(lsh->alphas);
	free(lsh->W);

	matrix_free(lsh->tmp);
	free(lsh->tmp2);

	if (lsh->recall.table != NULL) LSH_recall_cleanup(&lsh->recall);

	free(lsh);

	table->__private = NULL;
	return 0;
}

void LSH_hash2 (LSH_t *lsh, uint32_t **hash, uint32_t *hash2)
{
	int i, j;
	for (i = 0; i < lsh->L; i++)
	{
		hash2[i] = 0;
		for (j = 0; j < lsh->M; j++)
		{
			hash2[i] += lsh->rnd[i][j] * hash[i][j];
		}
		hash2[i] %= lsh->H;
	}
}

void LSH_hash (LSH_t *lsh, const float *pnt, uint32_t **hash)
{
	float s;
	int i, j, k, l;

	l = 0;
	for (i = 0; i < lsh->L; i++)
	{
		for (j = 0; j < lsh->M; j++)
		{
			s = lsh->betas[l];
			for (k = 0; k < lsh->D; k++)
			{
				s += pnt[k] * lsh->alphas[l][k];
			}

			s /= lsh->W[i];
			hash[i][j] = floor(s);
			l++;
		}
	}
}

static inline void LSH_hash_L (LSH_t *lsh, float *pnt, uint32_t **hash, int L)
{
	float s;
	int j, k, l;

	l = L * lsh->M;
	for (j = 0; j < lsh->M; j++)
	{
		s = lsh->betas[l];
		for (k = 0; k < lsh->D; k++)
		{
			s += pnt[k] * lsh->alphas[l][k];
		}

		s /= lsh->W[L];
		hash[L][j] = floor(s);
		l++;
	}
}

static inline void LSH_hash2_L (LSH_t *lsh, uint32_t **hash, uint32_t *hash2, int L)
{
	int j;
	hash2[L] = 0;
	for (j = 0; j < lsh->M; j++)
	{
		hash2[L] += lsh->rnd[L][j] * hash[L][j];
	}
	hash2[L] %= lsh->H;
}

static void LSH_insert (LSH_t *lsh, float *pnt, int tag)
{
	int l;
	LSH_hash(lsh, pnt, lsh->tmp);
	LSH_hash2(lsh, lsh->tmp, lsh->tmp2);

	for (l = 0; l < lsh->L; l++)
	{
		ohash_insert(&lsh->hash[l], lsh->tmp2[l], tag);
	}
	lsh->count++;
}

int LSH_batch_insert(cass_table_t *table, cass_dataset_t *parent, cass_vecset_id_t start, cass_vecset_id_t end)
{
	LSH_t *lsh = table->__private;
	uint32_t i, j, k;
	cass_vec_t *vec;
	for (i = start; i <= end; i++)
	{
		for (j = 0;  j < parent->vecset[i].num_regions; j++)
		{
			k =  j + parent->vecset[i].start_vecid;
			vec = DATASET_VEC(parent, k);
			LSH_insert(lsh, vec->u.float_data, k);
			lsh->count++;
		}
	}
	return 0;
}

/*
void LSH_mass_insert_parallel (LSH_t *lsh, float **pnt, int cnt)
{
	int l;
#pragma omp parallel for schedule(static, 1) shared(lsh, pnt, cnt)
	for (l =  0; l < lsh->L; l++)
	{
		printf("#TH = %d / %d\n", omp_get_thread_num(),
				omp_get_num_threads());
		int i;
		for (i = 0; i < cnt; i++)
		{
			LSH_hash_L(lsh, pnt[i], tmp, l);
			LSH_hash2_L(lsh, tmp, tmp2, l);
			ohash_insert(&lsh->hash[l], tmp2[l], i);
		}
	}
}
*/

int LSH_restore_private (cass_table_t *table, CASS_FILE *fin)
{
	int ret;
	cass_size_t L, M, D;
	int32_t recall;
	LSH_t *lsh;

	lsh = type_calloc(LSH_t, 1);

	ret = cass_read_size(&lsh->D, 1, fin);
	ret += cass_read_size(&lsh->M, 1, fin);
	ret += cass_read_size(&lsh->L, 1, fin);
	ret += cass_read_size(&lsh->H, 1, fin);
	ret += cass_read_size(&lsh->count, 1, fin);

	assert(ret == 5);

	lsh->W = type_calloc(float, lsh->L);
	ret  = cass_read_float(lsh->W, lsh->L, fin);
	assert(ret == lsh->L);

	L = lsh->L;
	M = lsh->M;
	D = lsh->D;

	lsh->alphas = type_matrix_alloc(float, L * M, D);
	lsh->betas = type_calloc(float, L * M);
	lsh->hash = NULL; //(ohash_t *)malloc(L * sizeof(ohash_t));

	lsh->tmp = type_matrix_alloc(uint32_t, L, M);
	lsh->tmp2 = malloc(sizeof(uint32_t) *L );

	ret = cass_read_float(lsh->alphas[0], L * M * D, fin);
	assert(ret == L * M * D);
	ret = cass_read_float(lsh->betas, L * M, fin);
	assert(ret == L * M);

	lsh->rnd = type_matrix_alloc(uint32_t, L, M);

	ret = cass_read_uint32(lsh->rnd[0], L * M, fin);
	assert(ret == L * M);

	ret = cass_read_int32(&recall, 1, fin);
	assert(ret == 1);

	if (recall)
	{
		if (LSH_recall_load(&lsh->recall, fin) != 0) assert(0);
	}

	table->__private = lsh;

	return 0;
}

int LSH_checkpoint_private (cass_table_t *table, CASS_FILE *fout)
{
	LSH_t *lsh = table->__private;
	int ret;
	cass_size_t L, M, D;
	int32_t recall;
	ret = cass_write_size(&lsh->D, 1, fout);
	ret += cass_write_size(&lsh->M, 1, fout);
	ret += cass_write_size(&lsh->L, 1, fout);
	ret += cass_write_size(&lsh->H, 1, fout);
	ret += cass_write_size(&lsh->count, 1, fout);
	assert(ret == 5);
	L = lsh->L;
	M = lsh->M;
	D = lsh->D;

	ret = cass_write_float(lsh->W, L, fout);
	assert(ret == L);

	ret = cass_write_float(lsh->alphas[0], L * M * D, fout);
	assert(ret == L * M * D);
	ret = cass_write_float(lsh->betas, L * M, fout);
	assert(ret == L * M);
	ret = cass_write_uint32(lsh->rnd[0], L * M, fout);
	assert(ret == L * M);

	recall = (lsh->recall.table != NULL);

	ret = cass_write_int32(&recall, 1, fout);
	assert(ret == 1);

	if (recall)
	{
		if (LSH_recall_dump(&lsh->recall, fout) != 0) assert(0);
	}

	if (lsh->hash != NULL && table->dirty) return LSH_dump(table);

	return 0;
}

int LSH_load (cass_table_t *table)
{
	LSH_t *lsh = table->__private;
	CASS_FILE *fin;
	uint32_t i;
	int ret;

	assert(lsh->hash == NULL);

	if (lsh->count == 0)
	{
		lsh->hash = type_calloc(ohash_t, lsh->L);
		for (i = 0; i < lsh->L; i++)
		{
			ohash_init(&(lsh->hash[i]), lsh->H);
		}
		return 0;
	}

	fin = cass_open(table->filename, "r");
	if (fin == NULL) return CASS_ERR_IO;

	lsh->hash = type_calloc(ohash_t, lsh->L);
	for (i = 0; i < lsh->L; i++)
	{
		ret = ohash_init_with_stream(&lsh->hash[i], fin);
		assert(ret == 0);
	}

	cass_close(fin);
	return 0;
}

int LSH_release (cass_table_t *table)
{
	int err = 0;
	uint32_t i;
	LSH_t *lsh = table->__private;
	if (table->loaded && table->dirty) err = LSH_dump(table);
	if (err != 0) return err;
	for (i = 0; i < lsh->L; i++)
	{
		ohash_cleanup(&(lsh->hash[i]));
	}
	free(lsh->hash);
	lsh->hash = NULL;
	return 0;
}

int LSH_dump (cass_table_t *table)
{
	uint32_t i;
	int ret;
	CASS_FILE *fout;
	LSH_t *lsh = table->__private;
	if (lsh->count == 0) return 0;
	assert(lsh->hash != NULL);
	fout = cass_open(table->filename, "w");
	assert(fout != NULL);
	for (i = 0; i < lsh->L; i++)
	{
		ret = ohash_dump_stream(&lsh->hash[i], fout);
		assert(ret == 0);
	}
	cass_close(fout);
	table->dirty = 0;
	return ret;
}

void LSH_stat (LSH_t *lsh)
{
	int i;
	for (i = 0; i < lsh->L; i++)
	{
		printf("%d:", i);
		ohash_stat(&lsh->hash[i]);
	}
}

/*
int LSH_est_dump_file (LSH_est_t *est, char *filename)
{
	int ret;
	CASS_FILE *fout = cass_open(filename, "wb");
	assert(fout != NULL);
	ret = 0;
	ret += cass_write_int(&est->a_step, 1, fout);
	ret += cass_write_double(&est->a_min, 1, fout);
	ret += cass_write_double(&est->a_max, 1, fout);
	ret += cass_write_int(&est->b_step, 1, fout);
	ret += cass_write_double(&est->b_min, 1, fout);
	ret += cass_write_double(&est->b_max, 1, fout);

	assert(ret == 6);

	est->a_delta = (est->a_max - est->a_min) / est->a_step;
	est->b_delta = (est->b_max - est->b_min) / est->b_step;

	type_matrix_dump_stream(double, fout, est->a_step, est->b_step, est->a_table);
	type_matrix_dump_stream(double, fout, est->a_step, est->b_step, est->b_table);

	cass_close(fout);

	return 0;
}

int LSH_est_init_with_file (LSH_est_t *est, char *filename)
{
	int ret;
	cass_size_t row, col;
	CASS_FILE *fin = cass_open(filename, "rb");
	assert(fin != NULL);
	ret = 0;
	ret += cass_read_int32(&est->a_step, 1, fin);
	ret += cass_read_double(&est->a_min, 1, fin);
	ret += cass_read_double(&est->a_max, 1, fin);
	ret += cass_read_int32(&est->b_step, 1, fin);
	ret += cass_read_double(&est->b_min, 1, fin);
	ret += cass_read_double(&est->b_max, 1, fin);

	assert(ret == 6);

	est->a_delta = (est->a_max - est->a_min) / est->a_step;
	est->b_delta = (est->b_max - est->b_min) / est->b_step;

	type_matrix_load_stream(double, fin, &row, &col, &est->a_table);
	assert(row == est->a_step);
	assert(col == est->b_step);
	type_matrix_load_stream(double, fin, &row, &col, &est->b_table);
	assert(row == est->a_step);
	assert(col == est->b_step);

	cass_close(fin);

	return 0;
}
*/

void LSH_est_cleanup (LSH_est_t *est)
{
	matrix_free(est->a_table);
	matrix_free(est->b_table);
}

double LSH_est (LSH_est_t *est, double a, double b, int K)
{
	int A, B;
	double alpha, beta;
	if (!(finite(a) && finite(b))) return 0.0;

	A = (a - est->a_min) / est->a_delta;
	B = (b - est->b_min) / est->b_delta;

	if (A < 0) return 0.0;
	if (B < 0) return 0.0;
	if (A >= est->a_step) A = est->a_step - 1;
	if (B >= est->a_step) B = est->a_step - 1;

	alpha = est->a_table[A][B];
	beta = est->b_table[A][B];

	return alpha * pow(K, beta);
}

int __LSH_query(cass_table_t *table,
                     cass_query_t *query, cass_result_t *result)
{
	LSH_t *lsh = (LSH_t *)table->__private;
	LSH_query_t query2;
	cass_table_t *parent;
	cass_dataset_t *ds;
	cass_vec_dist_t *vec_dist;
	cass_vecset_t *vecset;

	float recall;

	cass_size_t K, L, T;

	uint32_t i, j;

	assert(table->loaded);
	vec_dist = cass_reg_get(&table->env->vec_dist, query->vec_dist_id);
	assert(vec_dist != NULL);
	if (vec_dist->__class != &vec_dist_L2_float) debug("LSH only works for L2 distance for float.\n");

	if (table->parent_id == CASS_ID_INV) debug("LSH query requires parent.\n");
	parent = cass_reg_get(&table->env->table, table->parent_id);
	assert(parent != NULL);
	if (!parent->loaded) debug("LSH query requires parent to be loaded.\n");
	if (!(parent->opr->type & CASS_DATA)) debug("LSH only works on CASS_DATA.\n");

	ds = (cass_dataset_t *)parent->__private;

	K = query->topk;
	if (K == 0) debug("LSH only works for KNN queries.\n");

	L = param_get_int(query->extra_params, "-L", 1);
	T = param_get_int(query->extra_params, "-T", 0);

	recall = param_get_float(query->extra_params, "-recall", 0);

	LSH_query_init(&query2, lsh, ds, K, L, T);

	assert((query->flags & CASS_RESULT_BITMAPS) == 0);
	assert((query->flags & CASS_RESULT_BITMAP) == 0);
	assert((query->flags & CASS_RESULT_LIST) == 0);
	assert(query->flags & CASS_RESULT_LISTS);
	assert(query->candidate == NULL);

	result->flags = CASS_RESULT_LISTS;

	vecset = query->dataset->vecset + query->vecset_id;

	if (query->flags & CASS_RESULT_USERMEM)
	{
		assert(result->u.lists.len >= vecset->num_regions);
	}
	else
	{
		result->flags |= CASS_RESULT_MALLOC;
		ARRAY_INIT_SIZE(result->u.lists, vecset->num_regions);
		result->u.lists.len = vecset->num_regions;
	}

	for (i = 0; i < vecset->num_regions; i++)
	{
		if (query->flags & CASS_RESULT_USERMEM)
		{
			assert(result->u.lists.data[i].size >= query->topk);
		}
		else
		{
			ARRAY_INIT_SIZE(result->u.lists.data[i], query->topk);
		}
		result->u.lists.data[i].len = query->topk;

		if (recall == 0)
		{
			LSH_query(&query2, DATASET_VEC(query->dataset, vecset->start_vecid + i)->u.float_data);
		}
		else
		{
			LSH_query_recall(&query2, DATASET_VEC(query->dataset, vecset->start_vecid + i)->u.float_data, recall);
		}

		for (j = 0; j < K; j++)
		{
			result->u.lists.data[i].data[j] = query2.topk[j];
		}
	}

	if (query->flags & CASS_RESULT_SORT)
	{
		for (j = 0; j < vecset->num_regions; j++)
		{
			TOPK_SORT_MIN(result->u.lists.data[j].data, cass_list_entry_t, dist, query->topk);
		}
		result->flags |= CASS_RESULT_SORT;
	}

	LSH_query_cleanup(&query2);

	return 0;

}

int __LSH_batch_query(cass_table_t *table,
/* the dataset of all the queries should be the same */
	uint32_t count, cass_query_t **queries, cass_result_t **results)
{
	LSH_t *lsh = (LSH_t *)table->__private;
	cass_result_t *result;
	cass_query_t *query;

	LSH_query_t query2;

	const float **points;
	cass_list_entry_t **topks;

	cass_table_t *parent;
	cass_dataset_t *ds;
	cass_vec_dist_t *vec_dist;
	cass_vecset_t *vecset;

	cass_size_t K, L, T, vec_count;

	uint32_t i, j, k;

	query = queries[0];

	assert(table->loaded);
	vec_dist = cass_reg_get(&table->env->vec_dist, query->vec_dist_id);
	assert(vec_dist != NULL);
	if (vec_dist->__class != &vec_dist_L2_float) debug("LSH only works for L2 distance for float.\n");

	if (table->parent_id == CASS_ID_INV) debug("LSH query requires parent.\n");
	parent = cass_reg_get(&table->env->table, table->parent_id);
	assert(parent != NULL);
	if (!parent->loaded) debug("LSH query requires parent to be loaded.\n");
	if (!(parent->opr->type & CASS_DATA)) debug("LSH only works on CASS_DATA.\n");

	ds = (cass_dataset_t *)parent->__private;

	K = query->topk;
	if (K == 0) debug("LSH only works for KNN queries.\n");

	L = param_get_int(query->extra_params, "-L", 1);
	T = param_get_int(query->extra_params, "-T", 0);

	LSH_query_init(&query2, lsh, ds, K, L, T);

	assert((query->flags & CASS_RESULT_BITMAPS) == 0);
	assert((query->flags & CASS_RESULT_BITMAP) == 0);
	assert((query->flags & CASS_RESULT_LIST) == 0);
	assert(query->flags & CASS_RESULT_LISTS);
	assert(query->candidate == NULL);

	vec_count = 0;
	for (i = 0; i < count; i++)
	{
		result = results[i];
		query = queries[i];
		result->flags = CASS_RESULT_LISTS;
		vecset = query->dataset->vecset + query->vecset_id;

		if (query->flags & CASS_RESULT_USERMEM)
		{
			assert(result->u.lists.len >= vecset->num_regions);
		}
		else
		{
			result->flags |= CASS_RESULT_MALLOC;
			ARRAY_INIT_SIZE(result->u.lists, vecset->num_regions);
			result->u.lists.len = vecset->num_regions;
		}

		for (j = 0; j < vecset->num_regions; j++)
		{
			if (query->flags & CASS_RESULT_USERMEM)
			{
				assert(result->u.lists.data[j].size >= query->topk);
			}
			else
			{
				ARRAY_INIT_SIZE(result->u.lists.data[j], query->topk);
			}
			result->u.lists.data[j].len = K;
		}
		vec_count += vecset->num_regions;
	}

	points = type_calloc(const float *, vec_count);
	topks = type_calloc(cass_list_entry_t *, vec_count);

	k = 0;
	for (i = 0; i < count; i++)
	{
		result = results[i];
		query = queries[i];
		vecset = query->dataset->vecset + query->vecset_id;
		for (j = 0; j < vecset->num_regions; j++)
		{
			points[k] = DATASET_VEC(query->dataset, vecset->start_vecid + j)->u.float_data;
			topks[k] = result->u.lists.data[j].data;
			k++;
		}
	}

	assert(k == vec_count);

	if (strstr(query->extra_params, "-ca"))
	{
		LSH_query_batch_ca(&query2, vec_count, points, topks);
	}
	else
	{
		LSH_query_batch(&query2, vec_count, points, topks);
	}

	query = queries[0];

	if (query->flags & CASS_RESULT_SORT)
	{
		for (i = 0; i < vec_count; i++)
		{
			TOPK_SORT_MIN(topks[i], cass_list_entry_t, dist, K);
		}
		for (i = 0; i < count; i++)
		{
			results[i]->flags |= CASS_RESULT_SORT;
		}
	}

	LSH_query_cleanup(&query2);

	return 0;
}

cass_table_opr_t opr_lsh = {
	.name = "lsh",
	.type = CASS_VEC_INDEX,
	.vecset_type = CASS_ANY,
	.vec_type = CASS_VEC_FLOAT,
	.dist_vecset = CASS_ANY,
	.dist_vec = CASS_VEC_DIST_TYPE_L2,

	.cfg = NULL,
	.tune = NULL, //mplsh_tune,
	.init_private = LSH_init_private,
	.batch_insert = LSH_batch_insert,
	.query = __LSH_query, //mplsh_query,
	.batch_query = __LSH_batch_query,
	.load = LSH_load,
	.release = LSH_release,
	.checkpoint_private = LSH_checkpoint_private,
	.restore_private = LSH_restore_private,
	.free_private = LSH_free_private
};

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* rng/cmrg.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * 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_rng.h>

/* This is a combined multiple recursive generator. The sequence is,

   z_n = (x_n - y_n) mod m1

   where the two underlying generators x and y are,

   x_n = (a_{1} x_{n-1} + a_{2} x_{n-2} + a_{3} x_{n-3}) mod m1
   y_n = (b_{1} y_{n-1} + b_{2} y_{n-2} + b_{3} y_{n-3}) mod m2

   with coefficients a11 ... a23,

   a_{1} = 0,     a_{2} = 63308, a_{3} = -183326
   b_{1} = 86098, b_{2} = 0,     b_{3} = -539608

   and moduli m1, m2,

   m1 = 2^31 - 1 = 2147483647
   m2 = 2^31 - 2000169 = 2145483479

   We initialize the generator with

   x_1 = s_1 MOD m1, x_2 = s_2 MOD m1, x_3 = s_3 MOD m1
   y_1 = s_4 MOD m2, y_2 = s_5 MOD m2, y_3 = s_6 MOD m2

   where s_n = (69069 * s_{n-1}) mod 2^32 and s_0 = s is the
   user-supplied seed.

   NOTE: According to the paper the initial values for x_n must lie in
   the range 0 <= x_n <= (m1 - 1) and the initial values for y_n must
   lie in the range 0 <= y_n <= (m2 - 1), with at least one non-zero
   value -- our seeding procedure satisfies these constraints.

   We then use 7 iterations of the generator to "warm up" the internal
   state.

   The theoretical value of z_{10008} is 719452880. The subscript 10008
   means (1) seed the generator with s=1, (2) do the seven warm-up
   iterations that are part of the seeding process, (3) then do 10000
   actual iterations.

   The period of this generator is about 2^205.

   From: P. L'Ecuyer, "Combined Multiple Recursive Random Number
   Generators," Operations Research, 44, 5 (1996), 816--822.

   This is available on the net from L'Ecuyer's home page,

   http://www.iro.umontreal.ca/~lecuyer/myftp/papers/combmrg.ps
   ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/combmrg.ps */

static inline unsigned long int cmrg_get (void *vstate);
static double cmrg_get_double (void *vstate);
static void cmrg_set (void *state, unsigned long int s);

static const long int m1 = 2147483647, m2 = 2145483479;

static const long int a2 = 63308, qa2 = 33921, ra2 = 12979;
static const long int a3 = -183326, qa3 = 11714, ra3 = 2883;
static const long int b1 = 86098, qb1 = 24919, rb1 = 7417;
static const long int b3 = -539608, qb3 = 3976, rb3 = 2071;

typedef struct
  {
    long int x1, x2, x3;	/* first component */
    long int y1, y2, y3;	/* second component */
  }
cmrg_state_t;

static inline unsigned long int
cmrg_get (void *vstate)
{
  cmrg_state_t *state = (cmrg_state_t *) vstate;

  /* Component 1 */

  {
    long int h3 = state->x3 / qa3;
    long int p3 = -a3 * (state->x3 - h3 * qa3) - h3 * ra3;

    long int h2 = state->x2 / qa2;
    long int p2 = a2 * (state->x2 - h2 * qa2) - h2 * ra2;

    if (p3 < 0)
      p3 += m1;
    if (p2 < 0)
      p2 += m1;

    state->x3 = state->x2;
    state->x2 = state->x1;
    state->x1 = p2 - p3;
    if (state->x1 < 0)
      state->x1 += m1;
  }

  /* Component 2 */

  {
    long int h3 = state->y3 / qb3;
    long int p3 = -b3 * (state->y3 - h3 * qb3) - h3 * rb3;

    long int h1 = state->y1 / qb1;
    long int p1 = b1 * (state->y1 - h1 * qb1) - h1 * rb1;

    if (p3 < 0)
      p3 += m2;
    if (p1 < 0)
      p1 += m2;

    state->y3 = state->y2;
    state->y2 = state->y1;
    state->y1 = p1 - p3;
    if (state->y1 < 0)
      state->y1 += m2;
  }

  if (state->x1 < state->y1)
    return (state->x1 - state->y1 + m1);
  else
    return (state->x1 - state->y1);
}

static double
cmrg_get_double (void *vstate)
{
  return cmrg_get (vstate) / 2147483647.0 ;
}

static void
cmrg_set (void *vstate, unsigned long int s)
{
  /* An entirely adhoc way of seeding! This does **not** come from
     L'Ecuyer et al */

  cmrg_state_t *state = (cmrg_state_t *) vstate;

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

#define LCG(n) ((69069 * n) & 0xffffffffUL)
  s = LCG (s);
  state->x1 = s % m1;
  s = LCG (s);
  state->x2 = s % m1;
  s = LCG (s);
  state->x3 = s % m1;

  s = LCG (s);
  state->y1 = s % m2;
  s = LCG (s);
  state->y2 = s % m2;
  s = LCG (s);
  state->y3 = s % m2;

  /* "warm it up" */
  cmrg_get (state);
  cmrg_get (state);
  cmrg_get (state);
  cmrg_get (state);
  cmrg_get (state);
  cmrg_get (state);
  cmrg_get (state);
}

static const gsl_rng_type cmrg_type =
{"cmrg",			/* name */
 2147483646,			/* RAND_MAX */
 0,			        /* RAND_MIN */
 sizeof (cmrg_state_t),
 &cmrg_set,
 &cmrg_get,
 &cmrg_get_double};

const gsl_rng_type *gsl_rng_cmrg = &cmrg_type;


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

#undef __FUNCT__
#define __FUNCT__ "MatVecMPFR"
PetscErrorCode MatVecMPFR(PetscInt size_d, mpfr_t *Mat, mpfr_t *vec, mpfr_t *result) {

  PetscErrorCode ierr;
  mpfr_t temp1;
  mpfr_t *res;
  PetscFunctionBegin;

  mpfr_init(temp1);
  res = (mpfr_t*) malloc(sizeof(mpfr_t)*size_d);

  for(PetscInt k=0; k < size_d; ++k)
    mpfr_init(res[k]);

  for(PetscInt i=0; i<size_d; ++i) {
    mpfr_set_d(res[i], 0.0, MPFR_RNDN);
    for(PetscInt j=0; j<size_d; ++j) {
      mpfr_mul(temp1, Mat[i*size_d+j], vec[j], MPFR_RNDN);
      mpfr_add(res[i], res[i], temp1, MPFR_RNDN);
    }
  }
  for(PetscInt k=0; k < size_d; ++k) {
    mpfr_set(result[k], res[k], MPFR_RNDN);
    mpfr_clear(res[k]);
  }

  mpfr_clear(temp1);
  free(res);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "VecViewMPFR"
PetscErrorCode VecViewMPFR(PetscInt size_d, mpfr_t *vec) {
  PetscErrorCode ierr;
  double val;
  PetscFunctionBegin;
  printf("-------------------\n");
  for(PetscInt k=0; k < size_d; ++k) {
    val = mpfr_get_d(vec[k], MPFR_RNDN);
    printf("vec[%d] = %4.4f\n", k, val);
  }
  printf("-------------------\n");
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatViewMPFR"
PetscErrorCode MatViewMPFR(PetscInt m, PetscInt n, mpfr_t *mat) {
  PetscErrorCode ierr;
  double val;
  PetscFunctionBegin;
  printf("-------------------\n");
  for(PetscInt i=0; i < m; ++i) {
    for(PetscInt j=0; j < n; ++j) {
      val = mpfr_get_d(mat[i*n+j], MPFR_RNDN);
      printf("%3.3f ", val);
    }
    printf("\n");
  }
  printf("-------------------\n");

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "Vec2NormMPFR"
PetscErrorCode Vec2NormMPFR(PetscInt size_d, mpfr_t *vec, mpfr_t *norm)
{
  PetscErrorCode ierr;
  mpfr_t val;
  PetscFunctionBegin;
  mpfr_init(val);
  mpfr_set_d(*norm, 0.0, MPFR_RNDN);
  for(PetscInt k=0; k < size_d; ++k) {
    mpfr_mul(val, vec[k], vec[k], MPFR_RNDN);
    mpfr_add(*norm, *norm, val, MPFR_RNDN);
  }
  mpfr_sqrt(*norm, *norm, MPFR_RNDN);

  mpfr_clear(val);
  PetscFunctionReturn(0);
}

//returns (lambda/|v1|^2)*v1*v1^T*vec
#undef __FUNCT__
#define __FUNCT__ "OrthogMap"
PetscErrorCode OrthogMap(PetscInt size_d, mpfr_t *x, mpfr_t *v1, mpfr_t *lambda,  mpfr_t *result) {

  PetscErrorCode ierr;
  mpfr_t temp1, temp2, coef;
  PetscFunctionBegin;

  mpfr_inits(temp1, temp2, coef, NULL);

  ierr = Vec2NormMPFR(size_d, v1, &temp1);CHKERRQ(ierr);
  mpfr_mul(temp1, temp1, temp1, MPFR_RNDN);
  mpfr_div(coef, *lambda, temp1, MPFR_RNDN);

  //compute v1^T*vec
  mpfr_set_d(temp2, 0.0, MPFR_RNDN);
  for(PetscInt k=0; k < size_d; ++k) {
    mpfr_mul(temp1, v1[k], x[k], MPFR_RNDN);
    mpfr_add(temp2, temp2, temp1, MPFR_RNDN);
  }
  mpfr_mul(coef, coef, temp2, MPFR_RNDN);

  //multiply v1 by coef and store to result
  for(PetscInt k=0; k < size_d; ++k)
    mpfr_mul(result[k], v1[k], coef, MPFR_RNDN);

  mpfr_clears(temp1, temp2, coef, NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "eigsMPFR"
PetscErrorCode eigsMPFR(PetscInt size_d, mpfr_t *mat, PetscInt prec, mpfr_t *result) {

  PetscErrorCode ierr;
  mpfr_t temp1;
  mpfr_t norm;
  mpfr_t error, error2;
  mpfr_t *tempvec, *tempvec2;
  mpfr_t *eigvec;
  mpfr_t *prev;
  mpfr_t *eigvecs;
  mpfr_t *eigvals;
  PetscFunctionBegin;

  mpfr_inits(norm, error, error2, NULL);

  eigvec = (mpfr_t*) malloc(sizeof(mpfr_t)*size_d);
  prev = (mpfr_t*) malloc(sizeof(mpfr_t)*size_d);
  eigvals = (mpfr_t*) malloc(sizeof(mpfr_t)*size_d);
  tempvec = (mpfr_t*) malloc(sizeof(mpfr_t)*size_d);
  tempvec2 = (mpfr_t*) malloc(sizeof(mpfr_t)*size_d);

  for(PetscInt k=0; k < size_d; ++k) {
    mpfr_inits(tempvec2[k], tempvec[k], eigvals[k], eigvec[k], prev[k], NULL);
  }
  for(PetscInt i=0; i < size_d; ++i) {
    for(PetscInt l=0; l < size_d; ++l) {
      mpfr_set_d(eigvec[l], 0.5, MPFR_RNDN);
    }
    mpfr_set_d(error, 1.0, MPFR_RNDN);
    mpfr_set_d(error2, 1.0, MPFR_RNDN);
    mpfr_log10(error, error, MPFR_RNDN);
    mpfr_log10(error2, error2, MPFR_RNDN);
    PetscInt j=0;
    while(mpfr_get_d(error, MPFR_RNDN) > -prec && mpfr_get_d(error2, MPFR_RNDN) > -prec) {
      j++;
      for(PetscInt l=0; l < size_d; ++l)
	mpfr_set(prev[l], eigvec[l], MPFR_RNDN);
      ierr = MatVecMPFR(size_d, mat, eigvec, tempvec2);CHKERRQ(ierr);
      //subtract off already found eigenvectors
      for(PetscInt k=0; k<i; ++k) {
	//PetscErrorCode OrthogMap(PetscInt size_d, mpfr_t *x, mpfr_t *v1, mpfr_t *lambda,  mpfr_t *result) {
	OrthogMap(size_d, eigvec, result+k*size_d, eigvals+k, tempvec);
	for(PetscInt l=0; l<size_d; ++l) {
	  mpfr_sub(tempvec2[l], tempvec2[l], tempvec[l], MPFR_RNDN);
	}
      }
      //normalize eig vector
      ierr = Vec2NormMPFR(size_d, tempvec2, &norm);CHKERRQ(ierr);
      for(PetscInt k=0; k < size_d; ++k)
	mpfr_div(eigvec[k], tempvec2[k], norm, MPFR_RNDN);
      //compute error
      for(PetscInt k=0; k < size_d; ++k) {
	mpfr_sub(tempvec[k], eigvec[k], prev[k], MPFR_RNDN);
	mpfr_add(tempvec2[k], eigvec[k], prev[k], MPFR_RNDN);
      }
      ierr = Vec2NormMPFR(size_d, tempvec, &error);CHKERRQ(ierr);
      ierr = Vec2NormMPFR(size_d, tempvec2, &error2);CHKERRQ(ierr);
      ierr = mpfr_log10(error, error, MPFR_RNDN);
      ierr = mpfr_log10(error2, error2, MPFR_RNDN);
      //if(i==1) {
      //printf("e: %2.2f\n", mpfr_get_d(error, MPFR_RNDN));
      //printf("norm: %3.3e\n", mpfr_get_d(norm, MPFR_RNDN));
      //ierr = VecViewMPFR(size_d, eigvec);CHKERRQ(ierr);
      //}

    }
    //if(i==0) {
    //ierr = VecViewMPFR(size_d, eigvec);CHKERRQ(ierr);
    //}
    //add newly found eigvec to eigvecss
    for(PetscInt k=0; k < size_d; ++k)
      mpfr_set(result[i*size_d+k], eigvec[k], MPFR_RNDN);
    //add new eigval
    mpfr_set(eigvals[i], norm, MPFR_RNDN);

    //printf("EIGVAL IS %3.3f\n", mpfr_get_d(norm, MPFR_RNDN));

  }
  mpfr_clear(norm);
  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <assert.h>
#include <math.h>
#include <sys/time.h>
#include <mm_malloc.h>
#include <cblas.h>
#include <lapacke.h>

#include "random_matrix.h"
#include "driver.h"
#include "../../runtime/pcp.h"

static double compare_lower_triangular_matrices(int n, double *L1, int ldL1, double *L2, int ldL2)
{
#define L1(i,j) L1[(i) + (j) * ldL1]
#define L2(i,j) L2[(i) + (j) * ldL2]

    // Clear out the upper triangular part.
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < j; i++) {
            L1(i,j) = 0.0;
            L2(i,j) = 0.0;
        }
    }

    // Compute L1 := L1 - L2.
    for (int j = 0; j < n; j++) {
        for (int i = 0; i < n; i++) {
            L1(i,j) = L1(i,j) - L2(i,j);
        }
    }

    // Compute || L1 ||_F.
    double error = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', n, n, L1, ldL1);
    error /= n;
    return error;

#undef L1
#undef L2
}

static void usage(char *prog)
{
    printf("Usage: %s n b p q\n", prog);
    printf("\n");
    printf("n: The size of the matrix\n");
    printf("b: The tile size\n");
    printf("p: The number of threads (or 0 for one per core)\n");
    printf("q: The size of the reserved set (or 0 to use the regular mode)\n");
    exit(EXIT_FAILURE);
}

int main(int argc, char** argv)
{
    // Verify the number of command line options.
    if (argc != 5 && argc != 6) {
        usage(argv[0]);
    }

    // Parse command line.
    int n = atoi(argv[1]);
    int blksz = atoi(argv[2]);
    int num_threads = atoi(argv[3]);
    int reserved_set_size = atoi(argv[4]);
    int verify = 1;
    if (argc == 6) {
        verify = 0;
    }

    // Verify options.
    if (n < 1 ||
        blksz < 1 ||
        num_threads < 0 ||
        reserved_set_size < 0 ||
        (reserved_set_size >= num_threads && num_threads > 0)) {
        usage(argv[0]);
    }

    // Initialise random number generator.
    srand(time(NULL));

    // Generate random SPD matrix A.
    printf("Generating SPD matrix A...\n");
    int ldA = n;
    double *A = (double *) malloc(n * ldA * sizeof(double));
    generate_dense_spd_matrix(n, A, ldA);

    // Copy A to Ain.
    printf("Copying A to Ain...\n");
    double *Ain = (double *) malloc(n * ldA * sizeof(double));
    memcpy(Ain, A, n * ldA * sizeof(double));

    // Start the runtime system.
    pcp_start(num_threads);

    // Set the execution mode.
    if (reserved_set_size == 0) {
        pcp_set_mode(PCP_REGULAR);
    } else {
        pcp_set_mode(PCP_FIXED);
        pcp_set_reserved_set_size(reserved_set_size);
    }

    // Call the driver.
    printf("Calling the driver routine...\n");
    parallel_block_chol(n, blksz, A, ldA);

    if (verify) {
        // Verify the solution.
        printf("Verifying the solution... ");
        fflush(stdout);
        double tm = pcp_get_time();
        LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', n, Ain, ldA);
        tm = pcp_get_time() - tm;
        double error = compare_lower_triangular_matrices(n, Ain, ldA, A, ldA);
        if (error < 1e-14) {
            printf("PASSED\n");
        } else {
            printf("FAILED\n");
        }
        printf("LAPACK_dpotrf took %.6lf seconds\n", tm);
    } else {
        printf("Verification disabled\n");
    }

    // Save the trace.
    printf("Saving the trace...\n");
    pcp_view_trace_tikz();

    // View statistics.
    printf("Printing statistics...\n");
    pcp_view_statistics_stdout();

    // Stop the runtime system.
    pcp_stop();

    // Clean up.
    free(A);
    free(Ain);

    return EXIT_SUCCESS;
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
#ifndef __GSLPRINT_H__
#define __GSLPRINT_H__

#include <stdio.h>
#include <iostream>
#include <math.h>
#include <fstream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_integration.h>

#endif

void gsl_matrix_fprint(FILE * stream, gsl_matrix * target,int rows,int columns, char const * format);
void gsl_vector_fprint(FILE * stream, gsl_vector * target,int length, char const * format);
void gsl_matrix_printf(gsl_matrix * m, int rows, int columns, char const * format);
void gsl_vector_printf(gsl_vector * m,int length, char const * format);
void gsl_vector_complex_fprint(FILE * stream, gsl_vector_complex * target, int length, char const * format);
void gsl_matrix_complex_fprint(FILE * stream, gsl_matrix_complex * target, int rows, int columns, char const * format);
void gsl_matrix_complex_printf(gsl_matrix_complex * target, int rows, int columns, char const * format);
void gsl_vector_complex_fprint(gsl_vector_complex * target, int length, char const * format);

max_stars_repo_name: Walter-Feng/Hartree-Fock

max_issues_repo_name: Walter-Feng/HomebrewLib

max_forks_repo_name: Walter-Feng/HomebrewLib