text:
#ifndef ASSOCIATION_H
#define ASSOCIATION_H

#define TRUE 1
#define FALSE 0

#include <iostream>
// #include <omp.h>
#include <mkl.h>

extern "C" {
#include <gsl/gsl_rng.h>
#include <gsl/gsl_const_mksa.h>
}
#include "matrix.h"
#include "connectivity.h"
#include "potential.h"
/* #include "trajectory.h" */
/* #include "read_file_condition.h" */
#include "file_IO.h"
#include "geometry.h"

class ASSOCIATION
: public CONNECTIVITY
{
 public:
  MKL_LONG Nc; // number of chains per micelle
  MKL_LONG Tec; // tolerance of end chains
  MKL_LONG N_min; // 2Nc - Tec
  MKL_LONG N_max; // 2Nc + Tec
  MKL_LONG N_ASSOCIATION;
  bool MULTIPLE_CONNECTIONS;
  bool NUMBER_RESTRICTION;
  bool PROBABILITY_AREA;

  MATRIX N_ASSOCIATION_PARTICLE;
  MATRIX *weight;
  // related with suggestion probability (selecting a degenerated pair)
  MATRIX *CASE_SUGGESTION;
  MATRIX *dPDF_SUGGESTION;
  MATRIX *dCDF_SUGGESTION;
  double *Z_SUGGESTION;

  // related with association probability
  MATRIX *dCDF_ASSOCIATION;
  /* size_t **INDEX_ASSOCIATION; */
  MATRIX *INDEX_ASSOCIATION;
  MKL_LONG *TOKEN_ASSOCIATION;

  // member variables

  /* MKL_LONG flag_other, flag_new, flag_itself, flag_hash_other; */
  static const MKL_LONG flag_other = 0;
  static const MKL_LONG flag_new = 1;
  static const MKL_LONG flag_itself = 2;
  static const MKL_LONG flag_hash_other = 3;
  static const MKL_LONG flag_hash_backtrace = 4;

  bool
    (*ADD_ASSOCIATION)
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    (*OPP_DEL_ASSOCIATION)
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    (*CANCEL_ASSOCIATION)
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    (*CHECK_N_ADD_ASSOCIATION)
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    (*CHECK_N_MOV_ASSOCIATION)
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    (*CHECK_N_OPP_DEL_ASSOCIATION)
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);

  MKL_LONG
    N_TOTAL_ASSOCIATION
    ();

  MKL_LONG
  N_PARTICLE_ASSOCIATION
  (MKL_LONG index_particle);

  MKL_LONG
    N_CONNECTED_ENDS
    (MKL_LONG given_index);
  MKL_LONG
    N_TOTAL_CONNECTED_ENDS
    ();
  bool
    CHECK_EXIST
    (MKL_LONG index_A, MKL_LONG index_B);
  bool
    CHECK_EXIST_1D
    (MKL_LONG index_A, MKL_LONG index_B);
  MKL_LONG
    add_association
    (const MKL_LONG index_particle, MKL_LONG index_target);
  MKL_LONG
    add_association_INFO
    (POTENTIAL_SET& POTs, const MKL_LONG index_particle, MKL_LONG index_target, double distance);
  MKL_LONG
    opp_del_association
    (const MKL_LONG index_particle, MKL_LONG index_target);
  MKL_LONG
    opp_del_association_hash
    (const MKL_LONG index_particle, MKL_LONG hash_index_target);
  MKL_LONG
    opp_del_association_IK
    (MKL_LONG index_I, MKL_LONG hash_index_K);
  MKL_LONG
    opp_del_association_grab_IK
    (MKL_LONG index_I, MKL_LONG hash_index_K);

  MKL_LONG
    GET_INDEX_HASH_FROM_ROLL
    (const MKL_LONG index_particle, double rolled_p);
  MKL_LONG
    GET_HASH_FROM_ROLL
    (const MKL_LONG index_particle, double rolled_p);
  MKL_LONG
    FIND_HASH_INDEX
    (MKL_LONG index_A, MKL_LONG index_B);
  MKL_LONG
  update_N_association
  ();

  MKL_LONG
  get_N_association_particle
  (MKL_LONG index_particle);

  // At the moment, there is no namespace dependence between KINETICS::NORMALIZED and KINETICS::METROPOLIS.
  // The dependency is set with the w_function and transtion_probability that is developed with lib_potential.
  // For better performance, it would be better to tune the following with specific dependencies (future works).
  double
    update_CASE_SUGGESTION_particle_hash_target
    (POTENTIAL_SET& POTs, const MKL_LONG index_particle, MKL_LONG hash_index_target, double distance);

  double
    update_CASE_SUGGESTION_particle_target
    (POTENTIAL_SET& POTs, const MKL_LONG index_particle, MKL_LONG index_target, double distance);

  double
    update_Z_SUGGESTION_particle
    (const MKL_LONG index_particle);
  double
    update_dPDF_SUGGESTION_particle
    (const MKL_LONG index_particle);
  double
    update_dCDF_SUGGESTION_particle
    (const MKL_LONG index_particle);
  double
    update_CHAIN_SUGGESTION_MAP_particle
    (const MKL_LONG index_particle, POTENTIAL_SET& POTs, RDIST& R_boost);

  double
    update_ASSOCIATION_MAP_particle
    (const MKL_LONG index_particle, POTENTIAL_SET& POTs, RDIST& R_boost);

  MKL_LONG
    return_multiplicity
    (const MKL_LONG index_particle, const MKL_LONG target_particle)
  // this additional function is useful for tracking individual chain
  {
    for(MKL_LONG j=0; j<TOKEN[index_particle]; j++)
      {
        if(HASH[index_particle](j) == target_particle)
          return weight[index_particle](j);
      }
    return 0;
  }

 ASSOCIATION()
   : CONNECTIVITY(){}

  ASSOCIATION
    (COND& given_condition);

  ASSOCIATION
    (MKL_LONG number_of_particles, MKL_LONG number_of_chains_per_particles, MKL_LONG tolerance_connection, bool ALLOWING_MULTIPLE_CONNECTIONS);

  MKL_LONG
    dynamic_alloc
    ();
  MKL_LONG
    set_initial_condition
    ();
  MKL_LONG
    initial
    ();
  MKL_LONG
    initial_inheritance
    ();
  virtual
    ~ASSOCIATION()
    {
      delete[] weight;

      // related with suggestion probability
      delete[] CASE_SUGGESTION;
      delete[] dPDF_SUGGESTION;
      delete[] dCDF_SUGGESTION;
      delete[] Z_SUGGESTION;
      // related with association probability
      delete[] dCDF_ASSOCIATION;
      delete[] INDEX_ASSOCIATION;

      delete[] TOKEN_ASSOCIATION;
    }

  MKL_LONG
    read_exist_weight
  (const char* fn_weight, MKL_LONG N_steps);

};

MKL_LONG
index_set_val
(size_t* index, MKL_LONG given_val, MKL_LONG size_of_arr);

namespace TRUTH_MAP
{
  namespace MULTIPLE
  {
    bool
      CHECK_N_ADD_BOOST
      (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
    bool
      CHECK_N_OPP_DEL_BOOST
      (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
    bool
      CHECK_N_MOV_BOOST
      (ASSOCIATION& CONNECT, MKL_LONG index_set[]);

    namespace NO_RESTRICTION // this option is used when there are no restriction for number of attached chain ends
    {
      bool
	return_TRUE
	(ASSOCIATION& CONNECT, MKL_LONG index_set[]);

    }
  }

  namespace SINGLE
  {
    bool
      CHECK_N_ADD_BOOST
      (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
    bool
      CHECK_N_MOV_BOOST
      (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
    bool
      CHECK_N_OPP_DEL_BOOST
      (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  }
  bool
    CANCEL_ASSOCIATION_BOOST
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    ADD_ASSOCIATION_BOOST
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);
  bool
    OPP_DEL_ASSOCIATION_BOOST
    (ASSOCIATION& CONNECT, MKL_LONG index_set[]);

  bool
    OPP_DEL_ASSOCIATION_BASIC
    (ASSOCIATION& CONNECT, MKL_LONG index_itself, MKL_LONG index_target, MKL_LONG index_new);

  bool
    NEW_ASSOCIATION_BASIC
    (ASSOCIATION& CONNECT, MKL_LONG index_itself, MKL_LONG index_target, MKL_LONG index_new);
  bool
    NEW_ASSOCIATION_SINGLE
    (ASSOCIATION& CONNECT, MKL_LONG index_itself, MKL_LONG index_target, MKL_LONG index_new);

  bool
    MOV_ASSOCIATION_BASIC
    (ASSOCIATION& CONNECT, MKL_LONG index_itself, MKL_LONG index_target, MKL_LONG index_new);
  bool
    MOV_ASSOCIATION_SINGLE
    (ASSOCIATION& CONNECT, MKL_LONG index_itself, MKL_LONG index_target, MKL_LONG index_new);
}

namespace MICELLE
{
  double
  update_SOFT_REPULSION_P2(POTENTIAL_SET& given_POT, ASSOCIATION& CONNECT);

  double
  update_SOFT_REPULSION_P2_particle(POTENTIAL_SET& given_POT, ASSOCIATION& CONNECT, MKL_LONG index_particle);

}

namespace KINETICS
{
  double
    CONNECTIVITY_update_dCDF_SUGGESTION_particle
    (ASSOCIATION* CONNECT, const MKL_LONG index_particle);
  double
    CONNECTIVITY_update_CASE_SUGGESTION_particle_hash_target
    (ASSOCIATION* CONNECT, POTENTIAL_SET& POTs, const MKL_LONG index_particle, MKL_LONG hash_index_target, double distance);
  double
    CONNECTIVITY_update_CASE_SUGGESTION_particle_target
    (ASSOCIATION* CONNECT, POTENTIAL_SET& POTs, const MKL_LONG index_particle, MKL_LONG index_target, double distance);

  double
    CONNECTIVITY_update_Z_SUGGESTION_particle
    (ASSOCIATION* CONNECT, const MKL_LONG index_particle);
  double
    CONNECTIVITY_update_dPDF_SUGGESTION_particle
    (ASSOCIATION* CONNECT, const MKL_LONG index_particle);
}

class CHAIN_INFORMATION
{
  /*
    The CHAIN class is using CHAIN_NODE without inheritance.
    It contains dynamically allocated CHAIN_NODE.
  */
 public:

  CHAIN_NODE *CHAIN;
  MKL_LONG N_chains;
  MKL_LONG N_particles;

  bool INITIALIZATION;

  MKL_LONG&
    HEAD
    (MKL_LONG given_chain_index)
    {
      return CHAIN[given_chain_index].HEAD;
    }
  MKL_LONG&
    TAIL
    (MKL_LONG given_chain_index)
    {
      return CHAIN[given_chain_index].TAIL;
    }

  /*
    On the same reason, the ATTACHED will return pointer variable rather than reference variable. This break the existing interface.
  */
  MKL_LONG&
    ATTACHED
    (MKL_LONG given_chain_index, MKL_LONG flag_HEAD_TAIL)
    {
      return CHAIN[given_chain_index].index(flag_HEAD_TAIL);
    }

  MKL_LONG
    mov_attachment
    (MKL_LONG target_particle, MKL_LONG given_chain_index, MKL_LONG flag_HEAD_TAIL)
  {
    /*
      flag_HEAD_TAIL is 0 for HEAD while 1 for TAIL.
    */
    ATTACHED(given_chain_index, flag_HEAD_TAIL) = target_particle;

    return 0;
  }

  MKL_LONG
    add_attachment
    (const MKL_LONG index_particle, MKL_LONG given_chain_index, MKL_LONG flag_HEAD_TAIL)
  {
    // there is no distingushable between add and mov attachment for chain node
    mov_attachment(index_particle, given_chain_index, flag_HEAD_TAIL);
    return 0;
  }

  MKL_LONG
    del_attachment
    (MKL_LONG given_chain_index, MKL_LONG flag_HEAD_TAIL)
  {
    // it make detachment for subjected chain ends
    // then attach to particles which other chain ends are attached
    mov_attachment(ATTACHED(given_chain_index, (flag_HEAD_TAIL + 1)%2), given_chain_index, flag_HEAD_TAIL);
    return 0;
  }

  MKL_LONG
    initial
    (MKL_LONG number_of_chains, MKL_LONG number_of_particles)
  {
    printf("ST:initial:CHAIN_INFORMATION\n");
    N_chains = number_of_chains;
    N_particles = number_of_particles;
    CHAIN = new CHAIN_NODE [N_chains];

    for(MKL_LONG i=0; i<N_chains; i++)
      {
        /*
          it allocate the equivalent number of chain ends per particles
        */
        HEAD(i) = i%N_particles;
        TAIL(i) = i%N_particles;
      }
    INITIALIZATION = TRUE;
    return INITIALIZATION;
  }

  MKL_LONG
    initial
    (ASSOCIATION& CONNECT)
  {
    N_chains = CONNECT.Nc*CONNECT.Np;
    N_particles = CONNECT.Np;
    CHAIN = new CHAIN_NODE [N_chains];
    MKL_LONG cnt = 0;

    for(MKL_LONG i=0; i<N_particles; i++)
      {
        for(MKL_LONG j=i; j<N_particles; j++) // it includes i=j which is loop and prevents j<i means all the pair is once accounted.
          {
            MKL_LONG wij = CONNECT.return_multiplicity(i, j);
            if(i==j)
              {
                // note that weight count chain end rather than chain
                // when i != j, the duplication is removed because we only count the pair with condition j>i
                // when i == j, the duplication happens in the wij.
                wij /=2;
              }
            for(MKL_LONG w=0; w<wij; w++)
              {
                HEAD(cnt) = i;
                TAIL(cnt) = j;
                if(cnt++ > N_chains)
                  printf("Err: counted number of chain overflood\n");
              }
          }
      }
    printf("end_particle_allocation\n");
    INITIALIZATION = TRUE;
    return INITIALIZATION;
  }

  MKL_LONG
    inheritance_chain
    (COND& given_condition, ASSOCIATION& CONNECT)
  {
    N_chains = CONNECT.Nc*CONNECT.Np;
    N_particles = CONNECT.Np;
    /* CHAIN = (CHAIN_NODE*)mkl_malloc(N_chains*sizeof(CHAIN_NODE), BIT); */
    // the problem of the original is that 'malloc' will not call the constructor
    CHAIN = new CHAIN_NODE [N_chains];

    MKL_LONG cnt = 0;
    ifstream GIVEN_CHAIN_FILE;
    GIVEN_CHAIN_FILE.open(given_condition("CONTINUATION_CHAIN_FN").c_str());
    string line;
    while(getline(GIVEN_CHAIN_FILE, line))
      cnt ++;
    GIVEN_CHAIN_FILE.clear();

    GIVEN_CHAIN_FILE.seekg(0);
    MKL_LONG N_steps = atoi(given_condition("CONTINUATION_STEP").c_str());
    if(N_steps == -1)
      {
        for(MKL_LONG i=0; i<cnt-1; i++)
          getline(GIVEN_CHAIN_FILE, line);
      }
    else
      {
        for(MKL_LONG i=0; i<N_steps; i++)
          getline(GIVEN_CHAIN_FILE, line);
      }

    for(MKL_LONG k=0; k<N_chains; k++)
      GIVEN_CHAIN_FILE >> HEAD(k); // index from 0 to N_chains-1 : HEAD

    for(MKL_LONG k=0; k<N_chains; k++)
      GIVEN_CHAIN_FILE >> TAIL(k); // index from N_chains to 2*N_chains-1 : TAIL

    GIVEN_CHAIN_FILE.close();
    INITIALIZATION = TRUE;
    return INITIALIZATION;
  }

  CHAIN_INFORMATION()
    {
      std::cout << "ERR: CHAIN Class must have initialization argument\n";
    }
  CHAIN_INFORMATION(MKL_LONG number_of_chains, MKL_LONG number_of_particles)
    {
      initial(number_of_chains, number_of_particles);
    }
  CHAIN_INFORMATION(COND& given_condition)
    {
      if(given_condition("tracking_individual_chain") == "TRUE")
        {
          if(given_condition("CONTINUATION_CHAIN") == "TRUE")
            {
              printf("ERR: Without CONTINUATION TOPOLOGICAL INFORMATION (.hash, .weight), the tracking individual chain cannot inherit existing .chain file information\n");
            }
          else
            {
              printf("ST:Constructor:CHAIN_INFORMATION\n");
              initial(atoi(given_condition("N_chains_per_particle").c_str())*atoi(given_condition("Np").c_str()), atoi(given_condition("Np").c_str()));
              printf("END:Constructor:CHAIN_INFORMATION\n");
            }
        }
    }
  CHAIN_INFORMATION(COND& given_condition, ASSOCIATION& CONNECT)
    // it initialize based on the CONNECT information
    {
      if(given_condition("tracking_individual_chain") == "TRUE")
        {
          if(given_condition("CONTINUATION_CHAIN") == "TRUE")
            {
              inheritance_chain(given_condition, CONNECT);
            }
          else
            {
              printf("ST:Constructor:CHAIN_INFORMATION\n");
              initial(CONNECT);
              printf("END:Constructor:CHAIN_INFORMATION\n");
            }
        }
    }
  virtual ~CHAIN_INFORMATION()
    {
      if(INITIALIZATION)
        delete[] CHAIN;
    }
};

// double inline MICELLE::update_SOFT_REPULSION_P2_particle(POTENTIAL_SET& given_POT, ASSOCIATION& CONNECT, MKL_LONG index_particle)
// {
//   given_POT.force_variables.C_rep_arr(index_particle) = given_POT.force_variables.repulsion_coefficient_base*CONNECT.N_ASSOCIATION_PARTICLE(index_particle)*CONNECT.N_ASSOCIATION_PARTICLE(index_particle);
//   // C_rep(p) = C_0 * p^2

//   return given_POT.force_variables.C_rep_arr(index_particle);
// }

#endif


text:
//====---- Sudoku/Board_Section_iterator.h                            ----====//
//
// Location aware iterator implementation for Sections
//====--------------------------------------------------------------------====//
//
// The iterator can return a Location object.
// Included by Board_Section.h
//
//====--------------------------------------------------------------------====//
#pragma once

#include "Sudoku/Location.h"
#include "Sudoku/Location_Utilities.h"
#include "Sudoku/traits.h"
#include <gsl/gsl> // index

#include <iterator>
#include <type_traits>

#include "Board.fwd.h" // Forward declarations

#include <cassert>

namespace Sudoku::Board_Section
{
template<typename T, int N, Section S, bool is_const, bool is_reverse>
class Section_iterator;

//====--- Aliases --------------------------------------------------------====//
template<typename T, int N>
using Row_iterator = Section_iterator<T, N, Section::row, false, false>;
template<typename T, int N>
using Col_iterator = Section_iterator<T, N, Section::col, false, false>;
template<typename T, int N>
using Block_iterator = Section_iterator<T, N, Section::block, false, false>;
template<typename T, int N>
using const_Row_iterator = Section_iterator<T, N, Section::row, true, false>;
template<typename T, int N>
using const_Col_iterator = Section_iterator<T, N, Section::col, true, false>;
template<typename T, int N>
using const_Block_iterator =
	Section_iterator<T, N, Section::block, true, false>;
template<typename T, int N>
using reverse_Row_iterator = Section_iterator<T, N, Section::row, false, true>;
template<typename T, int N>
using reverse_Col_iterator = Section_iterator<T, N, Section::col, false, true>;
template<typename T, int N>
using reverse_Block_iterator =
	Section_iterator<T, N, Section::block, false, true>;
template<typename T, int N>
using const_reverse_Row_iterator =
	Section_iterator<T, N, Section::row, true, true>;
template<typename T, int N>
using const_reverse_Col_iterator =
	Section_iterator<T, N, Section::col, true, true>;
template<typename T, int N>
using const_reverse_Block_iterator =
	Section_iterator<T, N, Section::block, true, true>;

//====--------------------------------------------------------------------====//
template<typename T, int N, Section S, bool is_const, bool is_reverse>
class Section_iterator
{
	using owner_type =
		std::conditional_t<is_const, Board<T, N> const, Board<T, N>>;
	using Location       = ::Sudoku::Location<N>;
	using Location_Block = ::Sudoku::Location_Block<N>;

public:
	// member types
	using difference_type   = std::ptrdiff_t;
	using value_type        = T;
	using pointer           = std::conditional_t<is_const, T const*, T*>;
	using reference         = std::conditional_t<is_const, T const&, T&>;
	using iterator_category = std::random_access_iterator_tag;

	// Constructors
	constexpr Section_iterator() noexcept
	{ // defaults to [r]begin()
		if constexpr (is_reverse)
			elem_ = elem_size<N> - 1;
	}
	template<
		typename idT,
		typename elT,
		typename = std::enable_if_t<is_int_v<idT> && is_int_v<elT>>>
	constexpr Section_iterator(owner_type* owner, idT id, elT elem) noexcept
		: board_(owner), id_(id), elem_(elem)
	{
		assert(is_valid_size<N>(id));
	}
	// [[implicit]] Type Conversion to const_(reverse_)iterator
	// NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
	constexpr operator Section_iterator<T, N, S, true, is_reverse>()
		const noexcept // (guide format)
	{
		return Section_iterator<T, N, S, true, is_reverse>(board_, id_, elem_);
	}

	//====----------------------------------------------------------------====//
	[[nodiscard]] constexpr Location location() const noexcept
	{
		assert(is_valid_size<N>(elem_));
		using ::Sudoku::Board_Section::Section;
		switch (S)
		{
		case Section::row: return Location(id_, elem_); break;
		case Section::col: return Location(elem_, id_); break;
		case Section::block: return Location_Block(id_, elem_); break;
		}
	}
	// Allows for conversion to Location
	explicit constexpr operator Location() const noexcept { return location(); }

	//====----------------------------------------------------------------====//
	[[nodiscard]] constexpr reference operator*() const noexcept
	{ // Only valid in a dereferenceable location:
		assert(board_ != nullptr && is_valid_size<N>(elem_));
		return board_->operator[](location());
	}
	[[nodiscard]] constexpr pointer operator->() const noexcept
	{ // member access; equivalent to (*p).member
		return std::pointer_traits<pointer>::pointer_to(**this);
	}

	//====----------------------------------------------------------------====//
	constexpr Section_iterator& operator++() noexcept;
	constexpr Section_iterator operator++(int) noexcept;
	constexpr Section_iterator& operator--() noexcept;
	constexpr Section_iterator operator--(int) noexcept;
	constexpr Section_iterator&
		operator+=(const difference_type offset) noexcept;
	constexpr Section_iterator&
		operator-=(const difference_type offset) noexcept;

	[[nodiscard]] friend constexpr difference_type operator-(
		Section_iterator const& left, Section_iterator const& right) noexcept
	{ // difference
		assert(is_same_Section(left, right));
		if constexpr (is_reverse)
			return right.elem_ - left.elem_;
		else
			return left.elem_ - right.elem_;
	}

	[[nodiscard]] constexpr reference
		operator[](const difference_type offset) const noexcept
	{
		return (*(*this + offset));
	}

	//====----------------------------------------------------------------====//
	[[nodiscard]] friend constexpr bool operator==(
		Section_iterator const& left, Section_iterator const& right) noexcept
	{
		assert(is_same_Section(left, right));
		return is_same_Section(left, right) && left.elem_ == right.elem_;
	}
	[[nodiscard]] friend constexpr bool operator<(
		Section_iterator const& left, Section_iterator const& right) noexcept
	{
		assert(is_same_Section(left, right));
		if constexpr (is_reverse)
			return left.elem_ > right.elem_;
		else
			return left.elem_ < right.elem_;
	}

private:
	owner_type* board_{nullptr};
	gsl::index id_{0};
	gsl::index elem_{0};

	friend constexpr bool is_same_Section(
		Section_iterator const& A, Section_iterator const& B) noexcept
	{ // compare address of:
		return A.board_ == B.board_ && A.id_ == B.id_;
	}
};

//====--------------------------------------------------------------------====//
template<typename T, int N, Section S, bool C, bool is_reverse>
constexpr Section_iterator<T, N, S, C, is_reverse>&
	Section_iterator<T, N, S, C, is_reverse>::operator++() noexcept
{ // pre-increment
	assert(board_ != nullptr);
	if constexpr (is_reverse)
	{
		assert(elem_ >= 0);
		--elem_;
	}
	else
	{
		assert(elem_ < elem_size<N>);
		++elem_;
	}
	return (*this);
}
template<typename T, int N, Section S, bool C, bool R>
constexpr Section_iterator<T, N, S, C, R> // NOLINTNEXTLINE(readability/casting)
	Section_iterator<T, N, S, C, R>::operator++(int) noexcept
{ // post-increment
	const Section_iterator pre{*this};
	operator++();
	return pre;
}
template<typename T, int N, Section S, bool C, bool is_reverse>
constexpr Section_iterator<T, N, S, C, is_reverse>&
	Section_iterator<T, N, S, C, is_reverse>::operator--() noexcept
{ // pre-decrement
	assert(board_ != nullptr);
	if constexpr (is_reverse)
	{
		assert(elem_ < elem_size<N> - 1);
		++elem_;
	}
	else
	{
		assert(elem_ > 0);
		--elem_;
	}
	return (*this);
}
template<typename T, int N, Section S, bool C, bool R>
constexpr Section_iterator<T, N, S, C, R> // NOLINTNEXTLINE(readability/casting)
	Section_iterator<T, N, S, C, R>::operator--(int) noexcept
{ // post-decrement
	const Section_iterator pre{*this};
	operator--();
	return pre;
}
template<typename T, int N, Section S, bool C, bool is_reverse>
constexpr Section_iterator<T, N, S, C, is_reverse>&
	Section_iterator<T, N, S, C, is_reverse>::operator+=(
		const difference_type offset) noexcept
{
	assert(offset == 0 || board_ != nullptr);
	if constexpr (is_reverse)
	{
		elem_ -= offset;
		assert(elem_ >= -1);
		assert(elem_ < elem_size<N>);
	}
	else
	{
		elem_ += offset;
		assert(elem_ >= 0);
		assert(elem_ <= elem_size<N>);
	}
	return (*this);
}
template<typename T, int N, Section S, bool C, bool is_reverse>
constexpr Section_iterator<T, N, S, C, is_reverse>&
	Section_iterator<T, N, S, C, is_reverse>::operator-=(
		const difference_type offset) noexcept
{
	return operator+=(-offset);
}

//====--------------------------------------------------------------------====//
// Free functions (non-friend)
//====--------------------------------------------------------------------====//
template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] constexpr auto operator+(
	Section_iterator<T, N, S, is_const, is_reverse> left,
	typename Section_iterator<T, N, S, is_const, is_reverse>::
		difference_type const offset) noexcept
{ // itr + offset
	return (left += offset);
}
template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr auto operator+(
	typename Section_iterator<T, N, S, is_const, is_reverse>::
		difference_type const offset,
	Section_iterator<T, N, S, is_const, is_reverse> itr) noexcept
{ // offset + itr
	return (itr += offset);
}
template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] constexpr auto operator-(
	Section_iterator<T, N, S, is_const, is_reverse> left,
	typename Section_iterator<T, N, S, is_const, is_reverse>::
		difference_type const offset) noexcept
{ // itr - offset
	return (left += -offset);
}

template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator!=(
	Section_iterator<T, N, S, is_const, is_reverse> const& left,
	Section_iterator<T, N, S, is_const, is_reverse> const& right) noexcept
{
	return !(left == right);
}
template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator>(
	Section_iterator<T, N, S, is_const, is_reverse> const& left,
	Section_iterator<T, N, S, is_const, is_reverse> const& right) noexcept
{
	return (right < left);
}
template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator<=(
	Section_iterator<T, N, S, is_const, is_reverse> const& left,
	Section_iterator<T, N, S, is_const, is_reverse> const& right) noexcept
{
	return !(right < left);
}
template<typename T, int N, Section S, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator>=(
	Section_iterator<T, N, S, is_const, is_reverse> const& left,
	Section_iterator<T, N, S, is_const, is_reverse> const& right) noexcept
{
	return !(left < right);
}

} // namespace Sudoku::Board_Section

ext: h

max_stars_repo_name: Farwaykorse/fwkSudoku

max_issues_repo_name: Farwaykorse/fwkSudoku

max_forks_repo_name: Farwaykorse/fwkSudoku


text:
#include <stdio.h>
#include <cblas.h>
#include <sys/time.h>
#include <stdlib.h>
#include <omp.h>

#define NUM 64

#define GEMM_K 320
#define GEMM_N 4096
#define GEMM_M 256

void PACKA(float* A, float* Ac, long M, long K, long LK)
{
	long ii, jj, kk;
	for( ii = 0 ; ii < M; ii = ii + 8)
	{

		float *temp = A + ii * LK + kk;

		asm volatile(
				"	ldr		x0, %[Ac]							\n"
				"	ldr		x1, %[K]							\n"
				"	ldr		x2, %[temp]						\n"
				"	ldr		x30, %[LK]						\n"

				"	add		x3, x2, x30, lsl #2		\n"
				"	add		x4, x2, x30, lsl #3		\n"
				"	add		x5, x3, x30, lsl #3		\n"
				"	add		x6, x4, x30, lsl #3		\n"
				"	add		x7, x5, x30, lsl #3		\n"
				"	add		x8, x6, x30, lsl #3		\n"
				"	add		x9, x7, x30, lsl #3		\n"

				"	lsr		x21, x1, #3						\n"
				"	cmp		x21, #0								\n"
				"	beq		PACKA_END							\n"

				"PACKA:												\n"

				"	prfm	PLDL1KEEP, [x2, #128]	\n"
				"	prfm	PLDL1KEEP, [x3, #128]	\n"

				"	ldr		q0, [x2], #16					\n"
				"	ldr		q1, [x3], #16					\n"
				"	ldr		q2, [x4], #16					\n"
				"	ldr		q3, [x5], #16					\n"

				"	prfm	PLDL1KEEP, [x4, #128]	\n"
				"	prfm	PLDL1KEEP, [x5, #128]	\n"

				"	st4		{v0.s, v1.s, v2.s, v3.s}[0], [x0], #16	\n"

				"	ldr		q4, [x6], #16					\n"
				"	ldr		q5, [x7], #16					\n"
				"	ldr		q6, [x8], #16					\n"
				"	ldr		q7, [x9], #16					\n"

				"	st4		{v4.s, v5.s, v6.s, v7.s}[0], [x0], #16	\n"

				"	prfm	PLDL1KEEP, [x6, #128]	\n"
				"	prfm	PLDL1KEEP, [x7, #128]	\n"

				"	ldr		q8, [x2], #16					\n"
				"	st4		{v0.s, v1.s, v2.s, v3.s}[1], [x0], #16	\n"
				"	ldr		q9, [x3], #16					\n"
				"	st4		{v4.s, v5.s, v6.s, v7.s}[1], [x0], #16	\n"
				"	ldr		q10, [x4], #16				\n"
				"	st4		{v0.s, v1.s, v2.s, v3.s}[2], [x0], #16	\n"
				"	ldr		q11, [x5], #16				\n"
				"	st4		{v4.s, v5.s, v6.s, v7.s}[2], [x0], #16	\n"

				"	prfm	PLDL1KEEP, [x8, #128]	\n"
				"	prfm	PLDL1KEEP, [x9, #128]	\n"

				"	ldr		q12, [x6], #16				\n"
				"	st4		{v0.s, v1.s, v2.s, v3.s}[3], [x0], #16		\n"
				"	ldr		q13, [x7], #16				\n"
				"	st4		{v4.s, v5.s, v6.s, v7.s}[3], [x0], #16		\n"
				"	ldr		q14, [x8], #16				\n"
				"	st4		{v8.s, v9.s, v10.s, v11.s}[0], [x0], #16	\n"
				"	ldr		q15, [x9], #16				\n"
				"	st4		{v12.s, v13.s, v14.s, v15.s}[0], [x0], #16	\n"

				"	subs	x21, x21, #1					\n"

				"	st4		{v8.s, v9.s, v10.s, v11.s}[1], [x0], #16	\n"
				"	st4		{v12.s, v13.s, v14.s, v15.s}[1], [x0], #16	\n"
				"	st4		{v8.s, v9.s, v10.s, v11.s}[2], [x0], #16	\n"
				"	st4		{v12.s, v13.s, v14.s, v15.s}[2], [x0], #16	\n"
				"	st4		{v8.s, v9.s, v10.s, v11.s}[3], [x0], #16	\n"
				"	st4		{v12.s, v13.s, v14.s, v15.s}[3], [x0], #16	\n"

				"	bgt		PACKA 						\n"

				"	ands	x22, x1, #7				\n"
				"	beq		PACKA_END					\n"

				"	cmp		x22, #4						\n"
				"	blt 	K1_PACKA					\n"

				"K4_PACKA:								\n"

				"	ldr		q0, [x2], #16			\n"
				"	ldr		q1, [x3], #16			\n"
				"	ldr		q2, [x4], #16			\n"
				"	ldr		q3, [x5], #16			\n"

				"	st4		{v0.s, v1.s, v2.s, v3.s}[0], [x0], #16	\n"
				"	ldr		q4, [x6], #16			\n"
				"	st4		{v0.s, v1.s, v2.s, v3.s}[1], [x0], #16	\n"
				"	ldr		q5, [x7], #16			\n"
				"	st4		{v0.s, v1.s, v2.s, v3.s}[2], [x0], #16	\n"
				"	ldr		q6, [x8], #16			\n"
				"	st4		{v0.s, v1.s, v2.s, v3.s}[3], [x0], #16	\n"
				"	ldr		q7, [x9], #16			\n"

				"	st4		{v4.s, v5.s, v6.s, v7.s}[0], [x0], #16	\n"
				"	st4		{v4.s, v5.s, v6.s, v7.s}[1], [x0], #16	\n"
				"	st4		{v4.s, v5.s, v6.s, v7.s}[2], [x0], #16	\n"
				"	st4		{v4.s, v5.s, v6.s, v7.s}[3], [x0], #16	\n"

				"	subs	x22, x22, #4			\n"
				"	beq		PACKA_END					\n"

				"K1_PACKA:								\n"

				"	ldr		s0, [x2], #4			\n"
				"	ldr		s1, [x3], #4			\n"
				"	ldr		s2, [x4], #4			\n"
				"	ldr		s3, [x5], #4			\n"

				"	subs 	x22, x22, #1			\n"

				"	st4		{v0.s, v1.s, v2.s, v3.s}[0], [x0], #16	\n"
				"	ldr		s4, [x6], #4			\n"
				"	ldr		s5, [x7], #4			\n"
				"	ldr		s6, [x8], #4			\n"
				"	ldr		s7, [x9], #4			\n"

				"	st4		{v4.s, v5.s, v6.s, v7.s}[0], [x0], #16	\n"

				"	bgt		K1_PACKA				\n"

				"PACKA_END:							\n"

				:
				:
				[temp] "m" (temp),
				[Ac] "m" (Ac),
				[K] "m" (K),
				[LK] "m" (LK)
				:"x0", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8",
				"x9", "x10", "x11", "x12", "x13","x14", "x15", "x16",
				"x17", "x18", "x19", "x20", "x21", "x22", "x23", "x24","x25",
				"x26", "x27", "x28", "x29", "x30",
				"v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7",
				"v8", "v9", "v10", "v11", "v12", "v13", "v14", "v15",
				"v16", "v17", "v18", "v19", "v20", "v21", "v22", "v23",
				"v24", "v25", "v26", "v27", "v28", "v29", "v30", "v31"

		);

		Ac = Ac + K * 8;

	}

}

void Sin_PACK(float* A, float* Ac, long M, long K, long LK)
{
	long ii, jj, kk;
	long Kc;

	for( kk =0 ; kk < K; kk = kk + Kc)
	{

		Kc = GEMM_K;
		if(K - kk < GEMM_K)
			Kc= K - kk;

		for( ii = 0 ; ii < M; ii = ii + 8)
		{

			float *temp = A + ii * LK + kk;

			asm volatile(
					"	ldr		x0, %[Ac]							\n"
					"	ldr		x1, %[Kc]							\n"
					"	ldr		x2, %[temp]						\n"
					"	ldr		x30, %[LK]						\n"

					"	add		x3, x2, x30, lsl #2		\n"
					"	add		x4, x2, x30, lsl #3		\n"
					"	add		x5, x3, x30, lsl #3		\n"
					"	add		x6, x4, x30, lsl #3		\n"
					"	add		x7, x5, x30, lsl #3		\n"
					"	add		x8, x6, x30, lsl #3		\n"
					"	add		x9, x7, x30, lsl #3		\n"

					"	lsr		x21, x1, #3						\n"
					"	cmp		x21, #0								\n"
					"	beq		Sin_PACKA_END					\n"

					"Sin_PACKA:										\n"

					"	prfm	PLDL1KEEP, [x2, #128]	\n"
					"	prfm	PLDL1KEEP, [x3, #128]	\n"

					"	ldr		q0, [x2], #16					\n"
					"	ldr		q1, [x3], #16					\n"
					"	ldr		q2, [x4], #16					\n"
					"	ldr		q3, [x5], #16					\n"

					"	prfm	PLDL1KEEP, [x4, #128]	\n"
					"	prfm	PLDL1KEEP, [x5, #128]	\n"

					"	st4		{v0.s, v1.s, v2.s, v3.s}[0], [x0], #16	\n"

					"	ldr		q4, [x6], #16					\n"
					"	ldr		q5, [x7], #16					\n"
					"	ldr		q6, [x8], #16					\n"
					"	ldr		q7, [x9], #16					\n"

					"	st4		{v4.s, v5.s, v6.s, v7.s}[0], [x0], #16	\n"

					"	prfm	PLDL1KEEP, [x6, #128]	\n"
					"	prfm	PLDL1KEEP, [x7, #128]	\n"

					"	ldr		q8, [x2], #16					\n"
					"	st4		{v0.s, v1.s, v2.s, v3.s}[1], [x0], #16	\n"
					"	ldr		q9, [x3], #16					\n"
					"	st4		{v4.s, v5.s, v6.s, v7.s}[1], [x0], #16	\n"
					"	ldr		q10, [x4], #16				\n"
					"	st4		{v0.s, v1.s, v2.s, v3.s}[2], [x0], #16	\n"
					"	ldr		q11, [x5], #16				\n"
					"	st4		{v4.s, v5.s, v6.s, v7.s}[2], [x0], #16	\n"

					"	prfm	PLDL1KEEP, [x8, #128]	\n"
					"	prfm	PLDL1KEEP, [x9, #128]	\n"

					"	ldr		q12, [x6], #16				\n"
					"	st4		{v0.s, v1.s, v2.s, v3.s}[3], [x0], #16		\n"
					"	ldr		q13, [x7], #16				\n"
					"	st4		{v4.s, v5.s, v6.s, v7.s}[3], [x0], #16		\n"
					"	ldr		q14, [x8], #16				\n"
					"	st4		{v8.s, v9.s, v10.s, v11.s}[0], [x0], #16	\n"
					"	ldr		q15, [x9], #16				\n"
					"	st4		{v12.s, v13.s, v14.s, v15.s}[0], [x0], #16	\n"

					"	subs	x21, x21, #1					\n"

					"	st4		{v8.s, v9.s, v10.s, v11.s}[1], [x0], #16	\n"
					"	st4		{v12.s, v13.s, v14.s, v15.s}[1], [x0], #16	\n"
					"	st4		{v8.s, v9.s, v10.s, v11.s}[2], [x0], #16	\n"
					"	st4		{v12.s, v13.s, v14.s, v15.s}[2], [x0], #16	\n"
					"	st4		{v8.s, v9.s, v10.s, v11.s}[3], [x0], #16	\n"
					"	st4		{v12.s, v13.s, v14.s, v15.s}[3], [x0], #16	\n"

					"	bgt		Sin_PACKA 				\n"

					"	ands	x22, x1, #7				\n"
					"	beq		Sin_PACKA_END			\n"

					"	cmp		x22, #4						\n"
					"	blt 	Sin_K1_PACKA			\n"

					"Sin_K4_PACKA:						\n"

					"	ldr		q0, [x2], #16			\n"
					"	ldr		q1, [x3], #16			\n"
					"	ldr		q2, [x4], #16			\n"
					"	ldr		q3, [x5], #16			\n"

					"	st4		{v0.s, v1.s, v2.s, v3.s}[0], [x0], #16	\n"
					"	ldr		q4, [x6], #16			\n"
					"	st4		{v0.s, v1.s, v2.s, v3.s}[1], [x0], #16	\n"
					"	ldr		q5, [x7], #16			\n"
					"	st4		{v0.s, v1.s, v2.s, v3.s}[2], [x0], #16	\n"
					"	ldr		q6, [x8], #16			\n"
					"	st4		{v0.s, v1.s, v2.s, v3.s}[3], [x0], #16	\n"
					"	ldr		q7, [x9], #16			\n"

					"	st4		{v4.s, v5.s, v6.s, v7.s}[0], [x0], #16	\n"
					"	st4		{v4.s, v5.s, v6.s, v7.s}[1], [x0], #16	\n"
					"	st4		{v4.s, v5.s, v6.s, v7.s}[2], [x0], #16	\n"
					"	st4		{v4.s, v5.s, v6.s, v7.s}[3], [x0], #16	\n"

					"	subs	x22, x22, #4			\n"
					"	beq		Sin_PACKA_END			\n"

					"Sin_K1_PACKA:						\n"

					"	ldr		s0, [x2], #4			\n"
					"	ldr		s1, [x3], #4			\n"
					"	ldr		s2, [x4], #4			\n"
					"	ldr		s3, [x5], #4			\n"

					"	subs 	x22, x22, #1			\n"

					"	st4		{v0.s, v1.s, v2.s, v3.s}[0], [x0], #16	\n"
					"	ldr		s4, [x6], #4			\n"
					"	ldr		s5, [x7], #4			\n"
					"	ldr		s6, [x8], #4			\n"
					"	ldr		s7, [x9], #4			\n"

					"	st4		{v4.s, v5.s, v6.s, v7.s}[0], [x0], #16	\n"

					"	bgt		Sin_K1_PACKA			\n"

					"Sin_PACKA_END:						\n"

					:
	        :
	      		[temp] "m" (temp),
	      		[Ac] "m" (Ac),
	      		[Kc] "m" (Kc),
	      		[LK] "m" (LK)
	        :"x0", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8",
	        "x9", "x10", "x11", "x12", "x13","x14", "x15", "x16",
	        "x17", "x18", "x19", "x20", "x21", "x22", "x23", "x24","x25",
	        "x26", "x27", "x28", "x29", "x30",
	        "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7",
	        "v8", "v9", "v10", "v11", "v12", "v13", "v14", "v15",
	        "v16", "v17", "v18", "v19", "v20", "v21", "v22", "v23",
	        "v24", "v25", "v26", "v27", "v28", "v29", "v30", "v31"

			);

			Ac = Ac + Kc * 8;

		}
	}
}

max_stars_repo_name: AnonymousYWL/MYLIB

max_issues_repo_name: ProgrammerAnonymousWLY/MYLIB

max_forks_repo_name: ProgrammerAnonymousWLY/MYLIB


text:
/********************************************************
  Stanford Driving Software
  Copyright (c) 2011 Stanford University
  All rights reserved.

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

* Redistributions of source code must retain the above
  copyright notice, this list of conditions and the
  following disclaimer.
* Redistributions in binary form must reproduce the above
  copyright notice, this list of conditions and the
  following disclaimer in the documentation and/or other
  materials provided with the distribution.
* The names of the contributors may not be used to endorse
  or promote products derived from this software
  without specific prior written permission.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
  WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
  PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
  DAMAGE.
 ********************************************************/

#ifndef KALMAN_FILTER_H_
#define KALMAN_FILTER_H_

#include <vector>
#include <stdint.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit.h>
#include "gsl.h"

namespace vlr {

class KalmanFilter {

public:
  KalmanFilter();
  KalmanFilter(int state_dim);
  KalmanFilter(const KalmanFilter& kf);
  virtual ~KalmanFilter();

  void transitionUpdate();
  void observationUpdate(gsl_vector* innovation_vector);

  virtual gsl_vector* predictTransition(gsl_vector* orig);
  virtual void predictTransitions(std::vector<gsl_vector*>& orig_pts, std::vector<gsl_vector*>& next_pts);
  virtual gsl_vector* predictObservation(gsl_vector* pt);
  virtual void predictObservations(std::vector<gsl_vector*>& pts, std::vector<gsl_vector*>& obs);

protected:
  void calculateKalmanGain();

public:
  int getStateDim() const { return state_dim_; }
  void setStateDim(int state_dim);

  int getObservationDim() const { return observation_mat_->size1; }

  gsl_vector* getMean() const { return mean_; }
  gsl_matrix* getCovariance() const { return covariance_; }
  gsl_matrix* getKalmanGain() const { return kalman_gain_; }

  gsl_matrix* getTransitionMatrix() const { return transition_mat_; }
  void setTransitionMatrix(gsl_matrix* transition_mat) { transition_mat_ = transition_mat; }

  gsl_matrix* getTransitionCovariance() const { return transition_cov_; }
  void setTransitionCovariance(gsl_matrix* transition_cov) { transition_cov_ = transition_cov; }

  gsl_matrix* getObservationMatrix() const { return observation_mat_; }
  void setObservationMatrix(gsl_matrix* observation_mat) { observation_mat_ = observation_mat; }

  gsl_matrix* getObservationCovariance() const { return observation_cov_; }
  void setObservationCovariance(gsl_matrix* observation_cov) { observation_cov_ = observation_cov; }

  double getYaw() const { return yaw_; }
  void setYaw(double yaw) { yaw_ = yaw; }

  void print() const;

 protected:
  int state_dim_;

  gsl_vector* mean_;
  gsl_matrix* covariance_;

  // temp vars
  gsl_matrix* kalman_gain_;
  gsl_matrix* observation_mat_;
  gsl_matrix* observation_cov_;

  /*
   * TRANSITION MATRIX
   *
   * (for linear KF
   *  transition_mat_ is A)
   *
   *  A = [ 1 0 dt 0 ]
   *      [ 0 1 0 dt ]
   *      [ 0 0 1  0 ]
   *      [ 0 0 0  1 ]
   */
  gsl_matrix* transition_mat_;
  // do we own the pointer?
  uint32_t my_trans_mat_;
  /*
   * TRANSITION COVARIANCE
   *
   *  Q = [ 0 0  0  0  ]
   *      [ 0 0  0  0  ]
   *      [ 0 0 q^2 0  ]
   *      [ 0 0  0 q^2 ]
   */
  gsl_matrix* transition_cov_;
  // do we own the pointer?
  uint32_t my_trans_cov_;

  double yaw_;
};

} // namespace vlr

#endif /*KALMAN_FILTER_H_*/


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:04 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -notwiddleinv 5 */

/*
 * This function contains 32 FP additions, 12 FP multiplications,
 * (or, 26 additions, 6 multiplications, 6 fused multiply/add),
 * 16 stack variables, and 20 memory accesses
 */
static const fftw_real K250000000 = FFTW_KONST(+0.250000000000000000000000000000000000000000000);
static const fftw_real K951056516 = FFTW_KONST(+0.951056516295153572116439333379382143405698634);
static const fftw_real K587785252 = FFTW_KONST(+0.587785252292473129168705954639072768597652438);
static const fftw_real K559016994 = FFTW_KONST(+0.559016994374947424102293417182819058860154590);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftwi_no_twiddle_5(const fftw_complex *input, fftw_complex *output, int istride, int ostride)
{
     fftw_real tmp1;
     fftw_real tmp27;
     fftw_real tmp8;
     fftw_real tmp10;
     fftw_real tmp21;
     fftw_real tmp22;
     fftw_real tmp14;
     fftw_real tmp28;
     fftw_real tmp26;
     fftw_real tmp17;
     ASSERT_ALIGNED_DOUBLE;
     tmp1 = c_re(input[0]);
     tmp27 = c_im(input[0]);
     {
	  fftw_real tmp2;
	  fftw_real tmp3;
	  fftw_real tmp4;
	  fftw_real tmp5;
	  fftw_real tmp6;
	  fftw_real tmp7;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp2 = c_re(input[istride]);
	  tmp3 = c_re(input[4 * istride]);
	  tmp4 = tmp2 + tmp3;
	  tmp5 = c_re(input[2 * istride]);
	  tmp6 = c_re(input[3 * istride]);
	  tmp7 = tmp5 + tmp6;
	  tmp8 = tmp4 + tmp7;
	  tmp10 = K559016994 * (tmp4 - tmp7);
	  tmp21 = tmp2 - tmp3;
	  tmp22 = tmp5 - tmp6;
     }
     {
	  fftw_real tmp12;
	  fftw_real tmp13;
	  fftw_real tmp24;
	  fftw_real tmp15;
	  fftw_real tmp16;
	  fftw_real tmp25;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp12 = c_im(input[istride]);
	  tmp13 = c_im(input[4 * istride]);
	  tmp24 = tmp12 + tmp13;
	  tmp15 = c_im(input[2 * istride]);
	  tmp16 = c_im(input[3 * istride]);
	  tmp25 = tmp15 + tmp16;
	  tmp14 = tmp12 - tmp13;
	  tmp28 = tmp24 + tmp25;
	  tmp26 = K559016994 * (tmp24 - tmp25);
	  tmp17 = tmp15 - tmp16;
     }
     c_re(output[0]) = tmp1 + tmp8;
     {
	  fftw_real tmp18;
	  fftw_real tmp20;
	  fftw_real tmp11;
	  fftw_real tmp19;
	  fftw_real tmp9;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp18 = (K587785252 * tmp14) - (K951056516 * tmp17);
	  tmp20 = (K951056516 * tmp14) + (K587785252 * tmp17);
	  tmp9 = tmp1 - (K250000000 * tmp8);
	  tmp11 = tmp9 - tmp10;
	  tmp19 = tmp10 + tmp9;
	  c_re(output[2 * ostride]) = tmp11 - tmp18;
	  c_re(output[3 * ostride]) = tmp11 + tmp18;
	  c_re(output[ostride]) = tmp19 - tmp20;
	  c_re(output[4 * ostride]) = tmp19 + tmp20;
     }
     c_im(output[0]) = tmp27 + tmp28;
     {
	  fftw_real tmp23;
	  fftw_real tmp31;
	  fftw_real tmp30;
	  fftw_real tmp32;
	  fftw_real tmp29;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp23 = (K951056516 * tmp21) + (K587785252 * tmp22);
	  tmp31 = (K587785252 * tmp21) - (K951056516 * tmp22);
	  tmp29 = tmp27 - (K250000000 * tmp28);
	  tmp30 = tmp26 + tmp29;
	  tmp32 = tmp29 - tmp26;
	  c_im(output[ostride]) = tmp23 + tmp30;
	  c_im(output[4 * ostride]) = tmp30 - tmp23;
	  c_im(output[2 * ostride]) = tmp31 + tmp32;
	  c_im(output[3 * ostride]) = tmp32 - tmp31;
     }
}

fftw_codelet_desc fftwi_no_twiddle_5_desc =
{
     "fftwi_no_twiddle_5",
     (void (*)()) fftwi_no_twiddle_5,
     5,
     FFTW_BACKWARD,
     FFTW_NOTW,
     122,
     0,
     (const int *) 0,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <numpy/arrayobject.h>

#include <math.h>

#ifdef __INTEL_COMPILER
#include <mkl_cblas.h>
#else
#include <cblas.h>
#define kmp_set_blocktime(k)
#endif

static inline int min (int a, int b) { return a < b ? a : b; }
static inline int max (int a, int b) { return a < b ? b : a; }

static inline int square(int x) { return x*x; }

static void max_filter_1d(const float *vals, float *out_vals, int32_t *I,
                          int s, int step, int n, float a, float b) {
  int i;
  for (i = 0; i < n; i++) {
    float max_val = -INFINITY;
    int argmax     = 0;
    int first      = max(0, i-s);
    int last       = min(n-1, i+s);
    int j;
    for (j = first; j <= last; j++) {
      float val = *(vals + j*step) - a*square(i-j) - b*(i-j);
      if (val > max_val) {
        max_val = val;
        argmax  = j;
      }
    }
    *(out_vals + i*step) = max_val;
    *(I + i*step) = argmax;
  }
}

PyObject * deformation_cost (PyArrayObject * pydata, float ax, float bx, float ay, float by, int s) {
  npy_intp * dims = PyArray_DIMS(pydata);
  npy_intp * stride = PyArray_STRIDES(pydata);

  if (PyArray_NDIM(pydata) != 2) {
    PyErr_SetString(PyExc_TypeError, "data must be 2 dimensional.");
    return NULL;
  }

  if (PyArray_DESCR(pydata)->type_num != NPY_FLOAT) {
    PyErr_SetString(PyExc_TypeError, "data must be single precision floating point.");
    return NULL;
  }

  if (stride[0] != dims[1]*sizeof(float)) {
    PyErr_SetString(PyExc_TypeError, "Stride[0] must be sizeof(float).");
    return NULL;
  }

  if (stride[1] != sizeof(float)) {
    PyErr_SetString(PyExc_TypeError, "Stride[1] must be Dims[0]*sizeof(float).");
    return NULL;
  }

  PyArrayObject * pydeformed = (PyArrayObject*)PyArray_SimpleNew((npy_intp)2, dims, NPY_FLOAT);
  PyArrayObject * pyIx = (PyArrayObject*)PyArray_SimpleNew((npy_intp)2, dims, NPY_INT32);
  PyArrayObject * pyIy = (PyArrayObject*)PyArray_SimpleNew((npy_intp)2, dims, NPY_INT32);

  float *tmpM = (float *)calloc(dims[0]*dims[1], sizeof(float));
  int32_t *tmpIx = (int32_t*)calloc(dims[0]*dims[1], sizeof(int32_t));
  int32_t *tmpIy = (int32_t*)calloc(dims[0]*dims[1], sizeof(int32_t));

  int x, y;

  for (y = 0; y < dims[0]; y++)
    max_filter_1d((float*)PyArray_GETPTR2(pydata, y, 0), tmpM+y*dims[1], tmpIx+y*dims[1], s, 1, dims[1], ax, bx);

  for (x = 0; x < dims[1]; x++)
    max_filter_1d(tmpM+x, (float*)PyArray_GETPTR2(pydeformed, 0, x), tmpIy+x, s, dims[1], dims[0], ay, by);

  for (x = 0; x < dims[1]; ++x) {
    for (y = 0; y < dims[0]; ++y) {
      *(int32_t*)PyArray_GETPTR2(pyIy, y, x) = tmpIy[y*dims[1]+x];
      *(int32_t*)PyArray_GETPTR2(pyIx, y, x) = tmpIx[tmpIy[y*dims[1]+x]*dims[1]+x];
    }
  }

  free(tmpM);
  free(tmpIx);
  free(tmpIy);

  return Py_BuildValue("NNN", pydeformed, pyIx, pyIy);
}

PyObject * filter_image (PyArrayObject * pyfeatures, PyArrayObject * pyfilter, float bias, int width, int height) {
    npy_intp * features_dims = PyArray_DIMS(pyfeatures);
    npy_intp * filter_dims = PyArray_DIMS(pyfilter);
    int a, b, l;
    PyArrayObject * pyfiltered = NULL;
    npy_intp * features_stride = PyArray_STRIDES(pyfeatures);
    npy_intp * filtered_stride = NULL;
    npy_intp filtered_dims[2] = {0, 0};
    int tight_width;
    int tight_height;

    if (PyArray_NDIM(pyfeatures) != 3) {
        PyErr_SetString(PyExc_TypeError, "Features must be 3 dimensional.");
        return NULL;
    }

    if (PyArray_NDIM(pyfilter) != 3) {
        PyErr_SetString(PyExc_TypeError, "Filter must be 3 dimensional.");
        return NULL;
    }

    if (PyArray_DESCR(pyfeatures)->type_num != NPY_FLOAT) {
        PyErr_SetString(PyExc_TypeError, "Features must be single precision floating point.");
        return NULL;
    }

    if (PyArray_DESCR(pyfilter)->type_num != NPY_FLOAT) {
        PyErr_SetString(PyExc_TypeError, "Filter must be a single precision floating point.");
        return NULL;
    }

    if (features_dims[2] != 32) {
        PyErr_SetString(PyExc_TypeError, "features' feature dimsionality should be 32.");
        return NULL;
    }

    if (filter_dims[2] != 32) {
        PyErr_SetString(PyExc_TypeError, "filters' feature dimensionality should be 32.");
        return NULL;
    }

    tight_height = features_dims[0]-filter_dims[0]+1;
    tight_width = features_dims[1]-filter_dims[1]+1;

    filtered_dims[0] = height ? height : tight_height;
    filtered_dims[1] = width ? width : tight_width;

    if (filtered_dims[0] < 1 || filtered_dims[1] < 1) {
        PyErr_SetString(PyExc_TypeError, "Input features are too small for filter.");
        return NULL;
    }

    #pragma omp critical
    pyfiltered = (PyArrayObject*)PyArray_SimpleNew((npy_intp)2, filtered_dims, NPY_FLOAT);

    filtered_stride = PyArray_STRIDES(pyfiltered);

    /* zero out array */
    for (a = 0; a < tight_height; ++a) {
        for (b = 0; b < tight_width; ++b) {
            *(float*)PyArray_GETPTR2(pyfiltered, a, b) = -bias;
        }
    }

    /* iterate over filter which should be tiny compared to the image */
    int i;
    int stride_src = features_stride[1]/sizeof(float);
    int stride_dst = filtered_stride[1]/sizeof(float);
    for (i = 0; i < filter_dims[0]; ++i) {
        int j;
        for (j = 0; j < filter_dims[1]; ++j) {
            int k;
            for (k = 0; k < tight_height; ++k) {
                float * out = (float*)PyArray_GETPTR2(pyfiltered, k, 0);
                /* for each layer */
                for (l = 0; l < 32; ++l) {
                    float weight = *(float*)PyArray_GETPTR3(pyfilter, i, j, l);
                    float * in = (float*)PyArray_GETPTR3(pyfeatures, i+k, j, l);
                    cblas_saxpy(tight_width, weight, in, stride_src, out, stride_dst);
                }
            }
        }
    }

    for (a = tight_height; a < filtered_dims[0]; ++a) {
        for (b = 0; b < filtered_dims[1]; ++b) {
            *(float*)PyArray_GETPTR2(pyfiltered, a, b) = -INFINITY;
        }
    }

    for (a = 0; a < tight_height; ++a) {
        for (b = tight_width; b < filtered_dims[1]; ++b) {
            *(float*)PyArray_GETPTR2(pyfiltered, a, b) = -INFINITY;
        }
    }

    return Py_BuildValue("N", pyfiltered);
}

static PyObject * DeformationCost(PyObject * self, PyObject * args)
{
    PyArrayObject * pydata;
    float ax = 0.0f, bx = 0.0f, ay = 0.0f, by = 0.0f;
    int s = 0;
    if (!PyArg_ParseTuple(args, "O!ffffi", &PyArray_Type, &pydata, &ax, &bx, &ay, &by, &s))
        return NULL;
    return deformation_cost(pydata, ax, bx, ay, by, s);
}

static PyObject * FilterImage(PyObject * self, PyObject * args)
{
    PyArrayObject * pyfeatures;
    PyArrayObject * pyfilter;
    float bias = 0.0f;
    int width = 0;
    int height = 0;
    if (!PyArg_ParseTuple(args, "O!O!|fii", &PyArray_Type, &pyfeatures, &PyArray_Type, &pyfilter, &bias, &width, &height))
        return NULL;
    return filter_image(pyfeatures, pyfilter, bias, width, height);
}

static PyObject * FilterImages(PyObject * self, PyObject * args)
{
    PyObject * pyfeatures_list;
    PyObject * pydims_list = NULL;
    PyArrayObject * pyfilter;
    float bias = 0.0f;
    int numfilters;
    int numdims = 0;
    int i;
    PyObject ** objs = NULL;
    PyObject ** results = NULL;
    PyObject * pyresults_list;
    if (!PyArg_ParseTuple(args, "O!O!|fO!", &PyList_Type, &pyfeatures_list, &PyArray_Type, &pyfilter, &bias, &PyList_Type, &pydims_list))
        return NULL;

    numfilters = PyList_Size(pyfeatures_list);

    if (pydims_list) {
        numdims = PyList_Size(pydims_list);
    }

    if (numdims && numdims != numfilters) {
        PyErr_SetString(PyExc_TypeError, "If pad dims are specified, then it must be the same length as the features list.");
        return NULL;
    }

    objs = (PyObject**)calloc(numfilters, sizeof(PyObject*));
    int* widths = (int*)calloc(numfilters, sizeof(int));
    int* heights = (int*)calloc(numfilters, sizeof(int));
    results = (PyObject**)calloc(numfilters, sizeof(PyObject*));

    for (i = 0; i < numfilters; ++i) {
        objs[i] = PyList_GetItem(pyfeatures_list, i);
        if (!PyArray_Check(objs[i])) {
            free(objs);
            free(widths);
            free(heights);
            free(results);
            PyErr_SetString(PyExc_TypeError, "Must contain a list of numpy arrays.");
            return NULL;
        }

        if (pydims_list) {
            PyObject * dims = PyList_GetItem(pydims_list, i);
            if (!PyTuple_Check(dims) || 2 != PyTuple_Size(dims)) {
                free(objs);
                free(widths);
                free(heights);
                free(results);
                PyErr_SetString(PyExc_TypeError, "Must contain a list of tuples.");
                return NULL;
            }

            heights[i] = PyInt_AsLong(PyTuple_GetItem(dims, 0));
            widths[i] = PyInt_AsLong(PyTuple_GetItem(dims, 1));
        } else {
            widths[i] = 0;
            heights[i] = 0;
        }
    }

    kmp_set_blocktime(0);
    #pragma omp parallel for schedule(dynamic)
    for (i = 0; i < numfilters; ++i) {
        PyArrayObject * pyfeatures = (PyArrayObject*)objs[i];
        results[i] = filter_image(pyfeatures, pyfilter, bias, widths[i], heights[i]);
    }

    free(objs);
    free(widths);
    free(heights);

    pyresults_list = PyList_New(numfilters);

    for (i = 0; i < numfilters; ++i) {
        PyList_SetItem(pyresults_list, i, results[i]);
    }

    free(results);

    return Py_BuildValue("N", pyresults_list);
}

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef moduledef = {
    PyModuleDef_HEAD_INIT,
    "_detection",
    "Native convolution detection routine.",
    -1,
    NULL,
    NULL,
    NULL,
    NULL,
    NULL
};
#endif

#if PY_MAJOR_VERSION < 3
static PyMethodDef _detection_methods[] = {
    {"FilterImage", FilterImage, METH_VARARGS, "Compute a 2D cross correlation between a filter and image features.  Optionally add bias term."},
    {"FilterImages", FilterImages, METH_VARARGS, "Compute a 2D cross correlation between a filter and several image features in parallel.  Optionally add bias term."},
    {"DeformationCost", DeformationCost, METH_VARARGS, "Compute a fast bounded distance transform for the deformation cost."},
    {NULL}
};
#endif

#if PY_MAJOR_VERSION >= 3
PyMODINIT_FUNC
PyInit__detection(void)
#else
PyMODINIT_FUNC
init_detection(void)
#endif
{
    import_array();

#if PY_MAJOR_VERSION >= 3
    PyObject *m = PyModule_Create(&moduledef);
#else
    Py_InitModule3("_detection", _detection_methods, "Native convolution detection routine.");
#endif

#if PY_MAJOR_VERSION >= 3
    return m;
#endif
}

max_stars_repo_name: caomw/pydro

max_issues_repo_name: caomw/pydro

max_forks_repo_name: kmatzen/pydro


text:
/**
 *
 * @file testing_csymm.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:45:18 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_cmain.h"

static int check_solution(PLASMA_enum transA, PLASMA_enum transB, int M, int N,
                          PLASMA_Complex32_t alpha, PLASMA_Complex32_t *A, int LDA,
                          PLASMA_Complex32_t *B, int LDB,
                          PLASMA_Complex32_t beta, PLASMA_Complex32_t *Cref, PLASMA_Complex32_t *Cplasma, int LDC);

int testing_csymm(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 7 ){
        USAGE("SYMM", "alpha beta M N K LDA LDB LDC",
              "   - alpha : alpha coefficient \n"
              "   - beta : beta coefficient \n"
              "   - M : number of rows of matrices A and C \n"
              "   - N : number of columns of matrices B and C \n"
              "   - LDA : leading dimension of matrix A \n"
              "   - LDB : leading dimension of matrix B \n"
              "   - LDC : leading dimension of matrix C\n");
        return -1;
    }

    PLASMA_Complex32_t alpha = (PLASMA_Complex32_t) atol(argv[0]);
    PLASMA_Complex32_t beta  = (PLASMA_Complex32_t) atol(argv[1]);
    int M     = atoi(argv[2]);
    int N     = atoi(argv[3]);
    int LDA   = atoi(argv[4]);
    int LDB   = atoi(argv[5]);
    int LDC   = atoi(argv[6]);
    int MNmax = max(M, N);

    float eps;
    int info_solution;
    int i, j, s, u;
    int LDAxM = LDA*MNmax;
    int LDBxN = LDB*N;
    int LDCxN = LDC*N;

    PLASMA_Complex32_t *A      = (PLASMA_Complex32_t *)malloc(LDAxM*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *B      = (PLASMA_Complex32_t *)malloc(LDBxN*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *C      = (PLASMA_Complex32_t *)malloc(LDCxN*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *Cinit  = (PLASMA_Complex32_t *)malloc(LDCxN*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *Cfinal = (PLASMA_Complex32_t *)malloc(LDCxN*sizeof(PLASMA_Complex32_t));

    /* Check if unable to allocate memory */
    if ((!A)||(!B)||(!Cinit)||(!Cfinal)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_slamch_work('e');

    printf("\n");
    printf("------ TESTS FOR PLASMA CSYMM ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", M, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 10.\n");

    /*----------------------------------------------------------
    *  TESTING CSYMM
    */

    /* Initialize A */
    PLASMA_cplgsy( (float)0., MNmax, A, LDA, 51 );

    /* Initialize B */
    LAPACKE_clarnv_work(IONE, ISEED, LDBxN, B);

    /* Initialize C */
    LAPACKE_clarnv_work(IONE, ISEED, LDCxN, C);

    for (s=0; s<2; s++) {
        for (u=0; u<2; u++) {

            /* Initialize  Cinit / Cfinal */
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cinit[LDC*j+i] = C[LDC*j+i];
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cfinal[LDC*j+i] = C[LDC*j+i];

            /* PLASMA CSYMM */
            PLASMA_csymm(side[s], uplo[u], M, N, alpha, A, LDA, B, LDB, beta, Cfinal, LDC);

            /* Check the solution */
            info_solution = check_solution(side[s], uplo[u], M, N, alpha, A, LDA, B, LDB, beta, Cinit, Cfinal, LDC);

            if (info_solution == 0) {
                printf("***************************************************\n");
                printf(" ---- TESTING CSYMM (%5s, %5s) ....... PASSED !\n", sidestr[s], uplostr[u]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING CSYMM (%s, %s) ... FAILED !\n", sidestr[s], uplostr[u]);
                printf("************************************************\n");
            }
        }
    }

    free(A); free(B); free(C);
    free(Cinit); free(Cfinal);

    return 0;
}

/*--------------------------------------------------------------
 * Check the solution
 */

static int check_solution(PLASMA_enum side, PLASMA_enum uplo, int M, int N,
                   PLASMA_Complex32_t alpha, PLASMA_Complex32_t *A, int LDA,
                   PLASMA_Complex32_t *B, int LDB,
                   PLASMA_Complex32_t beta, PLASMA_Complex32_t *Cref, PLASMA_Complex32_t *Cplasma, int LDC)
{
    int info_solution, NrowA;
    float Anorm, Bnorm, Cinitnorm, Cplasmanorm, Clapacknorm, Rnorm;
    float eps;
    PLASMA_Complex32_t beta_const;
    float result;
    float *work = (float *)malloc(max(M, N)* sizeof(float));

    beta_const  = (PLASMA_Complex32_t)-1.0;

    NrowA = (side == PlasmaLeft) ? M : N;
    Anorm       = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), NrowA, NrowA, A,       LDA, work);
    Bnorm       = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,     N,     B,       LDB, work);
    Cinitnorm   = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,     N,     Cref,    LDC, work);
    Cplasmanorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,     N,     Cplasma, LDC, work);

    cblas_csymm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo, M, N, CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Cref, LDC);

    Clapacknorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Cref, LDC, work);

    cblas_caxpy(LDC * N, CBLAS_SADDR(beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Cref, LDC, work);

    eps = LAPACKE_slamch_work('e');

    printf("Rnorm %e, Anorm %e, Bnorm %e, Cinitnorm %e, Cplasmanorm %e, Clapacknorm %e\n",Rnorm,Anorm,Bnorm,Cinitnorm,Cplasmanorm,Clapacknorm);

    result = Rnorm / ((Anorm + Bnorm + Cinitnorm) * N * eps);

    printf("============\n");
    printf("Checking the norm of the difference against reference CSYMM \n");
    printf("-- ||Cplasma - Clapack||_oo/((||A||_oo+||B||_oo+||C||_oo).N.eps) = %e \n", result );

    if ( isinf(Clapacknorm) || isinf(Cplasmanorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
    }
    else {
        printf("-- The solution is CORRECT ! \n");
        info_solution= 0 ;
    }
    free(work);
    return info_solution;
}


text:
#ifndef QDM_IJK_H
#define QDM_IJK_H 1

#include <stddef.h>

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

typedef struct {
  size_t size1;
  size_t size2;
  size_t size3;

  double *data;
  bool owner;
} qdm_ijk;

typedef struct {
  qdm_ijk ijk;
} qdm_ijk_view;

qdm_ijk *
qdm_ijk_alloc(
    size_t size1,
    size_t size2,
    size_t size3
);

qdm_ijk *
qdm_ijk_calloc(
    size_t size1,
    size_t size2,
    size_t size3
);

void
qdm_ijk_free(qdm_ijk *t);

qdm_ijk_view
qdm_ijk_view_array(
    double *base,
    size_t size1,
    size_t size2,
    size_t size3
);

gsl_vector_view
qdm_ijk_get_k(
    qdm_ijk *t,
    size_t i,
    size_t j
);

gsl_matrix_view
qdm_ijk_get_ij(
    qdm_ijk *t,
    size_t k
);

double
qdm_ijk_get(
    qdm_ijk *t,
    size_t i,
    size_t j,
    size_t k
);

gsl_matrix *
qdm_ijk_cov(
    qdm_ijk *t
);

int
qdm_ijk_write(
    hid_t id,
    const char *name,
    const qdm_ijk *t
);

int
qdm_ijk_read(
    hid_t id,
    const char *name,
    qdm_ijk **t
);

#endif /* QDM_IJK_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_drot (const int N, double *X, const int incX, double *Y, const int incY,
	    const double c, const double s)
{
#define BASE double
#include "source_rot.h"
#undef BASE
}


text:
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include "cimple_polytope_library.h"

/**
 * "Constructor" Dynamically allocates the space a polytope needs
 */
struct polytope *polytope_alloc(size_t k,
                                size_t n)
{

    struct polytope *return_polytope = malloc (sizeof (struct polytope));

    return_polytope->H = gsl_matrix_alloc(k, n);
    if (return_polytope->H == NULL) {
        free (return_polytope);
        return NULL;
    }

    return_polytope->G = gsl_vector_alloc(k);
    if (return_polytope->H == NULL) {
        free (return_polytope);
        return NULL;
    }

    return_polytope->chebyshev_center = malloc(n* sizeof (double));
    if (return_polytope->chebyshev_center == NULL) {
        free (return_polytope);
        return NULL;
    }

    return return_polytope;
};

/**
 * "Destructor" Deallocates the dynamically allocated memory of the polytope
 */
void polytope_free(polytope *polytope)
{
    gsl_matrix_free(polytope->H);
    gsl_vector_free(polytope->G);
    free(polytope->chebyshev_center);
    free(polytope);
};

/**
 * "Constructor" Dynamically allocates the space a polytope needs
 */
struct cell *cell_alloc(size_t k,
                        size_t n,
                        int time_horizon)
{

    struct cell *return_cell = malloc (sizeof (struct cell));

    return_cell->safe_mode = malloc(sizeof(polytope)*time_horizon);
    if (return_cell->safe_mode == NULL) {
        free (return_cell);
        return NULL;
    }

    return_cell->polytope_description = polytope_alloc(k,n);
    if (return_cell->polytope_description == NULL) {
        free (return_cell);
        return NULL;
    }
    return return_cell;
};

/**
 * "Destructor" Deallocates the dynamically allocated memory of the region of polytopes
 */
void cell_free(cell *cell)
{
    polytope_free(cell->polytope_description);
    free(cell->safe_mode);
    free(cell);

};

/**
 * "Constructor" Dynamically allocates the space a region of polytope needs
 */
struct abstract_state *abstract_state_alloc(size_t *k,
                                            size_t k_hull,
                                            size_t n,
                                            int transitions_in_count,
                                            int transitions_out_count,
                                            int cells_count,
                                            int time_horizon)
{

    struct abstract_state *return_abstract_state = malloc (sizeof (struct abstract_state));

    /*
     * Default values: at initialization safe mode is not yet computed.
     * Thus it is assumed that the state does not contain an invariant set => invariant_set = NULL
     * next_state (next state the system has to transition to reach an invariant set) is unknown at initialization.
     */

    return_abstract_state->next_state = NULL;
    return_abstract_state->invariant_set = NULL;
//    return_abstract_state->distance_invariant_set = INFINITY;

    return_abstract_state->cells = malloc(sizeof(cell)*cells_count);
    if (return_abstract_state->cells == NULL) {
        free (return_abstract_state);
        return NULL;
    }

    for(int i = 0; i < cells_count; i++){
        return_abstract_state->cells[i] = cell_alloc(*(k+i), n, time_horizon);
        if (return_abstract_state->cells[i] == NULL) {
            free (return_abstract_state);
            return NULL;
        }
    }

    return_abstract_state->cells_count = cells_count;

    return_abstract_state->transitions_in = malloc(sizeof(struct abstract_state) * transitions_in_count);
    if (return_abstract_state->transitions_in == NULL) {
        free (return_abstract_state);
        return NULL;
    }

    return_abstract_state->transitions_in_count = transitions_in_count;

    return_abstract_state->transitions_out = malloc(sizeof(struct abstract_state) * transitions_out_count);
    if (return_abstract_state->transitions_out == NULL) {
        free (return_abstract_state);
        return NULL;
    }

    return_abstract_state->transitions_out_count = transitions_out_count;

    return_abstract_state->convex_hull = polytope_alloc(k_hull, n);
    if (return_abstract_state->convex_hull == NULL) {
        free (return_abstract_state);
        return NULL;
    }

    return return_abstract_state;
};

/**
 * "Destructor" Deallocates the dynamically allocated memory of the region of polytopes
 */
void abstract_state_free(abstract_state * abstract_state){
    polytope_free(abstract_state->convex_hull);
    for(int i = 0; i< abstract_state->cells_count; i++){
        cell_free(abstract_state->cells[i]);
    }
    free(abstract_state->transitions_out);
    free(abstract_state->transitions_in);
    free(abstract_state->cells);
    free(abstract_state);

};

/**
 * Converts two C arrays to a polytope consistent of a left side matrix (i.e. H) and right side vector (i.e. G)
 */
void polytope_from_arrays(polytope *polytope,
                          double *left_side,
                          double *right_side,
                          double *cheby,
                          char*name)
{

    gsl_matrix_from_array(polytope->H, left_side, name);
    gsl_vector_from_array(polytope->G, right_side, name);
    for(int i = 0; i < polytope->H->size2; i++){
        polytope->chebyshev_center[i] = cheby[i];
    }
};

/**
 * Converts a polytope in gsl form to cdd constraint form
 */
dd_PolyhedraPtr polytope_to_cdd(polytope *original,
                                dd_ErrorType *err)
{
    dd_PolyhedraPtr new;
    dd_MatrixPtr constraints;
    constraints = dd_CreateMatrix(original->H->size1, (original->H->size2+1));
    for (size_t k = 0; k < (original->H->size1); k++) {
        double value = gsl_vector_get(original->G, k);
        dd_set_d(constraints->matrix[k][0],value);
    }
    for(size_t i = 0; i<original->H->size1; i++){
        for (size_t j = 1; j < (original->H->size2+1); j++) {
            double value = gsl_matrix_get(original->H, i, j-1);
            dd_set_d(constraints->matrix[i][j],-1*value);
        }
    }
    constraints->representation=dd_Inequality;
    new = dd_DDMatrix2Poly(constraints, err);
    dd_FreeMatrix(constraints);
    return new;
};

/**
 * Converts a polytope in cdd constraint form to gsl form
 */
polytope * cdd_to_polytope(dd_PolyhedraPtr *original)
{

    dd_MatrixPtr constraints;
    constraints = dd_CopyInequalities(*original);
    polytope *new = polytope_alloc(constraints->rowsize, constraints->colsize-1);
    for (size_t k = 0; k < (constraints->rowsize); k++) {
        double value = dd_get_d(constraints->matrix[k][0]);
        gsl_vector_set(new->G, k, value);
    }
    for(size_t i = 0; i<constraints->rowsize; i++){
        for (size_t j = 0; j < (constraints->colsize-1); j++) {
            double value = dd_get_d(constraints->matrix[i][j+1]);
            gsl_matrix_set(new->H,i,j,-value);
        }
    }
    dd_FreeMatrix(constraints);
    return new;

};

/**
 * Generate a polytope representing a scaled unit cube
 */
polytope * polytope_scaled_unit_cube(double scale,
                                     int dimensions)
{

    dd_PolyhedraPtr cube = cdd_scaled_unit_cube(scale, dimensions);
    polytope * return_cube = cdd_to_polytope(&cube);
    dd_FreePolyhedra(cube);

    return return_cube;
};

/**
 * Checks whether a state is in a certain polytope
 */
bool polytope_check_state(polytope *polytope,
                         gsl_vector *x)
{
    gsl_vector * result = gsl_vector_alloc(polytope->G->size);
    gsl_blas_dgemv(CblasNoTrans, 1.0, polytope->H, x, 0.0, result);
    for(size_t i = 0; i< polytope->G->size; i++){
        if(gsl_vector_get(result, i) > gsl_vector_get(polytope->G, i)){
            gsl_vector_free(result);
            return false;
        }
    }
    gsl_vector_free(result);
    return true;
};

/**
 * Checks whether polytope P1 \ issubset P2
 */
bool polytope_is_subset(polytope *P1,
                        polytope *P2)
{
    dd_ErrorType err = dd_NoError;
    dd_PolyhedraPtr first = polytope_to_cdd(P1, &err);
    dd_MatrixPtr verticesFirst = dd_CopyGenerators(first);
    gsl_vector *vertex = gsl_vector_alloc(P1->H->size2);
    gsl_vector_set_zero(vertex);
    int is_included;
    bool is_subset = true;
    for(int i = 0; i<verticesFirst->rowsize;i++){
        double valueFirst0 = dd_get_d(verticesFirst->matrix[i][0]);
        if(valueFirst0 == 1){
            for(int j = 0; j<vertex->size; j++){
                double valueFirstj = dd_get_d(verticesFirst->matrix[i][j+1]);
                gsl_vector_set(vertex,(size_t)j, valueFirstj);
            }
            is_included = polytope_check_state(P2,vertex);
            if(!is_included){
                is_subset = false;
                break;
            }
        }
    }
    gsl_vector_free(vertex);
    dd_FreePolyhedra(first);
    dd_FreeMatrix(verticesFirst);
    return is_subset;
};

/**
 * Unite inequalities of P1 and P2 in new polytope and remove redundancies
 */
polytope * polytope_unite_inequalities(polytope *P1,
                                       polytope *P2)
{
    polytope *united = polytope_alloc(P1->H->size1+P2->H->size1,P1->H->size2);
    //Unite H
    gsl_matrix_set_zero(united->H);
    gsl_matrix_view H_P1 = gsl_matrix_submatrix(united->H, 0, 0, P1->H->size1, P1->H->size2);
    gsl_matrix_memcpy(&H_P1.matrix, P1->H);
    gsl_matrix_view H_P2 = gsl_matrix_submatrix(united->H, P1->H->size1, 0, P2->H->size1, P2->H->size2);
    gsl_matrix_memcpy(&H_P2.matrix, P2->H);

    //Unite G
    gsl_vector_set_zero(united->G);
    gsl_vector_view G_P1 = gsl_vector_subvector(united->G, 0, P1->G->size);
    gsl_vector_memcpy(&G_P1.vector, P1->G);
    gsl_vector_view G_P2 = gsl_vector_subvector(united->G, P1->G->size,P2->G->size);
    gsl_vector_memcpy(&G_P2.vector, P2->G);
    polytope * return_polytope = polytope_minimize(united);

    polytope_free(united);

    return return_polytope;

};

/**
 * Project original polytope (in cdd format) to the first n dimensions
 */
polytope * polytope_projection(polytope * original,
                               size_t n)
{

    dd_ErrorType err;
    dd_PolyhedraPtr orig_cdd = polytope_to_cdd(original, &err);
    dd_PolyhedraPtr new_cdd = NULL;

    cdd_projection(&orig_cdd, &new_cdd, n, &err);
    polytope * new = cdd_to_polytope(&new_cdd);
    dd_FreePolyhedra(orig_cdd);
    dd_FreePolyhedra(new_cdd);

    return new;

};

polytope * polytope_linear_transform(polytope *original,
                                     gsl_matrix *scale){

    dd_ErrorType err = dd_NoError;
    dd_PolyhedraPtr orig_cdd = polytope_to_cdd(original, &err);

    dd_MatrixPtr vertices = dd_CopyGenerators(orig_cdd);

    dd_MatrixPtr transformed_vertices = dd_CreateMatrix(1, scale->size1+1);
    dd_MatrixPtr transformed_vertex = dd_CreateMatrix(1, scale->size1+1);
    bool new_matrix_started = false;
    gsl_vector *vertex = gsl_vector_alloc((vertices->colsize - 1));
    gsl_vector *scaled_vertex = gsl_vector_alloc(scale->size1);
    gsl_vector_set_zero(vertex);
    gsl_vector_set_zero(scaled_vertex);

    for(int i = 0; i<vertices->rowsize; i++) {
        //Check whether row represents ray or vertex
        double is_vertex = dd_get_d(vertices->matrix[i][0]);
        if (is_vertex == 1) {
            for (size_t j = 1; j < vertices->colsize; j++) {
                double value = dd_get_d(vertices->matrix[i][j]);
                gsl_vector_set(vertex, j-1, value);
            }
            gsl_blas_dgemv(CblasNoTrans, 1.0, scale, vertex, 0.0, scaled_vertex);

            if (new_matrix_started) {
                dd_set_d(transformed_vertex->matrix[0][0], 1);
                for (size_t j = 0; j < scaled_vertex->size; j++) {
                    dd_set_d(transformed_vertex->matrix[0][j + 1], gsl_vector_get(scaled_vertex, j));
                }
                dd_MatrixAppendTo(&transformed_vertices, transformed_vertex);
            } else {
                dd_set_d(transformed_vertices->matrix[0][0], 1);
                for (size_t j = 0; j < scaled_vertex->size; j++) {
                    dd_set_d(transformed_vertices->matrix[0][j + 1], gsl_vector_get(scaled_vertex, j));
                }
                new_matrix_started = true;
            }

        }
    }

    transformed_vertices->representation = dd_Generator;
    dd_PolyhedraPtr transformed_cdd = dd_DDMatrix2Poly(transformed_vertices, &err);
    polytope *transformed = cdd_to_polytope(&transformed_cdd);

    //Clean up
    dd_FreeMatrix(transformed_vertex);
    dd_FreeMatrix(transformed_vertices);
    dd_FreePolyhedra(orig_cdd);
    dd_FreePolyhedra(transformed_cdd);
    dd_FreeMatrix(vertices);
    gsl_vector_free(vertex);
    gsl_vector_free(scaled_vertex);

    return transformed;

};
/**
 * Remove redundancies from gsl polytope inequalities
 */
polytope * polytope_minimize(polytope *original)
{

    dd_ErrorType err = dd_NoError;
    dd_MatrixPtr min_matrix;

    dd_PolyhedraPtr cdd_original = polytope_to_cdd(original, &err);

    dd_PolyhedraPtr cdd_minimized = cdd_minimize(&cdd_original, &err);
    min_matrix = dd_CopyInequalities(cdd_minimized);

    polytope * minimized = cdd_to_polytope(&cdd_minimized);
    dd_FreeMatrix(min_matrix);
    dd_FreePolyhedra(cdd_original);
    dd_FreePolyhedra(cdd_minimized);

    return minimized;

};

/**
 * Compute Minkowski sum of two polytopes
 */
polytope * polytope_minkowski(polytope *P1,
                              polytope *P2)
{
    dd_ErrorType err = dd_NoError;
    dd_PolyhedraPtr A = polytope_to_cdd(P1, &err);
    dd_PolyhedraPtr B = polytope_to_cdd(P2, &err);
    dd_PolyhedraPtr C = cdd_minkowski(A,B);
    polytope * returnPolytope = cdd_to_polytope(&C);
    dd_FreePolyhedra(A);
    dd_FreePolyhedra(B);
    dd_FreePolyhedra(C);
    return returnPolytope;
};

/**
 * Compute Pontryagin difference C = A-B s.t.:
 * A-B = {c \in A-B| c+b \in A, \forall b \in B}
 */
polytope * polytope_pontryagin(polytope* A,
                               polytope* B)
{

    dd_ErrorType err;
    dd_MatrixPtr verticesA, verticesB;
    //create cddPoly A,B
    dd_PolyhedraPtr cddA = polytope_to_cdd(A, &err);

    dd_PolyhedraPtr cddB = polytope_to_cdd(B, &err);

    polytope *C = NULL;

    verticesA = dd_CopyGenerators(cddA);
    verticesB = dd_CopyGenerators(cddB);
    for(int i = 0; i<verticesB->rowsize; i++){
        //Check whether row represents ray or vertex
        double value_B0 = dd_get_d(verticesB->matrix[i][0]);
        if(value_B0 == 1){
            dd_MatrixPtr tempA;
            //A-b (where b is the vertex)
            tempA = dd_CreateMatrix(verticesA->rowsize,verticesA->colsize);
            //each vertex of A displaced by b
            for(int j = 0; j<verticesA->rowsize; j++){
                double value_A0 = dd_get_d(verticesA->matrix[j][0]);
                if(value_A0 == 1){
                    dd_set_d(tempA->matrix[j][0], 1);
                    for(int k = 1; k<verticesA->colsize; k++){
                        double value_Ak = dd_get_d(verticesA->matrix[j][k]);
                        double value_Bk = dd_get_d(verticesB->matrix[i][k]);
                        dd_set_d(tempA->matrix[j][k],(value_Ak-value_Bk));
                    }
                }else{
                    dd_set_d(tempA->matrix[j][0], 0);
                    for(int k = 1; k<verticesA->colsize; k++){
                        double value_Ak = dd_get_d(verticesA->matrix[j][k]);
                        dd_set_d(tempA->matrix[j][k],(value_Ak));
                    }
                }
            }
            dd_PolyhedraPtr cdd_temp;
            tempA->representation = dd_Generator;
            cdd_temp = dd_DDMatrix2Poly(tempA, &err);
            polytope *tempC = cdd_to_polytope(&cdd_temp);
            dd_FreePolyhedra(cdd_temp);
            if(C == NULL){
                C = polytope_alloc(tempC->H->size1,tempC->H->size2);
                gsl_matrix_memcpy(C->H,tempC->H);
                gsl_vector_memcpy(C->G,tempC->G);
                polytope_free(tempC);
            }else{
                polytope *copyC = C;
                C = polytope_unite_inequalities(copyC, tempC);
                polytope_free(tempC);
                polytope_free(copyC);
            }
            dd_FreeMatrix(tempA);
        }
    }

    dd_FreeMatrix(verticesA);
    dd_FreeMatrix(verticesB);
    dd_FreePolyhedra(cddA);
    dd_FreePolyhedra(cddB);
    return C;
};

/**
 * Set up constraints in quadratic problem for GUROBI
 */
int polytope_to_constraints_gurobi(polytope *constraints,
                                   GRBmodel *model,
                                   size_t N)
{

    int error = 0;
    double constraint_val[N];
    int ind[N];
    for(size_t i = 0; i < constraints->H->size1;i++){
        char constraint_name[5];
        sprintf(constraint_name, "c_%d", (int)i);
        for(size_t j = 0; j < N;j++){
            ind[j] = (int)j;
            constraint_val[j] = gsl_matrix_get(constraints->H,i,j);
        }

        error = GRBaddconstr(model, (int)N, ind, constraint_val, GRB_LESS_EQUAL, gsl_vector_get(constraints->G,i), constraint_name);

        if(error){
            return error;
        }
    }
    return error;
};

/**
 * Generate a polytope representing a scaled unit cube
 */
dd_PolyhedraPtr cdd_scaled_unit_cube(double scale,
                                     int dimensions)
{

    dd_MatrixPtr constraints;
    dd_PolyhedraPtr cube = NULL;
    dd_ErrorType err = dd_NoError;
    constraints = dd_CreateMatrix(dimensions*2,dimensions+1);
    for(int i = 0; i<(dimensions); i++){

            dd_set_d(constraints->matrix[2*i][0],(scale*0.5));
            dd_set_d(constraints->matrix[2*i][i+1],-1);
            dd_set_d(constraints->matrix[2*i+1][0],(scale*0.5));
            dd_set_d(constraints->matrix[2*i+1][i+1],1);
    }
    constraints->representation=dd_Inequality;
    cube = dd_DDMatrix2Poly(constraints, &err);
    dd_FreeMatrix(constraints);

    return cube;
};

/**
 * Project original polytope (in cdd format) to the first n dimensions
 */
void cdd_projection(dd_PolyhedraPtr *original,
                    dd_PolyhedraPtr *new,
                    size_t n,
                    dd_ErrorType *err)
{
    dd_MatrixPtr full=NULL,projected=NULL;
    full = dd_CopyInequalities(*original);
    dd_colrange j,d;
    dd_rowset redset,impl_linset;
    dd_colset delset;
    dd_rowindex newpos;

    d=full->colsize;
    set_initialize(&delset, d);
    for (j=n+1; j<d; j++){
        set_addelem(delset, j+1);
    }

    projected=dd_BlockElimination(full, delset, err);

    dd_MatrixCanonicalize(&projected,&impl_linset,&redset,&newpos,err);

    dd_FreeMatrix(full);
    projected->representation = dd_Inequality;
    *new = dd_DDMatrix2Poly(projected, err);
    dd_FreeMatrix(projected);
    set_free(delset);
    set_free(redset);
    set_free(impl_linset);
    free(newpos);

};

/**
 * Remove redundancies from cdd polytope inequalities
 */
dd_PolyhedraPtr cdd_minimize(dd_PolyhedraPtr *original,
                             dd_ErrorType *err)
{

    dd_rowset redset,impl_linset;
    dd_rowindex newpos;
    dd_MatrixPtr full=NULL;
    full = dd_CopyInequalities(*original);
    dd_MatrixCanonicalize(&full,&impl_linset,&redset,&newpos,err);

    full->representation = dd_Inequality;
    dd_PolyhedraPtr new = dd_DDMatrix2Poly(full, err);
    dd_FreeMatrix(full);
    set_free(redset);
    set_free(impl_linset);
    free(newpos);
    return new;

};

max_stars_repo_name: shaesaert/TuLiPXML

max_issues_repo_name: shaesaert/TuLiPXML

max_forks_repo_name: shaesaert/TuLiPXML


text:
/* -*- Mode: C; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
#pragma once

#include <mcbp/protocol/status.h>
#include <stdint.h>
#include <gsl/gsl>

struct EngineIface;

/**
 * Callback for any function producing stats.
 *
 * @param key the stat's key
 * @param klen length of the key
 * @param val the stat's value in an ascii form (e.g. text form of a number)
 * @param vlen length of the value
 * @param cookie magic callback cookie
 */
typedef void (*ADD_STAT)(const char* key,
                         const uint16_t klen,
                         const char* val,
                         const uint32_t vlen,
                         gsl::not_null<const void*> cookie);

/**
 * Callback for adding a response backet
 * @param key The key to put in the response
 * @param keylen The length of the key
 * @param ext The data to put in the extended field in the response
 * @param extlen The number of bytes in the ext field
 * @param body The data body
 * @param bodylen The number of bytes in the body
 * @param datatype This is currently not used and should be set to 0
 * @param status The status code of the return packet (see in protocol_binary
 *               for the legal values)
 * @param cas The cas to put in the return packet
 * @param cookie The cookie provided by the frontend
 * @return true if return message was successfully created, false if an
 *              error occured that prevented the message from being sent
 */
typedef bool (*ADD_RESPONSE)(const void* key,
                             uint16_t keylen,
                             const void* ext,
                             uint8_t extlen,
                             const void* body,
                             uint32_t bodylen,
                             uint8_t datatype,
                             cb::mcbp::Status status,
                             uint64_t cas,
                             const void* cookie);


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

#ifndef _MCAS_NUPM_ARENA_
#define _MCAS_NUPM_ARENA_

#include <nupm/region_descriptor.h>
#include <common/logging.h>
#include <common/string_view.h>
#include <gsl/pointers> /* not_null */
#include <sys/uio.h> /* ::iovec */

#include <cstddef>
#include <vector>

namespace nupm
{
	struct registry_memory_mapped;
	struct space_registered;
}

struct arena
  : private common::log_source
{
  using region_descriptor = nupm::region_descriptor;
  using registry_memory_mapped = nupm::registry_memory_mapped;
  using space_registered = nupm::space_registered;
  using string_view = common::string_view;

  arena(const common::log_source &ls) : common::log_source(ls) {}
  virtual ~arena() {}
  virtual void debug_dump() const = 0;
  virtual region_descriptor region_get(const string_view &id) = 0;
  virtual region_descriptor region_create(const string_view &id, gsl::not_null<registry_memory_mapped *> mh, std::size_t size) = 0;
  virtual void region_resize(gsl::not_null<space_registered *> mh, std::size_t size) = 0;
  /* It is unknown whether region_erase may be used on an open region.
   * arena_fs assumes that it may, just as ::unlink can be used against
   * an open file.
   */
  virtual void region_erase(const string_view &id, gsl::not_null<registry_memory_mapped *> mh) = 0;
  virtual std::size_t get_max_available() = 0;
  virtual bool is_file_backed() const = 0;
protected:
  using common::log_source::debug_level;
};

#endif


text:
#include <stdlib.h>
#include <gsl/vector/gsl_vector.h>
#include <gsl/gsl_errno.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/vector/prop_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#ifndef _density_constructions_h
#define _density_constructions_h

#include "montecarlo.h"
#include "two_body_analytical.h"
#include <float.h>
#include <limits.h>
#include <gsl/gsl_sort.h>

int r_final_uniform(double x, double y, double vz, double mu, double t_peri ,int N, double a, char * directory);
int Vector_distribution(double x, double y, double vz, double mu, double t_peri,  int N, double a, char * directory);

typedef struct{
  double v[6];

} state_Vector;

int Save_Vector_Index(state_Vector *state_vector, int index, int N, char* directory);
int Chi_Square_Euler(double x, double y, double vz, double mu,  int N, char * directory, size_t Npoints);

#endif  // ! _density_constructions_h


text:
/**
 *
 * @file qwrapper_ztrtri.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_ztrtri(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, PLASMA_enum diag,
                       int n, int nb,
                       PLASMA_Complex64_t *A, int lda,
                       PLASMA_sequence *sequence, PLASMA_request *request,
                       int iinfo)
{
    QUARK_Insert_Task(
        quark, CORE_ztrtri_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,      VALUE,
        sizeof(PLASMA_enum),                &diag,      VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    A,                 INOUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        sizeof(int),                        &iinfo,     VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_ztrtri_quark = PCORE_ztrtri_quark
#define CORE_ztrtri_quark PCORE_ztrtri_quark
#endif
void CORE_ztrtri_quark(Quark *quark)
{
    PLASMA_enum uplo;
    PLASMA_enum diag;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int iinfo;

    int info;

    quark_unpack_args_8(quark, uplo, diag, N, A, LDA, sequence, request, iinfo);
    info = LAPACKE_ztrtri_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo), lapack_const(diag),
        N, A, LDA);
    if ((sequence->status == PLASMA_SUCCESS) && (info > 0))
        plasma_sequence_flush(quark, sequence, request, iinfo + info);
}


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

#include <gsl/gsl_monte.h>		//Relevant GSL functions for VEGAS numerical integration
#include <gsl/gsl_monte_vegas.h>

#include "timer.c"			//Timing Functions
#include "plotfun.h"			//Automatic Plotting Function

void timer_start (void);
double timer_stop (void);

void pfun(double x[], double y[], double u[], double v[], int points);

extern double g (double *t, size_t dim, void *params);
struct my_params {double r;};

double f (double x1, double x2, double y1, double y2, double z1, double z2, double r);

double dipole_approx (double r);

int main (void)
{
	int np = 20;
	double distance[20];
	double vegas_energy[20];
	double res, err;
        double error1[20];
	double error2[20];

	size_t dim = 6;
	double xl[] = {0., 0., 0., 0., 0., 0.};
	double xu[] = {1., 1., 1., 1., 1., 1.};

	gsl_rng *r = gsl_rng_alloc (gsl_rng_taus2);
	unsigned long seed = 1UL;

	gsl_rng_set (r, seed);

	size_t calls = 1000000;

	double dd = (4.-1.001)/np;

	timer_start();

	for(int q = 0; q < np; q++)
	{
		distance[q] = 1.001 + q*dd;
	}

	//VEGAS Integration

	for(int j = 0; j < np; j++)
	{
		struct my_params params = {distance[j]};
		gsl_monte_function G = { &g, dim , &params};

		gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (dim);

		gsl_monte_vegas_init (sv);

		gsl_monte_vegas_integrate (&G, xl, xu, dim, calls / 10, r, sv, &res, &err);

		do
		{
			gsl_monte_vegas_integrate (&G, xl, xu, dim, calls, r, sv, &res, &err);
			fflush (stdout);
		}
		while (fabs (gsl_monte_vegas_chisq (sv) - 1.0) > 0.2);
		vegas_energy[j] = fabs(res);
		error1[j] = fabs(err);

		gsl_monte_vegas_free (sv);
	}

	double t1 = timer_stop();

	gsl_rng_free (r);

	int iterations = 1000000;
	int tot = 20;
	double x1, x2, y1, y2, z1, z2;
	double handmade_energy[20];

	gsl_rng *zed = gsl_rng_alloc (gsl_rng_taus2);

	gsl_rng_set (zed, seed);

	timer_start();

	//Handmade Monte Carlo Integrator

	for(int j = 0; j < 20.; j++)
	{
	    double d = 1.001 + j * dd;
	    double ival = 0;
	    double i2val = 0;
	    for(int z = 0; z < 16; z++)
	    {
	      double value = 0;
       		for (int q = 0; q < iterations; q++)
        	{
            		x1 = gsl_rng_uniform(zed);
            		x2 = gsl_rng_uniform(zed);
            		y1 = gsl_rng_uniform(zed);
            		y2 = gsl_rng_uniform(zed);
            		z1 = gsl_rng_uniform(zed);
            		z2 = gsl_rng_uniform(zed);
            		value += f(x1, x2, y1, y2, z1, z2, d);
       		}
		ival += value;
		i2val += value*value;
	    }

	    double result = ival/iterations/16;
	    double result2 = i2val/iterations/iterations/16;

	    error2[j] = sqrt(fabs((result*result)-result2));

	    handmade_energy[j] = fabs(result);
	}

	double t2 = timer_stop();

	//Dipole Approximation

	double dipole_energy[20];
	for(int u = 0; u < tot; u++)
	{
		dipole_energy[u] = dipole_approx(distance[u]);
	}

	printf("Distance   VEGAS      Error      Handmade   Error      Dipole\n");

	for(int w = 0; w < tot; w++)
	{
	  printf("%.6f   %.6f   %.6f   %.6f   %.6f   %.6f\n", distance[w], vegas_energy[w], error1[w], handmade_energy[w], error2[w], dipole_energy[w]);
	}

	printf("Time for VEGAS:  %.6f\nTime for HANDMADE:  %.6f\n", t1, t2);

	//Plots the 3 methods vs. distance automatically (function in plotfun.h)

	pfun(distance, vegas_energy, handmade_energy, dipole_energy, 20);

	return 0;
}

double g (double *t, size_t dim, void *params)
{
    	double dist2, delta;
    	double rho1, rho2;
    	int sdim = ((int) dim)/2;
    	double r = *((double *) params);

    	dist2 = 0.;
    	for (int i = 0; i < sdim; i++)
    	{
        	delta = t[i] - t[i + sdim];
        	if (i == 0)
        	{
           		delta += r;
        	}
        	dist2 += delta*delta;
    	}

    	double norm = 512./(9.*(M_PI - 2.));
    	rho1 = norm * atan(2*(t[0]-.5))*pow(sin(M_PI*t[1]),4.)*pow(sin(M_PI*t[2]),4.);
    	rho2 = norm * atan(2*(t[3]-.5))*pow(sin(M_PI*t[4]),4.)*pow(sin(M_PI*t[5]),4.);

    	return rho1*rho2/sqrt(dist2);
}

double f (double x1, double x2, double y1, double y2, double z1, double z2, double r)
{
	double dist2;
    	double rho1, rho2;

    	dist2 = (x1 - x2 + r)*(x1 - x2 + r) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2);

    	double norm = 512./(9.*(M_PI - 2.));
    	rho1 = norm * atan(2*(x1-.5))*pow(sin(M_PI*y1),4.)*pow(sin(M_PI*z1),4.);
    	rho2 = norm * atan(2*(x2-.5))*pow(sin(M_PI*y2),4.)*pow(sin(M_PI*z2),4.);

    	return rho1*rho2/sqrt(dist2);
}

double dipole_approx (double r)
{
	return 2/(r*r*r);
}


text:
// Copyright (c) 2005 - 2011 Marc de Kamps
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
//
//    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
//    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation
//      and/or other materials provided with the distribution.
//    * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
//      without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
// USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
//      If you use this software in work leading to a scientific publication, you should cite
//      the 'currently valid reference', which can be found at http://miind.sourceforge.net
#ifndef _CODE_LIBS_NUMTOOLSLIB_RANDOMGENERATOR_INCLUDE_GUARD
#define _CODE_LIBS_NUMTOOLSLIB_RANDOMGENERATOR_INCLUDE_GUARD

#include <vector>
#include <iostream>
#include <gsl/gsl_rng.h>
#include "ChangeSeedException.h"

using std::vector;

namespace NumtoolsLib
{

	//! Global default seed for RandomGenerator

	const long GLOBAL_SEED = 9875987L;

	//! This class contains a specific uniform random generator that can be used to generate other distributions
	//! GaussianDistribution and UniformDistribution, for example, use this random generator, which gsl_rng_mt19937
	//! of the Gnu Scietific Library. This class can be used to modify the seed and to keep track of the number
	//! of draws made by the generator. Distributions use one copy of this class.

	class RandomGenerator
	{

	public:

		RandomGenerator(long lseed = GLOBAL_SEED);
		~RandomGenerator();

		double	NextSampleValue();

	private:

		RandomGenerator(const RandomGenerator&);
		RandomGenerator& operator=(const RandomGenerator&);

		double	Ran2(long*);

		// new, uniformly distributed value

		void    AddOne         (){ std::cout << "zopa" << std::endl; ++_number_of_draws; }

		long			_initial_seed;    // initial seed value
		unsigned long	_number_of_draws; // number of calls to NextSampleValue

		gsl_rng*		_p_generator;

	}; // end of RandomGenerator

	//! Global Random Generator
	extern RandomGenerator GLOBAL_RANDOM_GENERATOR;

	inline double RandomGenerator::NextSampleValue()
	{
		AddOne();
		return gsl_rng_uniform(_p_generator);
	}

} // end of Numtools

#endif // include guard

max_stars_repo_name: dekamps/miind

max_issues_repo_name: dekamps/miind

max_forks_repo_name: dekamps/miind