Merge branch 'master' into master

Author: dcwhite

src/Core/Algorithms/BrainStimulator/Tests/GenerateROIStatisticsAlgorithmTests.cc

@@ -75,7 +75,7 @@ namespace
    (*expected_result)(0, 1) = 0.2023192248404992;  (*expected_result)(1, 1) = std::numeric_limits<double>::quiet_NaN();  (*expected_result)(2, 1) = 0.270036948169627;  (*expected_result)(3, 1) = 0.2819928747450936;  (*expected_result)(4, 1) = 0.2464919130115238;  (*expected_result)(5, 1) = 0.348358459345021;  (*expected_result)(6, 1) = 0.3382031144990396;  (*expected_result)(7, 1) = 0.1699586184078796;
    (*expected_result)(0, 2) = 0.1868726045543786;  (*expected_result)(1, 2) = 0.7546866819823609; (*expected_result)(2, 2) = 0.1189976815583766; (*expected_result)(3, 2) = 0.223811939491137;  (*expected_result)(4, 2) = 0.2550951154592691;     (*expected_result)(5, 2) = 0.1386244428286791;  (*expected_result)(6, 2) = 0.2435249687249893;  (*expected_result)(7, 2) = 0.1965952504312082;
    (*expected_result)(0, 3) = 0.7951999011370632;  (*expected_result)(1, 3) = 0.7546866819823609;  (*expected_result)(2, 3) = 0.6797026768536748;  (*expected_result)(3, 3) = 0.9597439585160811;  (*expected_result)(4, 3) = 0.8909032525357985;  (*expected_result)(5, 3) = 0.9592914252054443;  (*expected_result)(6, 3) = 0.9292636231872278;  (*expected_result)(7, 3) = 0.6160446761466392;
-   (*expected_result)(0, 4) = 9;  (*expected_result)(1, 4) = 1;  (*expected_result)(2, 4) = 5;    (*expected_result)(3, 4) = 5;  (*expected_result)(4, 4) = 5;  (*expected_result)(5, 4) = 5;  (*expected_result)(6, 4) = 5;    (*expected_result)(7, 4) = 5; 
+   (*expected_result)(0, 4) = 9;  (*expected_result)(1, 4) = 1;  (*expected_result)(2, 4) = 5;    (*expected_result)(3, 4) = 5;  (*expected_result)(4, 4) = 5;  (*expected_result)(5, 4) = 5;  (*expected_result)(6, 4) = 5;    (*expected_result)(7, 4) = 5;
 
     return expected_result;
   }
@@ -153,7 +153,7 @@ TEST(GenerateROIStatisticsAlgorithm, ValidInput_TestCoordinateSpace)
   FieldHandle atlas = CreateAtlas();
   FieldHandle createcoordcubesshiftby10 = CreateCoordinatesCubesShiftedBy10();
   FieldHandle createmeshsolution = CreateGetMeshSolution();
-  const std::string& atlasMeshLabelsStr = std::string("Region1;Region2;Region3;Region4;Region5;Region6;Region7;Region8");
+  std::string atlasMeshLabelsStr("Region1;Region2;Region3;Region4;Region5;Region6;Region7;Region8");
   DenseMatrixHandle specROI(new DenseMatrix(5, 1));
   (*specROI)(0, 0) = 10.5;
   (*specROI)(1, 0) = 10.5;

src/Core/Algorithms/Math/SolveLinearSystemWithEigen.cc

@@ -34,45 +34,46 @@
 #include <Core/Datatypes/MatrixTypeConversions.h>
 #include <Eigen/Sparse>
 
+using namespace SCIRun;
 using namespace SCIRun::Core::Algorithms::Math;
 using namespace SCIRun::Core::Datatypes;
 using namespace SCIRun::Core::Algorithms;
 using namespace SCIRun::Core;
 
 namespace
 {
-  template <class ColumnMatrixType>
+  template <class ColumnMatrixType, template <typename> class SolverType>
   class SolveLinearSystemAlgorithmEigenCGImpl
   {
   public:
-    SolveLinearSystemAlgorithmEigenCGImpl(const ColumnMatrixType& rhs, double tolerance, int maxIterations) :
+    SolveLinearSystemAlgorithmEigenCGImpl(SharedPointer<ColumnMatrixType> rhs, double tolerance, int maxIterations) :
         tolerance_(tolerance), maxIterations_(maxIterations), rhs_(rhs) {}
 
     using SolutionType = ColumnMatrixType;
 
     template <class MatrixType>
     typename ColumnMatrixType::EigenBase solveWithEigen(const MatrixType& lhs)
     {
-      Eigen::ConjugateGradient<typename MatrixType::EigenBase> cg;
-      cg.compute(lhs);
+      SolverType<typename MatrixType::EigenBase> solver;
+      solver.compute(lhs);
 
-      if (cg.info() != Eigen::Success)
+      if (solver.info() != Eigen::Success)
         BOOST_THROW_EXCEPTION(AlgorithmInputException()
-          << LinearAlgebraErrorMessage("Conjugate gradient initialization was unsuccessful")
-          << EigenComputationInfo(cg.info()));
-
-      cg.setTolerance(tolerance_);
-      cg.setMaxIterations(maxIterations_);
-      auto solution = cg.solve(rhs_).eval();
-      tolerance_ = cg.error();
-      maxIterations_ = cg.iterations();
+          << LinearAlgebraErrorMessage("Eigen solver initialization was unsuccessful")
+          << EigenComputationInfo(solver.info()));
+
+      solver.setTolerance(tolerance_);
+      solver.setMaxIterations(maxIterations_);
+      auto solution = solver.solve(*rhs_).eval();
+      tolerance_ = solver.error();
+      maxIterations_ = solver.iterations();
       return solution;
     }
 
     double tolerance_;
     int maxIterations_;
   private:
-    const ColumnMatrixType& rhs_;
+    SharedPointer<ColumnMatrixType> rhs_;
   };
 }
 
@@ -86,6 +87,14 @@ SolveLinearSystemAlgorithm::ComplexOutputs SolveLinearSystemAlgorithm::run(const
   return runImpl<ComplexInputs, ComplexOutputs>(input, params);
 }
 
+template <typename T>
+using CG = Eigen::ConjugateGradient<T>;
+// Not available yet, need to upgrade Eigen
+// template <typename T>
+// using LSCG = Eigen::LeastSquaresConjugateGradient<T>;
+template <typename T>
+using BiCG = Eigen::BiCGSTAB<T>;
+
 template <typename In, typename Out>
 Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& params) const
 {
@@ -101,9 +110,32 @@ Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& param
   int maxIterations = std::get<1>(params);
   ENSURE_POSITIVE_INT(maxIterations, "Max iterations out of range!");
 
-  using SolutionType = DenseMatrixGeneric<typename std::tuple_element<0, In>::type::element_type::value_type>;
-  using SolverType = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType>;
-  SolverType impl(*b, tolerance, maxIterations);
+  auto method = std::get<2>(params);
+
+  using SolutionType = DenseColumnMatrixGeneric<typename std::tuple_element<0, In>::type::element_type::value_type>;
+  using AlgoTypeCG = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType, CG>;
+  using AlgoTypeBiCG = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType, BiCG>;
+
+  if ("cg" == method)
+    return solve<AlgoTypeCG, In, Out>(input, params);
+  else if ("bicg" == method)
+    return solve<AlgoTypeBiCG, In, Out>(input, params);
+  else
+  {
+    BOOST_THROW_EXCEPTION(AlgorithmProcessingException() << ErrorMessage("Need to upgrade Eigen for LSCG."));
+  }
+}
+
+template <typename SolverType, typename In, typename Out>
+Out SolveLinearSystemAlgorithm::solve(const In& input, const Parameters& params) const
+{
+  auto A = std::get<0>(input);
+  auto b = std::get<1>(input);
+  double tolerance = std::get<0>(params);
+  int maxIterations = std::get<1>(params);
+
+  SolverType impl(b, tolerance, maxIterations);
+
   typename SolverType::SolutionType x;
   if (matrixIs::dense(A))
   {
@@ -120,8 +152,7 @@ Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& param
 
   if (x.size() != 0)
   {
-    /// @todo: move ctor
-    auto solution(boost::make_shared<SolutionType>(x));
+    auto solution(boost::make_shared<typename SolverType::SolutionType>(x));
     return Out(solution, impl.tolerance_, impl.maxIterations_);
   }
   else

src/Core/Algorithms/Math/SolveLinearSystemWithEigen.h

@@ -45,9 +45,9 @@ namespace Math {
   public:
     typedef std::tuple<SCIRun::Core::Datatypes::MatrixHandle, SCIRun::Core::Datatypes::DenseColumnMatrixHandle> Inputs;
     typedef std::tuple<SCIRun::Core::Datatypes::ComplexMatrixHandle, SCIRun::Core::Datatypes::ComplexDenseColumnMatrixHandle> ComplexInputs;
-    typedef std::tuple<double, int> Parameters;
-    typedef std::tuple<SCIRun::Core::Datatypes::DenseMatrixHandle, double, int> Outputs;
-    typedef std::tuple<SCIRun::Core::Datatypes::ComplexDenseMatrixHandle, double, int> ComplexOutputs;
+    typedef std::tuple<double, int, std::string> Parameters;
+    typedef std::tuple<SCIRun::Core::Datatypes::DenseColumnMatrixHandle, double, int> Outputs;
+    typedef std::tuple<SCIRun::Core::Datatypes::ComplexDenseColumnMatrixHandle, double, int> ComplexOutputs;
 
     Outputs run(const Inputs& input, const Parameters& params) const;
     ComplexOutputs run(const ComplexInputs& input, const Parameters& params) const;
@@ -56,6 +56,8 @@ namespace Math {
   private:
     template <typename In, typename Out>
     Out runImpl(const In& input, const Parameters& params) const;
+    template <typename SolverType, typename In, typename Out>
+    Out solve(const In& input, const Parameters& params) const;
   };
 
   typedef boost::error_info<struct tag_eigen_computation, Eigen::ComputationInfo> EigenComputationInfo;

src/Core/Algorithms/Math/Tests/SolveLinearSystemWithEigenTests.cc

@@ -344,7 +344,7 @@ TEST(EigenSparseSolverTest, DISABLED_CanSolveBigSystem)
     ScopedTimer t("using algorithm object");
     SolveLinearSystemAlgorithm algo;
 
-    x = algo.run(std::make_tuple(A, bCol), std::make_tuple(1e-20, 4000));
+    x = algo.run(std::make_tuple(A, bCol), std::make_tuple(1e-20, 4000, "cg"));
     MatrixHandle solution = std::get<0>(x);
 
     ASSERT_TRUE(solution.get() != nullptr);

src/Dataflow/Engine/Controller/ProvenanceManager.h

@@ -44,11 +44,11 @@ namespace Engine {
   class ProvenanceManager : boost::noncopyable
   {
   public:
-    typedef ProvenanceItem<Memento> Item;
-    typedef typename Item::Handle ItemHandle;
-    typedef std::stack<ItemHandle> Stack;
-    typedef typename Stack::container_type List;
-    typedef Engine::NetworkIOInterface<Memento> IOType;
+    using Item = ProvenanceItem<Memento>;
+    using ItemHandle = typename Item::Handle;
+    using Stack = std::stack<ItemHandle>;
+    using List = typename Stack::container_type;
+    using IOType = Engine::NetworkIOInterface<Memento>;
 
     explicit ProvenanceManager(IOType* networkIO);
     void setInitialState(const Memento& initialState);

src/Dataflow/Engine/Scheduler/GraphNetworkAnalyzer.cc

@@ -154,6 +154,8 @@ bool ExecuteSingleModule::operator()(SCIRun::Dataflow::Networks::ModuleHandle mo
     THROW_INVALID_ARGUMENT("Module not found in component map");
 
   auto rootId = module_->get_id();
+  if (rootId.name_ == "Subnet")
+    return false;
   auto rootIdIter = components_.find(rootId);
   if (rootIdIter == components_.end())
     THROW_INVALID_ARGUMENT("Current module not found in component map");