Quelle est la meilleure façon de créer un tableau clairsemé en C++ ?
https://stackoverflow.com/questions/4306
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Je travaille sur un projet qui nécessite la manipulation d'énormes matrices, notamment une sommation pyramidale pour un calcul de copule.
En bref, je dois garder une trace d'un nombre relativement petit de valeurs (généralement une valeur de 1 et, dans de rares cas, supérieure à 1) dans une mer de zéros dans la matrice (tableau multidimensionnel).
Un tableau clairsemé permet à l'utilisateur de stocker un petit nombre de valeurs et de supposer que tous les enregistrements non définis sont une valeur prédéfinie.Puisqu'il n'est physiquement pas possible de stocker toutes les valeurs en mémoire, je dois stocker uniquement les quelques éléments non nuls.Cela pourrait représenter plusieurs millions d’entrées.
La vitesse est une priorité majeure et j'aimerais également choisir dynamiquement le nombre de variables dans la classe au moment de l'exécution.
Je travaille actuellement sur un système qui utilise un arbre de recherche binaire (b-tree) pour stocker les entrées.Quelqu'un connaît-il un meilleur système ?
La solution
Pour C++, une carte fonctionne bien.Plusieurs millions d'objets ne seront pas un problème.10 millions d'éléments ont pris environ 4,4 secondes et environ 57 Mo sur mon ordinateur.
Mon application de test est la suivante :
#include <stdio.h>
#include <stdlib.h>
#include <map>
class triple {
public:
int x;
int y;
int z;
bool operator<(const triple &other) const {
if (x < other.x) return true;
if (other.x < x) return false;
if (y < other.y) return true;
if (other.y < y) return false;
return z < other.z;
}
};
int main(int, char**)
{
std::map<triple,int> data;
triple point;
int i;
for (i = 0; i < 10000000; ++i) {
point.x = rand();
point.y = rand();
point.z = rand();
//printf("%d %d %d %d\n", i, point.x, point.y, point.z);
data[point] = i;
}
return 0;
}
Maintenant pour choisir dynamiquement le nombre de variables, la solution la plus simple est de représenter index sous forme de chaîne, puis utilisez la chaîne comme clé pour la carte.Par exemple, un élément situé en [23][55] peut être représenté via la chaîne "23,55".On peut également étendre cette solution à des dimensions supérieures ;comme pour les trois dimensions, un index arbitraire ressemblera à "34,45,56".Une mise en œuvre simple de cette technique est la suivante :
std::map data<string,int> data;
char ix[100];
sprintf(ix, "%d,%d", x, y); // 2 vars
data[ix] = i;
sprintf(ix, "%d,%d,%d", x, y, z); // 3 vars
data[ix] = i;
Autres conseils
À titre général, une méthode utilisant des chaînes comme indices est en fait très lent.Une solution beaucoup plus efficace mais par ailleurs équivalente serait d'utiliser des vecteurs/tableaux.Il n'est absolument pas nécessaire d'écrire les indices dans une chaîne.
typedef vector<size_t> index_t;
struct index_cmp_t : binary_function<index_t, index_t, bool> {
bool operator ()(index_t const& a, index_t const& b) const {
for (index_t::size_type i = 0; i < a.size(); ++i)
if (a[i] != b[i])
return a[i] < b[i];
return false;
}
};
map<index_t, int, index_cmp_t> data;
index_t i(dims);
i[0] = 1;
i[1] = 2;
// … etc.
data[i] = 42;
Cependant, en utilisant un map en pratique, ce n'est souvent pas très efficace en raison de la mise en œuvre en termes d'arbre de recherche binaire équilibré.Une structure de données plus performante dans ce cas serait une table de hachage, telle que fournie par std::unordered_map.
Boost a une implémentation basée sur un modèle de BLAS appelée uBLAS qui contient une matrice clairsemée.
http://www.boost.org/doc/libs/1_36_0/libs/numeric/ublas/doc/index.htm
Petit détail dans la comparaison des indices.Il faut faire une comparaison lexicographique, sinon :
a= (1, 2, 1); b= (2, 1, 2);
(a<b) == (b<a) is true, but b!=a
Modifier:La comparaison devrait donc probablement être :
return lhs.x<rhs.x
? true
: lhs.x==rhs.x
? lhs.y<rhs.y
? true
: lhs.y==rhs.y
? lhs.z<rhs.z
: false
: false
Les tables de hachage ont une insertion et une recherche rapides.Vous pouvez écrire une simple fonction de hachage puisque vous savez que vous n'aurez affaire qu'à des paires entières comme clés.
Propre est une bibliothèque d'algèbre linéaire C++ qui possède un mise en œuvre d'une matrice clairsemée.Il prend même en charge les opérations matricielles et les solveurs (factorisation LU, etc.) optimisés pour les matrices clairsemées.
La liste complète des solutions peut être trouvée sur Wikipédia.Pour plus de commodité, j’ai cité les sections pertinentes comme suit.
https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_.28DOK.29
Dictionnaire de clés (DOK)
Dok se compose d'un dictionnaire qui mappe (ligne, colonne) -paires à la valeur des éléments.Les éléments qui manquent du dictionnaire sont considérés comme nuls.Le format est bon pour construire progressivement une matrice clairsemée dans un ordre aléatoire, mais médiocre pour itération sur les valeurs non nulles dans l'ordre lexicographique.On construit généralement une matrice dans ce format, puis se convertit en un autre format plus efficace pour le traitement. [1
Liste des listes (LIL)
LIL stocke une liste par ligne, chaque entrée contenant l'index de colonne et la valeur.En règle générale, ces entrées sont maintenues par index de colonne pour une recherche plus rapide.Il s'agit d'un autre format bon pour la construction de matrice incrémentielle. [2
Liste de coordonnées (COO)
COO stocke une liste de tuples (ligne, colonne, valeur).Idéalement, les entrées sont triées (par index de ligne, puis index de colonne) pour améliorer les temps d'accès aléatoires.Il s'agit d'un autre format qui est bon pour la construction de matrice incrémentielle. [3
Ligne clairsemée compressée (format CSR, CRS ou Yale)
Le format de stockage de lignes de ligne (CRS) compressé compressé (CSR) ou comprimé (CRS) représente une matrice M de trois tableaux (unidimensionnels), qui contiennent respectivement des valeurs non nulles, des étendus des lignes et des indices de colonne.Il est similaire au COO, mais comprime les indices de ligne, d'où le nom.Ce format permet l'accès rapide à la ligne et les multiplications de vecteur matricielle (MX).
La meilleure façon d’implémenter des matrices clairsemées est de ne pas les implémenter – du moins pas par vous-même.Je suggérerais à BLAS (qui, je pense, fait partie de LAPACK) qui peut gérer des matrices vraiment énormes.
Puisque seules les valeurs avec [a][b][c]...[w][x][y][z] ont des conséquences, nous stockons uniquement l'indice lui-même, pas la valeur 1 qui est à peu près partout - toujours le même + pas moyen de le hacher.Notant que la malédiction de la dimensionnalité est présente, suggérez d'utiliser un outil établi NIST ou Boost, lisez au moins les sources pour éviter des erreurs inutiles.
Si le travail doit capturer les distributions de dépendance temporelle et les tendances paramétriques d'ensembles de données inconnus, alors une carte ou un arbre B avec une racine uni-valuée n'est probablement pas pratique.Nous ne pouvons stocker que l'indice lui-même, haché si l'ordre (sensibilité pour la présentation) peut être subordonné à la réduction du domaine temporel au moment de l'exécution, pour les 1 valeurs.Étant donné que les valeurs non nulles autres qu'une sont peu nombreuses, un candidat évident pour celles-ci est toute structure de données que vous pouvez trouver facilement et comprendre.Si l'ensemble de données est vraiment de la taille d'un vaste univers, je suggère une sorte de fenêtre coulissante qui gère vous-même les fichiers/disques/persistant-io, en déplaçant des parties des données dans la portée selon les besoins.(écrire du code que vous pouvez comprendre) Si vous vous engagez à fournir une solution réelle à un groupe de travail, ne pas le faire vous laisse à la merci de systèmes d'exploitation grand public qui ont pour seul objectif de vous enlever votre déjeuner.
Voici une implémentation relativement simple qui devrait fournir une recherche raisonnablement rapide (à l'aide d'une table de hachage) ainsi qu'une itération rapide sur des éléments non nuls dans une ligne/colonne.
// Copyright 2014 Leo Osvald
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_
#define UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_
#include <algorithm>
#include <limits>
#include <map>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>
// A simple time-efficient implementation of an immutable sparse matrix
// Provides efficient iteration of non-zero elements by rows/cols,
// e.g. to iterate over a range [row_from, row_to) x [col_from, col_to):
// for (int row = row_from; row < row_to; ++row) {
// for (auto col_range = sm.nonzero_col_range(row, col_from, col_to);
// col_range.first != col_range.second; ++col_range.first) {
// int col = *col_range.first;
// // use sm(row, col)
// ...
// }
template<typename T = double, class Coord = int>
class SparseMatrix {
struct PointHasher;
typedef std::map< Coord, std::vector<Coord> > NonZeroList;
typedef std::pair<Coord, Coord> Point;
public:
typedef T ValueType;
typedef Coord CoordType;
typedef typename NonZeroList::mapped_type::const_iterator CoordIter;
typedef std::pair<CoordIter, CoordIter> CoordIterRange;
SparseMatrix() = default;
// Reads a matrix stored in MatrixMarket-like format, i.e.:
// <num_rows> <num_cols> <num_entries>
// <row_1> <col_1> <val_1>
// ...
// Note: the header (lines starting with '%' are ignored).
template<class InputStream, size_t max_line_length = 1024>
void Init(InputStream& is) {
rows_.clear(), cols_.clear();
values_.clear();
// skip the header (lines beginning with '%', if any)
decltype(is.tellg()) offset = 0;
for (char buf[max_line_length + 1];
is.getline(buf, sizeof(buf)) && buf[0] == '%'; )
offset = is.tellg();
is.seekg(offset);
size_t n;
is >> row_count_ >> col_count_ >> n;
values_.reserve(n);
while (n--) {
Coord row, col;
typename std::remove_cv<T>::type val;
is >> row >> col >> val;
values_[Point(--row, --col)] = val;
rows_[col].push_back(row);
cols_[row].push_back(col);
}
SortAndShrink(rows_);
SortAndShrink(cols_);
}
const T& operator()(const Coord& row, const Coord& col) const {
static const T kZero = T();
auto it = values_.find(Point(row, col));
if (it != values_.end())
return it->second;
return kZero;
}
CoordIterRange
nonzero_col_range(Coord row, Coord col_from, Coord col_to) const {
CoordIterRange r;
GetRange(cols_, row, col_from, col_to, &r);
return r;
}
CoordIterRange
nonzero_row_range(Coord col, Coord row_from, Coord row_to) const {
CoordIterRange r;
GetRange(rows_, col, row_from, row_to, &r);
return r;
}
Coord row_count() const { return row_count_; }
Coord col_count() const { return col_count_; }
size_t nonzero_count() const { return values_.size(); }
size_t element_count() const { return size_t(row_count_) * col_count_; }
private:
typedef std::unordered_map<Point,
typename std::remove_cv<T>::type,
PointHasher> ValueMap;
struct PointHasher {
size_t operator()(const Point& p) const {
return p.first << (std::numeric_limits<Coord>::digits >> 1) ^ p.second;
}
};
static void SortAndShrink(NonZeroList& list) {
for (auto& it : list) {
auto& indices = it.second;
indices.shrink_to_fit();
std::sort(indices.begin(), indices.end());
}
// insert a sentinel vector to handle the case of all zeroes
if (list.empty())
list.emplace(Coord(), std::vector<Coord>(Coord()));
}
static void GetRange(const NonZeroList& list, Coord i, Coord from, Coord to,
CoordIterRange* r) {
auto lr = list.equal_range(i);
if (lr.first == lr.second) {
r->first = r->second = list.begin()->second.end();
return;
}
auto begin = lr.first->second.begin(), end = lr.first->second.end();
r->first = lower_bound(begin, end, from);
r->second = lower_bound(r->first, end, to);
}
ValueMap values_;
NonZeroList rows_, cols_;
Coord row_count_, col_count_;
};
#endif /* UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_ */
Pour simplifier, c'est immutable, mais vous pouvez le rendre mutable ;assurez-vous de changer std::vector à std::set si vous voulez des "insertions" raisonnablement efficaces (changer un zéro en un non nul).
Je suggérerais de faire quelque chose comme :
typedef std::tuple<int, int, int> coord_t;
typedef boost::hash<coord_t> coord_hash_t;
typedef std::unordered_map<coord_hash_t, int, c_hash_t> sparse_array_t;
sparse_array_t the_data;
the_data[ { x, y, z } ] = 1; /* list-initialization is cool */
for( const auto& element : the_data ) {
int xx, yy, zz, val;
std::tie( std::tie( xx, yy, zz ), val ) = element;
/* ... */
}
Pour vous aider à garder vos données clairsemées, vous souhaiterez peut-être écrire une sous-classe de unorderd_map, dont les itérateurs ignorent (et effacent) automatiquement tous les éléments ayant une valeur de 0.
