Adds complex number handling to morph::range and morph::Scale

morph/MathImpl.h

@@ -226,7 +226,7 @@ namespace morph {
      * exposed by morph::MathAlgo (client code should always select functions from
      * MathAlgo).
      *
-     * number_type::value will have been 1 - scalar (see number_type.h).
+     * number_type::value will have been 1 - scalar (see morph::number_type in trait_tests.h).
      */
     template<>
     struct MathImpl<1>
@@ -262,4 +262,36 @@ namespace morph {
 
     };
 
+    /*!
+     * Comlex scalar MathAlgo algorithm specializations
+     *
+     * This is a specialization of MathImpl with vtype set to 2. The templates are applied if the
+     * type T is a complex scalar such as std::complex<float> or std::complex<double>.
+     *
+     * This specialization contains complex scalar implementations of algorithms which are exposed
+     * by morph::MathAlgo (client code should always select functions from MathAlgo).
+     *
+     * number_type::value will have been 2 - complex scalar (see morph::number_type in trait_tests.h).
+     */
+    template<>
+    struct MathImpl<2>
+    {
+        //! Complex scalar maxmin implementation
+        template <typename Container, std::enable_if_t<morph::is_copyable_container<Container>::value, int> = 0>
+        static morph::range<typename Container::value_type> maxmin (const Container& values)
+        {
+            using T = typename Container::value_type;
+            using T_el = typename T::value_type; // If T is std::complex<float>, T_el will be float
+            // Note that there's no specialization of numeric_limits for std::complex, so set it up manually
+            morph::range<T> r ({std::numeric_limits<T_el>::max(), std::numeric_limits<T_el>::max() }, T{0, 0});
+            for (auto v : values) {
+                // comparison operations on complex numbers commonly consider their modulus - how
+                // far the number is from the origin.
+                r.max = std::abs(v) > std::abs(r.max) ? v : r.max;
+                r.min = std::abs(v) < std::abs(r.min) ? v : r.min;
+            }
+            return r;
+        }
+    };
+
 } // namespace morph

morph/Scale.h

@@ -114,14 +114,16 @@ namespace morph {
                << " as: " << this->transform_str()
                << ". ready()=" << (this->ready() ? "true" : "false")
                << ", do_autoscale=" << (this->do_autoscale ? "true" : "false")
-               << ", params=" << this->params_str();
+               << ", params=" << this->params_str() << ", output range: " << this->output_range_str();
             return ss.str();
         }
 
         //! Output the params vvec as a string
         virtual std::string params_str() const = 0;
         //! Describe the transformation in a text string
         virtual std::string transform_str() const = 0;
+        //! Describe the output range in a text string
+        virtual std::string output_range_str() const = 0;
 
         /*!
          * \brief Transform a container of scalars or vectors.
@@ -347,6 +349,13 @@ namespace morph {
             return ss.str();
         }
 
+        virtual std::string output_range_str() const
+        {
+            std::stringstream ss;
+            ss << this->output_range;
+            return ss.str();
+        }
+
         virtual T inverse_one (const S& datum) const
         {
             T rtn = T{};
@@ -509,6 +518,13 @@ namespace morph {
             return ss.str();
         }
 
+        virtual std::string output_range_str() const
+        {
+            std::stringstream ss;
+            ss << this->output_range;
+            return ss.str();
+        }
+
         virtual T inverse_one (const S& datum) const
         {
             T rtn = T{0};
@@ -569,7 +585,7 @@ namespace morph {
         //! Reset the Scaling by emptying params
         void reset() { this->params.clear(); }
 
-    private:
+    protected:
         //! Linear transform for scalar type; y = mx + c
         S transform_one_linear (const T& datum) const
         {
@@ -597,7 +613,7 @@ namespace morph {
             return (std::exp (res));
         }
 
-        void compute_scaling_linear (T input_min, T input_max)
+        virtual void compute_scaling_linear (T input_min, T input_max)
         {
             // Here, we need to use the output type for the computations. Does that mean
             // params is stored in the output type? I think it does.
@@ -609,14 +625,18 @@ namespace morph {
                 // m = rise/run
                 this->params[0] = (this->output_range.max - this->output_range.min) / static_cast<S>(input_max - input_min);
                 // c = y - mx => min = m * input_min + c => c = min - (m * input_min)
+                // FIXME: May need inspiration from the vector implementation of morph::Scale, above.
                 this->params[1] = this->output_range.min - (this->params[0] * static_cast<S>(input_min));
             }
         }
 
-        void compute_scaling_log (T input_min, T input_max)
+        virtual void compute_scaling_log (T input_min, T input_max)
         {
-            if (input_min <= T{0} || input_max <= T{0}) {
-                throw std::runtime_error ("Can't logarithmically autoscale a range which includes zeros or negatives");
+            // Have to check here as ScaleImpl<2, T, S> is built from ScaleImpl<1, T, S> but <= operator makes no sense
+            if constexpr (morph::number_type<T>::value == 1) {
+                if (input_min <= T{0} || input_max <= T{0}) {
+                    throw std::runtime_error ("Can't logarithmically autoscale a range which includes zeros or negatives");
+                }
             }
             T ln_imin = std::log(input_min);
             T ln_imax = std::log(input_max);
@@ -628,6 +648,70 @@ namespace morph {
         morph::vvec<S> params;
     };
 
+    /*!
+     * \brief Experimental ScaleImpl for complex scalars\a T
+     *
+     * A specialized implementation base class for the template class Scale, which is used when the
+     * number type of \a T is std::complex<>
+     *
+     * \tparam ntype The 'number type' as contained in number_type::value. 0 for vectors, 1 for
+     * scalars. This class is active only for ntype==2 (complex scalar).
+     *
+     * \tparam T The type of the number to be scaled. Should be some complex scalar type such as
+     * std::complex<float> or std::complex<double>.
+     *
+     * \tparam S The output type for the scaled number. For integer T, this might well be a floating
+     * point type. Does it make sense to scale from complex to non complex? Maybe, if you're scaling
+     * to magnitude.
+     *
+     * \sa ScaleImplBase
+     */
+    template<typename T, typename S>
+    class ScaleImpl<2, T, S> : public ScaleImpl<1, T, S>
+    {
+    public:
+        // Any public overrides go here
+    protected:
+        // Scaling a set of complex numbers means stretching/squashing the plane out so that the
+        // complex values have magnitudes between output_min and output_max. This is a little
+        // different from scaling of regular scalars.
+        void compute_scaling_linear (T input_min, T input_max)
+        {
+            // We enforce output_range.min == {0, 0} for complex number scaling. We simply shrink or
+            // stretch the complex plane so that the input numbers with the largest magnitude are
+            // given the magniude of output_range.max.
+            if (this->output_range.min != T{0}) {
+                throw std::runtime_error ("ScaleImpl<2, T, S>::compute_scaling_linear: "
+                                          "output_range minimum must be (0 + 0i) for complex scaling");
+            }
+
+            // Here, we need to use the output type for the computations. That means params is
+            // stored in the output type, which is (should be) complex. Only the magnitude of each
+            // param matters.
+            //
+            // Note that we ONLY scale values on the complex plane by multiplication. The origin
+            // remains at the origin. So we leave params[1] equal to S{0}
+            this->params.resize (2, S{0});
+
+            if (input_min != input_max) {
+                // m = rise/run but note that we ignore input_min
+                this->params[0] = (std::abs(this->output_range.max) - std::abs(this->output_range.min)) / static_cast<S>(std::abs(input_max) /*- std::abs(input_min)*/);
+            }
+        }
+
+        void compute_scaling_log (T input_min, T input_max)
+        {
+            // Log scaling complex range?
+            if (std::abs(input_min) == T{0} || std::abs(input_max) == T{0}) {
+                throw std::runtime_error ("Can't logarithmically autoscale a complex range which includes zeros");
+            }
+            T ln_imin = std::log(input_min);
+            T ln_imax = std::log(input_max);
+            // Now just scale linearly between ln_imin and ln_imax
+            this->compute_scaling_linear (ln_imin, ln_imax);
+        }
+    };
+
     /*!
      * A class for scaling and normalizing signals.
      *

morph/range.h

@@ -12,6 +12,7 @@
 
 #include <ostream>
 #include <limits>
+#include <morph/trait_tests.h>
 
 namespace morph {
 
@@ -50,21 +51,52 @@ namespace morph {
         // std::cout << "The range of values in data was: " << r << std::endl;
         constexpr void search_init()
         {
-            this->min = std::numeric_limits<T>::max();
-            this->max = std::numeric_limits<T>::lowest();
+            if constexpr (morph::number_type<T>::value == 2) { // range is complex
+                this->min = { std::numeric_limits<typename T::value_type>::max(), std::numeric_limits<typename T::value_type>::max() };
+                this->max = T{0};
+            } else { // assume numeric_limits will work with T
+                this->min = std::numeric_limits<T>::max();
+                this->max = std::numeric_limits<T>::lowest();
+            }
         }
 
         // Extend the range to include the given datum. Return true if the range changed.
         constexpr bool update (const T& d)
         {
             bool changed = false;
-            this->min = d < this->min ? changed = true, d : this->min;
-            this->max = d > this->max ? changed = true, d : this->max;
+            if constexpr (morph::number_type<T>::value == 2) { // range is complex
+                this->min = std::abs(d) < std::abs(this->min) ? changed = true, d : this->min;
+                this->max = std::abs(d) > std::abs(this->max) ? changed = true, d : this->max;
+
+            } else if constexpr (morph::number_type<T>::value == 0) { // range is vector
+#if __cplusplus >= 202002L
+                []<bool flag = false>() { static_assert(flag, "Vector ranges are not yet supported"); }();
+#else
+                throw std::runtime_error ("Vector ranges are not yet supported");
+#endif
+            } else {
+                this->min = d < this->min ? changed = true, d : this->min;
+                this->max = d > this->max ? changed = true, d : this->max;
+            }
             return changed;
         }
 
         // Does the range include v?
-        constexpr bool includes (const T& v) { return (v <= this->max && v >= this->min); }
+        constexpr bool includes (const T& v)
+        {
+            if constexpr (morph::number_type<T>::value == 2) { // range is complex
+                return (std::abs(v) <= std::abs(this->max) && std::abs(v) >= std::abs(this->min));
+
+            } else if constexpr (morph::number_type<T>::value == 0) { // range type is vector
+#if __cplusplus >= 202002L
+                []<bool flag = false>() { static_assert(flag, "Vector ranges are not yet supported"); }();
+#else
+                throw std::runtime_error ("Vector ranges are not yet supported");
+#endif
+            } else {
+                return (v <= this->max && v >= this->min);
+            }
+        }
 
         // What's the 'span of the range'?
         constexpr T span() const { return this->max - this->min; }

morph/trait_tests.h

@@ -142,17 +142,20 @@ namespace morph {
      * From the typename T, set a #value attribute which says whether T is a scalar (like
      * float, double), or vector (basically, anything else).
      *
-     * Query the attribute `value`, which will be 1 for scalar and 0 for anything else including
-     * vectors.  This really just wraps std::is_scalar.
+     * Query the attribute `value`, which will be 1 for scalars, 2 for complex scalars and 0
+     * for anything else, which includes vectors, arrays. As long as you know that your 'anything
+     * else' is some sort of vector type, you can use this in template classes like morph::Scale for
+     * scalar/vector implementations.
      *
      * \tparam T the type to distinguish
      */
     template <typename T>
     struct number_type {
         //! is_scalar test
         static constexpr bool const scalar = std::is_scalar<std::decay_t<T>>::value;
-        //! Set value simply from the is_scalar test. 0 for vector, 1 for scalar
-        static constexpr int const value = scalar ? 1 : 0;
+        static constexpr bool const cplx = morph::is_complex<std::decay_t<T>>::value;
+        //! Set value simply from the is_scalar test. 0 for vector, 1 for scalar, 2 for complex scalar
+        static constexpr int const value = scalar ? 1 : cplx ? 2 : 0;
     };
 
 } // morph::

tests/CMakeLists.txt

@@ -350,6 +350,10 @@ add_test(testvec11 testvec11)
 add_executable(testScaleVector testScaleVector.cpp)
 add_test(testScaleVector testScaleVector)
 
+# Test scaling complex numbers
+add_executable(testScale_complex testScale_complex.cpp)
+add_test(testScale_complex testScale_complex)
+
 # Test morph::TransformMatrix (4x4 matrix)
 add_executable(testTransformMatrix testTransformMatrix.cpp)
 add_test(testTransformMatrix testTransformMatrix)
@@ -394,6 +398,9 @@ add_test(testrange testrange)
 add_executable(testrange_constexpr testrange_constexpr.cpp)
 add_test(testrange_constexpr testrange_constexpr)
 
+add_executable(testrange_complex testrange_complex.cpp)
+add_test(testrange_complex testrange_complex)
+
 # Test the colour mapping
 add_executable(testColourMap testColourMap.cpp)
 add_test(testColourMap testColourMap)
@@ -550,3 +557,6 @@ add_test(testloadpng testloadpng)
 
 add_executable(test_histo test_histo.cpp)
 add_test(test_histo test_histo)
+
+add_executable(test_number_type test_number_type.cpp)
+add_test(test_number_type test_number_type)