matrix_math_functions.h 20 KB
Newer Older
1
// Copyright (C) 2006  Davis E. King (davis@dlib.net)
2
3
4
5
// License: Boost Software License   See LICENSE.txt for the full license.
#ifndef DLIB_MATRIx_MATH_FUNCTIONS
#define DLIB_MATRIx_MATH_FUNCTIONS 

6
#include "matrix_math_functions_abstract.h"
7
8
9
10
11
12
13
14
15
16
17
18
19
#include "matrix_utilities.h"
#include "matrix.h"
#include "../algs.h"
#include <cmath>
#include <complex>
#include <limits>


namespace dlib
{

// ----------------------------------------------------------------------------------------

20
#define DLIB_MATRIX_SIMPLE_STD_FUNCTION(name,extra_cost) struct op_##name {     \
21
22
    template <typename EXP>                                                     \
    struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>          \
23
    {                                                                           \
24
        const static long cost = EXP::cost+(extra_cost);                        \
25
26
27
28
        typedef typename EXP::type type;                                        \
        template <typename M>                                                   \
        static type apply ( const M& m, long r, long c)                         \
        { return static_cast<type>(std::name(m(r,c))); }                        \
29
    };};                                                                        \
30
    template < typename EXP >                                                   \
31
    const matrix_unary_exp<EXP,op_##name> name (                                \
32
33
        const matrix_exp<EXP>& m)                                               \
    {                                                                           \
34
        return matrix_unary_exp<EXP,op_##name>(m.ref());                        \
35
36
37
38
    }                                                                           

// ----------------------------------------------------------------------------------------

39
40
41
42
DLIB_MATRIX_SIMPLE_STD_FUNCTION(sqrt,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(log,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(log10,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(exp,7)
43

44
DLIB_MATRIX_SIMPLE_STD_FUNCTION(conj,1)
45

46
47
DLIB_MATRIX_SIMPLE_STD_FUNCTION(ceil,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(floor,7)
48

49
50
51
52
53
54
55
56
57
DLIB_MATRIX_SIMPLE_STD_FUNCTION(sin,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(cos,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(tan,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(sinh,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(cosh,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(tanh,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(asin,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(acos,7)
DLIB_MATRIX_SIMPLE_STD_FUNCTION(atan,7)
58
59
60

// ----------------------------------------------------------------------------------------

61
    struct op_sigmoid 
62
    {
63
64
65
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
66
            const static long cost = EXP::cost+7;
67
68
69
70
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
71
                return static_cast<type>(1/(1 + std::exp(-m(r,c))));
72
73
            }
        };
74
75
76
77
78
    };

    template <
        typename EXP
        >
79
    const matrix_unary_exp<EXP,op_sigmoid> sigmoid (
80
81
82
        const matrix_exp<EXP>& m
    )
    {
83
        return matrix_unary_exp<EXP,op_sigmoid>(m.ref());
84
85
86
87
    }

// ----------------------------------------------------------------------------------------

88
    struct op_round_zeros 
89
    {
90
91
92
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
93
            const static long cost = EXP::cost+7;
94
95
96
97
98
99
100
101
102
103
104
            typedef typename EXP::type type;
            template <typename M, typename T>
            static type apply ( const M& m, const T& eps, long r, long c)
            { 
                const type temp = m(r,c);
                if (temp >= eps || temp <= -eps)
                    return temp;
                else
                    return 0;
            }
        };
105
106
107
108
109
    };

    template <
        typename EXP
        >
110
    const matrix_scalar_binary_exp<EXP,typename EXP::type,op_round_zeros> round_zeros (
111
112
113
        const matrix_exp<EXP>& m
    )
    {
114
115
        // you can only round matrices that contain built in scalar types like double, long, float, etc...
        COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value);
116
117
        typedef matrix_scalar_binary_exp<EXP,typename EXP::type, op_round_zeros> exp;
        return exp(m.ref(),10*std::numeric_limits<typename EXP::type>::epsilon());
118
119
120
121
122
    }

    template <
        typename EXP
        >
123
    const matrix_scalar_binary_exp<EXP,typename EXP::type,op_round_zeros> round_zeros (
124
125
126
127
        const matrix_exp<EXP>& m,
        typename EXP::type eps 
    )
    {
128
129
        // you can only round matrices that contain built in scalar types like double, long, float, etc...
        COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value);
130
        return matrix_scalar_binary_exp<EXP,typename EXP::type, op_round_zeros>(m.ref(),eps);
131
132
133
134
    }

// ----------------------------------------------------------------------------------------

135
    struct op_cubed 
136
    {
137
138
139
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
140
            const static long cost = EXP::cost+7;
141
142
143
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
144
145
146
147
            { 
                const type temp = m(r,c);
                return temp*temp*temp; 
            }
148
        };
149
150
151
152
153
    };

    template <
        typename EXP
        >
154
    const matrix_unary_exp<EXP,op_cubed> cubed (
155
156
157
        const matrix_exp<EXP>& m
    )
    {
158
        return matrix_unary_exp<EXP,op_cubed>(m.ref());
159
160
161
162
    }

// ----------------------------------------------------------------------------------------

163
    struct op_squared
164
    {
165
166
167
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
168
            const static long cost = EXP::cost+6;
169
170
171
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
172
173
174
175
            { 
                const type temp = m(r,c);
                return temp*temp; 
            }
176
        };
177
178
179
180
181
    };

    template <
        typename EXP
        >
182
    const matrix_unary_exp<EXP,op_squared> squared (
183
184
185
        const matrix_exp<EXP>& m
    )
    {
186
        return matrix_unary_exp<EXP,op_squared>(m.ref());
187
188
189
190
    }

// ----------------------------------------------------------------------------------------

191
    struct op_pow
192
    {
193
194
195
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
196
            const static long cost = EXP::cost+7;
197
198
199
200
201
            typedef typename EXP::type type;
            template <typename M, typename S>
            static type apply ( const M& m, const S& s, long r, long c)
            { return static_cast<type>(std::pow(m(r,c),s)); }
        };
202
203
204
205
206
207
    };

    template <
        typename EXP,
        typename S
        >
208
    const matrix_scalar_binary_exp<EXP,typename EXP::type,op_pow> pow (
209
210
211
212
213
214
215
216
217
218
        const matrix_exp<EXP>& m,
        const S& s
    )
    {
        // you can only round matrices that contain floats, doubles or long doubles.
        COMPILE_TIME_ASSERT((
                is_same_type<typename EXP::type,float>::value == true || 
                is_same_type<typename EXP::type,double>::value == true || 
                is_same_type<typename EXP::type,long double>::value == true 
        ));
219
        return matrix_scalar_binary_exp<EXP,typename EXP::type,op_pow>(m.ref(),s);
220
221
    }

222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
// ----------------------------------------------------------------------------------------

    struct op_pow2
    {
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
            const static long cost = EXP::cost+7;
            typedef typename EXP::type type;
            template <typename M, typename S>
            static type apply ( const M& m, const S& s, long r, long c)
            { return static_cast<type>(std::pow(s,m(r,c))); }
        };
    };

    template <
        typename EXP,
        typename S
        >
    const matrix_scalar_binary_exp<EXP,typename EXP::type,op_pow2> pow (
        const S& s,
        const matrix_exp<EXP>& m
    )
    {
        // you can only round matrices that contain floats, doubles or long doubles.
        COMPILE_TIME_ASSERT((
                is_same_type<typename EXP::type,float>::value == true || 
                is_same_type<typename EXP::type,double>::value == true || 
                is_same_type<typename EXP::type,long double>::value == true 
        ));
        return matrix_scalar_binary_exp<EXP,typename EXP::type,op_pow2>(m.ref(),s);
    }

255
256
// ----------------------------------------------------------------------------------------

257
    struct op_reciprocal
258
    {
259
260
261
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
262
            const static long cost = EXP::cost+6;
263
264
265
266
267
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                const type temp = m(r,c);
268
                if (temp != static_cast<type>(0))
269
                    return static_cast<type>((type)1.0/temp);
270
271
272
273
                else
                    return 0;
            }
        };
274
275
276
277
278
    };

    template <
        typename EXP
        >
279
    const matrix_unary_exp<EXP,op_reciprocal> reciprocal (
280
281
282
283
284
285
286
        const matrix_exp<EXP>& m
    )
    {
        // you can only compute reciprocal matrices that contain floats, doubles or long doubles.
        COMPILE_TIME_ASSERT((
                is_same_type<typename EXP::type,float>::value == true || 
                is_same_type<typename EXP::type,double>::value == true || 
287
288
289
290
                is_same_type<typename EXP::type,long double>::value == true  ||
                is_same_type<typename EXP::type,std::complex<float> >::value == true || 
                is_same_type<typename EXP::type,std::complex<double> >::value == true || 
                is_same_type<typename EXP::type,std::complex<long double> >::value == true 
291
        ));
292
        return matrix_unary_exp<EXP,op_reciprocal>(m.ref());
293
294
    }

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
// ----------------------------------------------------------------------------------------

    struct op_reciprocal_max
    {
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
            const static long cost = EXP::cost+6;
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                const type temp = m(r,c);
                if (temp != static_cast<type>(0))
                    return static_cast<type>((type)1.0/temp);
                else
                    return std::numeric_limits<type>::max();
            }
        };
    };

    template <
        typename EXP
        >
    const matrix_unary_exp<EXP,op_reciprocal_max> reciprocal_max (
        const matrix_exp<EXP>& m
    )
    {
        // you can only compute reciprocal_max matrices that contain floats, doubles or long doubles.
        COMPILE_TIME_ASSERT((
                is_same_type<typename EXP::type,float>::value == true || 
                is_same_type<typename EXP::type,double>::value == true || 
                is_same_type<typename EXP::type,long double>::value == true 
        ));
        return matrix_unary_exp<EXP,op_reciprocal_max>(m.ref());
    }

332
333
// ----------------------------------------------------------------------------------------

334
    struct op_normalize
335
    {
336
337
338
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
339
            const static long cost = EXP::cost+5;
340
341
342
343
344
345
346
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, const type& s, long r, long c)
            { 
                return m(r,c)*s;
            }
        };
347
348
349
350
351
    };

    template <
        typename EXP
        >
352
    const matrix_scalar_binary_exp<EXP,typename EXP::type,op_normalize> normalize (
353
354
355
356
357
358
359
360
361
        const matrix_exp<EXP>& m
    )
    {
        // you can only compute normalized matrices that contain floats, doubles or long doubles.
        COMPILE_TIME_ASSERT((
                is_same_type<typename EXP::type,float>::value == true || 
                is_same_type<typename EXP::type,double>::value == true || 
                is_same_type<typename EXP::type,long double>::value == true 
        ));
362
        typedef matrix_scalar_binary_exp<EXP,typename EXP::type, op_normalize> exp;
363
364
365
366
367

        typename EXP::type temp = std::sqrt(sum(squared(m)));
        if (temp != 0.0)
            temp = 1.0/temp;

368
        return exp(m.ref(),temp);
369
370
371
372
    }

// ----------------------------------------------------------------------------------------

373
    struct op_round
374
    {
375
        template <typename EXP, typename enabled = void>
376
377
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
378
            const static long cost = EXP::cost+7;
379
380
381
382
383
384
385
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                return static_cast<type>(std::floor(m(r,c)+0.5)); 
            }
        };
386
387
388
389
390

        template <typename EXP>
        struct op<EXP,typename enable_if_c<std::numeric_limits<typename EXP::type>::is_integer>::type > 
                : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
391
            const static long cost = EXP::cost;
392
393
394
395
396
397
398
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                return m(r,c);
            }
        };
399
400
401
402
403
    };

    template <
        typename EXP
        >
404
    const matrix_unary_exp<EXP,op_round> round (
405
406
407
        const matrix_exp<EXP>& m
    )
    {
408
409
        // you can only round matrices that contain built in scalar types like double, long, float, etc...
        COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value);
410
411
        typedef matrix_unary_exp<EXP,op_round> exp;
        return exp(m.ref());
412
413
    }

414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
// ----------------------------------------------------------------------------------------

    struct op_abs
    {
        template <typename EXP, typename return_type = typename EXP::type>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
            const static long cost = EXP::cost+7;
            typedef typename EXP::type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                return static_cast<type>(std::abs(m(r,c))); 
            }
        };

        template <typename EXP, typename T>
        struct op<EXP, std::complex<T> > : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
            const static long cost = EXP::cost;
            typedef T type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                return static_cast<type>(std::abs(m(r,c))); 
            }
        };
    };

    template <
        typename EXP
        >
    const matrix_unary_exp<EXP,op_abs> abs (
        const matrix_exp<EXP>& m
    )
    {
        typedef matrix_unary_exp<EXP,op_abs> exp;
        return exp(m.ref());
    }

454
455
// ----------------------------------------------------------------------------------------

456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
    struct op_complex_matrix 
    {
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
            const static long cost = EXP::cost+1;
            typedef std::complex<typename EXP::type> type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { 
                return type(m(r,c));
            }
        };
    };

    template <
        typename EXP
        >
    const matrix_unary_exp<EXP,op_complex_matrix> complex_matrix (
        const matrix_exp<EXP>& m
    )
    {
        return matrix_unary_exp<EXP,op_complex_matrix>(m.ref());
    }

// ----------------------------------------------------------------------------------------

    struct op_complex_matrix2
484
    {
485
486
487
        template <typename EXP1, typename EXP2>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP1,EXP2>
        {
488
            const static long cost = EXP1::cost+EXP2::cost+1;
489
490
491
492
493
494
            typedef std::complex<typename EXP1::type> type;

            template <typename M1, typename M2>
            static type apply ( const M1& m1, const M2& m2 , long r, long c)
            { return type(m1(r,c),m2(r,c)); }
        };
495
496
497
498
499
500
    };

    template <
        typename EXP1,
        typename EXP2
        >
501
    const matrix_binary_exp<EXP1,EXP2,op_complex_matrix2> complex_matrix (
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
        const matrix_exp<EXP1>& real_part,
        const matrix_exp<EXP2>& imag_part 
    )
    {
        COMPILE_TIME_ASSERT((is_same_type<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(real_part.nr() == imag_part.nr() &&
               real_part.nc() == imag_part.nc(), 
            "\tconst matrix_exp::type complex_matrix(real_part, imag_part)"
            << "\n\tYou can only make a complex matrix from two equally sized matrices"
            << "\n\treal_part.nr(): " << real_part.nr()
            << "\n\treal_part.nc(): " << real_part.nc() 
            << "\n\timag_part.nr(): " << imag_part.nr()
            << "\n\timag_part.nc(): " << imag_part.nc() 
            );
519
        typedef matrix_binary_exp<EXP1,EXP2,op_complex_matrix2> exp;
520
        return exp(real_part.ref(),imag_part.ref());
521
522
523
524
    }

// ----------------------------------------------------------------------------------------

525
    struct op_norm
526
    {
527
528
529
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
530
            const static long cost = EXP::cost+6;
531
532
533
534
535
            typedef typename EXP::type::value_type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { return std::norm(m(r,c)); }
        };
536
537
538
539
540
    };

    template <
        typename EXP
        >
541
    const matrix_unary_exp<EXP,op_norm> norm (
542
543
544
        const matrix_exp<EXP>& m
    )
    {
545
546
        typedef matrix_unary_exp<EXP,op_norm> exp;
        return exp(m.ref());
547
548
549
550
    }

// ----------------------------------------------------------------------------------------

551
    struct op_real
552
    {
553
554
555
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
556
            const static long cost = EXP::cost;
557
558
559
560
561
            typedef typename EXP::type::value_type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { return std::real(m(r,c)); }
        };
562
563
564
565
566
    };

    template <
        typename EXP
        >
567
    const matrix_unary_exp<EXP,op_real> real (
568
569
570
        const matrix_exp<EXP>& m
    )
    {
571
572
        typedef matrix_unary_exp<EXP,op_real> exp;
        return exp(m.ref());
573
574
575
576
    }

// ----------------------------------------------------------------------------------------

577
    struct op_imag
578
    {
579
580
581
        template <typename EXP>
        struct op : has_nondestructive_aliasing, preserves_dimensions<EXP>
        {
582
            const static long cost = EXP::cost;
583
584
585
586
587
            typedef typename EXP::type::value_type type;
            template <typename M>
            static type apply ( const M& m, long r, long c)
            { return std::imag(m(r,c)); }
        };
588
589
590
591
592
    };

    template <
        typename EXP
        >
593
    const matrix_unary_exp<EXP,op_imag> imag (
594
595
596
        const matrix_exp<EXP>& m
    )
    {
597
598
        typedef matrix_unary_exp<EXP,op_imag> exp;
        return exp(m.ref());
599
600
601
602
603
604
605
606
    }

// ----------------------------------------------------------------------------------------

}

#endif // DLIB_MATRIx_MATH_FUNCTIONS