Add pointwise_divide function (#1586)

* Add pointwise_divide operator

* Add tests for pointwise_divide function

* Replace in affine layer pointwise_multiply for division by its equivalent pointwise_divide.

dlib/dnn/layers.h

@@ -2039,7 +2039,7 @@ namespace dlib
             auto sg = gamma(temp,0);
             auto sb = beta(temp,gamma.size());
 
-            g = pointwise_multiply(mat(sg), 1.0f/sqrt(mat(item.running_variances)+item.get_eps()));
+            g = pointwise_divide(mat(sg), sqrt(mat(item.running_variances)+item.get_eps()));
             b = mat(sb) - pointwise_multiply(mat(g), mat(item.running_means));
         }
 

dlib/matrix/matrix_utilities.h

@@ -2772,6 +2772,159 @@ namespace dlib
 
 // ----------------------------------------------------------------------------------------
 
+    template <typename M1, typename M2>
+    struct op_pointwise_divide : basic_op_mm<M1,M2>
+    {
+        op_pointwise_divide( const M1& m1_, const M2& m2_) : basic_op_mm<M1,M2>(m1_,m2_){}
+
+        typedef typename impl::compatible<typename M1::type, typename M2::type>::type type;
+        typedef const type const_ret_type;
+        const static long cost = M1::cost + M2::cost + 1;
+
+        const_ret_type apply ( long r, long c) const
+        { return this->m1(r,c)/this->m2(r,c); }
+    };
+
+    template <
+        typename EXP1,
+        typename EXP2
+        >
+    inline const matrix_op<op_pointwise_divide<EXP1,EXP2> > pointwise_divide (
+        const matrix_exp<EXP1>& a,
+        const matrix_exp<EXP2>& b
+    )
+    {
+        COMPILE_TIME_ASSERT((impl::compatible<typename EXP1::type,typename EXP2::type>::value == true));
+        COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0);
+        COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0);
+        DLIB_ASSERT(a.nr() == b.nr() &&
+               a.nc() == b.nc(),
+            "\tconst matrix_exp pointwise_divide(const matrix_exp& a, const matrix_exp& b)"
+            << "\n\tYou can only make a do a pointwise divide with two equally sized matrices"
+            << "\n\ta.nr(): " << a.nr()
+            << "\n\ta.nc(): " << a.nc()
+            << "\n\tb.nr(): " << b.nr()
+            << "\n\tb.nc(): " << b.nc()
+            );
+        typedef op_pointwise_divide<EXP1,EXP2> op;
+        return matrix_op<op>(op(a.ref(),b.ref()));
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename M1, typename M2, typename M3>
+    struct op_pointwise_divide3 : basic_op_mmm<M1,M2,M3>
+    {
+        op_pointwise_divide3( const M1& m1_, const M2& m2_, const M3& m3_) :
+            basic_op_mmm<M1,M2,M3>(m1_,m2_,m3_){}
+
+        typedef typename M1::type type;
+        typedef const typename M1::type const_ret_type;
+        const static long cost = M1::cost + M2::cost + M3::cost + 2;
+
+        const_ret_type apply (long r, long c) const
+        { return this->m1(r,c)/this->m2(r,c)/this->m3(r,c); }
+    };
+
+    template <
+        typename EXP1,
+        typename EXP2,
+        typename EXP3
+        >
+    inline const matrix_op<op_pointwise_divide3<EXP1,EXP2,EXP3> >
+        pointwise_divide (
+        const matrix_exp<EXP1>& a,
+        const matrix_exp<EXP2>& b,
+        const matrix_exp<EXP3>& c
+    )
+    {
+        COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true));
+        COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true));
+        COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0);
+        COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NR == 0 || EXP2::NC == 0);
+        COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0);
+        COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0);
+        DLIB_ASSERT(a.nr() == b.nr() &&
+               a.nc() == b.nc() &&
+               b.nr() == c.nr() &&
+               b.nc() == c.nc(),
+            "\tconst matrix_exp pointwise_divide(a,b,c)"
+            << "\n\tYou can only make a do a pointwise divide between equally sized matrices"
+            << "\n\ta.nr(): " << a.nr()
+            << "\n\ta.nc(): " << a.nc()
+            << "\n\tb.nr(): " << b.nr()
+            << "\n\tb.nc(): " << b.nc()
+            << "\n\tc.nr(): " << c.nr()
+            << "\n\tc.nc(): " << c.nc()
+            );
+
+        typedef op_pointwise_divide3<EXP1,EXP2,EXP3> op;
+        return matrix_op<op>(op(a.ref(),b.ref(),c.ref()));
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename M1, typename M2, typename M3, typename M4>
+    struct op_pointwise_divide4 : basic_op_mmmm<M1,M2,M3,M4>
+    {
+        op_pointwise_divide4( const M1& m1_, const M2& m2_, const M3& m3_, const M4& m4_) :
+            basic_op_mmmm<M1,M2,M3,M4>(m1_,m2_,m3_,m4_){}
+
+        typedef typename M1::type type;
+        typedef const typename M1::type const_ret_type;
+        const static long cost = M1::cost + M2::cost + M3::cost + M4::cost + 3;
+
+        const_ret_type apply (long r, long c) const
+        { return this->m1(r,c)/this->m2(r,c)/this->m3(r,c)/this->m4(r,c); }
+    };
+
+
+    template <
+        typename EXP1,
+        typename EXP2,
+        typename EXP3,
+        typename EXP4
+        >
+    inline const matrix_op<op_pointwise_divide4<EXP1,EXP2,EXP3,EXP4> > pointwise_divide (
+        const matrix_exp<EXP1>& a,
+        const matrix_exp<EXP2>& b,
+        const matrix_exp<EXP3>& c,
+        const matrix_exp<EXP4>& d
+    )
+    {
+        COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true));
+        COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true));
+        COMPILE_TIME_ASSERT((is_same_type<typename EXP3::type,typename EXP4::type>::value == true));
+        COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0);
+        COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0);
+        COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0);
+        COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0);
+        COMPILE_TIME_ASSERT(EXP3::NR == EXP4::NR || EXP3::NR == 0 || EXP4::NR == 0);
+        COMPILE_TIME_ASSERT(EXP3::NC == EXP4::NC || EXP3::NC == 0 || EXP4::NC == 0);
+        DLIB_ASSERT(a.nr() == b.nr() &&
+               a.nc() == b.nc() &&
+               b.nr() == c.nr() &&
+               b.nc() == c.nc() &&
+               c.nr() == d.nr() &&
+               c.nc() == d.nc(),
+            "\tconst matrix_exp pointwise_divide(a,b,c,d)"
+            << "\n\tYou can only make a do a pointwise divide between equally sized matrices"
+            << "\n\ta.nr(): " << a.nr()
+            << "\n\ta.nc(): " << a.nc()
+            << "\n\tb.nr(): " << b.nr()
+            << "\n\tb.nc(): " << b.nc()
+            << "\n\tc.nr(): " << c.nr()
+            << "\n\tc.nc(): " << c.nc()
+            << "\n\td.nr(): " << d.nr()
+            << "\n\td.nc(): " << d.nc()
+            );
+
+        typedef op_pointwise_divide4<EXP1,EXP2,EXP3,EXP4> op;
+        return matrix_op<op>(op(a.ref(),b.ref(),c.ref(),d.ref()));
+	}
+
+    // ----------------------------------------------------------------------------------------
+
     template <
         typename P,
         int type = static_switch<

dlib/matrix/matrix_utilities_abstract.h

@@ -800,6 +800,46 @@ namespace dlib
         performs pointwise_multiply(pointwise_multiply(a,b),pointwise_multiply(c,d));
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp pointwise_divide(
+        const matrix_exp& a,
+        const matrix_exp& b
+    );
+    /*!
+        requires
+            - a.nr() == b.nr()
+            - a.nc() == b.nc()
+            - a and b both contain the same type of element (one or both
+              can also be of type std::complex so long as the underlying type
+              in them is the same)
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R has the same dimensions as a and b.
+                - for all valid r and c:
+                  R(r,c) == a(r,c) / b(r,c)
+    !*/
+
+    const matrix_exp pointwise_divide(
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp& c
+    );
+    /*!
+        performs pointwise_divide(pointwise_divide(a,b),c);
+    !*/
+
+    const matrix_exp pointwise_divide(
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp& c,
+        const matrix_exp& d
+    );
+    /*!
+        performs pointwise_divide(pointwise_divide(pointwise_divide(pointwise_divide(a,b),c),d));
+    !*/
+
 // ----------------------------------------------------------------------------------------
 
     const matrix_exp join_rows (

dlib/test/matrix.cpp

@@ -1164,6 +1164,10 @@ namespace
             DLIB_TEST((complex_matrix(ones_matrix<double>(3,3), zeros_matrix<double>(3,3)) == complex_matrix(ones_matrix<double>(3,3))));
             DLIB_TEST((pointwise_multiply(complex_matrix(ones_matrix<double>(3,3)), ones_matrix<double>(3,3)*2) ==
                        complex_matrix(2*ones_matrix<double>(3,3))));
+            DLIB_TEST((pointwise_divide(complex_matrix(ones_matrix<double>(3,3)), ones_matrix<double>(3,3)) ==
+                       complex_matrix(ones_matrix<double>(3,3))));
+            DLIB_TEST((pointwise_divide(complex_matrix(zeros_matrix<double>(3,3)), ones_matrix<double>(3,3)) ==
+                       complex_matrix(zeros_matrix<double>(3,3))));
         }
 
         {

dlib/test/matrix2.cpp

@@ -134,6 +134,10 @@ namespace
 
         DLIB_TEST(squared(m4) == pointwise_multiply(m4,m4));
         DLIB_TEST(cubed(m4) == pointwise_multiply(m4,m4,m4));
+        DLIB_TEST(m4 == pointwise_divide(squared(m4),m4));
+        DLIB_TEST(m4 == pointwise_divide(cubed(m4),m4,m4));
+        DLIB_TEST(m4 == pointwise_divide(pointwise_multiply(cubed(m4),m4),m4,m4,m4));
+        DLIB_TEST(squared(m4) == pointwise_divide(cubed(m4),m4));
         DLIB_TEST(pow(matrix_cast<double>(m4),2) == squared(matrix_cast<double>(m4)));
         DLIB_TEST(pow(matrix_cast<double>(m4),3) == cubed(matrix_cast<double>(m4)));
 
@@ -323,7 +327,11 @@ namespace
 
 
         set_all_elements(v,2);
-        v2 =  pointwise_multiply(v, v*2);
+        v2 = pointwise_divide(v*2,v);
+        DLIB_TEST(v == v2);
+        DLIB_TEST(v == tmp(v2));
+
+        v2 = pointwise_multiply(v,v*2);
         set_all_elements(v,8);
         DLIB_TEST(v == v2);
         DLIB_TEST(v == tmp(v2));
@@ -336,6 +344,8 @@ namespace
         m5 = array2;
         DLIB_TEST((m5*2 == pointwise_multiply(m5,uniform_matrix<int,3,3,2>())));
         DLIB_TEST((tmp(m5*2) == tmp(pointwise_multiply(m5,uniform_matrix<int,3,3,2>()))));
+        DLIB_TEST((m5/2 == pointwise_divide(m5,uniform_matrix<int,3,3,2>())));
+        DLIB_TEST((tmp(m5/2) == tmp(pointwise_divide(m5,uniform_matrix<int,3,3,2>()))));
 
         v = tmp(v);
 
@@ -535,10 +545,10 @@ namespace
             set_all_elements(bt1,2);
             set_all_elements(bt2,3);
 
-			float val = trans(bt1)*bt2;
+            float val = trans(bt1)*bt2;
             DLIB_TEST((float)(trans(bt1)*bt2) == 18);
             DLIB_TEST((float)(trans(bt1)*bt2) != 19);
-			DLIB_TEST(val == 18);
+            DLIB_TEST(val == 18);
         }
         {
             matrix<float,3,1> bt1;
@@ -546,10 +556,10 @@ namespace
             set_all_elements(bt1,2);
             set_all_elements(bt2,3);
 
-			float val = trans(bt1)*bt2;
+            float val = trans(bt1)*bt2;
             DLIB_TEST((float)(trans(bt1)*bt2) == 18);
             DLIB_TEST((float)(trans(bt1)*bt2) != 19);
-			DLIB_TEST(val == 18);
+            DLIB_TEST(val == 18);
         }
         {
             matrix<float> bt1(3,1);
@@ -557,10 +567,10 @@ namespace
             set_all_elements(bt1,2);
             set_all_elements(bt2,3);
 
-			float val = trans(bt1)*bt2;
+            float val = trans(bt1)*bt2;
             DLIB_TEST((float)(trans(bt1)*bt2) == 18);
             DLIB_TEST((float)(trans(bt1)*bt2) != 19);
-			DLIB_TEST(val == 18);
+            DLIB_TEST(val == 18);
         }
         {
             matrix<float,3,1> bt1;
@@ -568,10 +578,10 @@ namespace
             set_all_elements(bt1,2);
             set_all_elements(bt2,3);
 
-			float val = trans(bt1)*bt2;
+            float val = trans(bt1)*bt2;
             DLIB_TEST((float)(trans(bt1)*bt2) == 18);
             DLIB_TEST((float)(trans(bt1)*bt2) != 19);
-			DLIB_TEST(val == 18);
+            DLIB_TEST(val == 18);
         }
 
 
@@ -968,7 +978,8 @@ namespace
             m = val1;
             m2 = val2;
 
-            DLIB_TEST(equal(reciprocal(m) , m2));
+            DLIB_TEST(equal(reciprocal(m),m2));
+            DLIB_TEST(equal(pointwise_multiply(m,m2),pointwise_divide(m,m)));
         }
         {
             matrix<complex<float> > m(2,2), m2(2,2);
@@ -976,7 +987,8 @@ namespace
             m = val1;
             m2 = val2;
 
-            DLIB_TEST(equal(reciprocal(m) , m2));
+            DLIB_TEST(equal(reciprocal(m),m2));
+            DLIB_TEST(equal(pointwise_multiply(m,m2),pointwise_divide(m,m)));
         }
 
         {