diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8b3c76b..2be014b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -19,7 +19,7 @@ ci: repos: # Standard hooks - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.4.0 + rev: v5.0.0 hooks: - id: check-added-large-files - id: check-case-conflict @@ -35,21 +35,21 @@ repos: # Black, the code formatter, natively supports pre-commit - repo: https://github.com/psf/black-pre-commit-mirror - rev: 23.7.0 + rev: 25.1.0 hooks: - id: black exclude: ^(docs) # Check linting and style issues - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.0.287" + rev: "v0.11.13" hooks: - id: ruff args: ["--fix", "--show-fixes"] # Changes tabs to spaces - repo: https://github.com/Lucas-C/pre-commit-hooks - rev: v1.5.4 + rev: v1.5.5 hooks: - id: remove-tabs exclude: ^(docs) diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json index 2150588..67eb602 100644 --- a/.vscode/c_cpp_properties.json +++ b/.vscode/c_cpp_properties.json @@ -24,4 +24,4 @@ } ], "version": 4 -} \ No newline at end of file +} diff --git a/.vscode/settings.json b/.vscode/settings.json index 5b6c867..f3e1b94 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -71,4 +71,4 @@ "valarray": "cpp", "variant": "cpp" } -} \ No newline at end of file +} diff --git a/CMakeLists.txt b/CMakeLists.txt index c6620e4..988beef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,8 @@ find_package(pybind11 CONFIG REQUIRED) # Include the include/barry library include_directories(include/barry) -python_add_library(_core MODULE src/main.cpp src/terms.cpp src/model.cpp WITH_SOABI) +python_add_library(_core MODULE src/main.cpp src/terms.cpp src/model.cpp + WITH_SOABI) target_link_libraries(_core PRIVATE pybind11::headers) target_compile_definitions(_core PRIVATE VERSION_INFO=${PROJECT_VERSION}) diff --git a/Makefile b/Makefile index 13f8a10..2b8237f 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,5 @@ build: - pip3 install . - + pip3 install . + update: - rsync -avz ../barry/include/barry include/ \ No newline at end of file + rsync -avz ../barry/include/barry include/ diff --git a/README.md b/README.md index 54abecd..c14f3f8 100644 --- a/README.md +++ b/README.md @@ -66,7 +66,7 @@ m.term_formula(obj, "{y0}") m.term_formula(obj, "{y1}") m.term_formula(obj, "{0y0, y1}") obj.init() -obj.print() +obj.print() ``` Num. of Arrays : 6 diff --git a/README.qmd b/README.qmd index 63b991c..09f4ce1 100644 --- a/README.qmd +++ b/README.qmd @@ -42,7 +42,7 @@ ERGMs. More information in [this preprint](https://arxiv.org/abs/2211.00627). - `pip install ./pydefm` -# Examples +# Examples ## Base setup @@ -64,7 +64,7 @@ x = np.array([1, 2.0, 3.4, 3, 1, 2]) id = np.array([1, 1, 1, 2, 2, 2]) ``` -Once we have the needed data -- the binary array `y`, the covariates `x` and the ids `id` -- we can create a `defm` object. +Once we have the needed data -- the binary array `y`, the covariates `x` and the ids `id` -- we can create a `defm` object. ```{python} # Creating the object @@ -83,7 +83,7 @@ m.term_formula(obj, "{y0}") m.term_formula(obj, "{y1}") m.term_formula(obj, "{0y0, y1}") obj.init() -obj.print() +obj.print() ``` ```{python} @@ -122,4 +122,4 @@ m.simulate(obj, par) # Acknowledgements -This port work was supported by the Office of the Assistant Secretary of Defense for Health Affairs through the Epilepsy Research Program under Award No. HT94252310221. Opinions, interpretations, conclusions, and recommendations are those of the author and are not necessarily endorsed by the Department of Defense. \ No newline at end of file +This port work was supported by the Office of the Assistant Secretary of Defense for Health Affairs through the Epilepsy Research Program under Award No. HT94252310221. Opinions, interpretations, conclusions, and recommendations are those of the author and are not necessarily endorsed by the Department of Defense. diff --git a/include/barry/barray-bones.hpp b/include/barry/barray-bones.hpp index 0fb6849..7635bbf 100644 --- a/include/barry/barray-bones.hpp +++ b/include/barry/barray-bones.hpp @@ -4,7 +4,7 @@ // #include "cell-bones.hpp" // #include "barraycell-bones.hpp" -#ifndef BARRAY_BONES_HPP +#ifndef BARRAY_BONES_HPP #define BARRAY_BONES_HPP 1 template @@ -15,7 +15,7 @@ class BArrayCell_const; /** * @brief Baseline class for binary arrays. - * + * * `BArray` class objects are arbitrary arrays * in which non-empty cells hold data of type `Cell_Type`. The non-empty cells * are stored by row and indexed using `unordered_map`s, i.e. @@ -43,8 +43,8 @@ class BArray { static const bool dense = false; public: - - /** + + /** * This is as a reference, if we need to iterate through the cells and we need * to keep track which were visited, we use this as a reference. So that if * cell.visited = true and visited = true, it means that we haven't been here @@ -52,11 +52,11 @@ class BArray { * beginning of the routine. */ bool visited = false; - + /** * @name Constructors - * + * * @param N_ Number of rows * @param M_ Number of columns * @param source An unsigned vector ranging from 0 to N_ @@ -64,13 +64,13 @@ class BArray { * @param target When `true` tries to add repeated observations. */ ///@{ - + /** @brief Zero-size array */ BArray() : N(0u), M(0u), NCells(0u), el_ij(0u), el_ji(0u) {}; - + /** @brief Empty array */ BArray (size_t N_, size_t M_) : N(N_), M(M_), NCells(0u), el_ij(N_), el_ji(M_) {}; - + /** @brief Edgelist with data */ BArray ( size_t N_, size_t M_, @@ -79,7 +79,7 @@ class BArray { const std::vector< Cell_Type > & value, bool add = true ); - + /** @brief Edgelist with no data (simpler) */ BArray ( size_t N_, size_t M_, @@ -87,10 +87,10 @@ class BArray { const std::vector< size_t > & target, bool add = true ); - + /** @brief Copy constructor */ BArray(const BArray & Array_, bool copy_data = false); - + /** @brief Assignment constructor */ BArray & operator=(const BArray & Array_); @@ -100,20 +100,20 @@ class BArray { /** @brief Move assignment */ BArray & operator=(BArray && x) noexcept; ///@} - + bool operator==(const BArray & Array_); ~BArray(); - + // In principle, copy can be faster by using openmp on the rows // since those are independent. // BArray(BArray & A); - + /** * @brief Set the data object - * - * @param data_ - * @param delete_data_ + * + * @param data_ + * @param delete_data_ */ ///@{ void set_data(Data_Type * data_, bool delete_data_ = false); @@ -123,11 +123,11 @@ class BArray { const Data_Type & D() const; void flush_data(); ///@} - + // Function to access the elements // bool check_cell void out_of_range(size_t i, size_t j) const; - Cell_Type get_cell(size_t i, size_t j, bool check_bounds = true) const; + Cell_Type get_cell(size_t i, size_t j, bool check_bounds = true) const; std::vector< Cell_Type > get_col_vec(size_t i, bool check_bounds = true) const; std::vector< Cell_Type > get_row_vec(size_t i, bool check_bounds = true) const; void get_col_vec(std::vector< Cell_Type > * x, size_t i, bool check_bounds = true) const; @@ -137,12 +137,12 @@ class BArray { /** * @brief Get the edgelist - * + * * `Entries` is a class with three objects: Two `std::vector` with the row and * column coordinates respectively, and one `std::vector` with the corresponding * value of the cell. - * - * @return Entries + * + * @return Entries */ Entries get_entries() const; @@ -170,37 +170,37 @@ class BArray { * delete/add), or, in the case of `swap_cells`, check if either of both * cells exists/don't exist. */ - ///@{ + ///@{ BArray & operator+=(const std::pair & coords); BArray & operator-=(const std::pair & coords); BArrayCell operator()(size_t i, size_t j, bool check_bounds = true); const Cell_Type operator()(size_t i, size_t j, bool check_bounds = true) const; - + void rm_cell(size_t i, size_t j, bool check_bounds = true, bool check_exists = true); - + void insert_cell(size_t i, size_t j, const Cell< Cell_Type > & v, bool check_bounds, bool check_exists); void insert_cell(size_t i, size_t j, Cell< Cell_Type > && v, bool check_bounds, bool check_exists); void insert_cell(size_t i, size_t j, Cell_Type v, bool check_bounds, bool check_exists); - + void swap_cells( size_t i0, size_t j0, size_t i1, size_t j1, bool check_bounds = true, int check_exists = CHECK::BOTH, int * report = nullptr ); - + void toggle_cell(size_t i, size_t j, bool check_bounds = true, int check_exists = EXISTS::UKNOWN); void toggle_lock(size_t i, size_t j, bool check_bounds = true); ///@} - + /**@name Column/row wise interchange*/ ///@{ void swap_rows(size_t i0, size_t i1, bool check_bounds = true); void swap_cols(size_t j0, size_t j1, bool check_bounds = true); - + void zero_row(size_t i, bool check_bounds = true); void zero_col(size_t j, bool check_bounds = true); ///@} - + void transpose(); void clear(bool hard = true); void resize(size_t N_, size_t M_); @@ -208,14 +208,14 @@ class BArray { // Advances operators // void toggle_iterator - + // Misc void print(const char * fmt = nullptr, ...) const; void print_n(size_t nrow, size_t ncol, const char * fmt = nullptr, ...) const; /** * @name Arithmetic operators - * + * */ ///@{ BArray& operator+=(const BArray& rhs); @@ -223,11 +223,11 @@ class BArray { BArray& operator-=(const BArray& rhs); BArray& operator-=(const Cell_Type & rhs); - + BArray& operator/=(const Cell_Type & rhs); BArray& operator*=(const Cell_Type & rhs); ///@} - + // /** // * @name Casting between types // */ diff --git a/include/barry/barray-iterator.hpp b/include/barry/barray-iterator.hpp index 610de85..799ddc7 100644 --- a/include/barry/barray-iterator.hpp +++ b/include/barry/barray-iterator.hpp @@ -3,33 +3,33 @@ // #include "typedefs.hpp" // #include "barray-bones.hpp" -#ifndef BARRAY_ITERATOR_HPP +#ifndef BARRAY_ITERATOR_HPP #define BARRAY_ITERATOR_HPP 1 template class ConstBArrayRowIter { public: - + size_t current_row, current_col; typename Row_type::const_iterator iter; const BArray * Array; - + ConstBArrayRowIter(const BArray * Array_) : Array(Array_) { - + // Finding the first entry of the iterator for (size_t i = 0u; i < Array->nrow(); ++i) if (A_ROW(i).size() != 0u) { iter = A_ROW(i).begin(); break; } - + return; }; ~ConstBArrayRowIter() {}; - + // operat - + }; #endif diff --git a/include/barry/barray-meat-operators.hpp b/include/barry/barray-meat-operators.hpp index c20168d..4b3deb2 100644 --- a/include/barry/barray-meat-operators.hpp +++ b/include/barry/barray-meat-operators.hpp @@ -35,7 +35,7 @@ BARRAY_TEMPLATE(BARRAY_TYPE()&, operator+=) ( // Must be compatible checkdim_(*this, rhs); - + for (size_t i = 0u; i < nrow(); ++i) for (size_t j = 0u; j < ncol(); ++j) this->operator()(i, j) += rhs.get_cell(i, j); @@ -62,7 +62,7 @@ BARRAY_TEMPLATE(BARRAY_TYPE()&, operator-=) ( // Must be compatible checkdim_(*this, rhs); - + for (size_t i = 0u; i < nrow(); ++i) { for (size_t j = 0u; j < ncol(); ++j) { this->operator()(i, j) -= rhs.get_cell(i, j); @@ -126,4 +126,4 @@ BARRAY_TEMPLATE(BARRAY_TYPE()&, operator/=) ( #undef ROW #undef COL -#endif \ No newline at end of file +#endif diff --git a/include/barry/barray-meat.hpp b/include/barry/barray-meat.hpp index ea35005..7b3deac 100644 --- a/include/barry/barray-meat.hpp +++ b/include/barry/barray-meat.hpp @@ -8,14 +8,14 @@ template class Cell_const; #ifndef BARRY_BARRAY_MEAT_HPP -#define BARRY_BARRAY_MEAT_HPP +#define BARRY_BARRAY_MEAT_HPP #define ROW(a) this->el_ij[a] #define COL(a) this->el_ji[a] template -Cell BArray::Cell_default = Cell(static_cast(1.0)); +Cell BArray::Cell_default = Cell(static_cast(1.0)); // Edgelist with data @@ -26,38 +26,38 @@ template inline BArray & value, bool add ) { - + if (source.size() != target.size()) throw std::length_error("-source- and -target- don't match on length."); if (source.size() != value.size()) throw std::length_error("-sorce- and -value- don't match on length."); - + // Initializing N = N_; M = M_; el_ij.resize(N); el_ji.resize(M); - - + + // Writing the data for (size_t i = 0u; i < source.size(); ++i) { - + // Checking range bool empty = this->is_empty(source[i], target[i], true); if (add && !empty) { ROW(source[i])[target[i]].add(value[i]); continue; - } - + } + if (!empty) throw std::logic_error("The value already exists. Use 'add = true'."); - + this->insert_cell(source[i], target[i], value[i], false, false); } - + return; - + } // Edgelist with data @@ -68,41 +68,41 @@ inline BArray::BArray ( const std::vector< size_t > & target, bool add ) { - + std::vector< Cell_Type > value(source.size(), (Cell_Type) 1.0); if (source.size() != target.size()) throw std::length_error("-source- and -target- don't match on length."); if (source.size() != value.size()) throw std::length_error("-sorce- and -value- don't match on length."); - + // Initializing N = N_; M = M_; - + el_ij.resize(N); el_ji.resize(M); - - + + // Writing the data for (size_t i = 0u; i < source.size(); ++i) { - + // Checking range if ((source[i] >= N_) || (target[i] >= M_)) throw std::range_error("Either source or target point to an element outside of the range by (N,M)."); - + // Checking if it exists auto search = ROW(source[i]).find(target[i]); if (search != ROW(source[i]).end()) { if (!add) throw std::logic_error("The value already exists. Use 'add = true'."); - + // Increasing the value (this will automatically update the // other value) ROW(source[i])[target[i]].add(value[i]); continue; } - + // Adding the value and creating a pointer to it ROW(source[i]).emplace( std::pair >( @@ -110,7 +110,7 @@ inline BArray::BArray ( Cell< Cell_Type >(value[i], visited) ) ); - + COL(target[i]).emplace( source[i], &ROW(source[i])[target[i]] @@ -119,9 +119,9 @@ inline BArray::BArray ( NCells++; } - + return; - + } template @@ -130,11 +130,11 @@ inline BArray::BArray ( bool copy_data ) : N(Array_.N), M(Array_.M) { - + // Dimensions // el_ij.resize(N); // el_ji.resize(M); - + std::copy(Array_.el_ij.begin(), Array_.el_ij.end(), std::back_inserter(el_ij)); std::copy(Array_.el_ji.begin(), Array_.el_ji.end(), std::back_inserter(el_ji)); @@ -149,7 +149,7 @@ inline BArray::BArray ( this->NCells = Array_.NCells; this->visited = Array_.visited; - + // Data if (Array_.data != nullptr) { @@ -168,42 +168,42 @@ inline BArray::BArray ( } } - + return; - + } template inline BArray & BArray:: operator= ( const BArray & Array_ ) { - + // Clearing if (this != &Array_) { - + this->clear(true); this->resize(Array_.N, Array_.M); - + // Entries for (size_t i = 0u; i < N; ++i) { - + if (Array_.nnozero() == nnozero()) break; - - for (auto& r : Array_.row(i, false)) + + for (auto& r : Array_.row(i, false)) this->insert_cell(i, r.first, r.second.value, false, false); - + } - + // Data if (data != nullptr) { if (delete_data) delete data; - + data = nullptr; delete_data = false; @@ -216,11 +216,11 @@ inline BArray & BArray:: operator= delete_data = true; } - + } - + return *this; - + } template inline BArray::BArray ( @@ -233,16 +233,16 @@ template inline BArrayclear(true); this->resize(x.N, x.M); - + // Entries for (size_t i = 0u; i < N; ++i) { - + if (x.nnozero() == nnozero()) break; - - for (auto& r : x.row(i, false)) + + for (auto& r : x.row(i, false)) this->insert_cell(i, r.first, r.second.value, false, false); - + } // Managing data @@ -268,24 +268,24 @@ template inline BArray inline BArray & BArray:: operator= ( BArray && x ) noexcept { - + // Clearing if (this != &x) { - + this->clear(true); this->resize(x.N, x.M); - + // Entries for (size_t i = 0u; i < N; ++i) { - + if (x.nnozero() == nnozero()) break; - - for (auto& r : x.row(i, false)) + + for (auto& r : x.row(i, false)) this->insert_cell(i, r.first, r.second.value, false, false); - + } - + // Data if (data != nullptr) { @@ -307,51 +307,51 @@ template inline BArray inline bool BArray:: operator== ( const BArray & Array_ ) { - + // Dimension and number of cells used if ((N != Array_.nrow()) | (M != Array_.ncol()) | (NCells != Array_.nnozero())) return false; - + // One holds, and the other doesn't. if ((!data & Array_.data) | (data & !Array_.data)) return false; - + if (this->el_ij != Array_.el_ij) return false; - + return true; } template inline BArray::~BArray () { - + if (delete_data && (data != nullptr)) delete data; - + return; } template inline void BArray:: set_data ( Data_Type * data_, bool delete_data_ -) { +) { if ((data != nullptr) && delete_data) delete data; - + data = data_; delete_data = delete_data_; - + return; - + } template inline Data_Type * BArray:: D_ptr () @@ -404,28 +404,28 @@ template inline void BArray inline Cell_Type BArray:: get_cell ( size_t i, size_t j, bool check_bounds ) const { - - // Checking boundaries + + // Checking boundaries if (check_bounds) out_of_range(i,j); - + if (ROW(i).size() == 0u) return (Cell_Type) 0.0; - + // If it is not empty, then find and return auto search = ROW(i).find(j); if (search != ROW(i).end()) return search->second.value; - + // This is if it is empty return (Cell_Type) 0.0; - + } template inline std::vector< Cell_Type > BArray:: get_row_vec ( @@ -433,14 +433,14 @@ template inline std::vector< Cell_Type > bool check_bounds ) const { - // Checking boundaries - if (check_bounds) + // Checking boundaries + if (check_bounds) out_of_range(i, 0u); std::vector< Cell_Type > ans(ncol(), (Cell_Type) false); - for (const auto & iter : row(i, false)) + for (const auto & iter : row(i, false)) ans[iter.first] = iter.second.value; //this->get_cell(i, iter->first, false); - + return ans; } @@ -451,13 +451,13 @@ template inline void BArrayat(iter.first) = iter.second.value; // this->get_cell(i, iter->first, false); - + } template inline std::vector< Cell_Type > BArray:: get_col_vec ( @@ -465,14 +465,14 @@ template inline std::vector< Cell_Type > bool check_bounds ) const { - // Checking boundaries - if (check_bounds) + // Checking boundaries + if (check_bounds) out_of_range(0u, i); std::vector< Cell_Type > ans(nrow(), (Cell_Type) false); - for (const auto iter : col(i, false)) + for (const auto iter : col(i, false)) ans[iter.first] = iter.second->value;//this->get_cell(iter->first, i, false); - + return ans; } @@ -483,13 +483,13 @@ template inline void BArrayat(iter.first) = iter.second->value;//this->get_cell(iter->first, i, false); - + } template inline const Row_type< Cell_Type > & BArray:: row ( @@ -513,25 +513,25 @@ template inline const Col_type< Cell_Typ out_of_range(0u, i); return this->el_ji[i]; - + } template inline Entries< Cell_Type > BArray:: get_entries () const { - + Entries res(NCells); - + for (size_t i = 0u; i < N; ++i) { - + if (ROW(i).size() == 0u) continue; - + for (auto col = ROW(i).begin(); col != ROW(i).end(); ++col) { res.source.push_back(i), res.target.push_back(col->first), res.val.push_back(col->second.value); } } - + return res; } @@ -540,20 +540,20 @@ template inline bool BArray inline Cell< Cell_Type > BArra template inline BArray & BArray:: operator+= ( const std::pair & coords ) { - + this->insert_cell( coords.first, coords.second, this->Cell_default, true, true ); - + return *this; - + } template inline BArray & BArray:: operator-= ( const std::pair & coords ) { - + this->rm_cell( coords.first, coords.second, true, true ); - + return *this; - + } template -inline BArrayCell BArray::operator()( +inline BArrayCell BArray::operator()( size_t i, size_t j, bool check_bounds ) { - + return BArrayCell(this, i, j, check_bounds); - + } template -inline const Cell_Type BArray::operator() ( +inline const Cell_Type BArray::operator() ( size_t i, size_t j, bool check_bounds ) const { - + return get_cell(i, j, check_bounds); - + } template inline void BArray:: rm_cell ( @@ -632,33 +632,33 @@ template inline void BArray inline void BArray:: insert_cell ( @@ -667,45 +667,45 @@ template inline void BArray & v, bool check_bounds, bool check_exists - ) { - + ) { + if (check_bounds) - out_of_range(i,j); - + out_of_range(i,j); + if (check_exists) { - + // Checking if nothing here, then we move along if (ROW(i).size() == 0u) { - + ROW(i).insert(std::pair< size_t, Cell>(j, v)); COL(j).emplace(i, &ROW(i)[j]); NCells++; return; - + } - - // In this case, the row exists, but we are checking that the value is empty + + // In this case, the row exists, but we are checking that the value is empty if (ROW(i).find(j) == ROW(i).end()) { - - ROW(i).insert(std::pair< size_t, Cell>(j, v)); + + ROW(i).insert(std::pair< size_t, Cell>(j, v)); COL(j).emplace(i, &ROW(i)[j]); NCells++; - + } else { throw std::logic_error("The cell already exists."); } - - + + } else { - + ROW(i).insert(std::pair< size_t, Cell>(j, v)); COL(j).emplace(i, &ROW(i)[j]); NCells++; - + } - + return; - + } template inline void BArray:: insert_cell ( @@ -714,45 +714,45 @@ template inline void BArray && v, bool check_bounds, bool check_exists - ) { - + ) { + if (check_bounds) - out_of_range(i,j); - + out_of_range(i,j); + if (check_exists) { - + // Checking if nothing here, then we move along if (ROW(i).size() == 0u) { - + ROW(i).insert(std::pair< size_t, Cell>(j, v)); COL(j).emplace(i, &ROW(i)[j]); NCells++; return; - + } - - // In this case, the row exists, but we are checking that the value is empty + + // In this case, the row exists, but we are checking that the value is empty if (ROW(i).find(j) == ROW(i).end()) { - - ROW(i).insert(std::pair< size_t, Cell>(j, v)); + + ROW(i).insert(std::pair< size_t, Cell>(j, v)); COL(j).emplace(i, &ROW(i)[j]); NCells++; - + } else { throw std::logic_error("The cell already exists."); } - - + + } else { - + ROW(i).insert(std::pair< size_t, Cell>(j, v)); COL(j).emplace(i, &ROW(i)[j]); NCells++; - + } - + return; - + } template inline void BArray:: insert_cell ( @@ -762,7 +762,7 @@ template inline void BArray(v, visited), check_bounds, check_exists); } @@ -774,24 +774,24 @@ template inline void BArray inline void BArray c1(std::move(ROW(i1)[j1])); rm_cell(i1, j1, false, false); - + // Inserting the cells by reference, these will be deleted afterwards insert_cell(i0, j0, c1, false, false); insert_cell(i1, j1, c0, false, false); - + return; - + } - + bool check0, check1; if (check_exists == CHECK::BOTH) { - + check0 = !is_empty(i0, j0, false); check1 = !is_empty(i1, j1, false); - + } else if (check_exists == CHECK::ONE) { - + check0 = !is_empty(i0, j0, false); check1 = true; - + } else if (check_exists == CHECK::TWO) { - + check0 = true; check1 = !is_empty(i1, j1, false); - + } - - if (report != nullptr) + + if (report != nullptr) (*report) = EXISTS::NONE; - + // If both cells exists if (check0 & check1) { - - if (report != nullptr) + + if (report != nullptr) (*report) = EXISTS::BOTH; - + // If source and target coincide, we do nothing - if ((i0 == i1) && (j0 == j1)) + if ((i0 == i1) && (j0 == j1)) return; - + Cell c0(std::move(ROW(i0)[j0])); rm_cell(i0, j0, false, false); Cell c1(std::move(ROW(i1)[j1])); rm_cell(i1, j1, false, false); - + insert_cell(i0, j0, c1, false, false); insert_cell(i1, j1, c0, false, false); - + } else if (!check0 & check1) { // If only the second exists - - if (report != nullptr) + + if (report != nullptr) (*report) = EXISTS::TWO; - + insert_cell(i0, j0, ROW(i1)[j1], false, false); rm_cell(i1, j1, false, false); - + } else if (check0 & !check1) { - - if (report != nullptr) + + if (report != nullptr) (*report) = EXISTS::ONE; - + insert_cell(i1, j1, ROW(i0)[j0], false, false); rm_cell(i0, j0, false, false); - + } - + return; } @@ -876,31 +876,31 @@ template inline void BArray::Cell_default, false, false); ROW(i)[j].visited = visited; } else rm_cell(i, j, false, false); - + } else if (check_exists == EXISTS::AS_ONE) { - + rm_cell(i, j, false, false); - + } else if (check_exists == EXISTS::AS_ZERO) { - + insert_cell(i, j, BArray::Cell_default, false, false); ROW(i)[j].visited = visited; - + } - + return; - + } template inline void BArray:: swap_rows ( @@ -908,42 +908,42 @@ template inline void BArray inline void BArray col_tmp = COL(j1); Col_type col1 = COL(j0); for (auto iter = col1.begin(); iter != col1.end(); ++iter) { - + // Swapping values (col-wise) swap_cells(iter->first, j0, iter->first, j1, false, CHECK::TWO, &status); - + // Need to remove it, so we don't swap that as well if (status == EXISTS::BOTH) col_tmp.erase(iter->first); } - + // If there's anything left to move, we start moving it, otherwise, we just // skip it if (col_tmp.size() != 0u) { - + for (auto iter = col_tmp.begin(); iter != col_tmp.end(); ++iter) { insert_cell(iter->first, j0, *iter->second, false, false); rm_cell(iter->first, j1); } - + } - + } else if (check0 && !check1) { - + // 1 is empty, so we just add new cells and remove the other ones for (auto iter = COL(j0).begin(); iter != COL(j0).begin(); ++iter) insert_cell(iter->first, j1, *iter->second, false, false); - + // Setting the column to be zero COL(j0).empty(); - + } else if (!check0 && check1) { - + // 1 is empty, so we just add new cells and remove the other ones for (auto iter = COL(j1).begin(); iter != COL(j1).begin(); ++iter) { - + // Swapping values (col-wise) insert_cell(iter->first, j0, *iter->second, false, false); } - + // Setting the column to be zero COL(j1).empty(); - + } - - + + return; } @@ -1024,101 +1024,101 @@ template inline void BArrayfirst, false, false); - + return; - + } template inline void BArray:: zero_col ( size_t j, bool check_bounds ) { - + if (check_bounds) out_of_range(0u, j); - + // Nothing to do if (COL(j).size() == 0u) return; - + // Else, remove all elements auto col0 = COL(j); - for (auto col = col0.begin(); col != col0.end(); ++col) + for (auto col = col0.begin(); col != col0.end(); ++col) rm_cell(col->first, j, false, false); - + return; - + } template inline void BArray:: transpose () { - - // Start by flipping the switch + + // Start by flipping the switch visited = !visited; - + // Do we need to resize (increase) either? if (N > M) el_ji.resize(N); else if (N < M) el_ij.resize(M); - + // size_t N0 = N, M0 = M; int status; for (size_t i = 0u; i < N; ++i) { - + // Do we need to move anything? if (ROW(i).size() == 0u) continue; - + // We now iterate changing rows Row_type row = ROW(i); for (auto col = row.begin(); col != row.end(); ++col) { - + // Skip if in the diagoal if (i == col->first) { ROW(i)[i].visited = visited; continue; } - + // We have not visited this yet, we need to change that if (ROW(i)[col->first].visited != visited) { - + // First, swap the contents swap_cells(i, col->first, col->first, i, false, CHECK::TWO, &status); - + // Changing the switch if (status == EXISTS::BOTH) ROW(i)[col->first].visited = visited; - + ROW(col->first)[i].visited = visited; - + } - + } - + } - + // Shreding. Note that no information should have been lost since, hence, no // change in NCells. if (N > M) el_ij.resize(M); else if (N < M) el_ji.resize(N); - + // Swapping the values std::swap(N, M); - + return; } @@ -1126,55 +1126,55 @@ template inline void BArray inline void BArray:: clear ( bool hard ) { - + if (hard) { - + el_ji.clear(); el_ij.clear(); - + el_ij.resize(N); el_ji.resize(M); NCells = 0u; - + } else { - + for (size_t i = 0u; i < N; ++i) zero_row(i, false); - + } - + return; - + } template inline void BArray:: resize ( size_t N_, size_t M_ ) { - + // Removing rows if (N_ < N) for (size_t i = N_; i < N; ++i) zero_row(i, false); - + // Removing cols if (M_ < M) for (size_t j = M_; j < M; ++j) zero_col(j, false); - + // Resizing will invalidate pointers and values out of range if (M_ != M) { el_ji.resize(M_); M = M_; } - + if (N_ != N) { el_ij.resize(N_); N = N_; } - - + + return; } @@ -1184,12 +1184,12 @@ inline void BArray:: reserve () { #ifdef BARRAY_USE_UNORDERED_MAP for (size_t i = 0u; i < N; i++) ROW(i).reserve(M); - + for (size_t i = 0u; i < M; i++) COL(i).reserve(N); #endif return; - + } template @@ -1197,13 +1197,13 @@ inline void BArray:: print ( const char * fmt, ... ) const { - + std::va_list args; va_start(args, fmt); print_n(N, M, fmt, args); - va_end(args); - + va_end(args); + return; } @@ -1240,9 +1240,9 @@ inline void BArray:: print_n ( for (size_t j = 0u; j < ncol; ++j) { if (this->is_empty(i, j, false)) printf_barry(" . "); - else + else printf_barry(" %.2f ", static_cast(this->get_cell(i, j, false))); - + } printf_barry("\n"); @@ -1257,8 +1257,8 @@ inline void BArray:: print_n ( if (nrow < N || ncol < M) printf_barry("\n"); - - + + return; } @@ -1267,4 +1267,3 @@ inline void BArray:: print_n ( #undef COL #endif - diff --git a/include/barry/barraycell-bones.hpp b/include/barry/barraycell-bones.hpp index bc6e274..6d91af0 100644 --- a/include/barry/barraycell-bones.hpp +++ b/include/barry/barraycell-bones.hpp @@ -6,14 +6,14 @@ template class BArrayCell { private: - + BArray * Array; size_t i; size_t j; - + public: - - BArrayCell(BArray * Array_, size_t i_, size_t j_, bool check_bounds = true) : + + BArrayCell(BArray * Array_, size_t i_, size_t j_, bool check_bounds = true) : Array(Array_), i(i_), j(j_) { if (check_bounds) @@ -37,7 +37,7 @@ class BArrayCell { operator Cell_Type() const; bool operator==(const Cell_Type & val) const; - + }; @@ -45,14 +45,14 @@ class BArrayCell { template class BArrayCell_const { private: - + const BArray * Array; size_t i; size_t j; - + public: - - BArrayCell_const(const BArray * Array_, size_t i_, size_t j_, bool check_bounds = true) : + + BArrayCell_const(const BArray * Array_, size_t i_, size_t j_, bool check_bounds = true) : Array(Array_), i(i_), j(j_) { if (check_bounds) { @@ -63,9 +63,9 @@ class BArrayCell_const { } }; - + ~BArrayCell_const(){}; - + operator Cell_Type() const; bool operator==(const Cell_Type & val) const; bool operator!=(const Cell_Type & val) const; @@ -73,7 +73,7 @@ class BArrayCell_const { bool operator>(const Cell_Type & val) const; bool operator<=(const Cell_Type & val) const; bool operator>=(const Cell_Type & val) const; - + }; -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraycell-meat.hpp b/include/barry/barraycell-meat.hpp index d4a50db..46c3363 100644 --- a/include/barry/barraycell-meat.hpp +++ b/include/barry/barraycell-meat.hpp @@ -5,7 +5,7 @@ template inline void BArrayCell::operator=(const Cell_Type & val) { - + if (Array->is_empty(i, j, false)) { Array->insert_cell(i, j, val, false, false); } else { @@ -16,7 +16,7 @@ inline void BArrayCell::operator=(const Cell_Type & val) { template inline void BArrayCell::operator+=(const Cell_Type & val) { - + if (Array->is_empty(i, j, false)) { Array->insert_cell(i, j, val, false, false); } else { @@ -27,7 +27,7 @@ inline void BArrayCell::operator+=(const Cell_Type & val) { template inline void BArrayCell::operator-=(const Cell_Type & val) { - + if (Array->is_empty(i, j, false)) { Array->insert_cell(i, j, -val, false, false); } else { @@ -38,7 +38,7 @@ inline void BArrayCell::operator-=(const Cell_Type & val) { template inline void BArrayCell::operator*=(const Cell_Type & val) { - + if (!Array->is_empty(i, j, false)) { Array->el_ij.at(i).at(j).value *= val; } @@ -47,7 +47,7 @@ inline void BArrayCell::operator*=(const Cell_Type & val) { template inline void BArrayCell::operator/=(const Cell_Type & val) { - + if (!Array->is_empty(i, j, false)) { Array->el_ij.at(i).at(j).value /= val; } @@ -61,7 +61,7 @@ inline BArrayCell::operator Cell_Type() const { template inline bool BArrayCell::operator==(const Cell_Type & val) const { - return Array->get_cell(i, j, false) == static_cast(val); + return Array->get_cell(i, j, false) == static_cast(val); } template @@ -71,7 +71,7 @@ inline BArrayCell_const::operator Cell_Type() const { template inline bool BArrayCell_const::operator==(const Cell_Type & val) const { - return Array->get_cell(i, j, false) == static_cast(val); + return Array->get_cell(i, j, false) == static_cast(val); } template @@ -81,22 +81,22 @@ inline bool BArrayCell_const::operator!=(const Cell_Type & template inline bool BArrayCell_const::operator<(const Cell_Type & val) const { - return Array->get_cell(i, j, false) < static_cast(val); + return Array->get_cell(i, j, false) < static_cast(val); } template inline bool BArrayCell_const::operator>(const Cell_Type & val) const { - return Array->get_cell(i, j, false) > static_cast(val); + return Array->get_cell(i, j, false) > static_cast(val); } template inline bool BArrayCell_const::operator<=(const Cell_Type & val) const { - return Array->get_cell(i, j, false) <= static_cast(val); + return Array->get_cell(i, j, false) <= static_cast(val); } template inline bool BArrayCell_const::operator>=(const Cell_Type & val) const { - return Array->get_cell(i, j, false) >= static_cast(val); + return Array->get_cell(i, j, false) >= static_cast(val); } -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraycol-bones.hpp.old b/include/barry/barraycol-bones.hpp.old index 50c0894..bac0a86 100644 --- a/include/barry/barraycol-bones.hpp.old +++ b/include/barry/barraycol-bones.hpp.old @@ -6,14 +6,14 @@ template class BArrayCol { private: - + BArray * Array; std::vector< Cell > iter_vec; size_t i; - + public: - - BArrayCol(BArray * Array_, size_t i_, bool check_bounds = true) : + + BArrayCol(BArray * Array_, size_t i_, bool check_bounds = true) : Array(Array_), i(i_) { if (check_bounds) { @@ -35,19 +35,19 @@ public: Col_type::iterator begin() noexcept; Col_type::iterator end() noexcept; - + }; template class BArrayCol_const { private: - + const BArray * Array; size_t i; - + public: - - BArrayCol(const BArray * Array_, size_t i_, bool check_bounds = true) : + + BArrayCol(const BArray * Array_, size_t i_, bool check_bounds = true) : Array(Array_), i(i_) { if (check_bounds) { @@ -58,7 +58,7 @@ public: }; ~BArrayCol_const(){}; - + // operator std::vector< Cell_Type() > const; bool operator==(const Cell_Type & val) const; bool operator!=(const Cell_Type & val) const; @@ -66,7 +66,7 @@ public: bool operator>(const Cell_Type & val) const; bool operator<=(const Cell_Type & val) const; bool operator>=(const Cell_Type & val) const; - + }; -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraycol-meat.hpp.old b/include/barry/barraycol-meat.hpp.old index 5036455..59b252f 100644 --- a/include/barry/barraycol-meat.hpp.old +++ b/include/barry/barraycol-meat.hpp.old @@ -5,7 +5,7 @@ template inline void BArrayCol::operator=(const Cell_Type & val) { - + if (Array->is_empty(i, j, false)) { Array->insert_cell(i, j, val, false, false); } else { @@ -16,7 +16,7 @@ inline void BArrayCol::operator=(const Cell_Type & val) { template inline void BArrayCol::operator+=(const Cell_Type & val) { - + if (Array->is_empty(i, j, false)) { Array->insert_cell(i, j, val, false, false); } else { @@ -27,7 +27,7 @@ inline void BArrayCol::operator+=(const Cell_Type & val) { template inline void BArrayCol::operator-=(const Cell_Type & val) { - + if (Array->is_empty(i, j, false)) { Array->insert_cell(i, j, -val, false, false); } else { @@ -38,7 +38,7 @@ inline void BArrayCol::operator-=(const Cell_Type & val) { template inline void BArrayCol::operator*=(const Cell_Type & val) { - + if (!Array->is_empty(i, j, false)) { Array->el_ij.at(i).at(j).value *= val; } @@ -47,7 +47,7 @@ inline void BArrayCol::operator*=(const Cell_Type & val) { template inline void BArrayCol::operator/=(const Cell_Type & val) { - + if (!Array->is_empty(i, j, false)) { Array->el_ij.at(i).at(j).value /= val; } @@ -61,7 +61,7 @@ inline BArrayCol::operator Cell_Type() const { template inline bool BArrayCol::operator==(const Cell_Type & val) const { - return Array->get_cell(i, j, false) == static_cast(val); + return Array->get_cell(i, j, false) == static_cast(val); } template @@ -75,7 +75,7 @@ inline Col_type::iterator BArrayCol::end() noexc } /******************************************************************************* - * Const Col + * Const Col * ****************************************************************************/ template @@ -85,7 +85,7 @@ inline BArrayCol_const::operator Cell_Type() const { template inline bool BArrayCol_const::operator==(const Cell_Type & val) const { - return Array->get_cell(i, j, false) == static_cast(val); + return Array->get_cell(i, j, false) == static_cast(val); } template @@ -95,22 +95,22 @@ inline bool BArrayCol_const::operator!=(const Cell_Type & v template inline bool BArrayCol_const::operator<(const Cell_Type & val) const { - return Array->get_cell(i, j, false) < static_cast(val); + return Array->get_cell(i, j, false) < static_cast(val); } template inline bool BArrayCol_const::operator>(const Cell_Type & val) const { - return Array->get_cell(i, j, false) > static_cast(val); + return Array->get_cell(i, j, false) > static_cast(val); } template inline bool BArrayCol_const::operator<=(const Cell_Type & val) const { - return Array->get_cell(i, j, false) <= static_cast(val); + return Array->get_cell(i, j, false) <= static_cast(val); } template inline bool BArrayCol_const::operator>=(const Cell_Type & val) const { - return Array->get_cell(i, j, false) >= static_cast(val); + return Array->get_cell(i, j, false) >= static_cast(val); } -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraydense-bones.hpp b/include/barry/barraydense-bones.hpp index 3da38d4..4344bd3 100644 --- a/include/barry/barraydense-bones.hpp +++ b/include/barry/barraydense-bones.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_BARRAYDENSE_BONES_HPP +#ifndef BARRY_BARRAYDENSE_BONES_HPP #define BARRY_BARRAYDENSE_BONES_HPP 1 template @@ -21,7 +21,7 @@ class BArrayDenseCell_const; /** * @brief Baseline class for binary arrays. - * + * * `BArrayDense` class objects are arbitrary dense-arrays. The data * is stored internally in the `el` member, which can be accessed * using the member function `get_data()`, by column. @@ -52,8 +52,8 @@ class BArrayDense { static const bool dense = true; public: - - /** + + /** * This is as a reference, if we need to iterate through the cells and we need * to keep track which were visited, we use this as a reference. So that if * cell.visited = true and visited = true, it means that we haven't been here @@ -61,11 +61,11 @@ class BArrayDense { * beginning of the routine. */ bool visited = false; - + /** * @name Constructors - * + * * @param N_ Number of rows * @param M_ Number of columns * @param source An unsigned vector ranging from 0 to N_ @@ -74,15 +74,15 @@ class BArrayDense { * @param value Cell_Type defaul fill-in value (zero, by default.) */ ///@{ - + /** @brief Zero-size array */ BArrayDense() : N(0u), M(0u), el(0u), el_rowsums(0u), el_colsums(0u) {}; - + /** @brief Empty array */ BArrayDense (size_t N_, size_t M_, Cell_Type value = static_cast(0)) : N(N_), M(M_), el(N_ * M_, value), el_rowsums(N_, static_cast(value * M_)), el_colsums(M_, static_cast(value * N_)) {}; - + /** @brief Edgelist with data */ BArrayDense ( size_t N_, @@ -92,7 +92,7 @@ class BArrayDense { const std::vector< Cell_Type > & value, bool add = true ); - + /** @brief Edgelist with no data (simpler) */ BArrayDense ( size_t N_, size_t M_, @@ -100,10 +100,10 @@ class BArrayDense { const std::vector< size_t > & target, bool add = true ); - + /** @brief Copy constructor */ BArrayDense(const BArrayDense & Array_, bool copy_data = false); - + /** @brief Assignment constructor */ BArrayDense & operator=(const BArrayDense & Array_); @@ -113,20 +113,20 @@ class BArrayDense { /** @brief Move assignment */ BArrayDense & operator=(BArrayDense && x) noexcept; ///@} - + bool operator==(const BArrayDense & Array_); ~BArrayDense(); - + // In principle, copy can be faster by using openmp on the rows // since those are independent. // BArrayDense(BArrayDense & A); - + /** * @brief Set the data object - * - * @param data_ - * @param delete_data_ + * + * @param data_ + * @param delete_data_ */ ///@{ void set_data(Data_Type * data_, bool delete_data_ = false); @@ -135,16 +135,16 @@ class BArrayDense { Data_Type & D(); const Data_Type & D() const; ///@} - + // Function to access the elements // bool check_cell void out_of_range(size_t i, size_t j) const; - Cell_Type get_cell(size_t i, size_t j, bool check_bounds = true) const; + Cell_Type get_cell(size_t i, size_t j, bool check_bounds = true) const; std::vector< Cell_Type > get_col_vec(size_t i, bool check_bounds = true) const; std::vector< Cell_Type > get_row_vec(size_t i, bool check_bounds = true) const; void get_col_vec(std::vector< Cell_Type > * x, size_t i, bool check_bounds = true) const; void get_row_vec(std::vector< Cell_Type > * x, size_t i, bool check_bounds = true) const; - + BArrayDenseRow & row(size_t i, bool check_bounds = true); const BArrayDenseRow_const row(size_t i, bool check_bounds = true) const; @@ -153,12 +153,12 @@ class BArrayDense { /** * @brief Get the edgelist - * + * * `Entries` is a class with three objects: Two `std::vector` with the row and * column coordinates respectively, and one `std::vector` with the corresponding * value of the cell. - * - * @return Entries + * + * @return Entries */ Entries get_entries() const; @@ -186,37 +186,37 @@ class BArrayDense { * delete/add), or, in the case of `swap_cells`, check if either of both * cells exists/don't exist. */ - ///@{ + ///@{ BArrayDense & operator+=(const std::pair & coords); BArrayDense & operator-=(const std::pair & coords); BArrayDenseCell operator()(size_t i, size_t j, bool check_bounds = true); const Cell_Type operator()(size_t i, size_t j, bool check_bounds = true) const; - + void rm_cell(size_t i, size_t j, bool check_bounds = true, bool check_exists = true); - + void insert_cell(size_t i, size_t j, const Cell< Cell_Type > & v, bool check_bounds, bool); // void insert_cell(size_t i, size_t j, Cell< Cell_Type > && v, bool check_bounds, bool check_exists); void insert_cell(size_t i, size_t j, Cell_Type v, bool check_bounds, bool); - + void swap_cells( size_t i0, size_t j0, size_t i1, size_t j1, bool check_bounds = true, int check_exists = CHECK::BOTH, int * report = nullptr ); - + void toggle_cell(size_t i, size_t j, bool check_bounds = true, int check_exists = EXISTS::UKNOWN); void toggle_lock(size_t i, size_t j, bool check_bounds = true); ///@} - + /**@name Column/row wise interchange*/ ///@{ void swap_rows(size_t i0, size_t i1, bool check_bounds = true); void swap_cols(size_t j0, size_t j1, bool check_bounds = true); - + void zero_row(size_t i, bool check_bounds = true); void zero_col(size_t j, bool check_bounds = true); ///@} - + void transpose(); void clear(bool hard = true); void resize(size_t N_, size_t M_); @@ -224,13 +224,13 @@ class BArrayDense { // Advances operators // void toggle_iterator - + // Misc void print(const char * fmt = nullptr, ...) const; /** * @name Arithmetic operators - * + * */ ///@{ BArrayDense& operator+=(const BArrayDense& rhs); @@ -238,11 +238,11 @@ class BArrayDense { BArrayDense& operator-=(const BArrayDense& rhs); BArrayDense& operator-=(const Cell_Type & rhs); - + BArrayDense& operator/=(const Cell_Type & rhs); BArrayDense& operator*=(const Cell_Type & rhs); ///@} - + // /** // * @name Casting between types // */ @@ -252,7 +252,7 @@ class BArrayDense { // operator BArrayDense() const; // operator BArrayDense() const; // ///@} - + bool is_dense() const noexcept {return dense;}; const std::vector< Cell_Type > & get_data() const; diff --git a/include/barry/barraydense-bones.hpp.old b/include/barry/barraydense-bones.hpp.old index 12f6029..c87bc80 100644 --- a/include/barry/barraydense-bones.hpp.old +++ b/include/barry/barraydense-bones.hpp.old @@ -4,10 +4,10 @@ /** * @brief Dense bi-dimensional array - * + * * @details elements is stored in a std::vector, in col-major order. - * - * @tparam Cell_Type + * + * @tparam Cell_Type */ template class BArrayDense { @@ -26,7 +26,7 @@ public: ~BArrayDense() {}; void fill(const Cell_Type & d); - + const std::vector< Cell_Type > & elements_raw() const noexcept; const std::vector< Cell_Type > * elements_ptr() const noexcept; size_t nrow() const noexcept; @@ -63,7 +63,7 @@ inline BArrayDense::BArrayDense( elements.resize(elements_.size()); std::copy(elements_.begin(), elements_.end(), elements.begin()); - + } template @@ -79,7 +79,7 @@ inline const std::vector< Cell_Type > & { return elements; } - + template inline const std::vector< Cell_Type > * BArrayDense::elements_ptr() const noexcept @@ -102,7 +102,7 @@ inline size_t BArrayDense::ncol() const noexcept template inline Cell_Type BArrayDense::operator()( size_t i, size_t j, bool check_bounds -) const +) const { if (check_bounds) { @@ -119,8 +119,8 @@ inline Cell_Type BArrayDense::operator()( } template -inline Cell_Type BArrayDense::operator[](size_t i) const -{ +inline Cell_Type BArrayDense::operator[](size_t i) const +{ return elements[i]; } @@ -139,4 +139,4 @@ inline void BArrayDense::print() const { printf_barry("\n"); } -} \ No newline at end of file +} diff --git a/include/barry/barraydense-meat-operators.hpp b/include/barry/barraydense-meat-operators.hpp index 287e12e..59b6ac6 100644 --- a/include/barry/barraydense-meat-operators.hpp +++ b/include/barry/barraydense-meat-operators.hpp @@ -37,7 +37,7 @@ BDENSE_TEMPLATE(BDENSE_TYPE()&, operator+=) ( // Must be compatible checkdim_(*this, rhs); - + for (size_t i = 0u; i < nrow(); ++i) for (size_t j = 0u; j < ncol(); ++j) this->operator()(i, j) += rhs.get_cell(i, j); @@ -64,7 +64,7 @@ BDENSE_TEMPLATE(BDENSE_TYPE()&, operator-=) ( // Must be compatible checkdim_(*this, rhs); - + for (size_t i = 0u; i < nrow(); ++i) { for (size_t j = 0u; j < ncol(); ++j) { this->operator()(i, j) -= rhs.get_cell(i, j); @@ -78,11 +78,11 @@ BDENSE_TEMPLATE(BDENSE_TYPE()&, operator-=) ( const Cell_Type& rhs ) { - for (size_t i = 0u; i < nrow(); ++i) - for (size_t j = 0u; j < ncol(); ++j) + for (size_t i = 0u; i < nrow(); ++i) + for (size_t j = 0u; j < ncol(); ++j) this->operator()(i, j) -= rhs; - - + + return *this; } @@ -91,7 +91,7 @@ BDENSE_TEMPLATE(BDENSE_TYPE()&, operator*=) ( const Cell_Type& rhs ) { - for (size_t i = 0u; i < nrow(); ++i) + for (size_t i = 0u; i < nrow(); ++i) for (size_t j = 0u; j < nrow(); ++j) el[POS(i, j)] *= rhs; @@ -102,7 +102,7 @@ BDENSE_TEMPLATE(BDENSE_TYPE()&, operator/=) ( const Cell_Type& rhs ) { - for (size_t i = 0u; i < nrow(); ++i) + for (size_t i = 0u; i < nrow(); ++i) for (size_t j = 0u; j < nrow(); ++j) el[POS(i, j)] /= rhs; @@ -118,4 +118,4 @@ BDENSE_TEMPLATE(BDENSE_TYPE()&, operator/=) ( #undef POS #undef POS_N -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraydense-meat.hpp b/include/barry/barraydense-meat.hpp index 3cd2a09..ad44fb7 100644 --- a/include/barry/barraydense-meat.hpp +++ b/include/barry/barraydense-meat.hpp @@ -2,7 +2,7 @@ // #include "barraydense-bones.hpp" #ifndef BARRY_BARRAYDENSE_MEAT_HPP -#define BARRY_BARRAYDENSE_MEAT_HPP +#define BARRY_BARRAYDENSE_MEAT_HPP template class BArrayDenseRow; @@ -26,7 +26,7 @@ class BArrayDenseCell; #define POS_N(a,b,c) (b)*(c) + (a) template -Cell_Type BArrayDense::Cell_default = static_cast< Cell_Type >(1.0); +Cell_Type BArrayDense::Cell_default = static_cast< Cell_Type >(1.0); #define ZERO_CELL static_cast(0.0) @@ -40,45 +40,45 @@ inline BArrayDense::BArrayDense( const std::vector< Cell_Type > & value, bool add ) : N(N_), M(M_), el(N_ * M_, ZERO_CELL), el_rowsums(N_, ZERO_CELL), el_colsums(M_, ZERO_CELL) { - + if (source.size() != target.size()) throw std::length_error("-source- and -target- don't match on length."); if (source.size() != value.size()) throw std::length_error("-sorce- and -value- don't match on length."); - + // Writing the data for (size_t i = 0u; i < source.size(); ++i) { - + // Checking range bool empty = is_empty(source[i], target[i], true); if (add && !empty) { Cell_Type tmp = el[POS(source[i], target[i])]; - + el_rowsums[source[i]] += (value[i] - tmp); el_colsums[target[i]] += (value[i] - tmp); el[POS(source[i], target[i])] += value[i]; - + continue; - } - + } + if (!empty) throw std::logic_error("The value already exists. Use 'add = true'."); - + el[POS(source[i], target[i])] = value[i]; el_rowsums[source[i]] += value[i]; el_colsums[target[i]] += value[i]; - + } - + return; - + } // Edgelist without data @@ -89,7 +89,7 @@ inline BArrayDense:: BArrayDense( const std::vector< size_t > & target, bool add ) : N(N_), M(M_), el(N_ * M_, ZERO_CELL), el_rowsums(N_, ZERO_CELL), el_colsums(M_, ZERO_CELL) { - + std::vector< Cell_Type > value(source.size(), static_cast(1.0)); if (source.size() != target.size()) @@ -97,38 +97,38 @@ inline BArrayDense:: BArrayDense( if (source.size() != value.size()) throw std::length_error("-sorce- and -value- don't match on length."); - + // Writing the data for (size_t i = 0u; i < source.size(); ++i) { - + // Checking range bool empty = is_empty(source[i], target[i], true); if (add && !empty) { Cell_Type tmp = el[POS(source[i], target[i])]; - + el_rowsums[source[i]] += (value[i] - tmp); el_colsums[target[i]] += (value[i] - tmp); el[POS(source[i], target[i])] += value[i]; - + continue; - } - + } + if (!empty) throw std::logic_error("The value already exists. Use 'add = true'."); - + el[POS(source[i], target[i])] = value[i]; el_rowsums[source[i]] += value[i]; el_colsums[target[i]] += value[i]; - + } - + } template @@ -136,7 +136,7 @@ inline BArrayDense:: BArrayDense( const BArrayDense & Array_, bool copy_data ) : N(Array_.N), M(Array_.M){ - + // Dimensions el = Array_.el; el_rowsums = Array_.el_rowsums; @@ -144,14 +144,14 @@ inline BArrayDense:: BArrayDense( // el.resize(0u); // el_rowsums.resize(0u); // el_colsums.resize(0u); - + // std::copy(Array_.el.begin(), Array_.el.end(), std::back_inserter(el)); // std::copy(Array_.el_rowsums.begin(), Array_.el_rowsums.end(), std::back_inserter(el_rowsums)); // std::copy(Array_.el_colsums.begin(), Array_.el_colsums.end(), std::back_inserter(el_colsums)); // this->NCells = Array_.NCells; this->visited = Array_.visited; - + // Data if (Array_.data != nullptr) { @@ -170,27 +170,27 @@ inline BArrayDense:: BArrayDense( } } - + return; - + } template inline BArrayDense & BArrayDense::operator=( const BArrayDense & Array_ ) { - + // Clearing if (this != &Array_) { - + el = Array_.el; el_rowsums = Array_.el_rowsums; el_colsums = Array_.el_colsums; // el.resize(0u); // el_rowsums.resize(0u); // el_colsums.resize(0u); - + // // Entries // std::copy(Array_.el.begin(), Array_.el.end(), std::back_inserter(el)); // std::copy(Array_.el_rowsums.begin(), Array_.el_rowsums.end(), std::back_inserter(el_rowsums)); @@ -200,7 +200,7 @@ inline BArrayDense & BArrayDense::ope // this->NCells = Array_.NCells; this->N = Array_.N; this->M = Array_.M; - + // Data if (data != nullptr) { @@ -218,11 +218,11 @@ inline BArrayDense & BArrayDense::ope delete_data = true; } - + } - + return *this; - + } template @@ -247,19 +247,19 @@ template inline BArrayDense & BArrayDense::operator=( BArrayDense && x ) noexcept { - + // Clearing if (this != &x) { - + N = x.N; M = x.M; // NCells = x.NCells; - + std::swap(el, x.el); std::swap(el_rowsums, x.el_rowsums); std::swap(el_colsums, x.el_colsums); - + // Data if (data != nullptr) { @@ -280,38 +280,38 @@ inline BArrayDense & BArrayDense::ope x.data = nullptr; } - + } - + return *this; - + } template inline bool BArrayDense::operator== ( const BArrayDense & Array_ ) { - + // Dimension and number of cells used if ( (N != Array_.nrow()) | (M != Array_.ncol()) ) return false; - + // One holds, and the other doesn't. if ((!data & Array_.data) | (data & !Array_.data)) return false; - + if (this->el != Array_.el) return false; - + return true; } template inline BArrayDense::~BArrayDense () { - + if (delete_data && (data != nullptr)) delete data; - + return; } @@ -319,16 +319,16 @@ template inline void BArrayDense::set_data ( Data_Type * data_, bool delete_data_ -) { +) { if ((data != nullptr) && delete_data) delete data; - + data = data_; delete_data = delete_data_; - + return; - + } template @@ -371,20 +371,20 @@ inline void BArrayDense::out_of_range ( return; } - + template inline Cell_Type BArrayDense::get_cell ( size_t i, size_t j, bool check_bounds ) const { - - // Checking boundaries + + // Checking boundaries if (check_bounds) out_of_range(i,j); - + return el[POS(i, j)]; - + } template @@ -393,15 +393,15 @@ inline std::vector< Cell_Type > BArrayDense::get_row_vec ( bool check_bounds ) const { - // Checking boundaries - if (check_bounds) + // Checking boundaries + if (check_bounds) out_of_range(i, 0u); std::vector< Cell_Type > ans; ans.reserve(ncol()); - for (size_t j = 0u; j < M; ++j) + for (size_t j = 0u; j < M; ++j) ans.push_back(el[POS(i, j)]); - + return ans; } @@ -412,13 +412,13 @@ template inline void BArrayDenseoperator[](j) = el[POS(i, j)]; - + } template inline std::vector< Cell_Type > BArrayDense:: get_col_vec( @@ -426,15 +426,15 @@ template inline std::vector< Cell_Type > bool check_bounds ) const { - // Checking boundaries - if (check_bounds) + // Checking boundaries + if (check_bounds) out_of_range(0u, i); std::vector< Cell_Type > ans; - ans.reserve(nrow()); - for (size_t j = 0u; j < N; ++j) + ans.reserve(nrow()); + for (size_t j = 0u; j < N; ++j) ans.push_back(el[POS(j, i)]); - + return ans; } @@ -445,8 +445,8 @@ template inline void BArrayDense inline void BArrayDenseoperator[](j) = el[POS(j, i)];//this->get_cell(iter->first, i, false); - + } template inline const BArrayDenseRow_const BArrayDense::row( @@ -486,7 +486,7 @@ inline BArrayDenseRow & BArrayDense:: } template -inline const BArrayDenseCol_const +inline const BArrayDenseCol_const BArrayDense::col( size_t j, bool check_bounds @@ -500,7 +500,7 @@ BArrayDense::col( } template -inline BArrayDenseCol & +inline BArrayDenseCol & BArrayDense::col( size_t j, bool check_bounds @@ -514,11 +514,11 @@ BArrayDense::col( } template inline Entries< Cell_Type > BArrayDense:: get_entries() const { - + size_t nzero = this->nnozero(); Entries res(nzero); - + for (size_t i = 0u; i < N; ++i) { for (size_t j = 0u; j < M; ++j) @@ -532,12 +532,12 @@ template inline Entries< Cell_Type > BAr res.val.push_back(el[POS(i, j)]); } - + } } - + return res; } @@ -547,12 +547,12 @@ template inline bool BArrayDense inline size_t BArrayDense:: nrow() const noexcept { @@ -582,7 +582,7 @@ template inline BArrayDense & BArrayDense::operator+=( const std::pair & coords ) { - + size_t i = coords.first; size_t j = coords.second; @@ -592,16 +592,16 @@ inline BArrayDense & BArrayDense::ope el[POS(i,j)] += 1; el_rowsums[i] += 1; el_colsums[j] += 1; - + return *this; - + } template inline BArrayDense & BArrayDense::operator-=( const std::pair & coords ) { - + size_t i = coords.first; size_t j = coords.second; @@ -612,24 +612,24 @@ inline BArrayDense & BArrayDense::ope el[POS(i,j)] = ZERO_CELL; el_rowsums[i] -= old; el_colsums[j] -= old; - + return *this; - + } template -inline BArrayDenseCell BArrayDense::operator()( +inline BArrayDenseCell BArrayDense::operator()( size_t i, size_t j, bool check_bounds ) { - + return BArrayDenseCell(this, i, j, check_bounds); - + } template -inline const Cell_Type BArrayDense::operator()( +inline const Cell_Type BArrayDense::operator()( size_t i, size_t j, bool check_bounds @@ -637,9 +637,9 @@ inline const Cell_Type BArrayDense::operator()( if (check_bounds) out_of_range(i, j); - + return el[POS(i,j)]; - + } template @@ -649,18 +649,18 @@ inline void BArrayDense::rm_cell ( bool check_bounds, bool check_exists ) { - + // Checking the boundaries if (check_bounds) out_of_range(i,j); // BARRY_UNUSED(check_exists) - + // Remove the pointer first (so it wont point to empty) el_rowsums[i] -= el[POS(i, j)]; - el_colsums[j] -= el[POS(i, j)]; + el_colsums[j] -= el[POS(i, j)]; el[POS(i, j)] = BARRY_ZERO_DENSE; - + return; } @@ -672,18 +672,18 @@ inline void BArrayDense::insert_cell ( const Cell< Cell_Type> & v, bool check_bounds, bool -) { - +) { + if (check_bounds) - out_of_range(i,j); + out_of_range(i,j); if (el[POS(i,j)] == BARRY_ZERO_DENSE) { el_rowsums[i] += v.value; el_colsums[j] += v.value; - - } + + } else { @@ -697,7 +697,7 @@ inline void BArrayDense::insert_cell ( return; - + } template inline void BArrayDense:: insert_cell( @@ -707,17 +707,17 @@ template inline void BArrayDense inline void BArrayDense inline void BArrayDense inline void BArrayDense:: swap_rows ( @@ -792,7 +792,7 @@ template inline void BArrayDense inline void BArrayDense inline void BArrayDense inline void BArrayDense inline void BArrayDense inline void BArrayDense inline void BArrayDense:: zero_col ( size_t j, bool check_bounds ) { - + if (check_bounds) out_of_range(0u, j); - + if (el_colsums[j] == ZERO_CELL) return; - + // Else, remove all elements - for (size_t row = 0u; row < N; row++) + for (size_t row = 0u; row < N; row++) rm_cell(row, j, false, false); - + return; - + } template inline void BArrayDense:: transpose () { - + // if (NCells == 0u) // { @@ -890,20 +890,20 @@ template inline void BArrayDense > tmp_el(std::move(el)); el.resize(N * M, ZERO_CELL); - for (size_t i = 0u; i < N; ++i) + for (size_t i = 0u; i < N; ++i) for (size_t j = 0u; j < M; ++j) std::swap(tmp_el[POS(i, j)], el[POS_N(j, i, M)]); - + // Swapping the values std::swap(N, M); std::swap(el_rowsums, el_colsums); - + return; } @@ -911,15 +911,15 @@ template inline void BArrayDense inline void BArrayDense:: clear ( bool hard ) { - + BARRY_UNUSED(hard) - + std::fill(el.begin(), el.end(), ZERO_CELL); std::fill(el_rowsums.begin(), el_rowsums.end(), ZERO_CELL); std::fill(el_colsums.begin(), el_colsums.end(), ZERO_CELL); - + return; - + } template inline void BArrayDense:: resize ( @@ -953,7 +953,7 @@ template inline void BArrayDense inline void BArrayDense inline void BArrayDense:: print ( @@ -972,7 +972,7 @@ template inline void BArrayDense inline void BArrayDenseis_empty(i, j, false)) printf_barry(" . "); - else + else printf_barry(" %.2f ", static_cast(this->get_cell(i, j, false))); - + } printf_barry("\n"); } - + return; - + } template inline const std::vector< Cell_Type > & BArrayDense:: get_data() const @@ -1024,4 +1024,3 @@ template inline const Cell_Type BArrayDe #undef ZERO_CELL #endif - diff --git a/include/barry/barraydensecell-bones.hpp b/include/barry/barraydensecell-bones.hpp index f59ccb7..050fa75 100644 --- a/include/barry/barraydensecell-bones.hpp +++ b/include/barry/barraydensecell-bones.hpp @@ -20,19 +20,19 @@ class BArrayDenseCell { friend class BArrayDenseCol; friend class BArrayDenseCol_const; private: - + BArrayDense * dat; size_t i; size_t j; - + public: - + BArrayDenseCell( BArrayDense * Array_, size_t i_, size_t j_, bool check_bounds = true - ) : + ) : i(i_), j(j_) { @@ -62,10 +62,10 @@ class BArrayDenseCell { operator Cell_Type() const; bool operator==(const Cell_Type & val) const; - + }; #undef POS -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraydensecell-meat.hpp b/include/barry/barraydensecell-meat.hpp index 137f31d..14abd5b 100644 --- a/include/barry/barraydensecell-meat.hpp +++ b/include/barry/barraydensecell-meat.hpp @@ -3,13 +3,13 @@ #ifndef BARRY_BARRAYDENSECELL_MEAT_HPP #define BARRY_BARRAYDENSECELL_MEAT_HPP 1 -#define POS(a, b) (a) + (b) * dat->N +#define POS(a, b) (a) + (b) * dat->N template inline BArrayDenseCell& BArrayDenseCell::operator=( const BArrayDenseCell & other ) { - + Cell_Type val = static_cast(other); #ifdef BARRY_DEBUG Cell_Type old = dat->el.at(POS(i,j)); @@ -41,12 +41,12 @@ inline void BArrayDenseCell::operator=(const Cell_Type & va dat->el_rowsums[i] += (val - old); dat->el_colsums[j] += (val - old); #endif - + } template inline void BArrayDenseCell::operator+=(const Cell_Type & val) { - + #ifdef BARRY_DEBUG dat->el.at(POS(i,j)) += val; dat->el_rowsums.at(i) += val; @@ -61,7 +61,7 @@ inline void BArrayDenseCell::operator+=(const Cell_Type & v template inline void BArrayDenseCell::operator-=(const Cell_Type & val) { - + #ifdef BARRY_DEBUG dat->el.at(POS(i,j)) -= val; dat->el_rowsums.at(i) -= val; @@ -76,7 +76,7 @@ inline void BArrayDenseCell::operator-=(const Cell_Type & v template inline void BArrayDenseCell::operator*=(const Cell_Type & val) { - + #ifdef BARRY_DEBUG Cell_Type old = dat->el.at(POS(i,j)); dat->el_colsums.at(j) += (old * val - old); @@ -93,7 +93,7 @@ inline void BArrayDenseCell::operator*=(const Cell_Type & v template inline void BArrayDenseCell::operator/=(const Cell_Type & val) { - + #ifdef BARRY_DEBUG Cell_Type old = dat->el.at(POS(i,j)); dat->el_rowsums.at(i) += (old/val - old); @@ -115,9 +115,9 @@ inline BArrayDenseCell::operator Cell_Type() const { template inline bool BArrayDenseCell::operator==(const Cell_Type & val) const { - return dat->el[POS(i,j)] == val; + return dat->el[POS(i,j)] == val; } #undef POS -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraydensecol-bones.hpp b/include/barry/barraydensecol-bones.hpp index f0476ec..e5b41a8 100644 --- a/include/barry/barraydensecol-bones.hpp +++ b/include/barry/barraydensecol-bones.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_BARRAYDENSECOL_BONES +#ifndef BARRY_BARRAYDENSECOL_BONES #define BARRY_BARRAYDENSECOL_BONES #define POS(a,b) (b)*N + (a) @@ -23,14 +23,14 @@ class BArrayDenseCol { for (size_t i = 0u; i < array->N; ++i) { - + if (array->el[POS_N(i, index, array->N)] != ZERO_CELL) col[i] = col[POS_N(i, index, array->N)]; - + } col_filled = true; - + } } @@ -85,10 +85,10 @@ class BArrayDenseCol_const { for (size_t i = 0u; i < array->N; ++i) { - + if (array->el[POS_N(i, index, array->N)] != ZERO_CELL) col[i] = col[POS_N(i, index, array->N)]; - + } }; @@ -120,4 +120,4 @@ class BArrayDenseCol_const { #undef POS_N #undef ZERO_CELL -#endif \ No newline at end of file +#endif diff --git a/include/barry/barraydenserow-bones.hpp b/include/barry/barraydenserow-bones.hpp index 1c48c2f..620fcb1 100644 --- a/include/barry/barraydenserow-bones.hpp +++ b/include/barry/barraydenserow-bones.hpp @@ -23,14 +23,14 @@ class BArrayDenseRow { for (size_t j = 0u; j < array->M; ++j) { - + if (array->el[POS_N(index, j, array->N)] != ZERO_CELL) row[j] = row[POS_N(index, j, array->N)]; - + } row_filled = true; - + } } @@ -94,10 +94,10 @@ class BArrayDenseRow_const { for (size_t j = 0u; j < array->M; ++j) { - + if (array->el[POS_N(index, j, array->M)] != ZERO_CELL) row[j] = row[POS_N(index, j, array->M)]; - + } return; @@ -131,4 +131,4 @@ class BArrayDenseRow_const { #undef POS_N #undef ZERO_CELL -#endif \ No newline at end of file +#endif diff --git a/include/barry/barrayrow-bones.hpp b/include/barry/barrayrow-bones.hpp index f193357..20d420d 100644 --- a/include/barry/barrayrow-bones.hpp +++ b/include/barry/barrayrow-bones.hpp @@ -4,13 +4,13 @@ template class BArrayRow { private: - + BArray * Array; size_t i; - + public: - - BArrayRow(BArray * Array_, size_t i_,, bool check_bounds = true) : + + BArrayRow(BArray * Array_, size_t i_,, bool check_bounds = true) : Array(Array_), i(i_), j(j_) { if (check_bounds) @@ -32,7 +32,7 @@ class BArrayRow { operator BArrayRow() const; bool operator==(const BArrayRow & val) const; - + }; @@ -40,13 +40,13 @@ class BArrayRow { template class BArrayRow_const { private: - + const BArray * Array; size_t i; - + public: - - BArrayRow_const(const BArray * Array_, size_t i_, bool check_bounds = true) : + + BArrayRow_const(const BArray * Array_, size_t i_, bool check_bounds = true) : Array(Array_), i(i_), { if (check_bounds) { @@ -55,9 +55,9 @@ class BArrayRow_const { } }; - + ~BArrayRow_const(){}; - + operator BArrayRow_const() const; bool operator==(const BArrayRow_const & val) const; bool operator!=(const BArrayRow_const & val) const; @@ -65,7 +65,7 @@ class BArrayRow_const { bool operator>(const BArrayRow_const & val) const; bool operator<=(const BArrayRow_const & val) const; bool operator>=(const BArrayRow_const & val) const; - + }; -#endif \ No newline at end of file +#endif diff --git a/include/barry/barrayrow-meat.hpp b/include/barry/barrayrow-meat.hpp index d9e3f2a..33ba9ee 100644 --- a/include/barry/barrayrow-meat.hpp +++ b/include/barry/barrayrow-meat.hpp @@ -9,7 +9,7 @@ template BROW_TEMPLATE_ARGS() inline a BROW_TYPE()::b BROW_TEMPLATE(void, operator=) (const BROW_TYPE() & val) { - + // First, zeroout the row this->Array->zero_row(j); @@ -23,10 +23,10 @@ BROW_TEMPLATE(void, operator=) (const BROW_TYPE() & val) { } BROW_TEMPLATE(void, operator+=) (const BROW_TYPE() & val) { - + for (auto & v : val) this->Array->operator(i, v.first) += v.second; - + return; } @@ -34,10 +34,10 @@ BROW_TEMPLATE(void, operator+=) (const BROW_TYPE() & val) { BROW_TEMPLATE(void, operator-=) ( const BROW_TYPE() & val ) { - + for (auto & v : val) this->Array->operator(i, v.first) -= v.second; - + return; } @@ -45,7 +45,7 @@ BROW_TEMPLATE(void, operator-=) ( BROW_TEMPLATE(void, operator*=) ( const BROW_TYPE() & val ) { - + if (!Array->is_empty(i, j, false)) { Array->el_ij.at(i).at(j).value *= val; } @@ -55,7 +55,7 @@ BROW_TEMPLATE(void, operator*=) ( BROW_TEMPLATE(void, operator/=) ( const BROW_TYPE() & val ) { - + if (!Array->is_empty(i, j, false)) { Array->el_ij.at(i).at(j).value /= val; } @@ -69,7 +69,7 @@ inline BArrayCell::operator Cell_Type() const { template inline bool BArrayCell::operator==(const Cell_Type & val) const { - return Array->get_cell(i, j, false) == static_cast(val); + return Array->get_cell(i, j, false) == static_cast(val); } template @@ -79,7 +79,7 @@ inline BArrayCell_const::operator Cell_Type() const { template inline bool BArrayCell_const::operator==(const Cell_Type & val) const { - return Array->get_cell(i, j, false) == static_cast(val); + return Array->get_cell(i, j, false) == static_cast(val); } template @@ -89,26 +89,26 @@ inline bool BArrayCell_const::operator!=(const Cell_Type & template inline bool BArrayCell_const::operator<(const Cell_Type & val) const { - return Array->get_cell(i, j, false) < static_cast(val); + return Array->get_cell(i, j, false) < static_cast(val); } template inline bool BArrayCell_const::operator>(const Cell_Type & val) const { - return Array->get_cell(i, j, false) > static_cast(val); + return Array->get_cell(i, j, false) > static_cast(val); } template inline bool BArrayCell_const::operator<=(const Cell_Type & val) const { - return Array->get_cell(i, j, false) <= static_cast(val); + return Array->get_cell(i, j, false) <= static_cast(val); } template inline bool BArrayCell_const::operator>=(const Cell_Type & val) const { - return Array->get_cell(i, j, false) >= static_cast(val); + return Array->get_cell(i, j, false) >= static_cast(val); } #undef BROW_TYPE #undef BROW_TEMPLATE_ARGS #undef BROW_TEMPLATE -#endif \ No newline at end of file +#endif diff --git a/include/barry/barrayvector-bones.hpp b/include/barry/barrayvector-bones.hpp index 61f2c60..b3781a9 100644 --- a/include/barry/barrayvector-bones.hpp +++ b/include/barry/barrayvector-bones.hpp @@ -3,14 +3,14 @@ /** * @brief Row or column of a `BArray` - * - * @tparam Cell_Type - * @tparam Data_Type + * + * @tparam Cell_Type + * @tparam Data_Type */ template class BArrayVector { private: - + BArray * Array; std::vector< std::pair< size_t, Cell_Type > > vec; size_t dim; @@ -18,12 +18,12 @@ class BArrayVector { void init_vec(); bool vec_initialized = false; - + public: - + /** * @brief Construct a new BArrayVector object - * + * * @param Array_ Pointer to a `BArray` object * @param dim_ Dimension. 0 means row and 1 means column. * @param i_ Element to point. @@ -34,7 +34,7 @@ class BArrayVector { size_t & dim_ size_t & i_, bool check_bounds = true - ) : + ) : Array(Array_), vec(0u), dim(dim_), i(i_) { if (dim > 1u) @@ -66,13 +66,13 @@ class BArrayVector { operator std::vector< Cell_Type >() const; bool operator==(const Cell_Type & val) const; - + }; template class BArrayVector_const { private: - + const BArray * Array; std::vector< std::pair< size_t, Cell_Type > > vec; size_t dim; @@ -80,15 +80,15 @@ class BArrayVector_const { void init_vec(); bool vec_initialized = false; - + public: - + BArrayVector_const( const BArray * Array_, size_t & dim_ size_t & i_, bool check_bounds = true - ) : + ) : Array(Array_), vec(0u), dim(dim_), i(i_) { if (dim > 1u) @@ -104,7 +104,7 @@ class BArrayVector_const { } }; - + ~BArrayVector_const() {}; bool is_row() const noexcept; @@ -112,7 +112,7 @@ class BArrayVector_const { size_t size() const noexcept; std::vector< Cell_Type >::const_iterator begin() noexcept; std::vector< Cell_Type >::const_iterator end() noexcept; - + operator std::vector() const; bool operator==(const Cell_Type & val) const; bool operator!=(const Cell_Type & val) const; @@ -120,7 +120,7 @@ class BArrayVector_const { bool operator>(const Cell_Type & val) const; bool operator<=(const Cell_Type & val) const; bool operator>=(const Cell_Type & val) const; - + }; -#endif \ No newline at end of file +#endif diff --git a/include/barry/barrayvector-meat.hpp b/include/barry/barrayvector-meat.hpp index 29528bb..74e594a 100644 --- a/include/barry/barrayvector-meat.hpp +++ b/include/barry/barrayvector-meat.hpp @@ -12,7 +12,7 @@ inline void BArrayVector::init_vec() { for (const auto& a : Array->el_ij[i]) vec.push_back(a); - + } else { for (const auto& a : Array->el_ji[i]) @@ -42,13 +42,13 @@ inline size_t BArrayVector::size() const noexcept { return Array->el_ij[i].size(); else return Array->el_ji[i].size(); - + } template inline std::vector< Cell_Type >::const_iterator BArrayVector::begin() noexcept { - + // For this, we will need the iterator init_vec(); @@ -67,10 +67,10 @@ inline std::vector< Cell_Type >::const_iterator BArrayVector inline void BArrayVector::operator=(const Cell_Type & val) { - + size_t k = 0u; size_t N_ = (dim == 0u) ? Array->nrow() : Array->ncol(); - + if (dim == 0u) { @@ -89,10 +89,10 @@ inline void BArrayVector::operator=(const Cell_Type & val) template inline void BArrayVector::operator+=(const Cell_Type & val) { - + size_t k = 0u; size_t N_ = (dim == 0u) ? Array->nrow() : Array->ncol(); - + if (dim == 0u) { @@ -110,10 +110,10 @@ inline void BArrayVector::operator+=(const Cell_Type & val) template inline void BArrayVector::operator-=(const Cell_Type & val) { - + size_t k = 0u; size_t N_ = (dim == 0u) ? Array->nrow() : Array->ncol(); - + if (dim == 0u) { @@ -131,10 +131,10 @@ inline void BArrayVector::operator-=(const Cell_Type & val) template inline void BArrayVector::operator*=(const Cell_Type & val) { - + size_t k = 0u; size_t N_ = (dim == 0u) ? Array->nrow() : Array->ncol(); - + if (dim == 0u) { @@ -152,10 +152,10 @@ inline void BArrayVector::operator*=(const Cell_Type & val) template inline void BArrayVector::operator/=(const Cell_Type & val) { - + size_t k = 0u; size_t N_ = (dim == 0u) ? Array->nrow() : Array->ncol(); - + if (dim == 0u) { @@ -178,7 +178,7 @@ inline BArrayVector::operator std::vector< Cell_Type >() co return Array->get_row_vec(i, false); else return Array->get_col_vec(i, false); - + } template @@ -205,7 +205,7 @@ inline bool BArrayVector::operator==(const Cell_Type & val) } return true; - + } template @@ -215,12 +215,12 @@ inline BArrayVector_const::operator std::vector< Cell_Type return Array->get_row_vec(i, false); else return Array->get_col_vec(i, false); - + } template inline bool BArrayVector_const::operator==(const Cell_Type & val) const { - + if (dim == 0u) { for (size_t j = 0u; j < Array->ncol(); ++j) @@ -252,7 +252,7 @@ inline bool BArrayVector_const::operator!=(const Cell_Type template inline bool BArrayVector_const::operator<(const Cell_Type & val) const { - + if (dim == 0u) { for (size_t j = 0u; j < Array->ncol(); ++j) @@ -279,7 +279,7 @@ inline bool BArrayVector_const::operator<(const Cell_Type & template inline bool BArrayVector_const::operator<=(const Cell_Type & val) const { - + if (dim == 0u) { for (size_t j = 0u; j < Array->ncol(); ++j) @@ -313,7 +313,7 @@ inline bool BArrayVector_const::operator>(const Cell_Type & template inline bool BArrayVector_const::operator>=(const Cell_Type & val) const { - return !(this->operator<(val)); -} + return !(this->operator<(val)); +} -#endif \ No newline at end of file +#endif diff --git a/include/barry/barry-configuration.hpp b/include/barry/barry-configuration.hpp index b85645b..e0024d5 100644 --- a/include/barry/barry-configuration.hpp +++ b/include/barry/barry-configuration.hpp @@ -4,18 +4,18 @@ /** * @name Configuration MACROS * @details These are mostly related to performance. The definitions follow: - * + * * - `BARRY_USE_UNORDERED_MAP` If specified, then barry is compiled using * `std::unordered_map`. Otherwise it will use `std::map` for the arrays. - * + * * - `BARRY_USE_SAFE_EXP` When specified, it will multiply all likelihoods * in `Model` by (1/-100)/(1/-100) so that numerical overflows are avoided. - * - * - `BARRY_USE_ISFINITE` When specified, it will introduce a macro that + * + * - `BARRY_USE_ISFINITE` When specified, it will introduce a macro that * checks whether the likelihood is finite or not. - * + * * - `printf_barry` If not specified, will be defined as `printf`. - * + * * - `BARRY_DEBUG_LEVEL`, when defined, will make things verbose. */ ///@{ @@ -28,7 +28,7 @@ #endif #ifdef BARRY_USE_SAFE_EXP - #define BARRY_SAFE_EXP + #define BARRY_SAFE_EXP #else #define BARRY_SAFE_EXP -100.0 #endif @@ -37,7 +37,7 @@ #define BARRY_ISFINITE(a) if (!std::isfinite( (a) )) \ throw std::overflow_error("The likelihood function has overflowed."); #else - #define BARRY_ISFINITE(a) + #define BARRY_ISFINITE(a) #endif #ifdef BARRAY_USE_CHECK_SUPPORT @@ -73,4 +73,4 @@ ///@} -#endif \ No newline at end of file +#endif diff --git a/include/barry/barry-debug.hpp b/include/barry/barry-debug.hpp index cdf499f..6112a21 100644 --- a/include/barry/barry-debug.hpp +++ b/include/barry/barry-debug.hpp @@ -13,7 +13,7 @@ template void BARRY_DEBUG_VEC_PRINT(const std::vector & a) { printf_barry("%s [", BARRY_DEBUG_HEADER); - for(const auto & iter : (a)) + for(const auto & iter : (a)) printf_barry("%.4f ", static_cast(iter)); printf_barry("]\n"); return; @@ -23,7 +23,7 @@ template<> inline void BARRY_DEBUG_VEC_PRINT(const std::vector< int > & a) { printf_barry("%s [", BARRY_DEBUG_HEADER); - for(const auto & iter : (a)) + for(const auto & iter : (a)) printf_barry("%i ", iter); printf_barry("]\n"); return; @@ -32,11 +32,11 @@ template<> inline void BARRY_DEBUG_VEC_PRINT(const std::vector< std::string > & a) { printf_barry("%s \n", BARRY_DEBUG_HEADER); - for(const auto & iter : (a)) + for(const auto & iter : (a)) printf_barry("%s %s\n", BARRY_DEBUG_HEADER, iter.c_str()); printf_barry("%s \n", BARRY_DEBUG_HEADER); return; } #endif -#endif \ No newline at end of file +#endif diff --git a/include/barry/barry-macros.hpp b/include/barry/barry-macros.hpp index c32f36c..bc8e7d3 100644 --- a/include/barry/barry-macros.hpp +++ b/include/barry/barry-macros.hpp @@ -11,9 +11,9 @@ #if defined(_OPENMP) || defined(__OPENMP) #define BARRY_NCORES_ARG(default) size_t ncores default -#else +#else #define BARRY_NCORES_ARG(default) size_t ncores default #endif -#endif \ No newline at end of file +#endif diff --git a/include/barry/barry.hpp b/include/barry/barry.hpp index 74b42a0..0735031 100644 --- a/include/barry/barry.hpp +++ b/include/barry/barry.hpp @@ -26,7 +26,7 @@ #endif #ifndef BARRY_HPP -#define BARRY_HPP +#define BARRY_HPP #define BARRY_VERSION_MAYOR 0 #define BARRY_VERSION_MINOR 1 @@ -36,7 +36,7 @@ * @brief barry: Your go-to motif accountant */ namespace barry { - + //! Tree class and TreeIterator class #include "typedefs.hpp" #include "barry-macros.hpp" @@ -60,7 +60,7 @@ namespace barry { #include "barraydense-meat.hpp" #include "barraydensecell-meat.hpp" #include "barraydense-meat-operators.hpp" - + #include "counters-bones.hpp" #include "counters-meat.hpp" @@ -75,16 +75,16 @@ namespace barry { #include "model-bones.hpp" #include "model-meat.hpp" - + #include "rules-bones.hpp" #include "rules-meat.hpp" - + namespace counters { namespace network { #include "counters/network.hpp" } } - + } namespace netcounters = barry::counters::network; @@ -103,4 +103,4 @@ namespace netcounters = barry::counters::network; Rule_fun_type a = \ [](const Array_Type & Array, size_t i, size_t j, Data_Type & data) -#endif \ No newline at end of file +#endif diff --git a/include/barry/cell-bones.hpp b/include/barry/cell-bones.hpp index 69c5bce..5d794f9 100644 --- a/include/barry/cell-bones.hpp +++ b/include/barry/cell-bones.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_CELL_BONES_HPP +#ifndef BARRY_CELL_BONES_HPP #define BARRY_CELL_BONES_HPP 1 /** @@ -16,33 +16,33 @@ template class Cell { Cell(Cell_Type value_, bool visited_ = false, bool active_ = true) : value(value_), visited(visited_), active(active_) {}; ~Cell() {}; - + // This is an explicit declaration since in other cases it seems // to try to use the move operator, which I do not intent to use. Cell(const Cell& arg) : value(arg.value), visited(arg.visited), active(arg.active) {}; - + // Copy by assignment Cell& operator=(const Cell& other); - + // Move constructor Cell(Cell&& arg) noexcept: value(std::move(arg.value)), visited(std::move(arg.visited)), active(std::move(arg.active)) {} ; - + // Move assign operator Cell& operator=(Cell&& other) noexcept; - + void add(Cell_Type x); - + // Casting operator (implicit and explicit) // int x = Cell(1); // returns 1 operator Cell_Type() const {return this->value;}; bool operator==(const Cell& rhs ) const; bool operator!=(const Cell& rhs ) const; - + }; #endif diff --git a/include/barry/cell-meat.hpp b/include/barry/cell-meat.hpp index 6d8b3c1..fe08b62 100644 --- a/include/barry/cell-meat.hpp +++ b/include/barry/cell-meat.hpp @@ -22,7 +22,7 @@ bool Cell::operator==(const Cell& rhs ) const { if (this == *rhs) return true; - + return this->value == rhs.value; } @@ -31,7 +31,7 @@ template bool Cell::operator!=(const Cell& rhs ) const { return !this->operator==(rhs); - + } @@ -64,4 +64,4 @@ template<> inline Cell< size_t >::Cell() : value(1u), visited(false), active(tru template<> inline Cell< int >::Cell() : value(1), visited(false), active(true) {} template<> inline Cell< bool >::Cell() : value(true), visited(false), active(true) {} -#endif \ No newline at end of file +#endif diff --git a/include/barry/col-bones.hpp b/include/barry/col-bones.hpp index ed95ed2..4dbb20b 100644 --- a/include/barry/col-bones.hpp +++ b/include/barry/col-bones.hpp @@ -54,4 +54,4 @@ // public: // }; -// #endif \ No newline at end of file +// #endif diff --git a/include/barry/counters-bones.hpp b/include/barry/counters-bones.hpp index 0233141..cf58937 100644 --- a/include/barry/counters-bones.hpp +++ b/include/barry/counters-bones.hpp @@ -6,35 +6,35 @@ * @details `barry` includes a flexible way to generate counters based on change * statistics. Since most of the time we are counting many motifs in a graph, * change statistics make a reasonable (and efficient) way to make such counts. - * + * * In particular, let the motif be defined as \f$s(y)\f$, with \f$y\f$ as the * binary array. The change statistic when adding cell \f$y_{ij}\f$, i.e. when * the cell moves from being emty to have a one, is defined as - * + * * \f[ * \delta(y_{ij}) = s^+_{ij}(y) - s^-_{ij}(y), * \f] - * + * * where \f$s^+_{ij}(y)\f$ and \f$s^-_{ij}(y)\f$ represent the motif statistic * with and without the ij-cell. For example, in the case of networks, the change - * statistic for the number of edges is always 1. - * - * To count statistics in an array, the [Counter] class will empty the array, + * statistic for the number of edges is always 1. + * + * To count statistics in an array, the [Counter] class will empty the array, * initialize the counters, and then start counting while adding at each step - * a single cell, until matching the original array. + * a single cell, until matching the original array. */ /** * @ingroup counting Implementation of motif counting * @brief A counter function based on change statistics. - * + * * This class is used by `CountStats` and `StatsCounter` as a way to count * statistics using change statistics. */ template , typename Data_Type = bool> class Counter { public: - + Counter_fun_type count_fun; Counter_fun_type init_fun; Hasher_fun_type hasher_fun; @@ -45,7 +45,7 @@ class Counter { /** * @name Creator passing a counter and an initializer - * + * * @param count_fun_ The main counter function. * @param init_fun_ The initializer function can also be used to check if the * `BArray` as the needed variables (see BArray::data). @@ -55,17 +55,17 @@ class Counter { */ ///@{ Counter() : count_fun(nullptr), init_fun(nullptr), hasher_fun(nullptr) {}; - + Counter( Counter_fun_type count_fun_, Counter_fun_type init_fun_, Hasher_fun_type hasher_fun_, Data_Type data_, - std::string name_ = "", + std::string name_ = "", std::string desc_ = "" ): count_fun(count_fun_), init_fun(init_fun_), hasher_fun(hasher_fun_), data(data_), name(name_), desc(desc_) {}; - + Counter(const Counter & counter_); ///< Copy constructor Counter(Counter && counter_) noexcept; ///< Move constructor Counter operator=(const Counter & counter_); ///< Copy assignment @@ -73,7 +73,7 @@ class Counter { ///@} ~Counter() {}; - + /*** * ! Main functions. */ @@ -84,87 +84,87 @@ class Counter { /** * @brief Get and set the hasher function - * + * * The hasher function is used to characterize the support of the array. * This way, if possible, the support enumeration is recycled. - * - * @param fun + * + * @param fun */ ///@{ void set_hasher(Hasher_fun_type fun); Hasher_fun_type get_hasher(); ///@} - + }; /** * @brief Vector of counters. - * + * * Various functions hold more than one counter, so this class is a helper class * that allows managing multiple counters efficiently. The main data is a vector * to pointers of counters. */ template , typename Data_Type = bool> class Counters { - + private: std::vector< Counter > data; Hasher_fun_type hasher; - -public: - + +public: + // Constructors Counters(); - + // Destructor needs to deal with the pointers ~Counters() {}; /** * @brief Copy constructor - * @param counter_ + * @param counter_ */ Counters(const Counters & counter_); - + /** * @brief Move constructor - * - * @param counters_ + * + * @param counters_ */ Counters(Counters && counters_) noexcept; /** * @brief Copy assignment constructor - * - * @param counter_ - * @return Counters + * + * @param counter_ + * @return Counters */ Counters operator=(const Counters & counter_); /** * @brief Move assignment constructor - * - * @param counter_ - * @return Counters& + * + * @param counter_ + * @return Counters& */ Counters & operator=(Counters && counter_) noexcept; - + /** * @brief Returns a pointer to a particular counter. - * + * * @param idx Id of the counter - * @return Counter* + * @return Counter* */ Counter & operator[](size_t idx); /** * @brief Number of counters in the set. - * - * @return size_t + * + * @return size_t */ std::size_t size() const noexcept { return data.size(); }; - + // Functions to add counters void add_counter(Counter counter); void add_counter( @@ -172,17 +172,17 @@ class Counters { Counter_fun_type init_fun_, Hasher_fun_type hasher_fun_, Data_Type data_, - std::string name_ = "", + std::string name_ = "", std::string desc_ = "" ); - + std::vector< std::string > get_names() const; std::vector< std::string > get_descriptions() const; /** * @brief Generates a hash for the given array according to the counters. - * - * @param array + * + * @param array * @param add_dims When `true` (default) the dimmension of the array will * be added to the hash. * @return std::vector< double > That can be hashed later. @@ -195,8 +195,7 @@ class Counters { void add_hash( Hasher_fun_type fun_ ); - + }; #endif - diff --git a/include/barry/counters-meat.hpp b/include/barry/counters-meat.hpp index d00976b..6598ca0 100644 --- a/include/barry/counters-meat.hpp +++ b/include/barry/counters-meat.hpp @@ -45,7 +45,7 @@ COUNTER_TEMPLATE(COUNTER_TYPE(),operator=)( this->init_fun = counter_.init_fun; this->hasher_fun = counter_.hasher_fun; - + this->data = counter_.data; this->name = counter_.name; this->desc = counter_.desc; @@ -155,7 +155,7 @@ COUNTERS_TEMPLATE(COUNTERS_TYPE(), operator=)(const Counters && counters_) noexcept +COUNTERS_TEMPLATE(COUNTERS_TYPE() &, operator=)(Counters && counters_) noexcept { if (this != &counters_) { @@ -169,9 +169,9 @@ COUNTERS_TEMPLATE(COUNTERS_TYPE() &, operator=)(Counters & COUNTERS_TEMPLATE(void, add_counter)(Counter counter) { - + data.push_back(counter); - + return; } @@ -184,7 +184,7 @@ COUNTERS_TEMPLATE(void, add_counter)( std::string desc_ ) { - + data.emplace_back(Counter( count_fun_, init_fun_, @@ -193,9 +193,9 @@ COUNTERS_TEMPLATE(void, add_counter)( name_, desc_ )); - + return; - + } COUNTERS_TEMPLATE(std::vector, get_names)() const @@ -212,7 +212,7 @@ COUNTERS_TEMPLATE(std::vector, get_names)() const COUNTERS_TEMPLATE(std::vector, get_descriptions)() const { - + std::vector< std::string > out; out.reserve(this->size()); for (size_t i = 0u; i < this->size(); ++i) @@ -228,7 +228,7 @@ COUNTERS_TEMPLATE(std::vector, gen_hash)( ) { std::vector res; - + // Iterating over the counters for (auto & c: data) { @@ -281,4 +281,4 @@ COUNTERS_TEMPLATE(void, add_hash)( #undef COUNTERS_TEMPLATE_ARGS #undef COUNTERS_TEMPLATE -#endif \ No newline at end of file +#endif diff --git a/include/barry/counters/network-css.hpp b/include/barry/counters/network-css.hpp index fa5f193..db6993b 100644 --- a/include/barry/counters/network-css.hpp +++ b/include/barry/counters/network-css.hpp @@ -1,7 +1,7 @@ #ifndef BARRY_CSS_COUNTERS #define BARRY_CSS_COUNTERS -// n: Net size, +// n: Net size, // s: Start of the i-th network // e: end of the i-th network // ego_id: Ego of the cell (i, j) @@ -29,7 +29,7 @@ }; // Variables in case that the current cell corresponds to the True -#define CSS_CASE_TRUTH() if ((i < n) && (j < n)) +#define CSS_CASE_TRUTH() if ((i < n) && (j < n)) // i_: i-th index of the cell // j_: j-th index of the cell @@ -71,7 +71,7 @@ /* The indices fall within the network */ \ if ((data.indices.at(0) > Array.ncol()) \ | (data.indices.at(2) > Array.ncol())) \ - throw std::range_error("The network does not match the prescribed size."); + throw std::range_error("The network does not match the prescribed size."); #define CSS_CHECK_SIZE() for (size_t i = 0u; i < end_.size(); ++i) {\ if (i == 0u) continue; \ @@ -93,13 +93,13 @@ /** - * @brief Counts errors of commission - * @param netsize Size of the reference (true) network + * @brief Counts errors of commission + * @param netsize Size of the reference (true) network * @param end_ Vector indicating one past the ending index of each network. (see details) * @param counter_type Size_t indicating the type of counter to use. Possible - * values are: 0: Count all, 1: Only count if perceiver is involved, and + * values are: 0: Count all, 1: Only count if perceiver is involved, and * 2: Only count if perceiver is not involved. - * @details + * @details * The `end_` parameter should be of length `N of networks` - 1. It is * assumed that the first network ends at `netsize`. */ @@ -110,9 +110,9 @@ inline void counter_css_partially_false_recip_commi( const std::vector< size_t > & end_, size_t counter_type = 0u ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + // Getting the network size CSS_SIZE() @@ -137,17 +137,17 @@ inline void counter_css_partially_false_recip_commi( } CSS_CASE_ELSE() return 0.0; - + }; CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("Partially false recip (comission)") return; - + } /** @brief Counts errors of omission */ @@ -158,12 +158,12 @@ inline void counter_css_partially_false_recip_omiss( const std::vector< size_t > & end_, size_t counter_type = 0u ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + // Getting the network size CSS_SIZE() - + // True network CSS_CASE_TRUTH() { @@ -176,7 +176,7 @@ inline void counter_css_partially_false_recip_omiss( } CSS_CASE_PERCEIVED() { CSS_PERCEIVED_CELLS() - return tji * tij * (1.0 - 2.0 * pji) - + return tji * tij * (1.0 - 2.0 * pji) - (1.0 - pji) * ((1.0 - tij) * tji + tij * (1.0 - tji)) ; @@ -184,15 +184,15 @@ inline void counter_css_partially_false_recip_omiss( return 0.0; }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("Partially false recip (omission)") - + return; - + } /** @brief Counts completely false reciprocity (comission) */ @@ -203,7 +203,7 @@ inline void counter_css_completely_false_recip_comiss( const std::vector< size_t > & end_, size_t counter_type = 0u ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Getting the network size @@ -225,15 +225,15 @@ inline void counter_css_completely_false_recip_comiss( return 0.0; }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("Completely false recip (comission)") - + return; - + } /** @brief Counts completely false reciprocity (omission) */ @@ -244,7 +244,7 @@ inline void counter_css_completely_false_recip_omiss( const std::vector< size_t > & end_, size_t counter_type = 0u ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Getting the network size @@ -264,17 +264,17 @@ inline void counter_css_completely_false_recip_omiss( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("Completely false recip (omission)") - + return; - + } /** @brief Counts mixed reciprocity errors */ @@ -285,7 +285,7 @@ inline void counter_css_mixed_recip( const std::vector< size_t > & end_, size_t counter_type = 0u ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Getting the network size @@ -305,17 +305,17 @@ inline void counter_css_mixed_recip( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("Mixed reciprocity errors") - + return; - + } /////////////////////////// CENSUS @@ -348,9 +348,9 @@ inline void counter_css_census01( } CSS_CASE_ELSE() return 0.0; - + }; - + // CSS_NET_COUNTER_LAMBDA_INIT() NETWORK_COUNTER_LAMBDA(tmp_init) { @@ -376,11 +376,11 @@ inline void counter_css_census01( return n_dbl * (n_dbl - 1.0) / 2.0; // / (Array.D().directed ? 1.0 : 2.0); }; - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(01) Accurate null") - + return; } @@ -412,15 +412,15 @@ inline void counter_css_census02( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(02) Partial false positive (null)") - + return; } @@ -452,15 +452,15 @@ inline void counter_css_census03( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(03) Complete false positive (null)") - + return; } @@ -492,15 +492,15 @@ inline void counter_css_census04( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(04) Partial false negative (assym)") - + return; } @@ -532,15 +532,15 @@ inline void counter_css_census05( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(05) Accurate assym") - + return; } @@ -572,15 +572,15 @@ inline void counter_css_census06( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(06) Mixed assym") - + return; } @@ -612,15 +612,15 @@ inline void counter_css_census07( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(07) Partial false positive (assym)") - + return; } @@ -652,15 +652,15 @@ inline void counter_css_census08( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(08) Complete false negative (full)") - + return; } @@ -692,15 +692,15 @@ inline void counter_css_census09( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(09) Partial false negative (full)") - + return; } @@ -732,15 +732,15 @@ inline void counter_css_census10( } CSS_CASE_ELSE() return 0.0; - + }; - + CSS_NET_COUNTER_LAMBDA_INIT() - + // checking sizes CSS_CHECK_SIZE() CSS_APPEND("(10) Accurate full") - + return; } diff --git a/include/barry/counters/network.hpp b/include/barry/counters/network.hpp index 782738f..99fdaec 100644 --- a/include/barry/counters/network.hpp +++ b/include/barry/counters/network.hpp @@ -2,28 +2,28 @@ #define BARRAY_NETWORK_H 1 /** - * @ingroup counting + * @ingroup counting * @details Details on the available counters for `NetworkData` can be found in * the \ref counters-network section. - * + * */ ///@{ /** * @brief Data class for Networks. - * + * * This holds information about whether the graph is directed or not, and, * if defined, vectors of node (vertex) attributes (`vertex_attr`). - * + * */ class NetworkData { public: - + bool directed = true; std::vector< std::vector< double > > vertex_attr; - + NetworkData() : vertex_attr(0u) {}; - + /** * @brief Constructor using a single attribute * @param vertex_attr_ Double vector of length equal to the number of vertices @@ -34,7 +34,7 @@ class NetworkData { std::vector< double > vertex_attr_, bool directed_ = true ) : directed(directed_), vertex_attr(1u, vertex_attr_) {}; - + /** * @brief Constructor using multiple attributes * @param vertex_attr_ Vector of double vectors. The size equals to the number @@ -46,8 +46,8 @@ class NetworkData { std::vector< std::vector< double > > vertex_attr_, bool directed_ = true ) : directed(directed_), vertex_attr(vertex_attr_) {}; - - + + ~NetworkData() {}; }; @@ -55,20 +55,20 @@ class NetworkData { * @brief Data class used to store arbitrary size_t or double vectors */ class NetCounterData { public: - + std::vector< size_t > indices; std::vector< double > numbers; - + NetCounterData() : indices(0u), numbers(0u) {}; NetCounterData( const std::vector< size_t > & indices_, const std::vector< double > & numbers_ ): indices(indices_), numbers(numbers_) {}; - + ~NetCounterData() {}; - + // const size_t get_size_t - + }; #define NET_C_DATA_IDX(i) (data.indices[i]) @@ -151,19 +151,19 @@ Rule_fun_type a = \ template inline void counter_edges(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(count_edges) { return 1.0; }; - + counters->add_counter( count_edges, nullptr, nullptr, - NetCounterData(), - "Edge counts", + NetCounterData(), + "Edge counts", "Number of edges" ); - + return; } @@ -174,33 +174,33 @@ inline void counter_edges(NetCounters * counters) template inline void counter_isolates(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { if (i == j) return 0.0; - + double res = 0.0; - + // i is sending its first tie if (Array.row(i).size() == 1u && Array.col(i).size() == 0u) res -= 1.0; - + // j is receiving its first tie, meaning that he // has no other tie but i's? if (Array.row(j).size() == 0u && Array.col(j).size() == 1u) res -= 1.0; - + return res; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { return static_cast(Array.nrow()); }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), @@ -214,15 +214,15 @@ inline void counter_isolates(NetCounters * counters) template<> inline void counter_isolates(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { if (i == j) return 0.0; - + double res = 0.0; - + // Checking the in and out degree if (Array.rowsum(i) == 1u && Array.colsum(i) == 0u) res -= 1.0; @@ -230,16 +230,16 @@ inline void counter_isolates(NetCounters * counters) // Now looking at j if (Array.rowsum(j) == 0u && Array.colsum(j) == 1u) res -= 1.0; - + return res; }; - + NETWORKDENSE_COUNTER_LAMBDA(tmp_init) { return static_cast(Array.nrow()); }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), @@ -255,16 +255,16 @@ inline void counter_isolates(NetCounters * counters) template inline void counter_mutual(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Is there any tie at ji? If not, then we have a new mutual! // but this only makes sence if the jth row and ith column exists - // if ((Array.nrow() > j) && (Array.ncol() > i)) + // if ((Array.nrow() > j) && (Array.ncol() > i)) if (i == j) return 0.0; - + // printf_barry("Checking if it is empty or not at (%i, %i)... ", i, j); if (!Array.is_empty(j, i, false)) { @@ -272,29 +272,29 @@ inline void counter_mutual(NetCounters * counters) return 1.0; } // printf_barry("No, no mutual.\n"); - + return 0.0; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { if (Array.nrow() != Array.ncol()) throw std::logic_error("The -mutual- counter only works on square arrays."); - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!Array.D_ptr()->directed) throw std::logic_error( "The -mutual- counter only works on directed (non-symmetric) arrays." ); - + return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), @@ -311,33 +311,33 @@ inline void counter_mutual(NetCounters * counters) template inline void counter_istar2(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Need to check the receiving, if he/she is getting a new set of stars // when looking at triads - + if (Array.col(j).size() == 1u) return 0.0; - + return static_cast(Array.col(j).size() - 1.0); }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), "Istar 2", "Indegree 2-star" ); - + return ; } template<> inline void counter_istar2(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { // Need to check the receiving, if he/she is getting a new set of stars @@ -354,19 +354,19 @@ inline void counter_istar2(NetCounters * counters) // if (indeg == 1) // return 0.0; - + // return static_cast(indeg - 1); return static_cast(Array.colsum(j) - 1); }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), "Istar 2", "Indegree 2-star" ); - + return ; } @@ -375,20 +375,20 @@ inline void counter_istar2(NetCounters * counters) template inline void counter_ostar2(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Need to check the receiving, if he/she is getting a new set of stars // when looking at triads - + if (Array.row(i).size() == 1u) return 0.0; - + return static_cast( Array.row(i).size() - 1.0); }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), @@ -397,13 +397,13 @@ inline void counter_ostar2(NetCounters * counters) ); return ; - + } template<> inline void counter_ostar2(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { @@ -418,12 +418,12 @@ inline void counter_ostar2(NetCounters * counters) // if (nties == 1u) // return 0.0; - + // return static_cast(nties - 1.0); return static_cast(Array.rowsum(i) - 1); }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), @@ -432,7 +432,7 @@ inline void counter_ostar2(NetCounters * counters) ); return ; - + } @@ -440,89 +440,89 @@ inline void counter_ostar2(NetCounters * counters) template inline void counter_ttriads(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { // Self ties do not count if (i == j) return 0.0; - + double ans = 0.0; - + // Case 1: i-j, i-k, j-k if (Array.row(j).size() < Array.row(i).size()) { - - for (auto j_row = Array.row(j).begin(); j_row != Array.row(j).end(); ++j_row) + + for (auto j_row = Array.row(j).begin(); j_row != Array.row(j).end(); ++j_row) if ((j != j_row->first) && (i != j_row->first) && !Array.is_empty(i, j_row->first, false)) ans += 1.0; - + } else { - - for (auto i_row = Array.row(i).begin(); i_row != Array.row(i).end(); ++i_row) + + for (auto i_row = Array.row(i).begin(); i_row != Array.row(i).end(); ++i_row) if ((i != i_row->first) && (i_row->first != j) && !Array.is_empty(j, i_row->first, false)) ans += 1.0; - + } - - // Case 2: i-j, i-k, k-j + + // Case 2: i-j, i-k, k-j if (Array.row(i).size() > Array.col(j).size()) { - + for (auto j_col = Array.col(j).begin(); j_col != Array.col(j).end(); ++j_col) if ((j != j_col->first) && (i != j_col->first) && !Array.is_empty(i, j_col->first, false)) ans += 1.0; - + } else { - - for (auto i_row = Array.row(i).begin(); i_row != Array.row(i).end(); ++i_row) + + for (auto i_row = Array.row(i).begin(); i_row != Array.row(i).end(); ++i_row) if ((i != i_row->first) && (j != i_row->first) && !Array.is_empty(i_row->first, j, false)) ans += 1.0; - + } - + // Case 3: i->j, k->j, k->i if (Array.col(i).size() > Array.col(j).size()) { - + for (auto j_col = Array.col(j).begin(); j_col != Array.col(j).end(); ++j_col) if ((j != j_col->first) && (i != j_col->first) && !Array.is_empty(j_col->first, i, false)) ans += 1.0; - + } else { - - for (auto i_col = Array.col(i).begin(); i_col != Array.col(i).end(); ++i_col) + + for (auto i_col = Array.col(i).begin(); i_col != Array.col(i).end(); ++i_col) if ((i != i_col->first) && (j != i_col->first) && !Array.is_empty(i_col->first, j, false)) ans += 1.0; - + } - + // The regular counter double counts return ans; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!(Array.D_ptr()->directed)) throw std::invalid_argument("The ttriads counter is only valid for directed networks. This is undirected."); return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), "Balance", "Number of directed triangles" ); - + return; } @@ -530,7 +530,7 @@ inline void counter_ttriads(NetCounters * counters) template<> inline void counter_ttriads(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { @@ -545,7 +545,7 @@ inline void counter_ttriads(NetCounters * counters) if (Array.rowsum(i) == BARRY_ZERO_NETWORK_DENSE) return 0.0; - + double ans = 0.0; for (size_t k = 0u; k < N; ++k) { @@ -563,43 +563,43 @@ inline void counter_ttriads(NetCounters * counters) if (dat[k * N + j]) ans += 1.0; - // Case 2: i-j, i-k, k-j + // Case 2: i-j, i-k, k-j if (dat[j * N + k] != BARRY_ZERO_NETWORK_DENSE) ans += 1.0; } - + // Case 3: i-j, k-i, k-j if ((dat[i * N + k] != BARRY_ZERO_NETWORK_DENSE) && (dat[j * N + k] != BARRY_ZERO_NETWORK_DENSE)) ans += 1.0; } } - + // The regular counter double counts return ans; }; - + NETWORKDENSE_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!(Array.D_ptr()->directed)) throw std::invalid_argument("The ttriads counter is only valid for directed networks. This is undirected."); return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), "Balance", "Number of directed triangles" ); - + return; } @@ -609,39 +609,39 @@ inline void counter_ttriads(NetCounters * counters) template inline void counter_ctriads(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { if (i == j) return 0.0; - + double ans = 0.0; if (Array.col(i).size() < Array.row(j).size()) { - - for (auto i_col = Array.col(i).begin(); i_col != Array.col(i).end(); ++i_col) + + for (auto i_col = Array.col(i).begin(); i_col != Array.col(i).end(); ++i_col) if ((i != i_col->first) && (j != i_col->first) && !Array.is_empty(j, i_col->first, false)) ans += 1.0; - + } else { - - for (auto j_row = Array.row(j).begin(); j_row != Array.row(j).end(); ++j_row) + + for (auto j_row = Array.row(j).begin(); j_row != Array.row(j).end(); ++j_row) if ((j != j_row->first) && (i != j_row->first) && !Array.is_empty(j_row->first, i, false)) ans += 1.0; - + } - + return ans; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!(Array.D_ptr()->directed)) throw std::invalid_argument( "The ctriads counter is only valid for directed networks. This is undirected." @@ -650,7 +650,7 @@ inline void counter_ctriads(NetCounters * counters) return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), @@ -658,22 +658,22 @@ inline void counter_ctriads(NetCounters * counters) ); return; - + } template<> inline void counter_ctriads(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { if (i == j) return 0.0; - + // i->j->k->i double ans = 0.0; - #if defined(__OPENMP) || defined(_OPENMP) + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:ans) #endif for (size_t k = 0u; k < Array.nrow(); ++k) @@ -694,17 +694,17 @@ inline void counter_ctriads(NetCounters * counters) } } - + return ans; }; - + NETWORKDENSE_COUNTER_LAMBDA(tmp_init) { if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!(Array.D_ptr()->directed)) throw std::invalid_argument( "The ctriads counter is only valid for directed networks. This is undirected." @@ -713,7 +713,7 @@ inline void counter_ctriads(NetCounters * counters) return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData(), @@ -721,25 +721,25 @@ inline void counter_ctriads(NetCounters * counters) ); return; - + } - + // Density -------------------------------------------------------------- template inline void counter_density(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + return 1.0/(Array.nrow() * (Array.ncol() - 1.0)) / ( (Array.D_ptr()->directed)? 1.0 : 2.0 ); - + }; - - // Preparing the counter data and returning. We make sure that the memory is + + // Preparing the counter data and returning. We make sure that the memory is // released so we set delete_data = true. counters->add_counter( tmp_count, nullptr, nullptr, @@ -749,28 +749,28 @@ inline void counter_density(NetCounters * counters) ); return ; - + } // idegree1.5 ------------------------------------------------------------- template inline void counter_idegree15(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + // In case of the first, we need to add if (Array.col(j).size() == 1u) return 1.0; - - return + + return pow(static_cast (Array.col(j).size()), 1.5) - pow(static_cast (Array.col(j).size() - 1), 1.5) ; - + }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), @@ -778,16 +778,16 @@ inline void counter_idegree15(NetCounters * counters) ); return; - + } template<> inline void counter_idegree15(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { - + // In case of the first, we need to add int ideg = 0; for (size_t k = 0u; k < Array.nrow(); ++k) @@ -802,7 +802,7 @@ inline void counter_idegree15(NetCounters * counters) if (ideg == 0) return 0.0; - + if (ideg == 1) return 1.0; @@ -811,14 +811,14 @@ inline void counter_idegree15(NetCounters * counters) if (std::isnan(res)) throw std::domain_error("Resulting indeg is undefined."); - - return + + return std::pow(static_cast (ideg), 1.5) - std::pow(static_cast (ideg - 1.0), 1.5) ; - + }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), @@ -826,28 +826,28 @@ inline void counter_idegree15(NetCounters * counters) ); return; - + } // odegree1.5 ------------------------------------------------------------- template inline void counter_odegree15(NetCounters * counters) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + // In case of the first, we need to add if (Array.row(i).size() == 1u) return 1.0; - - return + + return pow(static_cast(Array.row(i).size()), 1.5) - pow(static_cast(Array.row(i).size() - 1), 1.5) ; - + }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), @@ -855,16 +855,16 @@ inline void counter_odegree15(NetCounters * counters) ); return; - + } template<> inline void counter_odegree15(NetCounters * counters) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { - + // In case of the first, we need to add int odeg = 0; for (size_t k = 0u; k < Array.ncol(); ++k) @@ -883,14 +883,14 @@ inline void counter_odegree15(NetCounters * counters) if (odeg == 1) return 1.0; - - return + + return pow(static_cast(odeg), 1.5) - pow(static_cast(odeg - 1), 1.5) ; - + }; - + counters->add_counter( tmp_count, nullptr, nullptr, NetCounterData(), @@ -898,7 +898,7 @@ inline void counter_odegree15(NetCounters * counters) ); return; - + } @@ -910,43 +910,43 @@ inline void counter_absdiff( size_t attr_id, double alpha = 1.0 ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + return std::pow(std::fabs( - Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] - + Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] - Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][j] ), NET_C_DATA_NUM(0u)); - + }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (Array.D_ptr()->vertex_attr.size() == 0u) throw std::range_error("No attributes in the Array."); - + if ((NET_C_DATA_IDX(0u) != 0u) && (Array.D_ptr()->vertex_attr.size() <= (NET_C_DATA_IDX(0u) - 1u))) throw std::range_error("Attribute index out of range."); - + return 0.0; - + }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData({attr_id}, {alpha}), "Absdiff" ); - + return; - + } - + // ----------------------------------------------------------------------------- /**@brief Sum of attribute difference between ego and alter to pow(alpha)*/ template @@ -956,58 +956,58 @@ inline void counter_diff( double alpha = 1.0, double tail_head = true ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + return std::pow(NET_C_DATA_NUM(1u) * ( - Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] - + Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] - Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][j] ), NET_C_DATA_NUM(0u)); - + }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (Array.D_ptr()->vertex_attr.size() == 0u) throw std::range_error("No attributes in the Array."); - + if ((NET_C_DATA_IDX(0u) != 0u) && (Array.D_ptr()->vertex_attr.size() <= (NET_C_DATA_IDX(0u) - 1u))) throw std::range_error("Attribute index out of range."); - + return 0.0; - + }; - + counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData({attr_id}, {alpha, tail_head ? 1.0: -1.0}), "Absdiff^(" + std::to_string(alpha) + ")" ); - + return; - + } // Nodeicov, nodeocov, and Nodematch ------------------------------------------- NETWORK_COUNTER(init_single_attr) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (Array.D_ptr()->vertex_attr.size() == 0u) throw std::range_error("No attributes in the Array."); - + if ((NET_C_DATA_IDX(0u) != 0u) && (Array.D_ptr()->vertex_attr.size() <= (NET_C_DATA_IDX(0u) - 1u))) throw std::range_error("Attribute index out of range."); - + return 0.0; - + } // ----------------------------------------------------------------------------- @@ -1017,20 +1017,20 @@ inline void counter_nodeicov( NetCounters * counters, size_t attr_id ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + return Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][j]; - + }; - + counters->add_counter( tmp_count, init_single_attr, nullptr, NetCounterData({attr_id}, {}), "nodeicov", "Sum of ego attribute" ); - + return; } @@ -1042,20 +1042,20 @@ inline void counter_nodeocov( NetCounters * counters, size_t attr_id ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + return Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i]; - + }; - + counters->add_counter( tmp_count, init_single_attr, nullptr, NetCounterData({attr_id}, {}), "nodeocov", "Sum of alter attribute" ); - + return; } @@ -1067,21 +1067,21 @@ inline void counter_nodecov( NetCounters * counters, size_t attr_id ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + return Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] + Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][j]; - + }; - + counters->add_counter( tmp_count, init_single_attr, nullptr, NetCounterData({attr_id}, {}), "nodecov", "Sum of nodes covariates" ); - + return; } @@ -1092,19 +1092,19 @@ inline void counter_nodematch( NetCounters * counters, size_t attr_id ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - - return + + return ( - Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] == + Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][i] == Array.D_ptr()->vertex_attr[NET_C_DATA_IDX(0u)][j] )? 1.0 : 0.0; - + }; - - // Preparing the counter data and returning. We make sure that the memory is + + // Preparing the counter data and returning. We make sure that the memory is // released so we set delete_data = true. counters->add_counter( tmp_count, init_single_attr, nullptr, @@ -1112,9 +1112,9 @@ inline void counter_nodematch( "Homophily", "Number of homophilic ties" ); - + return ; - + } // ----------------------------------------------------------------------------- @@ -1127,33 +1127,33 @@ inline void counter_idegree( NETWORK_COUNTER_LAMBDA(tmp_count) { - + size_t d = Array.col(j).size(); if (d == NET_C_DATA_IDX(0u)) return 1.0; else if (d == (NET_C_DATA_IDX(0u) + 1)) return -1.0; - + return 0.0; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!Array.D_ptr()->directed) throw std::logic_error("-odegree- counter is only valid for directed graphs"); - + if (NET_C_DATA_IDX(0u) == 0u) return static_cast(Array.nrow()); - + return 0.0; }; - + for (auto iter = d.begin(); iter != d.end(); ++iter) counters->add_counter( tmp_count, tmp_init, nullptr, @@ -1161,8 +1161,8 @@ inline void counter_idegree( "Nodes indeg " + std::to_string(*iter), "Number of nodes with indigree " + std::to_string(*iter) ); - - return; + + return; } @@ -1174,7 +1174,7 @@ inline void counter_idegree( NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { - + size_t indeg = 0u; for (size_t k = 0u; k < Array.nrow(); ++k) if (Array(k, j) != BARRY_ZERO_NETWORK_DENSE) @@ -1184,27 +1184,27 @@ inline void counter_idegree( return 1.0; else if (indeg == (NET_C_DATA_IDX(0u) + 1)) return -1.0; - + return 0.0; }; - + NETWORKDENSE_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!Array.D_ptr()->directed) throw std::logic_error("-odegree- counter is only valid for directed graphs"); - + if (NET_C_DATA_IDX(0u) == 0u) return static_cast(Array.nrow()); - + return 0.0; }; - + for (auto iter = d.begin(); iter != d.end(); ++iter) counters->add_counter( tmp_count, tmp_init, nullptr, @@ -1212,8 +1212,8 @@ inline void counter_idegree( "Nodes indeg " + std::to_string(*iter), "Number of nodes with indigree " + std::to_string(*iter) ); - - return; + + return; } @@ -1224,99 +1224,99 @@ inline void counter_odegree( NetCounters * counters, std::vector d ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + size_t d = Array.row(i).size(); if (d == NET_C_DATA_IDX(0u)) return 1.0; else if (d == (NET_C_DATA_IDX(0u) + 1)) return -1.0; - + return 0.0; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!Array.D_ptr()->directed) throw std::logic_error("-odegree- counter is only valid for directed graphs"); - + if (NET_C_DATA_IDX(0u) == 0u) return static_cast(Array.nrow()); - + return 0.0; }; - - - for (auto iter = d.begin(); iter != d.end(); ++iter) + + + for (auto iter = d.begin(); iter != d.end(); ++iter) counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData({*iter}, {}), "Nodes w/ outdeg " + std::to_string(*iter), "Number of nodes with outdegree " + std::to_string(*iter) ); - - return; - + + return; + } - + template<> inline void counter_odegree( NetCounters * counters, std::vector d ) { - + NETWORKDENSE_COUNTER_LAMBDA(tmp_count) { - + size_t d = 0; for (size_t k = 0u; k < Array.ncol(); ++k) if (Array(i, k) != BARRY_ZERO_NETWORK_DENSE) d++; - + if (d == NET_C_DATA_IDX(0u)) return 1.0; else if (d == (NET_C_DATA_IDX(0u) + 1)) return -1.0; - + return 0.0; }; - + NETWORKDENSE_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (!Array.D_ptr()->directed) throw std::logic_error("-odegree- counter is only valid for directed graphs"); - + if (NET_C_DATA_IDX(0u) == 0u) return static_cast(Array.nrow()); - + return 0.0; }; - - - for (auto iter = d.begin(); iter != d.end(); ++iter) + + + for (auto iter = d.begin(); iter != d.end(); ++iter) counters->add_counter( tmp_count, tmp_init, nullptr, NetCounterData({*iter}, {}), "Nodes w/ outdeg " + std::to_string(*iter), "Number of nodes with outdegree " + std::to_string(*iter) ); - - return; - + + return; + } @@ -1327,33 +1327,33 @@ inline void counter_degree( NetCounters * counters, std::vector d ) { - + NETWORK_COUNTER_LAMBDA(tmp_count) { - + size_t d = Array.row(i).size(); if (d == NET_C_DATA_IDX(0u)) return 1.0; else if (d == (NET_C_DATA_IDX(0u) + 1)) return -1.0; - + return 0.0; }; - + NETWORK_COUNTER_LAMBDA(tmp_init) { - + if (Array.D_ptr() == nullptr) throw std::logic_error("The array data has not been initialized"); - + if (Array.D_ptr()->directed) throw std::logic_error("-degree- counter is only valid for undirected graphs"); - + if (NET_C_DATA_IDX(0u) == 0u) return static_cast(Array.nrow()); - + return 0.0; }; - - + + for (auto iter = d.begin(); iter != d.end(); ++iter) { counters->add_counter( @@ -1361,8 +1361,8 @@ inline void counter_degree( NetCounterData({*iter}, {}) ); } - - return; + + return; } #include "network-css.hpp" @@ -1379,17 +1379,17 @@ inline void counter_degree( /**@brief Number of edges */ template inline void rules_zerodiag(NetRules * rules) { - + NETWORK_RULE_LAMBDA(no_self_tie) { return i != j; }; - + rules->add_rule( no_self_tie, false, "No self-ties", "No self-ties" ); - + return; } diff --git a/include/barry/freqtable.hpp b/include/barry/freqtable.hpp index 17b14da..e5eeb09 100644 --- a/include/barry/freqtable.hpp +++ b/include/barry/freqtable.hpp @@ -1,24 +1,24 @@ -#ifndef BARRY_STATSDB_HPP +#ifndef BARRY_STATSDB_HPP #define BARRY_STATSDB_HPP 1 - + /** * @brief Frequency table of vectors - * + * * This is mostly used in `Support`. The main data is contained in the * `data` double vector. The matrix is stored in a row-wise fashion, where * the first element is the frequency with which the vector is observed. - * + * * For example, in a model with `k` terms the first k + 1 elements of * `data` would be: - * + * * - weights * - term 1 * - term 2 * - ... * - term k - * + * */ -template +template class FreqTable { private: @@ -28,18 +28,18 @@ class FreqTable { size_t n = 0u; typename std::unordered_map::iterator iter; - + public: // size_t ncols; FreqTable() {}; ~FreqTable() {}; - + size_t add(const std::vector< T > & x, size_t * h_precomp); - + Counts_type as_vector() const; const std::vector< double > & get_data() const {return data;}; const std::unordered_map & get_index() const {return index;}; - + void clear(); void reserve(size_t n, size_t k); void print() const; @@ -47,20 +47,20 @@ class FreqTable { /** * @brief Number of unique elements in the table. * ( - * @return size_t + * @return size_t */ size_t size() const noexcept; size_t make_hash(const std::vector< T > & x) const; - + }; -template +template inline size_t FreqTable::add( const std::vector< T > & x, size_t * h_precomp - ) { - + ) { + // The term exists, then we add it to the list and we initialize it // with a single count size_t h; @@ -91,13 +91,13 @@ inline size_t FreqTable::add( throw std::length_error( "The value you are trying to add doesn't have the same lenght used in the database." ); - + #if __cplusplus > 201700L auto iter2 = index.try_emplace(h, data.size()); - + if (!iter2.second) { - + data[(iter2.first)->second] += 1.0; } @@ -112,45 +112,45 @@ inline size_t FreqTable::add( if (iter == index.end()) { - + index.insert({h, data.size()}); data.push_back(1.0); data.insert(data.end(), x.begin(), x.end()); n++; - + return h; } data[(*iter).second] += 1.0; - + #endif - + } - + return h; } template inline Counts_type FreqTable::as_vector() const -{ - +{ + Counts_type ans; ans.reserve(index.size()); for (size_t i = 0u; i < n; ++i) { - + std::vector< double > tmp(k, 0.0); for (size_t j = 1u; j < (k + 1u); ++j) tmp[j - 1u] = data[i * (k + 1) + j]; - + ans.push_back( std::make_pair,size_t>( std::move(tmp), @@ -159,8 +159,8 @@ inline Counts_type FreqTable::as_vector() const ); } - - + + return ans; } @@ -241,14 +241,14 @@ inline size_t FreqTable::make_hash(const std::vector< T > & x) const std::hash< T > hasher; std::size_t hash = hasher(x[0u]); - + // ^ makes bitwise XOR // 0x9e3779b9 is a 32 bit constant (comes from the golden ratio) // << is a shift operator, something like lhs * 2^(rhs) if (x.size() > 1u) for (size_t i = 1u; i < x.size(); ++i) hash ^= hasher(x[i]) + 0x9e3779b9 + (hash<<6) + (hash>>2); - + return hash; } diff --git a/include/barry/model-bones.hpp b/include/barry/model-bones.hpp index 2adc79a..6f75929 100644 --- a/include/barry/model-bones.hpp +++ b/include/barry/model-bones.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_MODEL_BONES_HPP +#ifndef BARRY_MODEL_BONES_HPP #define BARRY_MODEL_BONES_HPP 1 /** @@ -13,14 +13,14 @@ * \sum_{A'\in \mathcal{A}}\exp{\left(\theta^{\mbox{t}}c(A')\right)} * } * \f] - * + * * This implementation aims to reduce the number of times that the support * needs to be computed. Models included here use more than a single array, and * thus allow the function to recycle support sets as needed. For example, * if we are looking at directed graphs all of the same size and without * vertex level features, i.e. a model that only counts edges, triangles, etc. * then the support needs to be fully computed only once. - * + * * @tparam Array_Type Class of `BArray` object. * @tparam Data_Counter_Type Any type. * @tparam Data_Rule_Type Any type. @@ -43,12 +43,12 @@ class Model { bool delete_rengine = false; /** - * @name Information about the arrays used in the model + * @name Information about the arrays used in the model * @details `stats_target` holds the observed sufficient statistics for each * array in the dataset. `array_frequency` contains the frequency with which - * each of the target stats_target (arrays) shows in the support. `array2support` + * each of the target stats_target (arrays) shows in the support. `array2support` * maps array indices (0, 1, ...) to the corresponding support. - * + * * Each vector of `stats_support` has the data stored in a row-wise order, * with each row starting with the weights, e.g., in a model with `k` terms * the first k + 1 elements of `stats_support` would be: @@ -87,7 +87,7 @@ class Model { std::vector< size_t > pset_sizes; ///< Number of vectors included in the support. std::vector< size_t > pset_locations; ///< Accumulated number of vectors included in the support. ///@} - + /** * @name Functions to compute statistics * @details Arguments are recycled to save memory and computation. @@ -99,7 +99,7 @@ class Model { Support support_fun; StatsCounter counter_fun; ///@} - + /**@brief Vector of the previously used parameters */ std::vector< std::vector > params_last; std::vector< std::vector > params_last_pset; @@ -112,16 +112,16 @@ class Model { /** * @brief Transformation of the model - * + * * @details When specified, this function will update the model by modifying * the linear equation. For example, if the user wanted to add interaction * terms, rescale, or apply other operations of the sorts, the user can do such * through this function. - * + * * The function should return `void` and receive the following arguments: * - `data` Pointer to the first element of the set of sufficient statistics * - `k` size_t indicating the number of sufficient statistics - * + * * @returns * Nothing, but it will modify the model data. */ @@ -129,7 +129,7 @@ class Model { transform_model_fun = nullptr; std::vector< std::string > transform_model_term_names; - + public: /** @@ -154,7 +154,7 @@ class Model { BARRY_NCORES_ARG(=1), int i = -1 ); - + void set_rengine(std::mt19937 * rengine_, bool delete_ = false) { if (delete_rengine) @@ -162,7 +162,7 @@ class Model { rengine = rengine_; delete_rengine = delete_; - + }; void set_seed(size_t s) { @@ -177,7 +177,7 @@ class Model { }; ///@} - + Model(); Model(size_t size_); Model(const Model & Model_); @@ -198,12 +198,12 @@ class Model { if (delete_rengine) delete rengine; }; - + void store_psets() noexcept; std::vector< double > gen_key(const Array_Type & Array_); - + /** - * @name Wrappers for the `Counters` member. + * @name Wrappers for the `Counters` member. * @details These will add counters to the model, which are shared by the * support and the actual counter function. */ @@ -217,9 +217,9 @@ class Model { void set_counters(Counters * counters_); void add_hasher(Hasher_fun_type fun_); ///@} - + /** - * @name Wrappers for the `Rules` member. + * @name Wrappers for the `Rules` member. * @details These will add rules to the model, which are shared by the * support and the actual counter function. */ @@ -229,7 +229,7 @@ class Model { Rule_fun_type count_fun_, Data_Rule_Type data_ ); - + void set_rules(Rules * rules_); void add_rule_dyn(Rule & rule); @@ -237,10 +237,10 @@ class Model { Rule_fun_type count_fun_, Data_Rule_Dyn_Type data_ ); - + void set_rules_dyn(Rules * rules_); ///@} - + /** * @brief Adds an array to the support of not already included. @@ -248,20 +248,20 @@ class Model { * @param force_new If `false`, it will use `keygen` to obtain a double vector * and create a hash of it. If the hash has been computed earlier, the support * is recycled. - * + * * @return The number of the array. */ size_t add_array(const Array_Type & Array_, bool force_new = false); - - + + /** * @name Likelihood functions. * @details Calculation of likelihood functions is done reusing normalizing - * constants. Before recalculating the normalizing constant, the function + * constants. Before recalculating the normalizing constant, the function * checks whether `params` matches the last set vector of parameters used * to compute it. - * - * + * + * * @param params Vector of parameters * @param as_log When `true`, the function returns the log-likelihood. */ @@ -272,7 +272,7 @@ class Model { bool as_log = false, bool no_update_normconst = false ); - + double likelihood( const std::vector & params, const Array_Type & Array_, @@ -280,7 +280,7 @@ class Model { bool as_log = false, bool no_update_normconst = false ); - + double likelihood( const std::vector & params, const std::vector & target_, @@ -296,7 +296,7 @@ class Model { bool as_log = false, bool no_update_normconst = false ); - + double likelihood_total( const std::vector & params, bool as_log = false, @@ -306,7 +306,7 @@ class Model { ///@} /** - * @name Extract elements by index + * @name Extract elements by index * @param i Index relative to the array in the model. * @param params A new vector of model parameters to compute the normalizing * constant. @@ -325,14 +325,14 @@ class Model { const size_t & i ); ///@} - + void print_stats(size_t i) const; /** * @brief Prints information about the model */ virtual void print() const; - + /** * @brief Sample a single array from the model * @param Array_ Baseline array to sample from. @@ -343,14 +343,14 @@ class Model { Array_Type sample(const Array_Type & Array_, const std::vector & params = {}); Array_Type sample(const size_t & i, const std::vector & params); ///@} - + /** * @brief Conditional probability ("Gibbs sampler") - * + * * @details Computes the conditional probability of observing * P{Y(i,j) = | Y^C, theta}, i.e., the probability of observing the entry Y(i,j) equal * to one given the rest of the array. - * + * * @param Array_ Array to check * @param params Vector of parameters * @param i Row entry @@ -363,14 +363,14 @@ class Model { size_t i, size_t j ); - + /** * @name Size of the model - * + * * @brief Number of different supports included in the model - * + * * This will return the size of `stats_target`. - * + * * @return `size()` returns the number of arrays in the model. * @return `size_unique()` returns the number of unique arrays (according to * the hasher) in the model. @@ -395,7 +395,7 @@ class Model { /** * @brief Raw pointers to the support and target statistics - * @details + * @details * The support of the model is stored as a vector of vector. Each * element of it contains the support for an specific type of array included. * It represents an array of size `(k + 1) x n unique elements`, with the data @@ -411,7 +411,7 @@ class Model { std::vector< size_t > * get_arrays2support(); std::vector< std::vector< Array_Type > > * get_pset_arrays(); std::vector< double > * get_pset_stats(); ///< Statistics of the support(s) - std::vector< double > * get_pset_probs(); + std::vector< double > * get_pset_probs(); std::vector< size_t > * get_pset_sizes(); std::vector< size_t > * get_pset_locations(); ///@} @@ -419,11 +419,11 @@ class Model { /** * @brief Set the transform_model_fun object * @details The transform_model function is used to transform the data - * - * @param data - * @param target - * @param n_arrays - * @param arrays2support + * + * @param data + * @param target + * @param n_arrays + * @param arrays2support */ ///@{ void set_transform_model( @@ -439,4 +439,4 @@ class Model { }; -#endif \ No newline at end of file +#endif diff --git a/include/barry/model-meat.hpp b/include/barry/model-meat.hpp index 720e1b5..d7ef6e2 100644 --- a/include/barry/model-meat.hpp +++ b/include/barry/model-meat.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_MODEL_MEAT_HPP +#ifndef BARRY_MODEL_MEAT_HPP #define BARRY_MODEL_MEAT_HPP 1 /** @@ -21,9 +21,9 @@ inline double update_normalizing_constant( { const double p = params[j]; - + #if defined(__OPENMP) || defined(_OPENMP) - #pragma omp simd + #pragma omp simd #elif defined(__GNUC__) && !defined(__clang__) #pragma GCC ivdep #endif @@ -32,7 +32,7 @@ inline double update_normalizing_constant( } - // Accumulate resv to a double res + // Accumulate resv to a double res #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:res) #elif defined(__GNUC__) && !defined(__clang__) @@ -65,7 +65,7 @@ inline double update_normalizing_constant( #endif return res; - + } inline double likelihood_( @@ -75,12 +75,12 @@ inline double likelihood_( size_t n_params, bool log_ = false ) { - + if (n_params != params.size()) throw std::length_error("-stats_target- and -params- should have the same length."); - + double numerator = 0.0; - + // Computing the numerator #ifdef __INTEL_LLVM_COMPILER #pragma code_align 32 @@ -124,7 +124,7 @@ inline double likelihood_( #endif return ans; - + } template < @@ -145,10 +145,10 @@ inline void Model 1u) && (n < 128u)) ncores = 1u; - + if (i >= 0) ncores = 1u; - + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp parallel for firstprivate(params) num_threads(ncores) \ shared(n, normalizing_constants, first_calc_done, \ @@ -174,7 +174,7 @@ inline void Model & params, size_t ncores ) { - + update_normalizing_constants(params, ncores); size_t n_params = params.size(); @@ -211,9 +211,9 @@ inline void Model -1) params_last[i] = params; @@ -241,7 +241,7 @@ inline void Model= 0) - ncores = 1u; + ncores = 1u; #if defined(__OPENMP) || defined(_OPENMP) #pragma omp parallel for num_threads(ncores) collapse(1) \ @@ -264,7 +264,7 @@ inline void Model::M delete_rules_dyn(true), transform_model_fun(nullptr), transform_model_term_names(0u) -{ +{ // Counters are shared support_fun.set_counters(counters); counter_fun.set_counters(counters); - + // Rules are shared support_fun.set_rules(rules); support_fun.set_rules_dyn(rules_dyn); - + return; - + } template < @@ -363,7 +363,7 @@ inline Model::M stats_support_n_arrays(0u), stats_target(0u), stats_likelihood(0u), - arrays2support(0u), keys2support(0u), + arrays2support(0u), keys2support(0u), pset_arrays(0u), pset_stats(0u), counters(new Counters()), rules(new Rules()), @@ -374,20 +374,20 @@ inline Model::M transform_model_fun(nullptr), transform_model_term_names(0u) { - + stats_target.reserve(size_); arrays2support.reserve(size_); // Counters are shared support_fun.set_counters(counters); counter_fun.set_counters(counters); - + // Rules are shared support_fun.set_rules(rules); support_fun.set_rules_dyn(rules_dyn); - + return; - + } template < @@ -398,7 +398,7 @@ template < > inline Model::Model( const Model & Model_ - ) : + ) : stats_support(Model_.stats_support), stats_support_sizes(Model_.stats_support_sizes), stats_support_sizes_acc(Model_.stats_support_sizes_acc), @@ -427,17 +427,17 @@ inline Model::M transform_model_fun(Model_.transform_model_fun), transform_model_term_names(Model_.transform_model_term_names) { - + // Counters are shared support_fun.set_counters(counters); counter_fun.set_counters(counters); - + // Rules are shared support_fun.set_rules(rules); support_fun.set_rules_dyn(rules_dyn); return; - + } template < @@ -446,11 +446,11 @@ template < typename Data_Rule_Type, typename Data_Rule_Dyn_Type > -inline Model & +inline Model & Model::operator=( const Model & Model_ ) { - + // Clearing if (this != &Model_) { @@ -459,10 +459,10 @@ inline Model & if (delete_rules) delete rules; - + if (delete_rules_dyn) delete rules_dyn; - + stats_support = Model_.stats_support; stats_support_sizes = Model_.stats_support_sizes; stats_support_sizes_acc = Model_.stats_support_sizes_acc; @@ -492,15 +492,15 @@ inline Model & // Counters are shared support_fun.set_counters(counters); counter_fun.set_counters(counters); - + // Rules are shared support_fun.set_rules(rules); support_fun.set_rules_dyn(rules_dyn); - + } - + return *this; - + } template @@ -513,14 +513,14 @@ template Model:: gen_key( const Array_Type & Array_ ) { - return this->counters->gen_hash(Array_); + return this->counters->gen_hash(Array_); } template inline void Model:: add_counter( Counter & counter ) { - + counters->add_counter(counter, Data_Counter_Type()); return; } @@ -531,15 +531,15 @@ inline void Model init_fun_, Data_Counter_Type data_ ) { - + counters->add_counter( count_fun_, init_fun_, data_ ); - + return; - + } template @@ -551,13 +551,13 @@ inline void Modelcounters = counters_; support_fun.set_counters(counters); counter_fun.set_counters(counters); - + return; - + } template @@ -575,7 +575,7 @@ template :: add_rule( Rule & rules ) { - + rules->add_rule(rules, Data_Rule_Type()); return; } @@ -603,7 +603,7 @@ template :: add_rule_dyn( Rule & rules_ ) { - + rules_dyn->add_rule(rules_, Data_Rule_Dyn_Type()); return; } @@ -613,14 +613,14 @@ inline void Model rule_fun_, Data_Rule_Dyn_Type data_ ) { - + rules_dyn->add_rule( rule_fun_, data_ ); - + return; - + } template @@ -646,13 +646,13 @@ Model::add_arr const Array_Type & Array_, bool force_new ) { - + // Array counts (target statistics) counter_fun.reset_array(&Array_); - + if (transform_model_fun) { - + auto tmpcounts = counter_fun.count_all(); stats_target.emplace_back( transform_model_fun(&tmpcounts[0u], tmpcounts.size()) @@ -660,7 +660,7 @@ Model::add_arr } else stats_target.push_back(counter_fun.count_all()); - + // If the data hasn't been analyzed earlier, then we need to compute // the support std::vector< double > key = counters->gen_hash(Array_); @@ -670,36 +670,36 @@ Model::add_arr // Current size of the support stats size_t stats_support_size = stats_support.size(); - + // Adding to the map keys2support[key] = stats_support_sizes.size(); stats_support_n_arrays.push_back(1u); // How many elements now arrays2support.push_back(stats_support_sizes.size()); // Map of the array id to the support - + // Computing support using the counters included in the model support_fun.reset_array(Array_); - + /** When computing with the powerset, we need to grow the corresponding * vectors on the fly */ if (with_pset) { - + // Making space for storing the support pset_arrays.resize(pset_arrays.size() + 1u); - + try { - + support_fun.calc( &(pset_arrays[pset_arrays.size() - 1u]), &pset_stats ); - + } catch (const std::exception& e) { - + printf_barry( "A problem ocurred while trying to add the array (and recording the powerset). " ); @@ -707,7 +707,7 @@ Model::add_arr printf_barry("Here is the array that generated the error.\n"); Array_.print(); throw std::logic_error(""); - + } // Recording the number of elements @@ -719,7 +719,7 @@ Model::add_arr pset_sizes.push_back(pset_arrays.back().size()); - + } else { @@ -727,7 +727,7 @@ Model::add_arr { support_fun.calc(); - + } catch (const std::exception& e) { @@ -742,14 +742,14 @@ Model::add_arr } } - + if (transform_model_fun) { auto tmpsupport = support_fun.get_counts(); size_t k = counter_fun.size(); size_t n = tmpsupport.size() / (k + 1); - std::vector< double > s_new(0u); + std::vector< double > s_new(0u); s_new.reserve(tmpsupport.size()); for (size_t i = 0u; i < n; ++i) @@ -766,13 +766,13 @@ Model::add_arr for (auto & s : s_new) stats_support.push_back(s); - + } else { for (const auto & s: support_fun.get_counts()) stats_support.push_back(s); } - + // Making room for the previous parameters. This will be used to check if // the normalizing constant has been updated or not. params_last.push_back(stats_target[0u]); @@ -783,32 +783,32 @@ Model::add_arr // Incrementing the size of the support set if (stats_support_sizes.size() == 0u) { - stats_support_sizes_acc.push_back(0u); + stats_support_sizes_acc.push_back(0u); } else { stats_support_sizes_acc.push_back( - stats_support_sizes.back() + + stats_support_sizes.back() + stats_support_sizes_acc.back() ); } stats_support_sizes.push_back( - + (stats_support.size() - stats_support_size)/ (counter_fun.size() + 1u) ); - + return arrays2support.size() - 1u; - + } - + // Increasing the number of arrays in that stat ++stats_support_n_arrays[locator->second]; - + // Adding the corresponding map arrays2support.push_back(locator->second); - + return arrays2support.size() - 1u; } @@ -820,7 +820,7 @@ inline double Model= arrays2support.size()) throw std::range_error("The requested support is out of range"); @@ -830,13 +830,13 @@ inline double Modelstats_support_sizes[idx] == 0u) return as_log ? -std::numeric_limits::infinity() : 0.0; - + // Checking if we have updated the normalizing constant or not if (!no_update_normconst && (!first_calc_done[idx] || !vec_equal_approx(params, params_last[idx]))) { - + first_calc_done[idx] = true; - + size_t k = params.size() + 1u; size_t n = stats_support_sizes[idx]; @@ -845,11 +845,11 @@ inline double Model @@ -868,7 +868,7 @@ inline double Model key = counters->gen_hash(Array_); MapVec_type< double, size_t >::const_iterator locator = keys2support.find(key); - if (locator == keys2support.end()) + if (locator == keys2support.end()) throw std::range_error( "This type of array has not been included in the model." ); @@ -900,12 +900,12 @@ inline double Modelstats_support_sizes[loc] == 0u) return as_log ? -std::numeric_limits::infinity() : 0.0; - + // Counting stats_target StatsCounter< Array_Type, Data_Counter_Type> tmpstats(&Array_); tmpstats.set_counters(this->counters); - + std::vector< double > target_ = tmpstats.count_all(); if (transform_model_fun) @@ -914,26 +914,26 @@ inline double Model::infinity() : 0.0; - + return likelihood_( &target_[0u], params, @@ -941,7 +941,7 @@ inline double Model @@ -952,7 +952,7 @@ inline double Model= arrays2support.size()) throw std::range_error("The requested support is out of range"); @@ -974,19 +974,19 @@ inline double Modelstats_support_sizes[loc] == 0u) { throw std::logic_error("The support set for this array is empty."); } - + // Checking if we have updated the normalizing constant or not if (!no_update_normconst && (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc])) ) { - + first_calc_done[loc] = true; - + size_t k = params.size() + 1u; size_t n = stats_support_sizes[loc]; @@ -995,11 +995,11 @@ inline double Model @@ -1018,7 +1018,7 @@ inline double Model= arrays2support.size()) throw std::range_error("The requested support is out of range"); @@ -1053,12 +1053,12 @@ inline double Model @@ -1089,7 +1089,7 @@ inline double Modelstats_support_n_arrays[i]); } else { - + res = 1.0; size_t stats_target_size = stats_target.size(); - #if defined(__OPENMP) || defined(_OPENMP) + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(*:res) #endif for (size_t i = 0; i < stats_target_size; ++i) @@ -1153,13 +1153,13 @@ inline double Model inline const std::vector< double > & Model:: get_normalizing_constants() const { - + return normalizing_constants; - + } template< @@ -1183,9 +1183,9 @@ template< > inline const std::vector< double > & Model:: get_likelihoods() const { - + return stats_likelihood; - + } template @@ -1218,7 +1218,7 @@ Model::get_pse template inline void Model:: print_stats(size_t i) const { - + if (i >= arrays2support.size()) throw std::range_error("The requested support is out of range"); @@ -1234,10 +1234,10 @@ inline void Model @@ -1272,11 +1272,11 @@ inline void Model(stat) > max_v) max_v = static_cast(stat); - + if (static_cast(stat) < min_v) min_v = static_cast(stat); - } + } // The vectors in the support reflec the size of nterms x entries // max_v /= static_cast(nterms() + 1); @@ -1288,13 +1288,13 @@ inline void Modelsize_unique()); printf_barry("Support size range : [%i, %i]\n", min_v, max_v); } - else + else { printf_barry("Num. of Arrays : 0\n"); printf_barry("Support size : -\n"); printf_barry("Support size range : -\n"); } - + if (with_pset) { @@ -1313,7 +1313,7 @@ inline void Modelnrules() > 0u) { printf_barry("Model rules (%li) :\n", this->nrules()); - + for (auto & rn : rules->get_names()) { printf_barry(" - %s\n", rn.c_str()); @@ -1323,7 +1323,7 @@ inline void Modelnrules_dyn() > 0u) { printf_barry("Model rules dyn (% 2li) :\n", this->nrules_dyn()); - + for (auto & rn : rules_dyn->get_names()) { printf_barry(" - %s\n", rn.c_str()); @@ -1349,12 +1349,12 @@ inline size_t Modelstats_support_sizes.size(); -} +} template inline size_t Model:: nterms() const noexcept { - + if (transform_model_fun) return transform_model_term_names.size(); else @@ -1365,7 +1365,7 @@ inline size_t Model inline size_t Model:: nrules() const noexcept { - + return this->rules->size(); } @@ -1373,7 +1373,7 @@ inline size_t Model inline size_t Model:: nrules_dyn() const noexcept { - + return this->rules_dyn->size(); } @@ -1395,14 +1395,14 @@ inline size_t Model inline std::vector< std::string > Model:: colnames() const { - + if (transform_model_fun) return transform_model_term_names; else return counters->get_names(); } - + template < typename Array_Type, typename Data_Counter_Type, @@ -1424,7 +1424,7 @@ Model::sample( // Getting the index size_t a = arrays2support[i]; - + // Generating a random std::uniform_real_distribution<> urand(0, 1); double r = urand(*rengine); @@ -1447,8 +1447,8 @@ Model::sample( if (j > 0u) j--; - } else { - + } else { + update_pset_probs(params, 1u, static_cast(a)); const double * probs = &pset_probs[pset_locations[a]]; @@ -1468,9 +1468,9 @@ Model::sample( std::string(" r: ") + std::to_string(r) ); #endif - + } - + return this->pset_arrays.at(a).at(j); } @@ -1499,37 +1499,37 @@ inline Array_Type Model s_new(0u); + std::vector< double > s_new(0u); s_new.reserve(tmpsupport.size()); for (size_t i = 0u; i < n; ++i) @@ -1582,7 +1582,7 @@ inline Array_Type Model urand(0, 1); // double r = urand(*rengine); @@ -1640,8 +1640,8 @@ inline Array_Type Model 0u) // j--; - // } else { - + // } else { + // // probs.resize(pset_arrays[a].size()); // std::vector< double > temp_stats(params.size()); // const double * stats = &pset_stats[pset_locations[a] * k]; @@ -1675,12 +1675,12 @@ inline Array_Type Modelpset_arrays.at(a).at(j); + // return this->pset_arrays.at(a).at(j); // #else - // return this->pset_arrays[a][j]; + // return this->pset_arrays[a][j]; // #endif } @@ -1715,7 +1715,7 @@ inline double Model @@ -1840,11 +1840,11 @@ Model::set_tra if (transform_model_fun) throw std::logic_error("A transformation function for the model has already been established."); - + transform_model_fun = fun; transform_model_term_names = names; - size_t k = counters->size(); + size_t k = counters->size(); auto stats_support_old = stats_support; @@ -1855,7 +1855,7 @@ Model::set_tra // How many observations in the support size_t n = stats_support_sizes[nsupport]; - // Iterating through each observation in the nsupport'th + // Iterating through each observation in the nsupport'th for (size_t i = 0; i < n; ++i) { @@ -1871,7 +1871,7 @@ Model::set_tra if (res.size() != transform_model_term_names.size()) throw std::length_error( std::string("The transform vector from -transform_model_fun- ") + - std::string(" does not match the size of ") + + std::string(" does not match the size of ") + std::string("-transform_model_term_names-.") ); @@ -1951,4 +1951,4 @@ Model::set_tra #undef MODEL_TEMPLATE_ARGS #undef MODEL_TYPE -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/defm/counters.hpp b/include/barry/models/defm/counters.hpp index 5401dcd..8dae52e 100644 --- a/include/barry/models/defm/counters.hpp +++ b/include/barry/models/defm/counters.hpp @@ -4,10 +4,10 @@ #include "formula.hpp" /** - * @ingroup counting + * @ingroup counting * @details Details on the available counters for `DEFMworkData` can be found in * the \ref counters-network section. - * + * */ ///@{ @@ -19,10 +19,10 @@ /** * @brief Data for the counters - * + * * @details This class is used to store the data for the counters. It is * used by the `Counters` class. - * + * */ #define MAKE_DEFM_HASHER(hasher,a,cov) \ barry::Hasher_fun_type \ @@ -82,7 +82,7 @@ barry::Rule_fun_type a = \ // ----------------------------------------------------------------------------- /** * @brief Prevalence of ones - * + * * @param counters Pointer ot a vector of counters * @param covar_index If >= than 0, then the interaction */ @@ -96,7 +96,7 @@ inline void counter_ones( // Weighted by a feature of the array if (covar_index >= 0) - { + { MAKE_DEFM_HASHER(hasher, array, covar_index) @@ -121,8 +121,8 @@ inline void counter_ones( counters->add_counter( counter_tmp, nullptr, hasher, - DEFMCounterData({static_cast(covar_index)}, {}, {}, true), - "Num. of ones x " + vname, + DEFMCounterData({static_cast(covar_index)}, {}, {}, true), + "Num. of ones x " + vname, "Overall number of ones" ); @@ -132,7 +132,7 @@ inline void counter_ones( DEFM_COUNTER_LAMBDA(count_ones) { - + // Only count the current if (i != (Array.nrow() - 1)) return 0.0; @@ -146,7 +146,7 @@ inline void counter_ones( counters->add_counter( count_ones, nullptr, nullptr, dat, // DEFMCounterData(), - "Num. of ones", + "Num. of ones", "Overall number of ones" ); } @@ -162,7 +162,7 @@ inline void counter_ones( * @param counters A pointer to the DEFMCounters object. * @param n_y The number of response variables. * @param which A vector of indices indicating which response variables to use. If empty, all response variables are used. - * @param covar_index The index of the covariate to use as the intercept. + * @param covar_index The index of the covariate to use as the intercept. * @param vname The name of the variable to use as the intercept. If empty, the intercept is set to zero. * @param x_names A pointer to a vector of strings containing the names of the covariates. * @param y_names A pointer to a vector of strings containing the names of the response variables. @@ -213,8 +213,8 @@ inline void counter_logit_intercept( counters->add_counter( tmp_counter, nullptr, nullptr, - DEFMCounterData({i}, {}, {}, false), - "Logit intercept " + vname, + DEFMCounterData({i}, {}, {}, false), + "Logit intercept " + vname, "Equal to one if the outcome " + vname + " is one. Equivalent to the logistic regression intercept." ); @@ -256,8 +256,8 @@ inline void counter_logit_intercept( if (hasher_added) counters->add_counter( tmp_counter, nullptr, nullptr, - DEFMCounterData({i, static_cast(covar_index)}, {}, {}, false), - "Logit intercept " + yname + " x " + vname, + DEFMCounterData({i, static_cast(covar_index)}, {}, {}, false), + "Logit intercept " + yname + " x " + vname, "Equal to one if the outcome " + yname + " is one. Equivalent to the logistic regression intercept." ); else { @@ -266,8 +266,8 @@ inline void counter_logit_intercept( counters->add_counter( tmp_counter, nullptr, hasher, - DEFMCounterData({i, static_cast(covar_index)}, {}, {}, false), - "Logit intercept " + yname + " x " + vname, + DEFMCounterData({i, static_cast(covar_index)}, {}, {}, false), + "Logit intercept " + yname + " x " + vname, "Equal to one if the outcome " + yname + " is one. Equivalent to the logistic regression intercept." ); @@ -276,13 +276,13 @@ inline void counter_logit_intercept( } } - + } /** * @brief Prevalence of ones - * + * * @param counters Pointer ot a vector of counters * @param covar_index If >= than 0, then the interaction */ @@ -335,20 +335,20 @@ inline void counter_transition( if (sgn[k] && (cellv != 1)) return 0.0; } - + // If nothing happens, then is one or the covaridx return (covaridx < 1000) ? Array.D()(Array.nrow() - 1u, covaridx) : 1.0; - + }; DEFM_COUNTER_LAMBDA(count_ones) { - + auto dat = data.indices; auto sgn = data.logical; int covaridx = dat[dat.size() - 1u]; - // Checking if the observation is in the stat. We + // Checking if the observation is in the stat. We const auto & array = Array.get_data(); size_t loc = i + j * Array.nrow(); size_t n_cells = dat.size() - 1u; @@ -370,13 +370,13 @@ inline void counter_transition( if ((sgn[e] && (array[dat[e]] == 1)) || (!sgn[e] && (array[dat[e]] == 0))) n_now++; - + } // If i in array still false, then no change if (!i_in_array) return 0.0; - + size_t n_prev = n_now; if (baseline_value) n_prev--; @@ -386,14 +386,14 @@ inline void counter_transition( // Computing stats if (covaridx < 1000) { - + double val = Array.D()(Array.nrow() - 1u, covaridx); double value_now = n_now == n_cells ? val : 0.0; double value_prev = n_prev == n_cells ? val : 0.0; return value_now - value_prev; - } + } else { @@ -424,12 +424,12 @@ inline void counter_transition( size_t c = std::floor(coords[d] / (m_order + 1u)); size_t r = coords[d] - c * (m_order + 1u); motif(r, c) = signs[d] ? 1 : -1; - + } // Checking if any prior to the event bool any_before_event = false; - + for (size_t i = 0u; i < m_order; ++i) { for (size_t j = 0u; j < n_y; ++j) @@ -442,7 +442,7 @@ inline void counter_transition( } } - + #ifdef BARRY_WITH_LATEX name += "$"; #endif @@ -581,8 +581,8 @@ inline void counter_transition( counters->add_counter( count_ones, count_init, hasher, - DEFMCounterData(coords, {}, signs, coords.size() > 1u ? true : false), - name + " x " + vname, + DEFMCounterData(coords, {}, signs, coords.size() > 1u ? true : false), + name + " x " + vname, "Motif weighted by single attribute" ); @@ -590,13 +590,13 @@ inline void counter_transition( counters->add_counter( count_ones, count_init, nullptr, - DEFMCounterData(coords, {}, signs, coords.size() > 1u ? true : false), - name, + DEFMCounterData(coords, {}, signs, coords.size() > 1u ? true : false), + name, "Motif" ); } - + return; @@ -604,7 +604,7 @@ inline void counter_transition( /** * @brief Prevalence of ones - * + * * @param counters Pointer ot a vector of counters * @param covar_index If >= than 0, then the interaction */ @@ -671,7 +671,7 @@ inline void counter_transition_formula( ); } - else + else { counter_transition( @@ -686,7 +686,7 @@ inline void counter_transition_formula( /** * @brief Prevalence of ones - * + * * @param counters Pointer ot a vector of counters * @param covar_index If >= than 0, then the interaction */ @@ -718,7 +718,7 @@ inline void counter_fixed_effect( counters->add_counter( count_tmp, count_init, hasher, - DEFMCounterData({static_cast(covar_index)}, {k}, {}), + DEFMCounterData({static_cast(covar_index)}, {k}, {}), "Fixed effect feature (" + vname + ")^" + std::to_string(k) ); @@ -737,32 +737,32 @@ inline void rules_markov_fixed( DEFMRules * rules, size_t markov_order ) { - + DEFM_RULE_LAMBDA(no_self_tie) { return i >= data.idx(0u); }; - + rules->add_rule( no_self_tie, DEFMRuleData({},{markov_order}), std::string("Markov model of order ") + std::to_string(markov_order), std::string("Blocks the first morder cells of the array.") ); - + return; } /** * @brief Blocks switching a one to zero. - * - * @param rules + * + * @param rules * @param ids Ids of the variables that will follow this rule. */ inline void rules_dont_become_zero( DEFMSupport * support, std::vector ids ) { - + DEFM_RULE_LAMBDA(rule) { if (!data.init) @@ -797,14 +797,14 @@ inline void rules_dont_become_zero( return (Array(i - 1, j) != 1) || (Array(i, j) != 1); }; - + support->get_rules()->add_rule( rule, DEFMRuleData({}, {ids}), std::string("Ones can't become zero"), std::string("Blocks cells that have became equal to one.") ); - + return; } @@ -814,7 +814,7 @@ inline void rules_dont_become_zero( * @param pos Position of the focal statistic. * @param lb Lower bound * @param ub Upper bound - * @details + * @details * @return (void) adds a rule limiting the support of the model. */ inline void rule_constrain_support( @@ -824,7 +824,7 @@ inline void rule_constrain_support( double ub ) { - + DEFM_RULEDYN_LAMBDA(tmp_rule) { @@ -834,10 +834,10 @@ inline void rule_constrain_support( return false; else return true; - + }; - + support->get_rules_dyn()->add_rule( tmp_rule, DEFMRuleDynData( @@ -847,14 +847,14 @@ inline void rule_constrain_support( support->get_counters()->get_names()[pos] + "' within [" + std::to_string(lb) + ", " + std::to_string(ub) + std::string("]"), - std::string("When the support is ennumerated, only states where the statistic '") + + std::string("When the support is ennumerated, only states where the statistic '") + support->get_counters()->get_names()[pos] + std::to_string(pos) + "' falls within [" + std::to_string(lb) + ", " + std::to_string(ub) + "] are included." ); - + return; - + } diff --git a/include/barry/models/defm/defm-bones.hpp b/include/barry/models/defm/defm-bones.hpp index abf0c38..8a41149 100644 --- a/include/barry/models/defm/defm-bones.hpp +++ b/include/barry/models/defm/defm-bones.hpp @@ -17,7 +17,7 @@ class DEFM : public DEFMModel { std::shared_ptr> Y_shared; ///< Outcome variable std::shared_ptr> ID_shared; ///< Individual ids std::shared_ptr> X_shared; ///< Covariates - + size_t N; ///< Number of agents/individuals size_t ID_length; ///< Length of the vector IDs size_t Y_ncol; ///< Number of columns in the response @@ -91,4 +91,3 @@ class DEFM : public DEFMModel { }; #endif - diff --git a/include/barry/models/defm/defm-meat.hpp b/include/barry/models/defm/defm-meat.hpp index 22b982a..b5dd94a 100644 --- a/include/barry/models/defm/defm-meat.hpp +++ b/include/barry/models/defm/defm-meat.hpp @@ -5,7 +5,7 @@ inline std::vector< double > keygen_defm( const DEFMArray & Array_, DEFMCounterData * data ) { - + size_t nrow = Array_.nrow(); size_t ncol = Array_.ncol(); @@ -41,7 +41,7 @@ inline void DEFM::simulate( int * y_out ) { - size_t model_num = 0u; + size_t model_num = 0u; size_t n_entry = M_order * Y_ncol; auto idx = this->get_arrays2support(); DEFMArray last_array; @@ -50,7 +50,7 @@ inline void DEFM::simulate( // Figuring out how many processes can we observe DEFM_RANGES(i) - + DEFM_LOOP_ARRAYS(proc_n) { @@ -96,7 +96,7 @@ inline void DEFM::simulate( } - + } n_entry += M_order * Y_ncol; @@ -196,7 +196,7 @@ inline DEFM::DEFM( start_end.push_back(row); } - + } start_end.push_back(id_length - 1u); @@ -210,12 +210,12 @@ inline DEFM::DEFM( for (auto i = 0u; i < X_ncol; ++i) X_names.emplace_back(std::string("X") + std::to_string(i)); - return; + return; } -inline void DEFM::init() +inline void DEFM::init() { // Adding the rule @@ -341,7 +341,7 @@ inline barry::FreqTable DEFM::motif_census( // Figuring out how many processes can we observe DEFM_RANGES(i) - + DEFM_LOOP_ARRAYS(proc_n) { @@ -367,7 +367,7 @@ inline std::vector< double > DEFM::logodds( size_t i_, size_t j_ ) { - + std::vector< double > res(ID_length, std::nan("")); @@ -376,7 +376,7 @@ inline std::vector< double > DEFM::logodds( // Figuring out how many processes can we observe DEFM_RANGES(i) - + DEFM_LOOP_ARRAYS(n_proc) { @@ -463,4 +463,4 @@ inline bool DEFM::get_column_major() const noexcept #undef DEFM_RANGES #undef DEFM_LOOP_ARRAYS -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/defm/defm-types.hpp b/include/barry/models/defm/defm-types.hpp index d37de1d..e70a7d9 100644 --- a/include/barry/models/defm/defm-types.hpp +++ b/include/barry/models/defm/defm-types.hpp @@ -6,16 +6,16 @@ typedef barry::BArrayDense DEFMArray; /** * @brief Data class for DEFM arrays. - * + * * This holds information pointing to the data array, including information * regarding the number of observations, the time slices of the observation, * and the number of covariates in the data. - * + * */ class DEFMData { public: - + DEFMArray * array; // Pointer to the owner of this data const double * covariates; ///< Vector of covariates (complete vector) size_t obs_start; ///< Index of the observation in the data. @@ -24,9 +24,9 @@ class DEFMData { std::vector< size_t > covar_sort; /// Value where the sorting of the covariates is stored. std::vector< size_t > covar_used; /// Vector indicating which covariates are included in the model bool column_major; - + DEFMData() {}; - + /** * @brief Constructor * @param covariates_ Pointer to the attribute data. @@ -42,21 +42,21 @@ class DEFMData { size_t X_nrow_, bool column_major_ ) : array(array_), covariates(covariates_), obs_start(obs_start_), - X_ncol(X_ncol_), X_nrow(X_nrow_), column_major(column_major_) {}; + X_ncol(X_ncol_), X_nrow(X_nrow_), column_major(column_major_) {}; /** * @brief Access to the row (i) colum (j) data - * - * @param i - * @param j - * @return double + * + * @param i + * @param j + * @return double */ double operator()(size_t i, size_t j) const; double at(size_t i, size_t j) const; size_t ncol() const; size_t nrow() const; void print() const; - + ~DEFMData() {}; }; @@ -70,22 +70,22 @@ class DEFMCounterData { std::vector< double > numbers; std::vector< bool > logical; bool is_motif; ///< If false, then is a logit intercept. - + DEFMCounterData() : indices(0u), numbers(0u), logical(0u), is_motif(true) {}; DEFMCounterData( const std::vector< size_t > indices_, const std::vector< double > numbers_, const std::vector< bool > logical_, bool is_motif_ = true - ): indices(indices_), numbers(numbers_), + ): indices(indices_), numbers(numbers_), logical(logical_), is_motif(is_motif_) {}; size_t idx(size_t i) const {return indices[i];}; double num(size_t i) const {return numbers[i];}; bool is_true(size_t i) const {return logical[i];}; - + ~DEFMCounterData() {}; - + }; class DEFMRuleData { @@ -144,7 +144,7 @@ inline void DEFMData::print() const { for (size_t j = 0u; j < X_ncol; ++j) printf_barry("% 5.2f, ", operator()(i, j)); printf_barry("\n"); - + } } @@ -162,21 +162,21 @@ class DEFMRuleDynData { size_t pos; size_t lb; size_t ub; - + DEFMRuleDynData( const std::vector< double > * counts_, size_t pos_, size_t lb_, size_t ub_ ) : counts(counts_), pos(pos_), lb(lb_), ub(ub_) {}; - + ~DEFMRuleDynData() {}; const double operator()() const { return (*counts)[pos]; } - + }; /** @@ -196,4 +196,4 @@ typedef barry::Rule DEFMRuleDyn; typedef barry::Rules DEFMRulesDyn; ///@} -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/defm/formula.hpp b/include/barry/models/defm/formula.hpp index 4c7af40..9ae7e5f 100644 --- a/include/barry/models/defm/formula.hpp +++ b/include/barry/models/defm/formula.hpp @@ -2,46 +2,46 @@ #define BARRY_DEFM_MOTIF_FORMULA_HPP /** * @brief Parses a motif formula - * + * * @details This function will take the formula and generate the corresponding * input for defm::counter_transition(). Formulas can be specified in the * following ways: - * + * * - Intercept effect: {...} No transition, only including the current state. * - Transition effect: {...} > {...} Includes current and previous states. - * + * * The general notation is `[0]y[column id]_[row id]`. A preceeding zero * means that the value of the cell is considered to be zero. The column * id goes between 0 and the number of columns in the array - 1 (so it * is indexed from 0,) and the row id goes from 0 to m_order. - * + * * ## Intercept effects - * + * * Intercept effects only involve a single set of curly brackets. Using the * 'greater-than' symbol (i.e., '>') is only for transition effects. When * specifying intercept effects, users can skip the `row_id`, e.g., * `y0_0` is equivalent to `y0`. If the passed `row id` is different from * the Markov order, i.e., `row_id != m_order`, then the function returns - * with an error. - * + * with an error. + * * Examples: - * + * * - `"{y0, 0y1}"` is equivalent to set a motif with the first element equal - * to one and the second to zero. - * + * to one and the second to zero. + * * ## Transition effects - * + * * Transition effects can be specified using two sets of curly brackets and * an greater-than symbol, i.e., `{...} > {...}`. The first set of brackets, * which we call LHS, can only hold `row id` that are less than `m_order`. - * - * - * - * @param formula - * @param locations - * @param signs - * @param m_order - * @param y_ncol + * + * + * + * @param formula + * @param locations + * @param signs + * @param m_order + * @param y_ncol */ inline void defm_motif_parser( std::string formula, @@ -100,7 +100,7 @@ inline void defm_motif_parser( // Will indicate where the arrow is located at size_t arrow_position = match.position(4u); - // This pattern will match + // This pattern will match std::regex pattern("(0?)y([0-9]+)(_([0-9]+))?"); auto iter = std::sregex_iterator(formula.begin(), formula.end(), pattern); @@ -126,7 +126,7 @@ inline void defm_motif_parser( std::string tmp_str = i->operator[](4u).str(); if (m_order > 1) { - // If missing, we replace with the location + // If missing, we replace with the location if (tmp_str == "") { @@ -144,7 +144,7 @@ inline void defm_motif_parser( } else { - // If missing, we replace with the location + // If missing, we replace with the location if (tmp_str != "") y_row = std::stoul(tmp_str); else @@ -168,14 +168,14 @@ inline void defm_motif_parser( locations.push_back(y_col * (m_order + 1) + y_row); signs.push_back(is_positive); - + } return; - } - + } + std::regex_match(formula, match, pattern_intercept); if (!match.empty()) { @@ -195,14 +195,14 @@ inline void defm_motif_parser( vname = vname.substr(1, vname.size() - 2); } - // This pattern will match + // This pattern will match std::regex pattern("(0?)y([0-9]+)(_([0-9]+))?"); auto iter = std::sregex_iterator(formula.begin(), formula.end(), pattern); for (auto i = iter; i != empty; ++i) { - + // First value true/false bool is_positive = true; if (i->operator[](1u).str() == "0") @@ -237,19 +237,18 @@ inline void defm_motif_parser( locations.push_back(y_col * (m_order + 1) + y_row); signs.push_back(is_positive); - + } return; - } - + } + throw std::logic_error( "The motif specified in the formula: " + formula + " has the wrong syntax." ); - + } #endif - diff --git a/include/barry/models/geese.hpp b/include/barry/models/geese.hpp index 7b4721e..bd0d89e 100644 --- a/include/barry/models/geese.hpp +++ b/include/barry/models/geese.hpp @@ -21,7 +21,7 @@ namespace geese { #include "geese/flock-bones.hpp" #include "geese/flock-meat.hpp" - + #include "geese/counters.hpp" } diff --git a/include/barry/models/geese/counters.hpp b/include/barry/models/geese/counters.hpp index 772ecca..ccdf095 100644 --- a/include/barry/models/geese/counters.hpp +++ b/include/barry/models/geese/counters.hpp @@ -25,13 +25,13 @@ /** * @brief Extension of a simple counter. - * + * * It allows specifying extra arguments, in particular, the corresponding * sets of rows to which this statistic may be relevant. This could be important * in the case of, for example, counting correlation type statistics between * function 1 and 2, and between function 1 and 3. - * - * + * + * */ #define PHYLO_RULE_LAMBDA(a) barry::Rule_fun_type a = \ [](const PhyloArray & Array, size_t i, size_t j, PhyloRuleData & data) @@ -44,7 +44,7 @@ #define PHYLO_CHECK_MISSING() if (Array.D_ptr() == nullptr) \ throw std::logic_error("The array data is nullptr."); \ - + inline std::string get_last_name(size_t d) {return ((d == 1u)? " at duplication" : ((d == 0u)? " at speciation" : ""));} /** @@ -63,12 +63,12 @@ inline void counter_overall_gains( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { PHYLO_CHECK_MISSING(); - + return 0.0; }; @@ -77,11 +77,11 @@ inline void counter_overall_gains( { IF_NOTMATCHES() return 0.0; - + return Array.D_ptr()->states[i] ? 0.0 : 1.0; - + }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), @@ -89,7 +89,7 @@ inline void counter_overall_gains( ); return; - + } // ----------------------------------------------------------------------------- @@ -102,7 +102,7 @@ inline void counter_gains( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -121,7 +121,7 @@ inline void counter_gains( if (Array(k, o) == 1u) ngains += 1.0; } - + return ngains; }; @@ -135,20 +135,20 @@ inline void counter_gains( IF_MATCHES() return (i == data[1u]) ? 1.0 : 0.0; - + return 0.0; }; - + for (auto& i : nfun) counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, i}), "Gains " + std::to_string(i) + get_last_name(duplication) ); - + return; - + } @@ -163,7 +163,7 @@ inline void counter_gains_k_offspring( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -204,13 +204,13 @@ inline void counter_gains_k_offspring( else if (diff == 0) { return 1.0; - } else + } else // (c) Otherwise, nothing happens return 0.0; - + }; - + for (auto& i : nfun) counters->add_counter( tmp_count, tmp_init, nullptr, @@ -218,9 +218,9 @@ inline void counter_gains_k_offspring( std::to_string(k) + " genes gain " + std::to_string(i) + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -233,10 +233,10 @@ inline void counter_genes_changing( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { - + PHYLO_CHECK_MISSING(); IF_NOTMATCHES() @@ -247,14 +247,14 @@ inline void counter_genes_changing( for (auto s : Array.D_ptr()->states) { - if (s) + if (s) // Yup, we are loosing a function, so break return static_cast(Array.ncol()); - + } return 0.0; - + }; @@ -272,25 +272,25 @@ inline void counter_genes_changing( // Nah, this gene was already different. if ((k != i) && (Array.D_ptr()->states[k] != (Array(k, j, false) == 1u))) return 0.0; - + } - // Nope, this gene is now matching its parent, so we need to + // Nope, this gene is now matching its parent, so we need to // take it out from the count of genes that have changed. return Array.D_ptr()->states[i] ? -1.0 : 1.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), "Num. of genes changing" + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -304,10 +304,10 @@ inline void counter_preserve_pseudogene( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { - + PHYLO_CHECK_MISSING(); IF_NOTMATCHES() @@ -320,7 +320,7 @@ inline void counter_preserve_pseudogene( double n = static_cast(Array.ncol()); return n * (n - 1.0) / 2.0; - + }; @@ -359,18 +359,18 @@ inline void counter_preserve_pseudogene( return res; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), - "Preserve pseudo gene (" + + "Preserve pseudo gene (" + std::to_string(nfunA) + ", " + std::to_string(nfunB) + ")" + get_last_name(duplication) ); - + return; - + } @@ -384,10 +384,10 @@ inline void counter_prop_genes_changing( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { - + PHYLO_CHECK_MISSING(); IF_NOTMATCHES() @@ -400,9 +400,9 @@ inline void counter_prop_genes_changing( if (s) return 1.0; } - + return 0.0; - + }; PHYLO_COUNTER_LAMBDA(tmp_count) @@ -411,7 +411,7 @@ inline void counter_prop_genes_changing( // Checking the type of event IF_NOTMATCHES() return 0.0; - + // Setup bool j_diverges = false; const std::vector< bool > & par_state = Array.D_ptr()->states; @@ -465,16 +465,16 @@ inline void counter_prop_genes_changing( return -1.0/Array.ncol(); }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), "Proportion of genes changing" + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -486,7 +486,7 @@ inline void counter_overall_loss( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -495,17 +495,17 @@ inline void counter_overall_loss( IF_MATCHES() return -1.0; - else + else return 0.0; - + }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { IF_NOTMATCHES() return 0.0; - + double res = 0.0; for (auto s : Array.D_ptr()->states) if (s) @@ -514,13 +514,13 @@ inline void counter_overall_loss( return res * static_cast(Array.ncol()); }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), "Overall loses" + get_last_name(duplication) ); - + return; } @@ -540,7 +540,7 @@ inline void counter_maxfuns( PHYLO_COUNTER_LAMBDA(tmp_init) { - PHYLO_CHECK_MISSING(); + PHYLO_CHECK_MISSING(); IF_NOTMATCHES() return 0.0; @@ -549,11 +549,11 @@ inline void counter_maxfuns( // bound is zero if (data[1u] == 0) return static_cast(Array.ncol()); - + return 0.0; }; - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -562,7 +562,7 @@ inline void counter_maxfuns( int count = Array.colsum(j); int ub = data[2u]; - + // It now matches if (count == static_cast(data[1u])) return 1.0; @@ -582,11 +582,11 @@ inline void counter_maxfuns( "Genes with [" + std::to_string(lb) + ", " + std::to_string(ub) + "] funs" + get_last_name(duplication) ); - + return; - + } - + // ----------------------------------------------------------------------------- /** * @brief Total count of losses for an specific function. @@ -597,7 +597,7 @@ inline void counter_loss( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -606,11 +606,11 @@ inline void counter_loss( if (!Array.D_ptr()->states[i]) return 0.0; - + return (i == data[1u]) ? -1.0 : 0.0; }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -618,25 +618,25 @@ inline void counter_loss( IF_NOTMATCHES() return 0.0; - + auto f = data[1u]; if (!Array.D_ptr()->states[f]) return 0.0; - + return static_cast(Array.ncol()); }; - + for (auto& i : nfun) counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, i}), "Loss " + std::to_string(i) + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -648,7 +648,7 @@ inline void counter_overall_changes( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -657,11 +657,11 @@ inline void counter_overall_changes( if (Array.D_ptr()->states[i]) return -1.0; - else + else return 1.0; }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -681,7 +681,7 @@ inline void counter_overall_changes( return counts; - + }; @@ -690,10 +690,10 @@ inline void counter_overall_changes( PhyloCounterData({duplication}), "Overall changes" + get_last_name(duplication) ); - - + + return; - + } @@ -709,7 +709,7 @@ inline void counter_subfun( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -719,50 +719,50 @@ inline void counter_subfun( auto funA = data[1u]; auto funB = data[2u]; - + // Are we looking at either of the relevant functions? if ((funA != i) && (funB != i)) return 0.0; - + // Are A and B existant? if not, no change if (!Array.D_ptr()->states[funA] || !Array.D_ptr()->states[funB]) return 0.0; - + // Figuring out which is the first (reference) function size_t other = (i == funA)? funB : funA; double res = 0.0; // There are 4 cases: (first x second) x (had the second function) if (Array(other, j, false) == 1u) - { - + { + for (size_t off = 0u; off < Array.ncol(); ++off) { - + // Not on self if (off == j) continue; - + if ((Array(i, off, false) == 1u) && (Array(other, off, false) == 0u)) res -= 1.0; - + } - + } else { - + for (size_t off = 0u; off < Array.ncol(); ++off) { - + // Not on self if (off == j) continue; - + if ((Array(i, off, false) == 0u) && (Array(other, off, false) == 1u)) res += 1.0; - + } - + } - + return res; }; @@ -774,16 +774,16 @@ inline void counter_subfun( return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, nfunA, nfunB}), "Subfun between " + std::to_string(nfunA) + " and " + std::to_string(nfunB) + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -798,7 +798,7 @@ inline void counter_cogain( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -807,11 +807,11 @@ inline void counter_cogain( auto d1 = data[1u]; auto d2 = data[2u]; - + // Is the function in scope relevant? if ((i != d1) && (i != d2)) return 0.0; - + // None should have it if (!Array.D_ptr()->states[d1] && !Array.D_ptr()->states[d2]) { @@ -823,7 +823,7 @@ inline void counter_cogain( else return 0.0; - } else + } else return 0.0; }; @@ -834,16 +834,16 @@ inline void counter_cogain( return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, nfunA, nfunB}), "Co-gains " + std::to_string(nfunA) + " & " + std::to_string(nfunB) + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -853,7 +853,7 @@ inline void counter_longest( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -867,7 +867,7 @@ inline void counter_longest( int nmutate_longest = 0; auto states = Array.D_ptr()->states; - + for (auto off = 0u; off < Array.ncol(); ++off) { @@ -880,7 +880,7 @@ inline void counter_longest( { if ((Array(f, off) == 1u) != states[f]) { - + // If it happens that j != off and is not longest // then return 0 (a not longest was mutating prev) if (is_longest[off] && (off != j)) @@ -942,32 +942,32 @@ inline void counter_longest( nmutate_prev++; } - + // Just compute the change statistic directly return ( ((nmutate == 0) & (nmutate_longest > 0)) ? 1.0 : 0.0 ) + ( ((nmutate_prev == 0) & (nmutate_longest_prev > 0)) ? 1.0 : 0.0 ); }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { PHYLO_CHECK_MISSING(); - + if (Array.D_ptr()->blengths.size() != Array.ncol()) throw std::logic_error( "longest should be initialized with a vec of size Array.ncol()." ); - + // Finding the longest branch (or branches) -- size_t longest_idx = 0u; double diff = 0.0; - data.reserve(Array.ncol()); + data.reserve(Array.ncol()); data.push_back(0u); for (size_t ii = 1u; ii < Array.ncol(); ++ii) { - + diff = Array.D_ptr()->blengths[longest_idx] - Array.D_ptr()->blengths[ii]; if (diff > 0.0) continue; @@ -981,41 +981,41 @@ inline void counter_longest( } else if (diff == 0.0) data.push_back(ii); - + } data.shrink_to_fit(); - + if (data.size() == 0u) throw std::logic_error("The data on the longest branch has size 0."); - + // Starting the counter, since all in zero, then this will be equal to // the number of functions in 1 x number of longest branches for (size_t ii = 0u; ii < Array.nrow(); ++ii) { - + if (Array.D_ptr()->states[ii]) return (1.0 * static_cast(data.size())); } - + return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), "Longest branch mutates" + get_last_name(duplication) ); - + return; - + } //------------------------------------------------------------------------------ /** - * @brief Total number of neofunctionalization events + * @brief Total number of neofunctionalization events * @details Needs to specify pairs of function. */ inline void counter_neofun( @@ -1025,34 +1025,34 @@ inline void counter_neofun( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { // Is this node duplication? IF_NOTMATCHES() return 0.0; - + auto funA = data[1u]; auto funB = data[2u]; // Is the function in scope relevant? if ((i != funA) && (i != funB)) return 0.0; - + // Checking if the parent has both functions size_t other = (i == funA)? funB : funA; bool parent_i = Array.D_ptr()->states[i]; bool parent_other = Array.D_ptr()->states[other]; - - if (!parent_i & !parent_other) + + if (!parent_i & !parent_other) return 0.0; - else if (parent_i & parent_other) + else if (parent_i & parent_other) return 0.0; - + // Figuring out which is the first (reference) function double res = 0.0; - + if (Array(other, j) == 0u) { @@ -1068,13 +1068,13 @@ inline void counter_neofun( for (auto off = 0u; off < Array.ncol(); ++off) if ((Array(i,off) == 1) && (Array(other,off) == 0)) res -= 1.0; - + } - + return res; }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { PHYLO_CHECK_MISSING(); @@ -1088,14 +1088,14 @@ inline void counter_neofun( "Neofun between " + std::to_string(nfunA) + " and " + std::to_string(nfunB) + get_last_name(duplication) ); - + return; - + } //------------------------------------------------------------------------------ /** - * @brief Total number of neofunctionalization events + * @brief Total number of neofunctionalization events * sum_u sum_{w < u} [x(u,a)*(1 - x(w,a)) + (1 - x(u,a)) * x(w,a)] * change stat: delta{x(u,a): 0->1} = 1 - 2 * x(w,a) */ @@ -1105,22 +1105,22 @@ inline void counter_pairwise_neofun_singlefun( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { // Is this node duplication? IF_NOTMATCHES() return 0.0; - + // Is the function in scope relevant? if (i != data[1u]) return 0.0; - + // Checking if the parent has the function if (Array.D_ptr()->states[i]) return 0.0; - + // Figuring out which is the first (reference) function double res = 0.0; for (auto off = 0u; off < Array.ncol(); ++off) @@ -1131,7 +1131,7 @@ inline void counter_pairwise_neofun_singlefun( if ((Array(i, off) == 0)) res += 1.0; - else + else res -= 1.0; } @@ -1139,7 +1139,7 @@ inline void counter_pairwise_neofun_singlefun( return res; }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { PHYLO_CHECK_MISSING(); @@ -1153,14 +1153,14 @@ inline void counter_pairwise_neofun_singlefun( "Pairwise neofun function " + std::to_string(nfunA) + get_last_name(duplication) ); - + return; - + } //------------------------------------------------------------------------------ /** - * @brief Total number of neofunctionalization events + * @brief Total number of neofunctionalization events * @details Needs to specify pairs of function. */ inline void counter_neofun_a2b( @@ -1170,23 +1170,23 @@ inline void counter_neofun_a2b( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { // Is this node duplication? IF_NOTMATCHES() return 0.0; - + const size_t & funA = data[1u]; const size_t & funB = data[2u]; // Checking scope if ((i != funA) && (i != funB)) return 0.0; - + // Checking the parent doesn't have funA or has funB - if (!Array.D_ptr()->states[funA] || Array.D_ptr()->states[funB]) + if (!Array.D_ptr()->states[funA] || Array.D_ptr()->states[funB]) return 0.0; double res = 0.0; @@ -1214,7 +1214,7 @@ inline void counter_neofun_a2b( for (auto off = 0u; off < Array.ncol(); ++off) { - + if (off == j) continue; @@ -1249,7 +1249,7 @@ inline void counter_neofun_a2b( for (auto off = 0u; off < Array.ncol(); ++off) { - + if (off == j) continue; @@ -1263,9 +1263,9 @@ inline void counter_neofun_a2b( } return res; - + }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -1280,9 +1280,9 @@ inline void counter_neofun_a2b( "Neofun from " + std::to_string(nfunA) + " to " + std::to_string(nfunB) + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -1299,17 +1299,17 @@ inline void counter_neofun_a2b( inline void counter_co_opt( PhyloCounters * counters, size_t nfunA, - size_t nfunB, + size_t nfunB, size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) - { + { // Checking whether this is for duplication or not IF_NOTMATCHES() return 0.0; - + const size_t funA = data[1u]; const size_t funB = data[2u]; @@ -1326,7 +1326,7 @@ inline void counter_co_opt( // What was the state of the other function? If B is present, then // nothing changes. - if (Array(funB, j, false) == 1u) + if (Array(funB, j, false) == 1u) return 0.0; // Iterating through the sibs @@ -1341,7 +1341,7 @@ inline void counter_co_opt( // What was the state of the other function? If A is not present, then // nothing changes. - if (Array(funA, j, false) == 0u) + if (Array(funA, j, false) == 0u) return 0.0; // Iterating through the sibs @@ -1354,10 +1354,10 @@ inline void counter_co_opt( } - + }; - + PHYLO_COUNTER_LAMBDA(tmp_init) { PHYLO_CHECK_MISSING(); @@ -1383,10 +1383,10 @@ inline void counter_co_opt( "Coopt of " + std::to_string(nfunA) + " by " + std::to_string(nfunB) + get_last_name(duplication) ); - - + + return; - + } // ----------------------------------------------------------------------------- @@ -1400,10 +1400,10 @@ inline void counter_k_genes_changing( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { - + PHYLO_CHECK_MISSING(); IF_NOTMATCHES() @@ -1416,7 +1416,7 @@ inline void counter_k_genes_changing( return Array.ncol() == data[1u] ? 1.0 : 0.0; return data[1u] == 0 ? 1.0 : 0.0; - + }; PHYLO_COUNTER_LAMBDA(tmp_count) @@ -1425,9 +1425,9 @@ inline void counter_k_genes_changing( // Checking the type of event IF_NOTMATCHES() return 0.0; - + // How many genes diverge the parent - int count = 0; + int count = 0; bool j_diverges = false; const auto & par_state = Array.D_ptr()->states; @@ -1500,13 +1500,13 @@ inline void counter_k_genes_changing( return (count == k ? 1.0 : 0.0) - (count_prev == k ? 1.0 : 0.0); }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, k}), std::to_string(k) + " genes changing" + get_last_name(duplication) ); - + } // ----------------------------------------------------------------------------- @@ -1520,10 +1520,10 @@ inline void counter_less_than_p_prop_genes_changing( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_init) { - + PHYLO_CHECK_MISSING(); IF_NOTMATCHES() @@ -1535,7 +1535,7 @@ inline void counter_less_than_p_prop_genes_changing( // Only one if it was specified it was zero return 1.0; - + }; PHYLO_COUNTER_LAMBDA(tmp_count) @@ -1544,7 +1544,7 @@ inline void counter_less_than_p_prop_genes_changing( // Checking the type of event IF_NOTMATCHES() return 0.0; - + // Setup double count = 0.0; ///< How many genes diverge the parent @@ -1616,13 +1616,13 @@ inline void counter_less_than_p_prop_genes_changing( return ((count/ncol) <= p ? 1.0 : 0.0) - ((count_prev/ncol) <= p ? 1.0 : 0.0); }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, static_cast(p * 100)}), std::to_string(p) + " prop genes changing" + get_last_name(duplication) ); - + } // ----------------------------------------------------------------------------- @@ -1636,7 +1636,7 @@ inline void counter_gains_from_0( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -1669,7 +1669,7 @@ inline void counter_gains_from_0( return res; - + }; PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -1678,7 +1678,7 @@ inline void counter_gains_from_0( return 0.0; }; - + for (auto& i : nfun) counters->add_counter( tmp_count, tmp_init, nullptr, @@ -1686,9 +1686,9 @@ inline void counter_gains_from_0( "First gain " + std::to_string(i) + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -1701,7 +1701,7 @@ inline void counter_overall_gains_from_0( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -1718,7 +1718,7 @@ inline void counter_overall_gains_from_0( } return 1.0; - + }; PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -1727,16 +1727,16 @@ inline void counter_overall_gains_from_0( return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), "Overall first gains" + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -1749,7 +1749,7 @@ inline void counter_pairwise_overall_change( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -1770,9 +1770,9 @@ inline void counter_pairwise_overall_change( else if (funpar < Array(i, off)) res += 1.0; } - + return res; - + }; PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -1791,16 +1791,16 @@ inline void counter_pairwise_overall_change( return res; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication}), "Pairs of genes changing" + get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -1816,7 +1816,7 @@ inline void counter_pairwise_preserving( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -1843,7 +1843,7 @@ inline void counter_pairwise_preserving( { if (Array(k, j) == 1u) - return 0.0; + return 0.0; for (auto off = 0u; off < Array.ncol(); ++off) { @@ -1861,7 +1861,7 @@ inline void counter_pairwise_preserving( { if (Array(k, j) == 1u) - return 0.0; + return 0.0; for (auto off = 0u; off < Array.ncol(); ++off) { @@ -1879,7 +1879,7 @@ inline void counter_pairwise_preserving( { if (Array(k, j) == 0u) - return 0.0; + return 0.0; for (auto off = 0u; off < Array.ncol(); ++off) { @@ -1897,7 +1897,7 @@ inline void counter_pairwise_preserving( { if (Array(k, j) == 0u) - return 0.0; + return 0.0; for (auto off = 0u; off < Array.ncol(); ++off) { @@ -1912,7 +1912,7 @@ inline void counter_pairwise_preserving( } return res; - + }; PHYLO_COUNTER_LAMBDA(tmp_init) { @@ -1920,9 +1920,9 @@ inline void counter_pairwise_preserving( IF_NOTMATCHES() return 0.0; - + PHYLO_CHECK_MISSING(); - + double n = static_cast< double >(Array.ncol()); if (!Array.D_ptr()->states[data[1u]] && !Array.D_ptr()->states[data[2u]]) return n * (n - 1.0) / 2.0; @@ -1930,16 +1930,16 @@ inline void counter_pairwise_preserving( return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, nfunA, nfunB}), "Pariwise preserve (" + std::to_string(nfunA) + ", " + std::to_string(nfunB) + ")" +get_last_name(duplication) ); - + return; - + } // ----------------------------------------------------------------------------- @@ -1955,7 +1955,7 @@ inline void counter_pairwise_first_gain( size_t duplication = Geese::etype_default ) { - + PHYLO_COUNTER_LAMBDA(tmp_count) { @@ -2004,7 +2004,7 @@ inline void counter_pairwise_first_gain( else { - if (Array(k, off) == 1u) + if (Array(k, off) == 1u) // j: (0,0)\(0,1) -> (1,0)\(0,1), so less 1 res -= 1.0; else @@ -2016,29 +2016,29 @@ inline void counter_pairwise_first_gain( } } - + return res; - + }; PHYLO_COUNTER_LAMBDA(tmp_init) { PHYLO_CHECK_MISSING(); - + return 0.0; }; - + counters->add_counter( tmp_count, tmp_init, nullptr, PhyloCounterData({duplication, nfunA, nfunB}), "First gain (either " + std::to_string(nfunA) + " or " + std::to_string(nfunB) + ")" +get_last_name(duplication) ); - + return; - + } ///@} @@ -2081,7 +2081,7 @@ inline void rule_leafs( * @param pos Position of the focal statistic. * @param lb Lower bound * @param ub Upper bound - * @details + * @details * @return (void) adds a rule limiting the support of the model. */ inline void rule_dyn_limit_changes( @@ -2092,7 +2092,7 @@ inline void rule_dyn_limit_changes( size_t duplication = Geese::etype_default ) { - + PHYLO_RULE_DYN_LAMBDA(tmp_rule) { @@ -2104,7 +2104,7 @@ inline void rule_dyn_limit_changes( return true; else if (!Array.D_ptr()->duplication & (rule_type != Geese::etype_speciation)) return true; - + } if (data() < data.lb) @@ -2113,7 +2113,7 @@ inline void rule_dyn_limit_changes( return false; else return true; - + }; // Checking whether the rule makes sense (dupl) @@ -2128,27 +2128,27 @@ inline void rule_dyn_limit_changes( ) ); } - + support->get_rules_dyn()->add_rule( tmp_rule, PhyloRuleDynData( support->get_current_stats(), pos, lb, ub, duplication ), - std::string("Limiting changes in '") + + std::string("Limiting changes in '") + support->get_counters()->get_names()[pos] + "' to [" + std::to_string(lb) + ", " + std::to_string(ub) + std::string("]") + - get_last_name(duplication), - std::string("When the support is ennumerated, the number of changes in '") + + get_last_name(duplication), + std::string("When the support is ennumerated, the number of changes in '") + support->get_counters()->get_names()[pos] + std::to_string(pos) + "' is limited to [" + std::to_string(lb) + ", " + std::to_string(ub) + "]" + get_last_name(duplication) ); - + return; - + } ///@} @@ -2161,4 +2161,4 @@ inline void rule_dyn_limit_changes( #undef IF_NOTMATCHES -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/flock-bones.hpp b/include/barry/models/geese/flock-bones.hpp index fb06d92..fb242d8 100644 --- a/include/barry/models/geese/flock-bones.hpp +++ b/include/barry/models/geese/flock-bones.hpp @@ -9,7 +9,7 @@ class Geese; * @details This object buils a model with multiple trees (Geese objects), * with all of these using the same PhyloModel object. Available counters * (terms) can be found in \ref counter-phylo. - * + * */ class Flock { public: @@ -17,7 +17,7 @@ class Flock { std::vector< Geese > dat; size_t nfunctions = 0u; bool initialized = false; - + // Common components std::mt19937 rengine; PhyloModel model = PhyloModel(); @@ -27,7 +27,7 @@ class Flock { /** * @brief Add a tree to the flock - * + * * @param annotations see Geese::Geese. * @param geneid see Geese. * @param parent see Geese. @@ -43,13 +43,13 @@ class Flock { /** * @brief Set the seed of the model - * + * * @param s Passed to the `rengine.seed()` member object. */ void set_seed(const size_t & s); void init(size_t bar_width = BARRY_PROGRESS_BAR_WIDTH); - + // void add_geese(Geese x); PhyloCounters * get_counters(); PhyloSupport * get_support_fun(); @@ -59,12 +59,12 @@ class Flock { /** * @brief Returns the joint likelihood of the model - * + * * @param par Vector of model parameters. * @param as_log When `true` it will return the value as log. * @param use_reduced_sequence When `true` (default) will compute the * likelihood using the reduced sequence, which is faster. - * @return double + * @return double */ double likelihood_joint( const std::vector< double > & par, @@ -74,7 +74,7 @@ class Flock { ); /** - * @name Information about the model + * @name Information about the model */ ///@{ size_t nfuns() const noexcept; @@ -93,7 +93,7 @@ class Flock { /** * @brief Access the i-th geese element - * + * * @param i Element to access * @param check_bounds When true, it will check bounds. * @return Geese* @@ -102,4 +102,4 @@ class Flock { }; -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/flock-meat.hpp b/include/barry/models/geese/flock-meat.hpp index 670e6ab..08f4659 100644 --- a/include/barry/models/geese/flock-meat.hpp +++ b/include/barry/models/geese/flock-meat.hpp @@ -1,4 +1,4 @@ -#ifndef GEESE_FLOCK_MEET_HPP +#ifndef GEESE_FLOCK_MEET_HPP #define GEESE_FLOCK_MEET_HPP 1 // #include "flock-bones.hpp" @@ -17,7 +17,7 @@ inline size_t Flock::add_data( model.set_rengine(&this->rengine, false); model.add_hasher(keygen_full); - + model.store_psets(); } @@ -34,7 +34,7 @@ inline size_t Flock::add_data( if (dat.size() == 1u) this->nfunctions = dat[0].nfuns(); - + return dat.size() - 1u; } @@ -65,7 +65,7 @@ inline void Flock::init(size_t bar_width) a.rengine = &rengine; a.delete_rengine = false; - + } // Initializing the models. @@ -94,7 +94,7 @@ inline void Flock::init(size_t bar_width) } this->initialized = true; - + } inline PhyloCounters * Flock::get_counters() @@ -157,10 +157,10 @@ inline double Flock::likelihood_joint( #if defined(_OPENMP) || defined(__OPENMP) #pragma omp parallel for reduction(+:ans) num_threads(ncores) #endif - for (auto& d : this->dat) + for (auto& d : this->dat) ans += d.likelihood(par, as_log, use_reduced_sequence, 1u, true); } else { - for (auto& d : this->dat) + for (auto& d : this->dat) ans += d.likelihood(par, as_log, use_reduced_sequence, 1u, true); } @@ -168,20 +168,20 @@ inline double Flock::likelihood_joint( else { - if (ncores > 1u) + if (ncores > 1u) { #if defined(_OPENMP) || defined(__OPENMP) #pragma omp parallel for reduction(*:ans) num_threads(ncores) #endif - for (auto& d : this->dat) + for (auto& d : this->dat) ans *= d.likelihood(par, as_log, use_reduced_sequence, 1u, true); } else { - for (auto& d : this->dat) + for (auto& d : this->dat) ans *= d.likelihood(par, as_log, use_reduced_sequence, 1u, true); } - + } - + return ans; } @@ -277,7 +277,7 @@ inline size_t Flock::parse_polytomies( } -inline void Flock::print() const +inline void Flock::print() const { // Information relevant to print: @@ -300,23 +300,23 @@ inline void Flock::print() const nones += tree.n_ones; ndpl += tree.n_dupl_events; nspe += tree.n_spec_events; - + } printf_barry("FLOCK (GROUP OF GEESE)\nINFO ABOUT THE PHYLOGENIES\n"); - + printf_barry("# of phylogenies : %li\n", ntrees()); - + printf_barry("# of functions : %li\n", nfuns()); - + printf_barry("# of ann. [zeros; ones] : [%li; %li]\n", nzeros, nones); - + printf_barry("# of events [dupl; spec] : [%li; %li]\n", ndpl, nspe); - + printf_barry("Largest polytomy : %li\n", parse_polytomies(false)); - + printf_barry("\nINFO ABOUT THE SUPPORT\n"); - + return this->model.print(); } @@ -331,4 +331,4 @@ inline Geese* Flock::operator()(size_t i, bool check_bounds) } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-bones.hpp b/include/barry/models/geese/geese-bones.hpp index 2cee3ae..518bbb7 100644 --- a/include/barry/models/geese/geese-bones.hpp +++ b/include/barry/models/geese/geese-bones.hpp @@ -16,7 +16,7 @@ inline std::vector< Ta > vector_caster(const std::vector< Tb > & x) { ans.push_back(static_cast< Ta >(*i)); return ans; - + } #define INITIALIZED() if (!this->initialized) \ @@ -26,7 +26,7 @@ inline std::vector< Ta > vector_caster(const std::vector< Tb > & x) { RULE_FUNCTION(rule_empty_free) { return Array(i, j) == 9u; - + } @@ -43,11 +43,11 @@ inline std::vector< double > keygen_full( static_cast(array.nrow()) * 100000 + static_cast(array.ncol()), // state of the parent - 1000000.0, + 1000000.0, // type of the parent - array.D_ptr()->duplication ? 1.0 : 0.0, + array.D_ptr()->duplication ? 1.0 : 0.0, // Annotations with zeros - 0.0 + 0.0 }; // State of the parent @@ -65,7 +65,7 @@ inline std::vector< double > keygen_full( // } return dat; - + } inline bool vec_diff( @@ -91,24 +91,24 @@ class Flock; */ /** * @brief Class representing a phylogenetic tree model with annotations. - * + * * The `Geese` class represents a phylogenetic tree model with annotations. It * includes a total of `N + 1` nodes, the `+ 1` being the root node. The class * provides methods for initializing the model, calculating the likelihood, - * simulating trees, and making predictions. - * + * simulating trees, and making predictions. + * * The class includes shared objects within a `Geese` object, such as `rengine`, * `model`, `states`, `n_zeros`, `n_ones`, `n_dupl_events`, and `n_spec_events`. * It also includes information about the type of event, such as `etype_default`, * `etype_speciation`, `etype_duplication`, and `etype_either`. - * + * * The class provides constructors, a destructor, and methods for initializing * the model, inheriting support, calculating the sequence, calculating the * reduced sequence, calculating the likelihood, calculating the likelihood * exhaustively, getting probabilities, setting the seed, simulating trees, * parsing polytomies, getting observed counts, printing observed counts, * printing information about the GEESE, and making predictions. - * + * * @see Flock */ class Geese { @@ -120,10 +120,10 @@ class Geese { * @details * Since users may start adding counters before initializing the PhyloModel * object, the object `counter` is initialized first. - * + * * While the member `model` has an `rengine`, since `Geese` can sample trees, * we have the option to keep it separate. - * + * */ ///@{ std::mt19937 * rengine = nullptr; @@ -140,13 +140,13 @@ class Geese { // Data size_t nfunctions; std::map< size_t, Node > nodes; - + barry::MapVec_type< size_t > map_to_state_id; std::vector< std::vector< std::vector< size_t > > > pset_loc; ///< Locations of columns // Tree-traversal sequence std::vector< size_t > sequence; - std::vector< size_t > reduced_sequence; + std::vector< size_t > reduced_sequence; // Admin-related objects bool initialized = false; @@ -154,14 +154,14 @@ class Geese { bool delete_support = false; // Information about the type of event - + /*** * @name Information about the type of event * @details * The type of event is stored in the `etype` member. The possible values * are `etype_default`, `etype_speciation`, `etype_duplication`, and * `etype_either`. - * + * */ ///@{ static const size_t etype_default = 1ul; @@ -184,8 +184,8 @@ class Geese { * @param parent Id of the parent gene. Also of length `N` * @param duplication Logical scalar indicating the type of event (true: * duplication, false: speciation.) - * - * @details + * + * @details * The ordering of the entries does not matter. Passing the nodes in post * order or not makes no difference to the constructor. */ @@ -201,7 +201,7 @@ class Geese { // Copy constructor Geese(const Geese & model_, bool copy_data = true); - + // Constructor move Geese(Geese && x) noexcept; @@ -241,7 +241,7 @@ class Geese { ); /** - * @name Information about the model + * @name Information about the model * @param verb When `true` it will print out information about the encountered * polytomies. */ @@ -273,7 +273,7 @@ class Geese { /** * @name Geese prediction * @brief Calculate the conditional probability - * + * * @param par Vector of parameters (terms + root). * @param res_prob Vector indicating each nodes' state probability. * @param leave_one_out When `true`, it will compute the predictions using @@ -282,12 +282,12 @@ class Geese { * on the induced sub-tree with annotated leafs. * @param use_reduced_sequence Passed to the `likelihood` method. * @param preorder For the tree traversal. - * + * * @details When `res_prob` is specified, the function will attach * the member vector `probabilities` from the `Node`s objects. This * contains the probability that the ith node has either of the * possible states. - * + * * @return std::vector< double > Returns the posterior probability */ ///@{ @@ -298,7 +298,7 @@ class Geese { bool only_annotated = false, bool use_reduced_sequence = true ); - + std::vector< std::vector > predict_backend( const std::vector< double > & par, bool use_reduced_sequence, @@ -329,10 +329,10 @@ class Geese { /** * @name Non-const pointers to shared objects in `Geese` - * + * * @details These functions provide direct access to some member * objects that are shared by the nodes within `Geese`. - * + * * @return `get_rengine()` returns the Pseudo-RNG engine used. * @return `get_counters()` returns the vector of counters used. * @return `get_model()` returns the `Model` object used. @@ -344,16 +344,16 @@ class Geese { PhyloModel * get_model(); PhyloSupport * get_support_fun(); ///@} - + /** * @brief Powerset of a gene's possible states * @details This list of vectors is used throughout `Geese`. It lists * all possible combinations of functional states for any gene. Thus, * for `P` functions, there will be `2^P` possible combinations. - * + * * @return std::vector< std::vector< bool > > of length `2^P`. */ - std::vector< std::vector< bool > > get_states() const; + std::vector< std::vector< bool > > get_states() const; std::vector< size_t > get_annotated_nodes() const; ///< Returns the ids of the nodes with at least one annotation std::vector< size_t > get_annotations() const; ///< Returns the annotations of the nodes with at least one annotation diff --git a/include/barry/models/geese/geese-meat-constructors.hpp b/include/barry/models/geese/geese-meat-constructors.hpp index 5e75a2e..f9f168e 100644 --- a/include/barry/models/geese/geese-meat-constructors.hpp +++ b/include/barry/models/geese/geese-meat-constructors.hpp @@ -227,7 +227,7 @@ inline Geese::Geese( } -inline Geese::Geese(const Geese & model_, bool copy_data) : +inline Geese::Geese(const Geese & model_, bool copy_data) : states(model_.states), n_zeros(model_.n_zeros), n_ones(model_.n_ones), @@ -241,7 +241,7 @@ inline Geese::Geese(const Geese & model_, bool copy_data) : reduced_sequence(model_.reduced_sequence), initialized(model_.initialized) { - + // Replicating ------------------------------------------------------------- if (copy_data) { @@ -259,7 +259,7 @@ inline Geese::Geese(const Geese & model_, bool copy_data) : } } else { - + if (model_.rengine != nullptr) { rengine = model_.rengine; @@ -292,7 +292,7 @@ inline Geese::Geese(const Geese & model_, bool copy_data) : // Clearing offspring this->nodes[i].offspring.clear(); - // I cannot directly access the node since, if non existent, it will + // I cannot directly access the node since, if non existent, it will // create an entry with it (alegedly). auto n = model_.nodes.find(i); @@ -302,7 +302,7 @@ inline Geese::Geese(const Geese & model_, bool copy_data) : } return; - + } // Constructor move @@ -346,7 +346,7 @@ inline Geese::Geese(Geese && x) noexcept : model = x.model; delete_support = false; - + } // Figuring out if model needs to be updated @@ -359,4 +359,4 @@ inline Geese::Geese(Geese && x) noexcept : -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat-likelihood.hpp b/include/barry/models/geese/geese-meat-likelihood.hpp index af0f25a..b4d8728 100644 --- a/include/barry/models/geese/geese-meat-likelihood.hpp +++ b/include/barry/models/geese/geese-meat-likelihood.hpp @@ -16,7 +16,7 @@ inline void pset_loop( const std::vector< std::vector< size_t > > & locations, const std::vector & node_offspring, const double * psetprobs -) +) { // Retrieving the pset const auto & x = psets[n]; @@ -34,7 +34,7 @@ inline void pset_loop( // Setting the node const Node * n_off = node_offspring[o]; - + // In the case that the offspring is a leaf, then we need to // check whether the state makes sense. if (n_off->is_leaf()) @@ -51,7 +51,7 @@ inline void pset_loop( break; } - + } } @@ -59,7 +59,7 @@ inline void pset_loop( // Going out if (off_mult < 0) break; - + continue; } @@ -92,7 +92,7 @@ inline void pset_loop( err; throw std::runtime_error(err); - + } // Adding to the total probabilities @@ -165,18 +165,18 @@ inline double Geese::likelihood( // Making sure parallelization makes sense if (psets.size() < 128) ncores = 1u; - + // Summation over all possible values of X const auto & node_offspring = node.offspring; std::vector< double > totprob_n(psets.size(), 0.0); - for (size_t n = 0u; n < psets.size(); ++n) + for (size_t n = 0u; n < psets.size(); ++n) { pset_loop( n, s, nfunctions, node_id, array_id, totprob_n, - par0, states, psets, locations, + par0, states, psets, locations, node_offspring, psetsprobs_s ); - } + } // Setting the probability at the node node.subtree_prob[s] = 0.0; @@ -223,4 +223,4 @@ inline double Geese::likelihood( return as_log ? std::log(ll) : ll; } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat-likelihood_exhaust.hpp b/include/barry/models/geese/geese-meat-likelihood_exhaust.hpp index a5be5d1..b933975 100644 --- a/include/barry/models/geese/geese-meat-likelihood_exhaust.hpp +++ b/include/barry/models/geese/geese-meat-likelihood_exhaust.hpp @@ -2,7 +2,7 @@ #ifndef GEESE_MEAT_LIKELIHOOD_EXHAUST_HPP #define GEESE_MEAT_LIKELIHOOD_EXHAUST_HPP 1 // #include "../../barry.hpp" -// #include "geese-bones.hpp" +// #include "geese-bones.hpp" inline double Geese::likelihood_exhaust( const std::vector< double > & par @@ -33,7 +33,7 @@ inline double Geese::likelihood_exhaust( for (size_t i = 0u; i < nfuns(); ++i) base(i, n.second.ord) = n.second.annotations[i]; - + } PhyloPowerSet pset(base);//this->nfuns(), this->nnodes()); @@ -48,15 +48,15 @@ inline double Geese::likelihood_exhaust( std::reverse(preorder.begin(), preorder.end()); double totprob = 0.0; - - // This vector says whether the probability has to be included in + + // This vector says whether the probability has to be included in // the final likelihood or not. for (size_t p = 0u; p < pset.size(); ++p) { - + // ith state const PhyloArray * s = &pset[p]; - + // Following the sequence double prob = 1.0; std::vector< size_t > tmpstates(this->nfuns()); @@ -71,7 +71,7 @@ inline double Geese::likelihood_exhaust( // Root node first if (node->parent == nullptr) - { + { // Since it is the root, the first probability is computed using // the root only for (auto k = 0u; k < this->nfuns(); ++k) @@ -122,4 +122,4 @@ inline double Geese::likelihood_exhaust( return totprob; } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat-predict.hpp b/include/barry/models/geese/geese-meat-predict.hpp index ed768b9..6191cbc 100644 --- a/include/barry/models/geese/geese-meat-predict.hpp +++ b/include/barry/models/geese/geese-meat-predict.hpp @@ -26,16 +26,16 @@ inline std::vector< std::vector > Geese::predict_backend( for (size_t f = 0u; f < nfuns(); ++f) rootp[s] *= states[s][f] ? par_root[f] : (1.0 - par_root[f]); - + } - // Making room + // Making room std::vector< std::vector > res( nnodes(), std::vector(nfuns()) ); // Step 1: Computing the probability at the root node - std::vector< double > tmp_prob(nfuns(), 0.0); + std::vector< double > tmp_prob(nfuns(), 0.0); size_t root_id = preorder[0u]; Node * tmp_node = &nodes[root_id]; tmp_node->probability.resize(states.size(), 0.0); @@ -57,7 +57,7 @@ inline std::vector< std::vector > Geese::predict_backend( { throw std::runtime_error("Probability is not finite"); } - + // Marginalizing the probabilities P(x_sf | D) for (size_t f = 0u; f < nfuns(); ++f) { @@ -74,7 +74,7 @@ inline std::vector< std::vector > Geese::predict_backend( // Storing the final prob res[nodes[root_id].ord] = tmp_prob; - + // Retrieving the powersets probabilities const auto & pset_probs = *(model->get_pset_probs()); const auto & arrays2support = *(model->get_arrays2support()); @@ -95,7 +95,7 @@ inline std::vector< std::vector > Geese::predict_backend( everything_below.reserve(states.size()); everything_above.reserve(states.size()); - + // All combinations of the the parent states // So psets[s] = combinations of offspring given state s. // psets[s][i] = The ith combination of offspring given state s. @@ -133,7 +133,7 @@ inline std::vector< std::vector > Geese::predict_backend( // Adding to the map, we only do this during the first run, // afterwards, we need to actually look for the array. bool in_the_set = true; /// < True if the array belongs to the set - + // Everything below just need to be computed only once // and thus, if already added, no need to go through all of this! double everything_below_p = 1.0; @@ -155,7 +155,7 @@ inline std::vector< std::vector > Geese::predict_backend( in_the_set = false; break; } - + } if (!in_the_set) @@ -181,7 +181,7 @@ inline std::vector< std::vector > Geese::predict_backend( continue; pset.push_back(array_p); // Generating a copy - + // - With focal node, conditioning on it beening status s. // - But the offspring probabilities are the central ones here. // - So the saved values are for computing P(x_offspring | Data) @@ -204,7 +204,7 @@ inline std::vector< std::vector > Geese::predict_backend( everything_below.push_back(std::move(below)); everything_above.push_back(std::move(above)); - + } // end for states // Marginalizing at the state level for each offspring @@ -234,7 +234,7 @@ inline std::vector< std::vector > Geese::predict_backend( if (cvec[f] == 1u) res[parent.offspring[off]->ord][f] += everything_above[s][p]; } - + } @@ -242,7 +242,7 @@ inline std::vector< std::vector > Geese::predict_backend( } } - // Finally, we can marginalize the values at the + // Finally, we can marginalize the values at the // gene function level. for (const auto & off : parent.offspring) { @@ -258,8 +258,8 @@ inline std::vector< std::vector > Geese::predict_backend( std::to_string(res[off->ord][f]); throw std::logic_error(msg); - - } + + } if (res[off->ord][f] > 1.0) res[off->ord][f] = 1.0; @@ -271,7 +271,7 @@ inline std::vector< std::vector > Geese::predict_backend( } } // end for over preorder - + return res; } @@ -316,7 +316,7 @@ inline std::vector< std::vector > Geese::predict( } // In this case, we need to update the predictions, mostly of the annotated - // leaf nodes. Because of + // leaf nodes. Because of if (leave_one_out) { @@ -339,7 +339,7 @@ inline std::vector< std::vector > Geese::predict( break; } - + if (!use_it) continue; @@ -366,9 +366,9 @@ inline std::vector< std::vector > Geese::predict( } } - + return res; } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat-predict_exhaust.hpp b/include/barry/models/geese/geese-meat-predict_exhaust.hpp index db86387..4c89a8a 100644 --- a/include/barry/models/geese/geese-meat-predict_exhaust.hpp +++ b/include/barry/models/geese/geese-meat-predict_exhaust.hpp @@ -77,19 +77,19 @@ inline std::vector< std::vector > Geese::predict_exhaust_backend( PhyloRuleData() ); pset.calc(); - + // Making space for the expected values std::vector< double > expected(nnodes() * nfuns(), 0.0); - - // This vector says whether the probability has to be included in + + // This vector says whether the probability has to be included in // the final likelihood or not. for (size_t p = 0u; p < pset.size(); ++p) { - + // ith state const PhyloArray * s = &pset[p]; - - // Computing the likelihood of the state s + + // Computing the likelihood of the state s double current_prob = 1.0; for (auto & o: preorder) { @@ -111,7 +111,7 @@ inline std::vector< std::vector > Geese::predict_exhaust_backend( current_prob *= par_state[f] ? par_root[f] : (1.0 - par_root[f]); } - + // Generating a copy of the observed array // (data is copied so that we can chage the state of the parent) PhyloArray tmparray(n.array, true); @@ -123,7 +123,7 @@ inline std::vector< std::vector > Geese::predict_exhaust_backend( // Updating offspring annotations int loc = 0; for (auto & off : n.offspring) { - + for (size_t f = 0u; f < nfuns(); ++f) { @@ -138,19 +138,19 @@ inline std::vector< std::vector > Geese::predict_exhaust_backend( ++loc; } - + // Computing the likelihood current_prob *= model->likelihood(par_terms, tmparray, -1, false); } // this->update_annotations(n.second.id, s->get_col_vec(n.second.ord)); - + // Adding to the overall probability for (auto & n: nodes) for (size_t j = 0u; j < nfuns(); ++j) expected[n.second.ord + j * nnodes()] += s->operator()(j, n.second.ord) * current_prob/ baseline_likelihood; - + } // Coercing expected to a list vector @@ -166,4 +166,4 @@ inline std::vector< std::vector > Geese::predict_exhaust_backend( return res; } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat-predict_sim.hpp b/include/barry/models/geese/geese-meat-predict_sim.hpp index b2d44c7..4e98368 100644 --- a/include/barry/models/geese/geese-meat-predict_sim.hpp +++ b/include/barry/models/geese/geese-meat-predict_sim.hpp @@ -44,7 +44,7 @@ inline std::vector< std::vector > Geese::predict_sim( if (n.id == id) continue; - const auto & ord = nodes[id].ord; + const auto & ord = nodes[id].ord; const auto & n_w_ann = nodes[id].annotations; for (size_t f = 0u; f < nfuns(); ++f) { @@ -94,10 +94,10 @@ inline std::vector< std::vector > Geese::predict_sim( for (size_t f = 0u; f < nfuns(); ++f) res_vec[i][f] /= (static_cast< double >(counts[i]) + 1e-10); } - + return res_vec; } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat-simulate.hpp b/include/barry/models/geese/geese-meat-simulate.hpp index b7e6f9c..9698321 100644 --- a/include/barry/models/geese/geese-meat-simulate.hpp +++ b/include/barry/models/geese/geese-meat-simulate.hpp @@ -20,7 +20,7 @@ inline std::vector< std::vector< size_t > > Geese::simulate( p = std::exp(p)/(std::exp(p) + 1); } - // Making room + // Making room std::vector< std::vector< size_t > > res(nodes.size()); // Inverse sequence @@ -48,14 +48,14 @@ inline std::vector< std::vector< size_t > > Geese::simulate( } #ifdef BARRY_DEBUG - + // auto totprob = std::accumulate(rootp.begin(), rootp.end(), 0.0); // if (totprob < 0.9999999999999999 || totprob > 1.0000000000000001) // throw std::runtime_error("Root probabilities do not sum to 1!" // " (totprob = " + std::to_string(totprob) + ")"); - + #endif - + // We now know the state of the root res[nodes[preorder[0u]].ord] = vector_caster< size_t, bool>(states[idx]); @@ -87,4 +87,4 @@ inline std::vector< std::vector< size_t > > Geese::simulate( } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-meat.hpp b/include/barry/models/geese/geese-meat.hpp index 67e9c97..bcda0f5 100644 --- a/include/barry/models/geese/geese-meat.hpp +++ b/include/barry/models/geese/geese-meat.hpp @@ -65,7 +65,7 @@ inline void Geese::init_node(Node & n) n.narray.resize(states.size()); } - + // Here we have an issue: Some transitions may not be right // under the dynamic rules. So not all states can be valid. // The arrays and narrays need to be updated once the model @@ -99,7 +99,7 @@ inline void Geese::init_node(Node & n) " cannot be added to the model with error:\n" + err + "\n. This is likely due to a dynamic rule. " + "The array to be added was in the following state:"; - + std::string state_str = ""; for (const auto & ss : states[s]) state_str += std::to_string(ss) + " "; @@ -107,7 +107,7 @@ inline void Geese::init_node(Node & n) err += state_str + "\n"; throw std::runtime_error(err); - + } } @@ -148,7 +148,7 @@ inline void Geese::init(size_t bar_width) { // Checking rseed, this is relevant when dealing with a flock. In the case of // flock, both model and rengine are shared. - if (this->model->get_rengine() == nullptr) + if (this->model->get_rengine() == nullptr) this->model->set_rengine(this->rengine, false); // All combinations of the function @@ -164,7 +164,7 @@ inline void Geese::init(size_t bar_width) { { states.emplace_back(std::vector< bool >(nfunctions, false)); - + for (auto j = 0u; j < nfunctions; ++j) { @@ -192,10 +192,10 @@ inline void Geese::init(size_t bar_width) { // Only parents get a node if (!iter.second.is_leaf()) - this->init_node(iter.second); - + this->init_node(iter.second); + prog_bar.next(); - + } prog_bar.end(); @@ -211,8 +211,8 @@ inline void Geese::init(size_t bar_width) { // Only parents get a node if (!iter.second.is_leaf()) - this->init_node(iter.second); - + this->init_node(iter.second); + } } @@ -243,13 +243,13 @@ inline void Geese::init(size_t bar_width) { sup_array[a].get_col_vec(&tmpstate, o, false); pset_loc[s][a].push_back(map_to_state_id[tmpstate]); - - } + + } } } - + // So that others now know it was initialized initialized = true; @@ -259,7 +259,7 @@ inline void Geese::init(size_t bar_width) { inline void Geese::inherit_support(const Geese & model_, bool delete_support_) { - + if (this->model != nullptr) throw std::logic_error( "There is already a -model- in this Geese. Cannot set a -model- after one is present." @@ -278,9 +278,9 @@ inline void Geese::inherit_support(const Geese & model_, bool delete_support_) this->delete_rengine = false; } - + this->rengine = model_.rengine; - + return; } @@ -304,7 +304,7 @@ inline void Geese::update_annotations( // parent node nodes[nodeid].annotations = newann; - // This only makes sense (for now) if it is a tip + // This only makes sense (for now) if it is a tip if (!nodes[nodeid].is_leaf()) return; @@ -389,11 +389,11 @@ inline void Geese::calc_reduced_sequence() { // Checking, am I including any of my offspring? - for (auto& o : n.offspring) + for (auto& o : n.offspring) if (includeit[o->ord]) { - + includeit[n.ord] = true; reduced_sequence.push_back(i); break; @@ -414,7 +414,7 @@ inline std::vector< double > Geese::get_probabilities() const res.reserve( this->states.size() * nodes.size() ); - + for (auto& i : sequence) { @@ -424,7 +424,7 @@ inline std::vector< double > Geese::get_probabilities() const } return res; - + } inline size_t Geese::nfuns() const noexcept @@ -467,7 +467,7 @@ inline size_t Geese::support_size() const noexcept return 0u; return model->support_size(); - + } inline std::vector< size_t > Geese::nannotations() const noexcept @@ -509,7 +509,7 @@ inline size_t Geese::parse_polytomies( if (verb) printf_barry("Node id: %li has polytomy size %li\n", n.second.id, noff); - + } if (noff > largest) @@ -679,17 +679,17 @@ inline void Geese::print_nodes() const printf_barry("GEESE\nINFO ABOUT NODES\n"); for (const auto & n: nodes) - { + { printf_barry("% 4li - Id: %li -- ", n.second.ord, n.second.id); // Node type printf_barry( "node type: %s -- ", - n.second.is_leaf() ? + n.second.is_leaf() ? std::string("leaf").c_str() : std::string("internal").c_str() ); - + // Event type printf_barry( "event type: %s -- ", @@ -792,13 +792,13 @@ inline std::vector< size_t > Geese::get_annotated_nodes() const { } inline std::vector< size_t > Geese::get_annotations() const { - + // Makeing space for the annotations std::vector< size_t > ann(this->nfuns() * this->nnodes(), 9u); size_t nrows = this->nnodes(); for (auto & n : nodes) { - + // Getting the location size_t row = n.second.ord; @@ -812,11 +812,11 @@ inline std::vector< size_t > Geese::get_annotations() const { } } - + } - + return ann; - + } diff --git a/include/barry/models/geese/geese-node-bones.hpp b/include/barry/models/geese/geese-node-bones.hpp index 7ecf3e7..c605915 100644 --- a/include/barry/models/geese/geese-node-bones.hpp +++ b/include/barry/models/geese/geese-node-bones.hpp @@ -27,17 +27,17 @@ class Node { std::vector< double > subtree_prob; ///< Induced subtree probabilities std::vector< double > probability; ///< The probability of observing each state - + /** * @name Construct a new Node object - * + * */ ///@{ - + Node() : ord(std::numeric_limits< size_t >::max()) {}; Node(size_t id_, size_t ord_, bool duplication_); Node(size_t id_, size_t ord_, std::vector< size_t > annotations_, bool duplication_); - + // Move constructor Node(Node && x) noexcept; @@ -79,11 +79,11 @@ inline Node::Node(Node && x) noexcept : probability(std::move(x.probability)) { return; - + } inline Node::Node(const Node & x) : - id(x.id), ord(x.ord), array(x.array), + id(x.id), ord(x.ord), array(x.array), annotations(x.annotations), duplication(x.duplication), arrays(x.arrays), parent(x.parent), @@ -94,7 +94,7 @@ inline Node::Node(const Node & x) : probability(x.probability) { return; - + } inline int Node::get_parent() const { @@ -113,4 +113,4 @@ inline bool Node::is_leaf() const noexcept { return offspring.size() == 0u; } -#endif \ No newline at end of file +#endif diff --git a/include/barry/models/geese/geese-types.hpp b/include/barry/models/geese/geese-types.hpp index 445b83f..c946cf0 100644 --- a/include/barry/models/geese/geese-types.hpp +++ b/include/barry/models/geese/geese-types.hpp @@ -8,27 +8,27 @@ * */ /** * @brief Data definition for the `PhyloArray` class. - * + * * This holds basic information about a given node. - * + * */ class NodeData { public: - + /** * Branch length. */ std::vector< double > blengths = {}; - + /** * State of the parent node. */ std::vector< bool > states = {}; - + bool duplication = true; ///< Whether the node is a duplication. bool has_leaf = false; ///< Whether the node has a leaf as offspring. - - + + NodeData( const std::vector< double > & blengths_, const std::vector< bool > & states_, @@ -36,9 +36,9 @@ class NodeData { bool has_leaf_ = false ) : blengths(blengths_), states(states_), duplication(duplication_), has_leaf(has_leaf_) {}; - + // ~NodeData() {}; - + }; class PhyloCounterData { @@ -91,9 +91,9 @@ class PhyloRuleDynData { { return (*counts)[pos]; } - + ~PhyloRuleDynData() {}; - + }; @@ -116,7 +116,7 @@ typedef barry::Model PhyloPowerSet; ///@} -// template<> +// template<> // inline void PhyloArray::insert_cell( // size_t i, // size_t j, @@ -126,7 +126,7 @@ typedef barry::PowerSet PhyloPowerSet; // ) { // if (check_bounds) -// out_of_range(i,j); +// out_of_range(i,j); // auto & elptr = el[POS(i,j)]; @@ -135,8 +135,8 @@ typedef barry::PowerSet PhyloPowerSet; // el_rowsums[i] += v.value; // el_colsums[j] += v.value; - -// } + +// } // else // { @@ -153,4 +153,4 @@ typedef barry::PowerSet PhyloPowerSet; #undef POS -#endif \ No newline at end of file +#endif diff --git a/include/barry/powerset-bones.hpp b/include/barry/powerset-bones.hpp index 50ed1f9..7ba11d7 100644 --- a/include/barry/powerset-bones.hpp +++ b/include/barry/powerset-bones.hpp @@ -1,18 +1,18 @@ -#ifndef BARRY_POWERSET_BONES_HPP +#ifndef BARRY_POWERSET_BONES_HPP #define BARRY_POWERSET_BONES_HPP 1 /** * @brief Powerset of a binary array - * - * @tparam Array_Type - * @tparam Data_Rule_Type + * + * @tparam Array_Type + * @tparam Data_Rule_Type */ -template , typename Data_Rule_Type = bool> +template , typename Data_Rule_Type = bool> class PowerSet { - + private: - void calc_backend_sparse(size_t pos = 0u); - void calc_backend_dense(size_t pos = 0u); + void calc_backend_sparse(size_t pos = 0u); + void calc_backend_dense(size_t pos = 0u); public: Array_Type EmptyArray; @@ -27,13 +27,13 @@ class PowerSet { std::vector< size_t > coordinates_locked; size_t n_free; size_t n_locked; - + /** * @name Construct and destroy a PowerSet object - * + * */ ///@{ - PowerSet() : + PowerSet() : EmptyArray(), data(0u), rules(new Rules()), N(0u), M(0u) {}; PowerSet(size_t N_, size_t M_) : EmptyArray(N_, M_), data(0u), @@ -42,13 +42,13 @@ class PowerSet { ~PowerSet(); ///@} - + void init_support(); void calc(); void reset(size_t N_, size_t M_); - + /** - * @name Wrappers for the `Rules` member. + * @name Wrappers for the `Rules` member. * @details These will add rules to the model, which are shared by the * support and the actual counter function. */ @@ -59,7 +59,7 @@ class PowerSet { Data_Rule_Type data_ ); ///@} - + /** @name Getter functions */ ///@{ @@ -70,7 +70,7 @@ class PowerSet { std::size_t size() const noexcept {return data.size();}; const Array_Type& operator[](const size_t & i) const {return data.at(i);}; ///@} - + }; #endif diff --git a/include/barry/powerset-meat.hpp b/include/barry/powerset-meat.hpp index bf03d9b..551d000 100644 --- a/include/barry/powerset-meat.hpp +++ b/include/barry/powerset-meat.hpp @@ -18,7 +18,7 @@ inline PowerSet::~PowerSet() { template inline void PowerSet::init_support() { - + // Computing the locations coordinates_free.clear(); coordinates_locked.clear(); @@ -26,7 +26,7 @@ inline void PowerSet::init_support() n_free = coordinates_free.size() / 2u; n_locked = coordinates_locked.size() / 2u; - + // Computing initial statistics if (EmptyArray.nnozero() > 0u) { @@ -34,7 +34,7 @@ inline void PowerSet::init_support() if (EmptyArray.is_dense()) { - for (size_t i = 0u; i < n_free; ++i) + for (size_t i = 0u; i < n_free; ++i) EmptyArray( coordinates_free[i * 2u], coordinates_free[i * 2u + 1u] @@ -44,7 +44,7 @@ inline void PowerSet::init_support() else { - for (size_t i = 0u; i < n_free; ++i) + for (size_t i = 0u; i < n_free; ++i) EmptyArray.rm_cell( coordinates_free[i * 2u], coordinates_free[i * 2u + 1u], @@ -54,18 +54,18 @@ inline void PowerSet::init_support() } - + } // EmptyArray.clear(true); // EmptyArray.reserve(); - + // Resizing support - data.reserve(pow(2.0, n_free)); + data.reserve(pow(2.0, n_free)); // Adding the empty array to the set data.push_back(EmptyArray); - + return; } @@ -74,14 +74,14 @@ inline void PowerSet::calc_backend_sparse( size_t pos ) { - + // Did we reached the end?? if (pos >= n_free) return; - + // We will pass it to the next step, if the iteration makes sense. calc_backend_sparse(pos + 1u); - + // Toggle the cell (we will toggle it back after calling the counter) EmptyArray.insert_cell( coordinates_free[pos * 2u], @@ -98,20 +98,20 @@ inline void PowerSet::calc_backend_sparse( BARRY_USER_INTERRUPT } #endif - + // Again, we only pass it to the next level iff the next level is not // passed the last step. calc_backend_sparse(pos + 1u); - + // We need to restore the state of the cell EmptyArray.rm_cell( coordinates_free[pos * 2u], coordinates_free[pos * 2u + 1u], false, false - ); - + ); + return; - + } template @@ -119,33 +119,33 @@ inline void PowerSet::calc_backend_dense( size_t pos ) { - + // Did we reached the end?? if (pos >= n_free) return; - + // We will pass it to the next step, if the iteration makes sense. calc_backend_dense(pos + 1u); - + // Toggle the cell (we will toggle it back after calling the counter) EmptyArray(coordinates_free[pos * 2u], coordinates_free[pos * 2u + 1u]) = 1; data.push_back(EmptyArray); - + // Again, we only pass it to the next level iff the next level is not // passed the last step. calc_backend_dense(pos + 1u); - + // We need to restore the state of the cell EmptyArray(coordinates_free[pos * 2u], coordinates_free[pos * 2u + 1u]) = 0; - + return; - + } /*** - * Function to generate the powerset of the + * Function to generate the powerset of the */ template inline void PowerSet::calc() { @@ -160,7 +160,7 @@ inline void PowerSet::calc() { calc_backend_sparse(0u); return; - + } template @@ -168,10 +168,10 @@ inline void PowerSet::reset( size_t N_, size_t M_ ) { - + data.empty(); N = N_, M = M_; - + return; } @@ -180,7 +180,7 @@ template inline void PowerSet::add_rule( Rule rule ) { - + rules->add_rule(rule); return; } @@ -190,14 +190,14 @@ inline void PowerSet::add_rule( Rule_fun_type rule_fun_, Data_Rule_Type data_ ) { - + rules->add_rule( rule_fun_, data_ ); - + return; - + } -#endif \ No newline at end of file +#endif diff --git a/include/barry/progress.hpp b/include/barry/progress.hpp index d4daf17..6ec4933 100644 --- a/include/barry/progress.hpp +++ b/include/barry/progress.hpp @@ -16,7 +16,7 @@ class Progress { int last_loc; ///< Last location of the bar int cur_loc; ///< Last location of the bar int i; ///< Current iteration step - + public: Progress(int n_, int width_); @@ -55,4 +55,4 @@ inline void Progress::end() { } -#endif \ No newline at end of file +#endif diff --git a/include/barry/rules-bones.hpp b/include/barry/rules-bones.hpp index fecead4..dd11a7c 100644 --- a/include/barry/rules-bones.hpp +++ b/include/barry/rules-bones.hpp @@ -18,20 +18,20 @@ bool rule_fun_default(const Array_Type * array, size_t i, size_t j, Data_Type * */ template, typename Data_Type = bool> class Rule { - + private: Rule_fun_type fun; Data_Type dat; - + std::string name = ""; std::string desc = ""; - + public: /** * @name Construct a new Rule object * @brief Construct a new Rule object - * + * * @param fun_ A function of type `Rule_fun_type`. * @param dat_ Data pointer to be passed to `fun_` * @param delete_dat_ When `true`, the `Rule` destructor will delete the @@ -42,7 +42,7 @@ class Rule { Rule( Rule_fun_type fun_, Data_Type dat_, - std::string name_ = "", + std::string name_ = "", std::string desc_ = "" ) : fun(fun_), dat(dat_), name(name_), desc(desc_) {}; ///@} @@ -50,7 +50,7 @@ class Rule { ~Rule() {}; Data_Type & D(); ///< Read/Write access to the data. - + bool operator()(const Array_Type & a, size_t i, size_t j); std::string & get_name(); @@ -58,12 +58,12 @@ class Rule { std::string get_name() const; std::string get_description() const; - + }; /** * @brief Vector of objects of class Rule - * + * * @tparam Array_Type An object of class `BArray` * @tparam Data_Type Any type. */ @@ -72,7 +72,7 @@ class Rules { private: std::vector< Rule > data; - + public: Rules() {}; @@ -84,11 +84,11 @@ class Rules { size_t size() const noexcept { return data.size(); }; - + /** * @name Rule adding - * - * @param rule + * + * @param rule */ ///@{ void add_rule(Rule rule); @@ -102,7 +102,7 @@ class Rules { /** * @brief Check whether a given cell is free or locked - * + * * @param a A `BArray` object * @param i row position * @param j col position @@ -111,10 +111,10 @@ class Rules { */ bool operator()(const Array_Type & a, size_t i, size_t j); - + /** * @brief Computes the sequence of free and locked cells in an BArray - * + * * @param a An object of class `BArray`. * @param free Pointer to a vector of pairs (i, j) listing the free cells. * @param locked (optional) Pointer to a vector of pairs (i, j) listing the @@ -137,7 +137,7 @@ class Rules { typename std::vector< Rule >::iterator end() { return data.end(); }; - + }; diff --git a/include/barry/rules-meat.hpp b/include/barry/rules-meat.hpp index 1e707ab..1c54917 100644 --- a/include/barry/rules-meat.hpp +++ b/include/barry/rules-meat.hpp @@ -6,7 +6,7 @@ inline Rules::Rules( const Rules & rules_ ) { - // Copy all rules, if a rule is tagged as + // Copy all rules, if a rule is tagged as // to be deleted, then copy the value for (auto i = 0u; i != rules_.size(); ++i) this->add_rule(rules_.data[i]); @@ -22,7 +22,7 @@ Rules Rules::operator=( if (this != &rules_) { - // Copy all rules, if a rule is tagged as + // Copy all rules, if a rule is tagged as // to be deleted, then copy the value for (auto i = 0u; i != rules_.size(); ++i) this->add_rule(rules_.data[i]); @@ -72,9 +72,9 @@ template inline void Rules::add_rule( Rule rule ) { - + data.push_back(rule); - + return; } @@ -85,32 +85,32 @@ inline void Rules::add_rule( std::string name_, std::string description_ ) { - + data.push_back(Rule( rule_, data_, name_, description_ )); - + return; - + } template inline bool Rules::operator()( const Array_Type & a, size_t i, size_t j ) { - + if (data.size()==0u) return true; - + for (auto & f: data) if (!f.operator()(a, i, j)) return false; - + return true; - + } template @@ -120,14 +120,14 @@ inline void Rules::get_seq( std::vector< size_t > * locked ) { - + size_t N = a.nrow(); size_t K = a.ncol(); - + // Reserving some space (void) free->empty(); (void) free->reserve(2u * N * K); - + for (size_t i = 0u; i < N; ++i) { @@ -152,11 +152,11 @@ inline void Rules::get_seq( free->push_back(i); free->push_back(j); - + } } - + free->shrink_to_fit(); return; @@ -179,7 +179,7 @@ inline std::vector Rules::get_names() const template inline std::vector Rules::get_descriptions() const { - + std::vector< std::string > out; out.reserve(this->size()); for (size_t i = 0u; i < out.size(); ++i) diff --git a/include/barry/statscounter-bones.hpp b/include/barry/statscounter-bones.hpp index af0092e..1c162a3 100644 --- a/include/barry/statscounter-bones.hpp +++ b/include/barry/statscounter-bones.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_STATSCOUNTER_BONES_HPP +#ifndef BARRY_STATSCOUNTER_BONES_HPP #define BARRY_STATSCOUNTER_BONES_HPP 1 class NetworkDense; @@ -6,10 +6,10 @@ class NetCounterData; /** * @brief Count stats for a single Array. - * + * * Users can a list of functions that can be used with this. The baseline set of * arguments is a pointer to a binary array and a dataset to add the counts to. - */ + */ template class StatsCounter { @@ -19,7 +19,7 @@ class StatsCounter { const Array_Type * Array; Array_Type EmptyArray; std::vector< double > current_stats; - + // We will save the data here Counters * counters; bool counter_deleted = false; @@ -28,48 +28,48 @@ class StatsCounter { std::vector< double > count_all_sparse(); public: - + /** * @brief Creator of a `StatsCounter` - * + * * @param Array_ A const pointer to a `BArray`. */ StatsCounter(const Array_Type * Array_) : Array(Array_), EmptyArray(*Array_), counters(new Counters()) { - + // We are removing the entries without freeing the memory. This should // make the insertion faster. EmptyArray.clear(false); - + return; } /** * @brief Copy constructor - * - * @param counter + * + * @param counter */ StatsCounter(const StatsCounter & counter); - + /** * @brief Can be created without setting the array. - * + * */ StatsCounter() : Array(nullptr), EmptyArray(0u,0u), counters(new Counters()) {}; ~StatsCounter(); - + /** * @brief Changes the reference array for the counting. - * + * * @param Array_ A pointer to an array of class `Array_Type`. */ void reset_array(const Array_Type * Array_); - + void add_counter(Counter f_); void set_counters(Counters * counters_); - + /** * @brief Counter functions * This function recurses through the entries of `Array` and at each step of @@ -84,7 +84,7 @@ class StatsCounter { std::vector< std::string > get_descriptions() const; size_t size() const {return counters->size();}; - + }; #endif diff --git a/include/barry/statscounter-meat.hpp b/include/barry/statscounter-meat.hpp index ab05210..ef7455b 100644 --- a/include/barry/statscounter-meat.hpp +++ b/include/barry/statscounter-meat.hpp @@ -17,7 +17,7 @@ STATSCOUNTER_TEMPLATE(,StatsCounter)( EmptyArray = *Array; EmptyArray.clear(); current_stats = counter.current_stats; - + // We will save the data here counters = new Counters((*counter.counters)); counter_deleted = false; @@ -33,56 +33,56 @@ STATSCOUNTER_TEMPLATE(,~StatsCounter)() STATSCOUNTER_TEMPLATE(void, reset_array)(const Array_Type * Array_) { - + Array = Array_; EmptyArray = *Array_; EmptyArray.clear(); - + return; } STATSCOUNTER_TEMPLATE(void, add_counter)(Counter f_) { - + counters->add_counter(f_); - + return; - + } STATSCOUNTER_TEMPLATE(void, set_counters)(Counters * counters_) { - + // Cleaning up before replacing the memory if (!counter_deleted) delete counters; counter_deleted = true; counters = counters_; - + return; - + } STATSCOUNTER_TEMPLATE(void, count_init)(size_t i,size_t j) { - + // Do we have any counter? if (counters->size() == 0u) throw std::logic_error("No counters added: Cannot count without knowning what to count!"); - + // Iterating through the functions, and updating the set of // statistics. current_stats.resize(counters->size(), 0.0); // change_stats.resize(counters->size(), 0.0); - for (size_t n = 0u; n < counters->size(); ++n) + for (size_t n = 0u; n < counters->size(); ++n) current_stats[n] = counters->operator[](n).init(EmptyArray, i, j); - + return; } STATSCOUNTER_TEMPLATE(void, count_current)(size_t i, size_t j) { - + // Iterating through the functions, and updating the set of // statistics. for (size_t n = 0u; n < counters->size(); ++n) { @@ -92,7 +92,7 @@ STATSCOUNTER_TEMPLATE(void, count_current)(size_t i, size_t j) } return; - + } template @@ -101,7 +101,7 @@ inline std::vector< double > StatsCounter::count_all() if (Array->is_dense()) { - return count_all_dense(); + return count_all_dense(); } else { @@ -113,10 +113,10 @@ inline std::vector< double > StatsCounter::count_all() template inline std::vector< double > StatsCounter::count_all_sparse() { - + // Initializing the counter on the empty array count_init(0u, 0u); - + // Setting it to zero. EmptyArray.clear(false); @@ -126,17 +126,17 @@ inline std::vector< double > StatsCounter::count_all_spars BARRY_DEBUG_VEC_PRINT(this->get_names()); #endif #endif - + // Start iterating through the data for (size_t i = 0; i < Array->nrow(); ++i) { - + const auto & row = Array->row(i, false); // Any element? if (row.size() == 0u) continue; - + // If there's one, then update the statistic, by iterating for (const auto& col: row) { @@ -144,7 +144,7 @@ inline std::vector< double > StatsCounter::count_all_spars // We only insert if it is different from zero if (static_cast(col.second.value) == 0) continue; - + // Adding a cell EmptyArray.insert_cell(i, col.first, col.second, false, false); @@ -161,7 +161,7 @@ inline std::vector< double > StatsCounter::count_all_spars EmptyArray.print(); #endif #endif - #endif + #endif // Computing the change statistics count_current(i, col.first); @@ -171,23 +171,23 @@ inline std::vector< double > StatsCounter::count_all_spars BARRY_DEBUG_VEC_PRINT(current_stats); #endif #endif - - } - + + } + } - + // Adding to the sufficient statistics return current_stats; - + } template inline std::vector< double > StatsCounter::count_all_dense() { - + // Initializing the counter on the empty array count_init(0u, 0u); - + // Setting it to zero. EmptyArray.clear(false); @@ -197,7 +197,7 @@ inline std::vector< double > StatsCounter::count_all_dense BARRY_DEBUG_VEC_PRINT(this->get_names()); #endif #endif - + // Start iterating through the data for (size_t i = 0u; i < Array->nrow(); ++i) { @@ -207,7 +207,7 @@ inline std::vector< double > StatsCounter::count_all_dense // We only insert if it is different from zero if (Array->is_empty(i,j)) continue; - + // Adding a cell EmptyArray.insert_cell(i, j, 1, false, false); @@ -224,7 +224,7 @@ inline std::vector< double > StatsCounter::count_all_dense EmptyArray.print(); #endif #endif - #endif + #endif // Computing the change statistics count_current(i, j); @@ -235,12 +235,12 @@ inline std::vector< double > StatsCounter::count_all_dense #endif #endif } - + } - + // Adding to the sufficient statistics return current_stats; - + } template STATSCOUNTER_TEMPLATE_ARGS() @@ -262,4 +262,4 @@ STATSCOUNTER_TEMPLATE(std::vector< std::string >, get_descriptions)() const #undef STATSCOUNTER_TEMPLATE_ARGS #undef STATSCOUNTER_TEMPLATE -#endif \ No newline at end of file +#endif diff --git a/include/barry/support-bones.hpp b/include/barry/support-bones.hpp index e4e984b..11aae3e 100644 --- a/include/barry/support-bones.hpp +++ b/include/barry/support-bones.hpp @@ -1,4 +1,4 @@ -#ifndef BARRY_SUPPORT_BONES_HPP +#ifndef BARRY_SUPPORT_BONES_HPP #define BARRY_SUPPORT_BONES_HPP 1 template @@ -18,29 +18,29 @@ class Rule; /** * @brief Compute the support of sufficient statistics - * + * * Given an array and a set of counters, this object iterates throughout the * support set of the Array while at the same time computing the support of * the sufficient statitics. - * + * * The members `rule` and `rule_dyn` allow constraining the support. The first * will establish which cells of the array will be used to iterate, for example, * in the case of social networks, self-loops are not allowed, so the entire * diagonal would be fixed to zero, reducing the size of the support. - * + * * In the case of `rule_dyn`, the function will stablish dynamically whether * the current state will be included in the counts or not. For example, this * set of rules can be used to constrain the support to networks that have a - * prescribed degree sequence. - */ + * prescribed degree sequence. + */ template < typename Array_Type = BArray, typename Data_Counter_Type = bool, typename Data_Rule_Type = bool, - typename Data_Rule_Dyn_Type = bool + typename Data_Rule_Dyn_Type = bool > class Support { - + private: void calc_backend_sparse( size_t pos = 0u, @@ -62,17 +62,17 @@ class Support { Counters * counters; ///< Vector of couter functions. Rules * rules; ///< Vector of static rules (cells to iterate). Rules * rules_dyn; ///< Vector of dynamic rules (to include/exclude a realizaton). - + size_t iter_counter = 0u; ///< Number of iterations (internal use only). public: - + size_t N, M; bool delete_counters = true; bool delete_rules = true; bool delete_rules_dyn = true; size_t max_num_elements = BARRY_MAX_NUM_ELEMENTS; - + // Temp variables to reduce memory allocation std::vector< double > current_stats; std::vector< size_t > coordinates_free; @@ -83,7 +83,7 @@ class Support { std::vector< size_t > hashes; std::vector< bool > hashes_initialized; size_t n_counters; - + /**@brief Constructor passing a reference Array. */ Support(const Array_Type & Array_) : @@ -92,7 +92,7 @@ class Support { rules(new Rules()), rules_dyn(new Rules()), N(Array_.nrow()), M(Array_.ncol()), current_stats() {}; - + /**@brief Constructor specifying the dimensions of the array (empty). */ Support(size_t N_, size_t M_) : @@ -101,16 +101,16 @@ class Support { rules(new Rules()), rules_dyn(new Rules()), N(N_), M(M_), current_stats() {}; - + Support() : EmptyArray(0u, 0u), counters(new Counters()), rules(new Rules()), rules_dyn(new Rules()), N(0u), M(0u), current_stats() {}; - + ~Support() { - + if (delete_counters) delete counters; if (delete_rules) @@ -119,27 +119,27 @@ class Support { delete rules_dyn; }; - + void init_support( std::vector< Array_Type > * array_bank = nullptr, std::vector< double > * stats_bank = nullptr ); - + /** * @name Resets the support calculator - * + * * If needed, the counters of a support object can be reused. - * + * * @param Array_ New array over which the support will be computed. */ ///@{ void reset_array(); void reset_array(const Array_Type & Array_); ///@} - + /** - * @name Manage counters - * + * @name Manage counters + * * @param f_ A counter to be added. * @param counters_ A vector of counters to be added. */ @@ -147,10 +147,10 @@ class Support { void add_counter(Counter f_); void set_counters(Counters * counters_); ///@} - + /** - * @name Manage rules - * + * @name Manage rules + * * @param f_ A rule to be added. * @param counters_ A vector of rules to be added. */ @@ -166,32 +166,32 @@ class Support { /** * @brief Computes the entire support - * + * * Not to be used by the user. Sets the starting point in the array * (column-major). - * - * @param array_bank If specified, the counter will add to the vector each + * + * @param array_bank If specified, the counter will add to the vector each * possible state of the array, as it counts. - * + * * @param stats_bank If specified, the counter will add to the vector each * possible set of statistics, as it counts. - * + * */ void calc( std::vector< Array_Type > * array_bank = nullptr, std::vector< double > * stats_bank = nullptr, size_t max_num_elements_ = 0u ); - + const std::vector< double > & get_counts() const; std::vector< double > * get_current_stats(); ///< List current statistics. void print() const; - + const FreqTable< double > & get_data() const; Counters * get_counters(); ///< Vector of couter functions. Rules * get_rules(); ///< Vector of static rules (cells to iterate). Rules * get_rules_dyn(); ///< Vector of dynamic rules (to include/exclude a realizaton). - + }; diff --git a/include/barry/support-meat.hpp b/include/barry/support-meat.hpp index 61bda9b..ab5cd40 100644 --- a/include/barry/support-meat.hpp +++ b/include/barry/support-meat.hpp @@ -9,7 +9,7 @@ inline void Supportiter_counter = 0u; - + // Computing the locations coordinates_free.clear(); coordinates_locked.clear(); @@ -21,7 +21,7 @@ inline void Support 0u) { @@ -32,7 +32,7 @@ inline void Support(coordiantes_n_free)), counters->size() - ); + ); // Adding to the overall count bool include_it = rules_dyn->operator()(EmptyArray, 0u, 0u); @@ -96,10 +96,10 @@ inline void Supportpush_back(EmptyArray); - + if (include_it && (stats_bank != nullptr)) std::copy(current_stats.begin(), current_stats.end(), std::back_inserter(*stats_bank)); @@ -109,19 +109,19 @@ inline void Support inline void Support::reset_array() { - + data.clear(); - + } template inline void Support::reset_array(const Array_Type & Array_) { - + data.clear(); EmptyArray = Array_; N = Array_.nrow(); M = Array_.ncol(); - + } template @@ -130,7 +130,7 @@ inline void Support * array_bank, std::vector< double > * stats_bank ) { - + #ifdef BARRY_USER_INTERRUPT if (++iter_counter % 1000u == 0u) { @@ -141,10 +141,10 @@ inline void Support= coordiantes_n_free) return; - + // We will pass it to the next step, if the iteration makes sense. calc_backend_sparse(pos + 1u, array_bank, stats_bank); - + // Once we have returned, everything will be back as it used to be, so we // treat the data as if nothing has changed. const size_t & coord_i = coordinates_free[pos * 2u]; @@ -170,7 +170,7 @@ inline void Support -DBL_MIN)) { @@ -187,12 +187,12 @@ inline void Supportsize() > 0u) { - + if (rules_dyn->operator()( EmptyArray, coord_i, @@ -208,12 +208,12 @@ inline void Supportpush_back(EmptyArray); - + if (stats_bank != nullptr) std::copy(current_stats.begin(), current_stats.end(), std::back_inserter(*stats_bank)); } - + } else { @@ -225,35 +225,35 @@ inline void Supportpush_back(EmptyArray); - + if (stats_bank != nullptr) std::copy(current_stats.begin(), current_stats.end(), std::back_inserter(*stats_bank)); } - + // Again, we only pass it to the next level iff the next level is not // passed the last step. calc_backend_sparse(pos + 1u, array_bank, stats_bank); - + // We need to restore the state of the cell EmptyArray.rm_cell( coord_i, coord_j, false, false ); - + if (change_stats_different > 0u) { #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd #endif - for (size_t n = 0u; n < n_counters; ++n) + for (size_t n = 0u; n < n_counters; ++n) current_stats[n] -= change_stats[pos * n_counters + n]; } - - + + return; - + } template @@ -269,14 +269,14 @@ inline void Support= coordiantes_n_free) return; - + // We will pass it to the next step, if the iteration makes sense. calc_backend_dense(pos + 1u, array_bank, stats_bank); - + // Once we have returned, everything will be back as it used to be, so we // treat the data as if nothing has changed. const size_t & coord_i = coordinates_free[pos * 2u]; @@ -316,12 +316,12 @@ inline void Supportsize() > 0u) { - + if (rules_dyn->operator()(EmptyArray, coord_i, coord_j)) { @@ -333,12 +333,12 @@ inline void Supportpush_back(EmptyArray); - + if (stats_bank != nullptr) std::copy(current_stats.begin(), current_stats.end(), std::back_inserter(*stats_bank)); } - + } else @@ -352,30 +352,30 @@ inline void Supportpush_back(EmptyArray); - + if (stats_bank != nullptr) std::copy(current_stats.begin(), current_stats.end(), std::back_inserter(*stats_bank)); } - + // Again, we only pass it to the next level iff the next level is not // passed the last step. calc_backend_dense(pos + 1u, array_bank, stats_bank); - + // We need to restore the state of the cell EmptyArray.rm_cell(coord_i, coord_j, false, false); - + if (change_stats_different > 0u) { #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd #endif - for (size_t n = 0u; n < n_counters; ++n) + for (size_t n = 0u; n < n_counters; ++n) current_stats[n] -= change_stats[pos * n_counters + n]; } - + return; - + } template @@ -410,32 +410,32 @@ Support::calc( return; - + } template inline void Support::add_counter( Counter f_ ) { - + counters->add_counter(f_); return; - + } template inline void Support::set_counters( Counters * counters_ ) { - + // Cleaning up before replacing the memory if (delete_counters) delete counters; delete_counters = false; counters = counters_; - + return; - + } ///////////////////////////// @@ -444,70 +444,70 @@ template ::add_rule( Rule * f_ ) { - + rules->add_rule(f_); return; - + } template inline void Support::add_rule( Rule f_ ) { - + rules->add_rule(f_); return; - + } template inline void Support::set_rules( Rules * rules_ ) { - + // Cleaning up before replacing the memory if (delete_rules) delete rules; delete_rules = false; rules = rules_; - + return; - + } template inline void Support::add_rule_dyn( Rule * f_ ) { - + rules_dyn->add_rule(f_); return; - + } template inline void Support::add_rule_dyn( Rule f_ ) { - + rules_dyn->add_rule(f_); return; - + } template inline void Support::set_rules_dyn( Rules * rules_ ) { - + // Cleaning up before replacing the memory if (delete_rules_dyn) delete rules_dyn; delete_rules_dyn = false; rules_dyn = rules_; - + return; - + } template @@ -555,16 +555,16 @@ inline bool Support inline const std::vector< double > & Support::get_counts() const { - - return data.get_data(); - + + return data.get_data(); + } // template // inline const MapVec_type<> * Support::get_counts_ptr() const { - + // return data.get_data_ptr(); - + // } template @@ -592,8 +592,8 @@ inline const FreqTable & Support inline Counters * Support::get_counters() { return this->counters; -} - +} + template inline Rules * Support::get_rules() { return this->rules; @@ -605,4 +605,4 @@ inline Rules * Support* >; /** * @brief A wrapper class to store `source`, `target`, `val` from a `BArray` object. - * + * * @tparam Cell_Type Any type */ template @@ -80,7 +80,7 @@ class Entries { std::vector< size_t > source; std::vector< size_t > target; std::vector< Cell_Type > val; - + Entries() : source(0u), target(0u), val(0u) {}; Entries(size_t n) { source.reserve(n); @@ -88,16 +88,16 @@ class Entries { val.reserve(n); return; }; - + ~Entries() {}; - + void resize(size_t n) { source.resize(n); target.resize(n); val.resize(n); return; } - + }; // Relevant for anything using vecHasher function ------------------------------ @@ -107,35 +107,35 @@ struct vecHasher std::size_t operator()(std::vector< T > const& dat) const noexcept { - + std::hash< T > hasher; std::size_t hash = hasher(dat[0u]); - + // ^ makes bitwise XOR // 0x9e3779b9 is a 32 bit constant (comes from the golden ratio) // << is a shift operator, something like lhs * 2^(rhs) if (dat.size() > 1u) for (size_t i = 1u; i < dat.size(); ++i) hash ^= hasher(dat[i]) + 0x9e3779b9 + (hash<<6) + (hash>>2); - + return hash; - + } }; -template +template using MapVec_type = std::unordered_map< std::vector< Ta >, Tb, vecHasher>; /** * @brief Ascending sorting an array - * + * * It will sort an array solving ties using the next column. Data is * stored column-wise. - * - * @tparam T - * @param v - * @param nrows + * + * @tparam T + * @param v + * @param nrows * @return std::vector The sorting index. */ inline std::vector< size_t > sort_array( @@ -155,8 +155,8 @@ inline std::vector< size_t > sort_array( for (size_t j = 0u; j < ncols; ++j) { if (*(v + (nrows * j + i1+start)) == *(v + (nrows * j + i2 + start))) - continue; - else + continue; + else return *(v + (nrows * j + i1+start)) < *(v + (nrows * j + i2 + start)); } @@ -165,7 +165,7 @@ inline std::vector< size_t > sort_array( return idx; -} +} // Mostly relevant in the case of the stats count functions ------------------- @@ -193,8 +193,8 @@ using Rule_fun_type = std::function using Hasher_fun_type = std::function(const Array_Type &, Data_Type *)>; @@ -211,23 +211,23 @@ inline bool vec_equal( const std::vector< T > & a, const std::vector< T > & b ) { - + if (a.size() != b.size()) { - + std::string err = "-a- and -b- should have the same length. length(a) = " + std::to_string(a.size()) + " and length(b) = " + std::to_string(b.size()) + std::string("."); throw std::length_error(err); } - + size_t i = 0; while (a[i] == b[i]) { if (++i == a.size()) return true; } - + return false; } @@ -237,7 +237,7 @@ inline bool vec_equal_approx( const std::vector< T > & b, double eps = 1e-100 ) { - + if (a.size() != b.size()) { std::string err = "-a- and -b- should have the same length. length(a) = " + @@ -245,13 +245,13 @@ inline bool vec_equal_approx( std::string("."); throw std::length_error(err); } - + size_t i = 0; while (static_cast(std::fabs(a[i] - b[i])) < eps) { if (++i == a.size()) return true; } - + return false; } ///@} @@ -265,16 +265,16 @@ inline T vec_inner_prod( const T * b, size_t n ) { - + double res = 0.0; - #if defined(__OPENMP) || defined(_OPENMP) + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:res) #elif defined(__GNUC__) && !defined(__clang__) #pragma GCC ivdep #endif for (size_t i = 0u; i < n; ++i) res += (*(a + i) * *(b + i)); - + return res; } @@ -288,7 +288,7 @@ inline double vec_inner_prod( const double * b, size_t n ) { - + double res = 0.0; #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:res) @@ -297,10 +297,9 @@ inline double vec_inner_prod( #endif for (size_t i = 0u; i < n; ++i) res += (*(a + i) * *(b + i)); - + return res; } #endif - diff --git a/src/defm-common.hpp b/src/defm-common.hpp index 8fe18fc..7bab35e 100644 --- a/src/defm-common.hpp +++ b/src/defm-common.hpp @@ -3,7 +3,7 @@ using namespace pybind11::literals; -inline void pyprinter(const char * fmt, ...) +inline void pyprinter(const char * fmt, ...) { // Creating a buffer @@ -66,7 +66,7 @@ inline void check_covar( element_access = [](size_t i, size_t j, size_t, size_t ncol) -> size_t { \ return j + i * ncol; \ }; \ - } + } /** @@ -84,4 +84,3 @@ inline void check_covar( #endif - diff --git a/src/main.cpp b/src/main.cpp index ac1859a..dc0e7f0 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,12 +14,12 @@ namespace py = pybind11; std::shared_ptr< defm::DEFM > new_defm( - py::array_t< int > id, + py::array_t< int > id, /* python uses row-major. defm uses col-major. Need to address this with: https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html?highlight=mutable_data#arrays */ - py::array_t< int > y, + py::array_t< int > y, py::array_t< double > x, size_t order = 0u, bool copy = false, @@ -84,7 +84,7 @@ py::array_t< int > simulate( * @param object The DEFM object */ void print_y(const std::shared_ptr< defm::DEFM > & object) { - + auto Y = object->get_Y(); for (size_t i = 0u; i < object->get_n_y(); ++i) std::cout << (*(Y + i)) << " "; diff --git a/src/model.cpp b/src/model.cpp index 0d9466b..d358269 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1,40 +1,40 @@ #include -#include +#include #include "defm-common.hpp" namespace py = pybind11; //' Get sufficient statistics counts -//' +//' //' This function computes the individual counts of the sufficient statistics -//' included in the model. +//' included in the model. //' @param m An object of class [DEFM]. //' @export //' @return A matrix with the counts of the sufficient statistics. //' @examples //' data(valentesnsList) -//' +//' //' mymodel <- new_defm( //' id = valentesnsList$id, //' Y = valentesnsList$Y, //' X = valentesnsList$X, //' order = 1 //' ) -//' +//' //' # Adding the intercept terms and a motif from tobacco to mj //' term_defm_logit_intercept(mymodel) //' term_defm_transition_formula(mymodel, "{y1, 0y2} > {y1, y2}") -//' +//' //' # Initialize the model //' init_defm(mymodel) -//' +//' //' # Get the counts //' head(get_stats(mymodel)) // [[Rcpp::export(rng = false)]] py::array_t get_stats(std::shared_ptr< defm::DEFM > m) { - + // Getting sizes size_t nrows = m->get_n_rows(); size_t ncols = m->nterms(); @@ -43,12 +43,12 @@ py::array_t get_stats(std::shared_ptr< defm::DEFM > m) const int * ID = m->get_ID(); DEFM_WRAP_NUMPY(res, res_ptr, nrows, ncols, double) - + auto target = m->get_stats_target(); // Figure out wether is column or row major DEFM_DEFINE_ACCESS(m); - + size_t i_effective = 0u; size_t n_obs_i = 0u; for (size_t i = 0u; i < nrows; ++i) @@ -98,5 +98,3 @@ void init_get_stats(py::module_ &m) )pbdoc"); } - - diff --git a/src/pydefm/__init__.py b/src/pydefm/__init__.py index 38b343d..306b3ce 100644 --- a/src/pydefm/__init__.py +++ b/src/pydefm/__init__.py @@ -1,19 +1,26 @@ from __future__ import annotations -from ._core import __doc__, __version__, new_defm, print_y, term_formula, get_stats, simulate - +from ._core import ( + __doc__, + __version__, + get_stats, + new_defm, + print_y, + simulate, + term_formula, +) from .defm_mle import defm_mle # from scipy.optimize import minimize __all__ = [ - "__doc__", - "__version__", - "new_defm", - "print_y", - "term_formula", - "get_stats", - "simulate", - "defm_mle", - "defm_mle_fit" - ] + "__doc__", + "__version__", + "defm_mle", + "defm_mle_fit", + "get_stats", + "new_defm", + "print_y", + "simulate", + "term_formula", +] diff --git a/src/pydefm/defm_mle.py b/src/pydefm/defm_mle.py index b8ddcf9..56ff415 100644 --- a/src/pydefm/defm_mle.py +++ b/src/pydefm/defm_mle.py @@ -1,57 +1,65 @@ -from scipy.optimize import minimize +from __future__ import annotations + import numpy as np +from scipy.optimize import minimize + # Function to print an array using two digits def print_array(arr): - return "[" + ", ".join([f"{'' if x < 0 else ' '}{x:.2f}" for x in arr]) + "]" + return "[" + ", ".join([f"{'' if x < 0 else ' '}{x:.2f}" for x in arr]) + "]" + class defm_mle_fit: - """ - Represents the result of a Maximum Likelihood Estimation (MLE) fit in the DEFM model. - - Attributes: - par (list): The estimated parameter values. - or_ (list): The odds ratios corresponding to the estimated parameter values. - se (list): The standard errors of the estimated parameter values. - ll (float): The log-likelihood value of the fit. - optimres (object): The optimization result object containing additional information. - """ - - def __init__(self, par, or_, se, ll, optimres): - self.par = par - self.or_ = or_ - self.se = se - self.ll = ll - self.optimres = optimres - - def __str__(self): - return f"defm_mle_fit\n\tpar =" + \ - print_array(self.par) + \ - ",\n\tor =" + \ - print_array(self.or_) + \ - ",\n\tse ="+ \ - print_array(self.se) + \ - f"\n\tll ={self.ll: .2f}" + """ + Represents the result of a Maximum Likelihood Estimation (MLE) fit in the DEFM model. + + Attributes: + par (list): The estimated parameter values. + or_ (list): The odds ratios corresponding to the estimated parameter values. + se (list): The standard errors of the estimated parameter values. + ll (float): The log-likelihood value of the fit. + optimres (object): The optimization result object containing additional information. + """ + + def __init__(self, par, or_, se, ll, optimres): + self.par = par + self.or_ = or_ + self.se = se + self.ll = ll + self.optimres = optimres + + def __str__(self): + return ( + "defm_mle_fit\n\tpar =" + + print_array(self.par) + + ",\n\tor =" + + print_array(self.or_) + + ",\n\tse =" + + print_array(self.se) + + f"\n\tll ={self.ll: .2f}" + ) + def defm_mle(obj, par): - """ - Maximum Likelihood Estimation using the defm method. - - Parameters: - obj (object): The object containing the likelihood function. - par (array-like): The initial parameter values. - - Returns: - defm_mle_fit: An instance of the defm_mle_fit class containing the estimated parameters, odds ratios, standard errors, log-likelihood value, and optimization results. - """ - def ll(par): - return -obj.likelihood(par, as_log=True) - - res = minimize(ll, par) - return defm_mle_fit( - par=res.x, - or_=np.exp(res.x), - se=np.sqrt(np.diag(res.hess_inv)), - ll=-res.fun, - optimres=res - ) + """ + Maximum Likelihood Estimation using the defm method. + + Parameters: + obj (object): The object containing the likelihood function. + par (array-like): The initial parameter values. + + Returns: + defm_mle_fit: An instance of the defm_mle_fit class containing the estimated parameters, odds ratios, standard errors, log-likelihood value, and optimization results. + """ + + def ll(par): + return -obj.likelihood(par, as_log=True) + + res = minimize(ll, par) + return defm_mle_fit( + par=res.x, + or_=np.exp(res.x), + se=np.sqrt(np.diag(res.hess_inv)), + ll=-res.fun, + optimres=res, + ) diff --git a/tests/test_basic.py b/tests/test_basic.py index c0c1f1b..edbee4c 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,17 +1,18 @@ from __future__ import annotations + import warnings warnings.filterwarnings("ignore", category=ImportWarning) -import pydefm as m import numpy as np +import pydefm as m y = np.array([0, 10, 3]) x = np.array([1, 2.0, 3.4]) id = np.array([11, 2, 3]) -obj = m.new_defm(id, y, x, column_major = False) +obj = m.new_defm(id, y, x, column_major=False) obj.print() @@ -22,4 +23,5 @@ def test_version(): assert m.__version__ == "0.0.1" -print("Everything passed") \ No newline at end of file + +print("Everything passed")