dir* command modifications for fixed direction handling

cmd/dirflip.cpp

@@ -49,6 +49,9 @@ void usage ()
     + Option ("permutations", "number of permutations to try (default: " + str(default_permutations) + ")")
     +   Argument ("num").type_integer (1)
 
+    + Option ("preserve", "preserve the sign of some number of directions at the start of the set")
+    +   Argument ("num").type_integer(1)
+
     + Option ("cartesian", "Output the directions in Cartesian coordinates [x y z] instead of [az el].");
 }
 
@@ -57,10 +60,10 @@ using value_type = double;
 using vector3_type = Eigen::Vector3d;
 
 
-class Shared { 
+class Shared {
   public:
-    Shared (const Eigen::MatrixXd& directions, size_t target_num_permutations) :
-      directions (directions), target_num_permutations (target_num_permutations), num_permutations(0),
+    Shared (const Eigen::MatrixXd& directions, const size_t target_num_permutations, const size_t preserve) :
+      directions (directions), target_num_permutations (target_num_permutations), preserve (preserve), num_permutations(0),
       progress ("optimising directions for eddy-currents", target_num_permutations),
       best_signs (directions.rows(), 1), best_eddy (std::numeric_limits<value_type>::max()) { }
 
@@ -90,11 +93,13 @@ class Shared {
 
     vector<int> get_init_signs () const { return vector<int> (directions.rows(), 1); }
     const vector<int>& get_best_signs () const { return best_signs; }
+    size_t get_preserve() const { return preserve; }
 
 
   protected:
     const Eigen::MatrixXd& directions;
     const size_t target_num_permutations;
+    const size_t preserve;
     size_t num_permutations;
     ProgressBar progress;
     vector<int> best_signs;
@@ -107,12 +112,12 @@ class Shared {
 
 
 
-class Processor { 
+class Processor {
   public:
     Processor (Shared& shared) :
       shared (shared),
       signs (shared.get_init_signs()),
-      uniform (0, signs.size()-1) { }
+      uniform (shared.get_preserve(), signs.size() - shared.get_preserve() - 1) { }
 
     void execute () {
       while (eval());
@@ -155,10 +160,11 @@ void run ()
   auto directions = DWI::Directions::load_cartesian (argument[0]);
 
   size_t num_permutations = get_option_value<size_t> ("permutations", default_permutations);
+  size_t preserve = get_option_value<size_t> ("preserve", 0);
 
   vector<int> signs;
   {
-    Shared eddy_shared (directions, num_permutations);
+    Shared eddy_shared (directions, num_permutations, preserve);
     Thread::run (Thread::multi (Processor (eddy_shared)), "eval thread");
     signs = eddy_shared.get_best_signs();
   }

cmd/dirgen.cpp

@@ -68,6 +68,9 @@ void usage ()
     + Option ("restarts", "specify the number of restarts to perform (default: " + str(DEFAULT_RESTARTS) + ").")
     +   Argument ("num").type_integer (1, std::numeric_limits<int>::max())
 
+    + Option ("fixed", "specify a fixed direction (comm-separateed floats) that will always be included at the start of the scheme").allow_multiple()
+    +   Argument ("direction").type_sequence_float()
+
     + Option ("unipolar", "optimise assuming a unipolar electrostatic repulsion model rather than the bipolar model normally assumed in DWI")
 
     + Option ("cartesian", "Output the directions in Cartesian coordinates [x y z] instead of [az el].");
@@ -78,7 +81,7 @@ void usage ()
 
 
 // constrain directions to remain unit length:
-class ProjectedUpdate { 
+class ProjectedUpdate {
   public:
     bool operator() (
         Eigen::VectorXd& newx,
@@ -98,11 +101,11 @@ class ProjectedUpdate {
 
 
 
-class Energy { 
+class Energy {
   public:
-    Energy (ProgressBar& progress) :
+    Energy (ProgressBar& progress, const size_t ndirs) :
       progress (progress),
-      ndirs (to<int> (argument[0])),
+      ndirs (ndirs),
       bipolar (!(get_options ("unipolar").size())),
       power (0),
       directions (3 * ndirs) { }
@@ -125,7 +128,9 @@ FORCE_INLINE
     double init (Eigen::VectorXd& x)
     {
       Math::RNG::Normal<double> rng;
-      for (size_t n = 0; n < ndirs; ++n) {
+      for (size_t n = 0; n != fixed_directions.size(); ++n)
+        x.segment (3*n, 3) = fixed_directions[n];
+      for (size_t n = fixed_directions.size(); n < ndirs; ++n) {
         auto d = x.segment (3*n,3);
         d[0] = rng();
         d[1] = rng();
@@ -170,8 +175,10 @@ FORCE_INLINE
         }
       }
 
+      // don't move those directions that are to remain fixed
+      g.segment(0,3*fixed_directions.size()).setZero();
       // constrain gradients to lie tangent to unit sphere:
-      for (size_t n = 0; n < ndirs; ++n)
+      for (size_t n = fixed_directions.size(); n < ndirs; ++n)
         g.segment(3*n,3) -= x.segment(3*n,3).dot (g.segment(3*n,3)) * x.segment(3*n,3);
 
       return E;
@@ -209,8 +216,6 @@ FORCE_INLINE
           E = optim.value();
         }
 
-
-
         std::lock_guard<std::mutex> lock (mutex);
         if (E < best_E) {
           best_E = E;
@@ -225,6 +230,7 @@ FORCE_INLINE
     static size_t niter;
     static double best_E;
     static Eigen::VectorXd best_directions;
+    static vector<Eigen::Vector3d> fixed_directions;
 
   protected:
     ProgressBar& progress;
@@ -246,6 +252,7 @@ std::mutex Energy::mutex;
 std::atomic<size_t> Energy::current_start (0);
 double Energy::best_E = std::numeric_limits<double>::infinity();
 Eigen::VectorXd Energy::best_directions;
+vector<Eigen::Vector3d> Energy::fixed_directions;
 
 
 
@@ -255,15 +262,26 @@ void run ()
   Energy::restarts = get_option_value ("restarts", DEFAULT_RESTARTS);
   Energy::target_power = get_option_value ("power", DEFAULT_POWER);
   Energy::niter = get_option_value ("niter", DEFAULT_NITER);
+  auto opt = get_options ("fixed");
+  for (size_t i = 0; i != opt.size(); ++i) {
+    auto value = parse_floats(opt[i][0]);
+    switch (value.size()) {
+      case 2: { Eigen::Vector3d xyz; Math::Sphere::spherical2cartesian(Eigen::Map<Eigen::Vector2d> (value.data()), xyz); Energy::fixed_directions.push_back (xyz); } break;
+      case 3: Energy::fixed_directions.push_back (Eigen::Map<Eigen::Vector3d>(value.data()).normalized()); break;
+      default: throw Exception ("Fixed directions must be either spherical or cartesian directions (comma-separated 2- or 3-vectors)");
+    }
+  }
+  const size_t ndirs = argument[0];
+  if (Energy::fixed_directions.size() >= ndirs)
+    throw Exception ("No directions left to optimise after fixed directions specified");
 
   {
     ProgressBar progress ("Optimising directions up to power " + str(Energy::target_power) + " (" + str(Energy::restarts) + " restarts)");
-    Energy energy_functor (progress);
+    Energy energy_functor (progress, ndirs);
     auto threads = Thread::run (Thread::multi (energy_functor), "energy function");
   }
 
   CONSOLE ("final energy = " + str(Energy::best_E));
-  size_t ndirs = Energy::best_directions.size()/3;
   Eigen::MatrixXd directions_matrix (ndirs, 3);
   for (size_t n = 0; n < ndirs; ++n)
     directions_matrix.row (n) = Energy::best_directions.segment (3*n, 3);

cmd/dirorder.cpp

@@ -43,6 +43,8 @@ void usage ()
     + Argument ("output", "the output directions file").type_file_out();
 
   OPTIONS
+    + Option ("skip", "preserve some number of directions in their position at the start of the set")
+      + Argument ("num").type_integer(1)
     + Option ("cartesian", "Output the directions in Cartesian coordinates [x y z] instead of [az el].");
 }
 
@@ -52,12 +54,18 @@ using value_type = double;
 
 
 
-vector<size_t> optimise (const Eigen::MatrixXd& directions, const size_t first_volume)
+vector<size_t> optimise (const Eigen::MatrixXd& directions, const size_t skip, const size_t first_volume)
 {
-  vector<size_t> indices (1, first_volume);
+  assert (first_volume >= skip);
+  vector<size_t> indices;
+  indices.reserve (directions.rows());
+  for (size_t n = 0; n != skip; ++n)
+    indices.push_back(n);
+  indices.push_back (first_volume);
   vector<size_t> remaining;
-  for (size_t n = 0; n < size_t(directions.rows()); ++n)
-    if (n != indices[0])
+  remaining.reserve (directions.rows() - skip);
+  for (size_t n = skip; n < size_t(directions.rows()); ++n)
+    if (n != first_volume)
       remaining.push_back (n);
 
   while (remaining.size()) {
@@ -117,20 +125,24 @@ void run ()
 {
   auto directions = DWI::Directions::load_cartesian (argument[0]);
 
+  const size_t skip = get_option_value<size_t> ("skip", 0);
+
   size_t last_candidate_first_volume = directions.rows();
   if (size_t(directions.rows()) <= Math::SH::NforL (2)) {
     WARN ("Very few directions in input ("
           + str(directions.rows())
           + "); selection of first direction cannot be optimised"
-          + " (first direction in input will be first direction in output)");
-    last_candidate_first_volume = 1;
+          + (skip ? " (direction #" + str(skip+1) + " will be first direction in output "
+                    "as that is the first direction after those to be skipped)"
+                  : " (first direction in input will be first direction in output)"));
+    last_candidate_first_volume = skip + 1;
   }
 
   value_type min_cost = std::numeric_limits<value_type>::infinity();
   vector<size_t> best_order;
   ProgressBar progress ("Determining best reordering", directions.rows());
-  for (size_t first_volume = 0; first_volume != last_candidate_first_volume; ++first_volume) {
-    const vector<size_t> order = optimise (directions, first_volume);
+  for (size_t first_volume = skip; first_volume != last_candidate_first_volume; ++first_volume) {
+    const vector<size_t> order = optimise (directions, skip, first_volume);
     const value_type cost = calc_cost (directions, order);
     if (cost < min_cost) {
       min_cost = cost;

docs/reference/commands/dirflip.rst

@@ -28,6 +28,8 @@ Options
 
 -  **-permutations num** number of permutations to try (default: 100000000)
 
+-  **-preserve num** preserve the sign of some number of directions at the start of the set
+
 -  **-cartesian** Output the directions in Cartesian coordinates [x y z] instead of [az el].
 
 Standard options

docs/reference/commands/dirgen.rst

@@ -32,6 +32,8 @@ Options
 
 -  **-restarts num** specify the number of restarts to perform (default: 10).
 
+-  **-fixed direction** *(multiple uses permitted)* specify a fixed direction (comm-separateed floats) that will always be included at the start of the scheme
+
 -  **-unipolar** optimise assuming a unipolar electrostatic repulsion model rather than the bipolar model normally assumed in DWI
 
 -  **-cartesian** Output the directions in Cartesian coordinates [x y z] instead of [az el].

docs/reference/commands/dirorder.rst

@@ -26,6 +26,8 @@ The intent of this command is to reorder a set of gradient directions such that
 Options
 -------
 
+-  **-skip num** preserve some number of directions in their position at the start of the set
+
 -  **-cartesian** Output the directions in Cartesian coordinates [x y z] instead of [az el].
 
 Standard options