diff --git a/src/database/database.cpp b/src/database/database.cpp
index 3b9cf57..6d147f8 100644
--- a/src/database/database.cpp
+++ b/src/database/database.cpp
@@ -1,77 +1,84 @@
/*-------------------------------------------------------
- Module : database
- File : database.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "database.hpp"
#include <iostream>
#include "selector.hpp"
#include "canonicalizer.hpp"
#include "switch_basis.hpp"
#include "mineral_selector.hpp"
#include "reader.hpp"
namespace specmicp {
namespace database {
void Database::make_canonical()
{
canonicalize_database(data);
}
void Database::change_basis(std::vector<int>& new_basis)
{
BasisSwitcher(data).change_basis(new_basis);
}
// this is a proxy to a DatabaseSelector
void Database::remove_components(const std::vector<int>& id_components_to_remove)
{
DatabaseSelector selector(data);
selector.remove_component(id_components_to_remove);
}
+// this is a proxy to a DatabaseSelector
+void Database::keep_only_components(const std::vector<int>& id_components_to_keep)
+{
+ DatabaseSelector selector(data);
+ selector.keep_only_component(id_components_to_keep);
+}
+
void Database::minerals_keep_only(const std::vector<std::string>& minerals_to_keep)
{
MineralSelector(data).keep_only(minerals_to_keep);
}
void Database::swap_components(std::map<std::string, std::string>& swap_to_make)
{
std::vector<int> id_swap;
id_swap.reserve(data->nb_component);
for (int i=0; i<data->nb_component; ++i)
{
id_swap.push_back(i);
}
for (auto it=swap_to_make.begin(); it!=swap_to_make.end(); ++it)
{
if (it->first == "H2O") {
throw std::invalid_argument("Basis switch : water cannot be swapped from basis");
}
const int idc = safe_label_to_id(it->first, data->labels_basis);
const int idaq = data->nb_component + label_to_id(it->second, data->labels_aqueous);
id_swap[idc] = idaq;
}
change_basis(id_swap);
}
void Database::parse_database(std::string filepath)
{
DataReader reader(filepath);
reader.parse();
- data = reader.get_database();
+ set_database(reader.get_database());
}
} // end namespace database
} // end namespace specmicp
diff --git a/src/database/database.hpp b/src/database/database.hpp
index 7fdb0ca..2d1fab3 100644
--- a/src/database/database.hpp
+++ b/src/database/database.hpp
@@ -1,86 +1,87 @@
/*-------------------------------------------------------
- Module : database
- File : database.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_DATABASE_DATABASE_HPP
#define SPECMICP_DATABASE_DATABASE_HPP
//! \file database.hpp Database management
#include "module.hpp"
namespace specmicp {
//! \namespace database Database management
namespace database {
//! \brief Database management
//!
//! Prepare the database for the computation.
//! This class is a proxy for one-tasked modules.
class Database: public DatabaseModule
{
public:
//! \brief Default constructor
//!
//! the method parse_database must be called to parse the json database
Database() {}
//! \brief Initialise the database py parsing 'filepath'
Database(std::string filepath) {
parse_database(filepath);
make_canonical();
}
//! \brief initialise the database with the raw database
Database(std::shared_ptr<DataContainer> raw_data):
- data(raw_data)
+ DatabaseModule(raw_data)
{}
//! \brief Parse the database 'filepath'
void parse_database(std::string filepath);
//! \brief Transform the database into canonical form
//!
//! This process is necessary before any computation
void make_canonical();
//! \brief Return the database
//!
//! Return a smart pointer of a DataCotnainer instance
//! Note : this is a read write access, be careful
std::shared_ptr<DataContainer> get_database() {return data;}
//! \brief Change the basis
//!
//! @param new_basis list of id of the new basis
//!
//! The new basis is a list of id, id = id_component for no swapping
//! or id = id_aqueous + nb_component for swapping a secondary species
void change_basis(std::vector<int>& new_basis);
//! \brief Swap some component in the basis
void swap_components(std::map<std::string, std::string>& swap_to_make);
- //! \brief Remove components non present in the system
+ //! \brief Remove components not present in the system
void remove_components(const std::vector<int>& id_components_to_remove);
+ //! \brief Keep only components in the id_components to keep
+ void keep_only_components(const std::vector<int>& id_components_to_keep);
+
//! \brief Keep only some minerals at equilibrium
//!
//! The effect is to flag all the other minerals as "kinetic"
void minerals_keep_only(const std::vector<std::string>& minerals_to_keep);
-private:
- std::shared_ptr<DataContainer> data; //!< The data
};
} // end namespace database
} // end namespace specmicp
#endif // SPECMICP_DATABASE_DATABASE_HPP
diff --git a/src/database/module.hpp b/src/database/module.hpp
index 43649ce..ea06fca 100644
--- a/src/database/module.hpp
+++ b/src/database/module.hpp
@@ -1,102 +1,105 @@
/*-------------------------------------------------------
- Module : database
- File : module
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_DATAABASE_MODULE_HPP
#define SPECMICP_DATAABASE_MODULE_HPP
//! \file module.hpp base class for a database module
#include <memory>
#include "data_container.hpp"
namespace specmicp {
namespace database {
//! \brief Base class for a database module
class DatabaseModule
{
public:
//! The type of a smart pointer to a raw database
using RawDatabasePtr = std::shared_ptr<DataContainer>;
//! Index that indicates no known species
static const int no_species = -1;
//! Create a blank database
DatabaseModule() {
reset_database();
}
//! Use the database thedata
DatabaseModule(RawDatabasePtr thedata): data(thedata) {}
+ //! \brief set the database
+ void set_database(RawDatabasePtr thedata) {data = thedata;}
+
//! \brief (Re)set the database - create a blank database
void reset_database() {data = std::make_shared<DataContainer>();}
//! \brief Return true if the database is canonical
bool is_database_canonical() {return data->is_canonical;}
//! \brief Return the id of species "label" in list_labels
//!
//! return no_species if the species cannot be found in the database
static int label_to_id(const std::string& label, const std::vector<std::string>& list_labels)
{
auto search = std::find(list_labels.begin(), list_labels.end(), label);
if (search == list_labels.end()) return DatabaseModule::no_species;
else return (search - list_labels.begin());
}
//! \brief Return the id of species "label" in list_labels
//!
//! Throw std::invalid_argument if the species cannot be found in the list
static int safe_label_to_id(const std::string& label, const std::vector<std::string>& list_labels)
{
auto search = std::find(list_labels.begin(), list_labels.end(), label);
if (search == list_labels.end()) throw std::invalid_argument("Unknown species : '"+label+"'.");
else return (search - list_labels.begin());
}
//! \brief Return the id of the mineral "label"
//!
//! Return no_species if label is not a mineral (at equilibrium)
int mineral_label_to_id(const std::string& label)
{
return label_to_id(label, data->labels_minerals);
}
//! \brief Return the id of the mineral "label"
//!
//! Return no_species if label is not a mineral (at equilibrium)
int mineral_kinetic_label_to_id(const std::string& label)
{
return label_to_id(label, data->labels_minerals_kinetic);
}
//! \brief Return the id of the component "label"
//!
//! Return no_species if "label" is not a component
int component_label_to_id(const std::string& label)
{
return label_to_id(label, data->labels_basis);
}
//! \brief Return the id of the secondary aqueous species "label"
//!
//! Return no_species if "label" is not an aqueous species.
//! A component (of the basis) is not considered as a secondary aqueous species.
int aqueous_label_to_it(const std::string& label)
{
return label_to_id(label, data->labels_aqueous);
}
protected:
RawDatabasePtr data; //! The database
};
} // end namespace database
} // end namespace specmicp
#endif // SPECMICP_DATAABASE_MODULE_HPP
diff --git a/src/database/selector.cpp b/src/database/selector.cpp
index 8b236c7..7ac5d84 100644
--- a/src/database/selector.cpp
+++ b/src/database/selector.cpp
@@ -1,181 +1,197 @@
/*-------------------------------------------------------
- Module : database
- File : selector.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "selector.hpp"
namespace specmicp {
namespace database {
void DatabaseSelector::analyse_component(const std::vector<int> &id_component_to_remove)
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
m_is_component_to_remove[*it] = true;
}
m_nb_component_to_keep = data->nb_component - id_component_to_remove.size();
int new_i = 0;
for (int i=0; i<data->nb_component; ++i)
{
if (is_component_to_remove(i)) continue;
data->param_aq.row(new_i) = data->param_aq.row(i);
++new_i;
}
assert(new_i == nb_component_to_keep());
}
void DatabaseSelector::select_aqueous(const std::vector<int>& id_component_to_remove)
{
// first we select which aqueous species we should remove
std::vector<bool> aqueous_to_remove(data->nb_aqueous, false);
for (int j=0; j<data->nb_aqueous; ++j)
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_aqueous(j, *it) != 0)
{
aqueous_to_remove[j] = true;
break;
}
}
}
// then we remove them, in two steps
const int nb_aq_to_keep = std::count(aqueous_to_remove.begin(), aqueous_to_remove.end(), false);
int new_j = 0;
// 1) first we copy data in the beginning of the arrays
for (int j=0; j<data->nb_aqueous; ++j)
{
if (aqueous_to_remove[j]) continue;
data->logk_aqueous(new_j) = data->logk_aqueous(j);
data->param_aq.row(nb_component_to_keep()+new_j) = data->param_aq.row(data->nb_component+j);
data->labels_aqueous[new_j] = data->labels_aqueous[j];
int new_i = 0;
for (int i=0; i<data->nb_component; ++i)
{
if (is_component_to_remove(i)) continue;
data->nu_aqueous(new_j, new_i) = data->nu_aqueous(j, i);
++new_i;
}
assert(new_i == nb_component_to_keep());
++new_j;
}
assert(new_j == nb_aq_to_keep);
// 2) then we resize the arrays
data->nu_aqueous.conservativeResize(nb_aq_to_keep, nb_component_to_keep());
data->logk_aqueous.conservativeResize(nb_aq_to_keep);
data->labels_aqueous.resize(nb_aq_to_keep);
data->param_aq.conservativeResize(nb_component_to_keep()+nb_aq_to_keep, Eigen::NoChange);
data->nb_aqueous = nb_aq_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_minerals(const std::vector<int>& id_component_to_remove)
{
// first we select which minerals we should remove
std::vector<bool> minerals_to_remove(data->nb_mineral, false);
for (int m=0; m<data->nb_mineral; ++m)
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_mineral(m, *it) != 0)
{
minerals_to_remove[m] = true;
break;
}
}
}
// then we remove them, in two steps
const int nb_min_to_keep = std::count(minerals_to_remove.begin(), minerals_to_remove.end(), false);
int new_m = 0;
// 1) first we copy data in the beginning of the arrays
for (int m=0; m<data->nb_mineral; ++m)
{
if (minerals_to_remove[m]) continue;
data->logk_mineral(new_m) = data->logk_mineral(m);
data->labels_minerals[new_m] = data->labels_minerals[m];
int new_i = 0;
for (int i=0; i<data->nb_component; ++i)
{
if (is_component_to_remove(i)) continue;
data->nu_mineral(new_m, new_i) = data->nu_mineral(m, i);
++new_i;
}
assert(new_i == nb_component_to_keep());
++new_m;
}
assert(new_m == nb_min_to_keep);
// 2) then we resize the arrays
data->nu_mineral.conservativeResize(nb_min_to_keep, nb_component_to_keep());
data->logk_mineral.conservativeResize(nb_min_to_keep);
data->labels_minerals.resize(nb_min_to_keep);
data->nb_mineral = nb_min_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_minerals_kinetic(const std::vector<int>& id_component_to_remove)
{
// first we select which minerals we should remove
std::vector<bool> minerals_to_remove(data->nb_mineral_kinetic, false);
for (int m=0; m<data->nb_mineral_kinetic; ++m)
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_mineral_kinetic(m, *it) != 0)
{
minerals_to_remove[m] = true;
break;
}
}
}
// then we remove them, in two steps
const int nb_min_to_keep = std::count(minerals_to_remove.begin(), minerals_to_remove.end(), false);
int new_m = 0;
// 1) first we copy data in the beginning of the arrays
for (int m=0; m<data->nb_mineral_kinetic; ++m)
{
if (minerals_to_remove[m]) continue;
data->logk_mineral_kinetic(new_m) = data->logk_mineral_kinetic(m);
data->labels_minerals_kinetic[new_m] = data->labels_minerals_kinetic[m];
int new_i = 0;
for (int i=0; i<data->nb_component; ++i)
{
if (is_component_to_remove(i)) continue;
data->nu_mineral_kinetic(new_m, new_i) = data->nu_mineral_kinetic(m, i);
++new_i;
}
assert(new_i == nb_component_to_keep());
++new_m;
}
assert(new_m == nb_min_to_keep);
// 2) then we resize the arrays
data->nu_mineral_kinetic.conservativeResize(nb_min_to_keep, nb_component_to_keep());
data->logk_mineral_kinetic.conservativeResize(nb_min_to_keep);
data->labels_minerals_kinetic.resize(nb_min_to_keep);
data->nb_mineral_kinetic = nb_min_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_components()
{
int new_i = 0;
for (int i=0; i<data->nb_component; ++i)
{
if (is_component_to_remove(i)) continue;
data->labels_basis[new_i] = data->labels_basis[i];
++new_i;
}
assert(new_i == nb_component_to_keep());
data->labels_basis.resize(nb_component_to_keep());
data->nb_component = nb_component_to_keep();
}
+void DatabaseSelector::keep_only_component(const std::vector<int>& id_to_keep)
+{
+ // First build the list of components to remove
+ std::vector<int> id_to_remove;
+ id_to_remove.reserve(data->nb_component - id_to_keep.size());
+ for (int id=0; id<data->nb_component; ++id)
+ {
+ auto search = std::find(id_to_keep.begin(), id_to_keep.end(), id);
+ if (search == id_to_keep.end())
+ id_to_remove.push_back(id);
+ }
+ // Then remove them
+ remove_component(id_to_remove);
+}
+
+
} // end namespace database
} // end namespace specmicp
diff --git a/src/database/selector.hpp b/src/database/selector.hpp
index 50d3233..55c6bcd 100644
--- a/src/database/selector.hpp
+++ b/src/database/selector.hpp
@@ -1,71 +1,75 @@
/*-------------------------------------------------------
- Module : database
- File : selector.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_DATABASE_SELECTOR_HPP
#define SPECMICP_DATABASE_SELECTOR_HPP
//! \file selector.hpp select species int the database
#include <vector>
#include "module.hpp"
namespace specmicp {
namespace database {
//! \brief Select species in the database
//!
//! This is a special class with the only function is to remove some components from the database
class DatabaseSelector: public DatabaseModule
{
public:
DatabaseSelector(RawDatabasePtr thedata):
DatabaseModule(thedata),
m_is_component_to_remove(thedata->nb_component, false)
{}
+ //! Remove compoment in the list "id_component_to_remove"
void remove_component(const std::vector<int>& id_component_to_remove)
{
analyse_component(id_component_to_remove);
select_aqueous(id_component_to_remove);
select_minerals(id_component_to_remove);
select_minerals_kinetic(id_component_to_remove);
select_components();
}
+ //! keep only components in the "id_to_keep" list
+ void keep_only_component(const std::vector<int>& id_to_keep);
+
private:
//! \brief Analyse wich component should be removed
void analyse_component(const std::vector<int>& id_component_to_remove);
//! \brief Select the correct aqueous species
void select_aqueous(const std::vector<int>& id_component_to_remove);
//! \brief Select the correct minerals
void select_minerals(const std::vector<int>& id_component_to_remove);
//! \brief Select the correct minerals governed by kinetics
void select_minerals_kinetic(const std::vector<int>& id_component_to_remove);
//! \brief Select components
//!
//! This should be the last operation
void select_components();
//! \brief Return true if the component must be removed
bool is_component_to_remove(int id) {return m_is_component_to_remove[id];}
int nb_component_to_keep() {return m_nb_component_to_keep;}
std::vector<bool> m_is_component_to_remove;
int m_nb_component_to_keep;
};
} // end namespace database
} // end namespace specmicp
#endif // SPECMICP_DATABASE_SELECTOR_HPP