diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 0157edd8b889d01de253f8545b06cbf91b9b8416..156103100b8885036931a2cbb148adc8371a05a3 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -42,9 +42,9 @@ flags = [
'-I',
'.',
'-isystem',
- '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1'
+ '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
'-isystem',
- '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1/bits'
+ '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1/bits',
'-I',
'src',
# '-I',
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 8ddb2938f45b3cf7740becfcc94356b938a99ca9..e72bd5e82cf0098c22fc5d92d20c61d1a76e0ba9 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -30,6 +30,10 @@ class fvm_cell {
using index_view = typename index_type::view_type;
fvm_cell(nest::mc::cell const& cell);
+
+ private:
+
+ index_type parent_index_;
};
////////////////////////////////////////////////////////////
@@ -39,7 +43,7 @@ class fvm_cell {
template <typename T, typename I>
fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
{
- auto const& parent_index = cell.parent_index();
+ parent_index_ = cell.parent_index();
auto const& segment_index = cell.segment_index();
auto num_fv = segment_index.back();
@@ -65,7 +69,7 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
for(auto c : cable->compartments()) {
auto rhs = segment_index[seg_idx] + c.index;
- auto lhs = parent_index[rhs];
+ auto lhs = parent_index_[rhs];
auto rad_C = math::mean(left(c.radius), right(c.radius));
auto len = c.length/2;
@@ -85,9 +89,8 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
++seg_idx;
}
- //std::cout << "volumes " << fv_volumes << " : " << sum(fv_volumes) << "\n";
- std::cout << "areas " << sum(fv_areas) << ", " << cell.volume() << "\n";
- std::cout << "volumes " << sum(fv_volumes) << ", " << cell.area() << "\n";
+ //std::cout << "areas " << algorithms::sum(fv_volumes) << ", " << cell.volume() << "\n";
+ //std::cout << "volumes " << algorithms::sum(fv_areas) << ", " << cell.area() << "\n";
}
} // namespace fvm
diff --git a/src/math.hpp b/src/math.hpp
index af0c402c40b603c8397fbeeb631695e68fa8a58b..ec737fa2e9dd3c2b0ce49f6708da3e235c1172ba 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -29,6 +29,14 @@ namespace math {
return a*a*a;
}
+ // area of a circle
+ // pi r-squared
+ template <typename T>
+ T constexpr area_circle(T r)
+ {
+ return pi<T>() * square(r);
+ }
+
// volume of a conic frustrum
template <typename T>
T constexpr area_frustrum(T L, T r1, T r2)
diff --git a/src/matrix.hpp b/src/matrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1b1daad2fc17af2d8a3ad87ced5fc03435deca96
--- /dev/null
+++ b/src/matrix.hpp
@@ -0,0 +1,87 @@
+#pragma once
+
+#include "../vector/include/Vector.hpp"
+
+namespace nest {
+namespace mc {
+
+/// matrix storage structure
+template<typename T, typename I>
+class matrix {
+ public:
+
+ // define basic types
+ using value_type = T;
+ using size_type = I;
+
+ // define storage types
+ using vector_type = memory::HostVector<value_type>;
+ using view_type = typename vector_type::view_type;
+ using index_type = memory::HostVector<size_type>;
+ using index_view = typename index_type::view_type;
+
+
+ /// the dimension of the matrix (i.e. the number of rows or colums)
+ size_type size() const
+ {
+ return parent_indices_.size();
+ }
+
+ /// the total memory used to store the matrix
+ std::size_t memory() const
+ {
+ auto s = sizeof(value_type) * data_.size();
+ s += sizeof(size_type) * parent_indices_.size();
+ return s + sizeof(matrix);
+ }
+
+ /// the number of cell matrices that have been packed together
+ size_type num_cells() const
+ {
+ return cell_index_.size();
+ }
+
+ private:
+
+ /// calculate the number of values per sub-array plus padding
+ size_type reservation() const
+ {
+ constexpr auto alignment = vector_type::coordinator_type::alignment();
+ auto const padding = memory::impl::get_padding<value_type>(alignment, size());
+ return size() + padding;
+ }
+
+ /// subdivide the memory in data_ into the sub-arrays that store the matrix
+ /// diagonals and rhs/solution vectors
+ void setup()
+ {
+ const auto n = size();
+ const auto n_alloc = reservation();
+
+ // point sub-vectors into storage
+ a_ = data_( 0, n);
+ d_ = data_( n_alloc, n_alloc + n);
+ b_ = data_(2*n_alloc, 2*n_alloc + n);
+ rhs_ = data_(3*n_alloc, 3*n_alloc + n);
+ v_ = data_(4*n_alloc, 4*n_alloc + n);
+ }
+
+ // the parent indices are stored as an index array
+ index_type parent_indices_;
+ index_type cell_index_;
+
+ size_type num_cells_;
+
+ // all the data fields point into a single block of storage
+ vector_type data_;
+
+ // the individual data files that point into block storage
+ view_type a_;
+ view_type d_;
+ view_type b_;
+ view_type rhs_;
+ view_type v_;
+};
+
+} // namespace nest
+} // namespace mc
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 28e189a6f5a19df31c9ee961ffe0fe79be1a7676..f6a328956af2469a278b3be357f15bb239b1bd5e 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -6,10 +6,9 @@
#include <swcio.hpp>
-namespace nestmc
-{
-namespace io
-{
+namespace nest {
+namespace mc {
+namespace io {
//
// cell_record implementation
@@ -217,5 +216,6 @@ cell_record_range_clean::cell_record_range_clean(std::istream &is)
}
}
-} // end of nestmc::io
-} // end of nestmc
+} // namespace io
+} // namespace mc
+} // namespace nest
diff --git a/src/swcio.hpp b/src/swcio.hpp
index e267c184b4e462b2015fb8c38866e2133d60b8b1..45e83f741cb30455e9607bf82df5d97ced37f083 100644
--- a/src/swcio.hpp
+++ b/src/swcio.hpp
@@ -6,12 +6,9 @@
#include <string>
#include <vector>
-namespace nestmc
-{
-
-namespace io
-{
-
+namespace nest {
+namespace mc {
+namespace io {
class cell_record
{
@@ -399,5 +396,7 @@ template<typename T = swc_io_clean>
return typename T::cell_range_type(is);
}
-} // end of nestmc::io
-} // end of nestmc
+} // namespace io
+} // namespace mc
+} // namespace nest
+
diff --git a/src/tree.hpp b/src/tree.hpp
index 27822daad87a1ed8480223e784bbe8ee6b0e8cf1..8a296cece8c44f750b2a6798780449231e7cf878 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -296,7 +296,6 @@ class tree {
// check for the senitel that indicates that the old root has
// been processed
if(old_node==-1) {
- //assert(parent_node<0); // sanity check
return new_node;
}
@@ -378,8 +377,17 @@ std::vector<int> make_parent_index(tree const& t, C const& counts)
// taking care for the case where the root node has -1 as its parent
auto parent = t.parent(i);
parent = parent>=0 ? parent : 0;
- parent_index[pos++] = index[parent];
- while(pos<index[i]) {
+
+ // the index of the first compartment in the segment
+ // is calculated differently for the root (i.e when i==parent)
+ if(i!=parent) {
+ parent_index[pos++] = index[parent+1]-1;
+ }
+ else {
+ parent_index[pos++] = parent;
+ }
+ // number the remaining compartments in the segment consecutively
+ while(pos<index[i+1]) {
parent_index[pos] = pos-1;
pos++;
}
diff --git a/src/util.hpp b/src/util.hpp
index 738c3d0f6302cdc467a28a467427c50b8a96a593..cb68c2d1f7aa05c5ad4da7a37a5a1abae11c19dc 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -29,6 +29,17 @@ operator << (std::ostream &o, std::vector<T>const& v)
return o;
}
+template <typename T>
+std::ostream& print(std::ostream &o, std::vector<T>const& v)
+{
+ o << "[";
+ for(auto const& i: v) {
+ o << i << ", ";
+ }
+ o << "]";
+ return o;
+}
+
namespace nest {
namespace mc {
namespace util {
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 9f35f70ca5f9b32c045cb9ff9dd8bd46d4aebb6b..c489429926c39d8ff424f46bb3d6e24dc2a4b591 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -16,14 +16,6 @@ TEST(run, init)
EXPECT_EQ(cell.tree().num_segments(), 2u);
- /*
- for(auto &s : cell.segments()) {
- std::cout << "volume : " << s->volume()
- << " area : " << s->area()
- << " ratio : " << s->volume()/s->area() << std::endl;
- }
- */
-
// in this context (i.e. attached to a segment on a high-level cell)
// a mechanism is essentially a set of parameters
// - the only "state" is that used to define parameters
@@ -43,12 +35,11 @@ TEST(run, init)
EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
- cell.segment(1)->set_compartments(2);
+ cell.segment(1)->set_compartments(200);
using fvm_cell = fvm::fvm_cell<double, int>;
fvm_cell fvcell(cell);
- // print out the parameters if you want...
- //std::cout << soma_hh << "\n";
+
}
// test out the parameter infrastructure
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 7ecd5fc013d98e66428d7a85cb658360e4ca6f2d..96405e4f155ef210b53af6cee757aa7110051fc3 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -10,8 +10,8 @@
#include <swcio.hpp>
// SWC tests
-void expect_cell_equals(const nestmc::io::cell_record &expected,
- const nestmc::io::cell_record &actual)
+void expect_cell_equals(const nest::mc::io::cell_record &expected,
+ const nest::mc::io::cell_record &actual)
{
EXPECT_EQ(expected.id(), actual.id());
EXPECT_EQ(expected.type(), actual.type());
@@ -24,7 +24,7 @@ void expect_cell_equals(const nestmc::io::cell_record &expected,
TEST(cell_record, construction)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// force an invalid type
@@ -91,7 +91,7 @@ TEST(cell_record, construction)
TEST(cell_record, comparison)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// check comparison operators
@@ -107,7 +107,7 @@ TEST(cell_record, comparison)
TEST(swc_parser, invalid_input)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// check incomplete lines; missing parent
@@ -136,7 +136,7 @@ TEST(swc_parser, invalid_input)
TEST(swc_parser, valid_input)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// check empty file; no cell may be parsed
@@ -198,7 +198,7 @@ TEST(swc_parser, valid_input)
TEST(swc_parser, from_allen_db)
{
- using namespace nestmc;
+ using namespace nest::mc;
auto fname = "../data/example.swc";
std::ifstream fid(fname);
@@ -219,7 +219,7 @@ TEST(swc_parser, from_allen_db)
TEST(swc_parser, input_cleaning)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// Check duplicates
@@ -296,7 +296,7 @@ TEST(swc_parser, input_cleaning)
TEST(cell_record_ranges, raw)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// Check valid usage