matrix.h 58.8 KB
Newer Older
1
2
3
4
5
// Copyright (C) 2006  Davis E. King (davisking@users.sourceforge.net)
// License: Boost Software License   See LICENSE.txt for the full license.
#ifndef DLIB_MATRIx_
#define DLIB_MATRIx_

6
#include "matrix_fwd.h"
7
8
9
10
11
12
13
#include "matrix_abstract.h"
#include "../algs.h"
#include "../serialize.h"
#include "../enable_if.h"
#include <sstream>
#include <algorithm>
#include "../memory_manager.h"
14
#include "../is_kind.h"
15
#include "matrix_data_layout.h"
16
#include "matrix_assign_fwd.h"
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

#ifdef _MSC_VER
// Disable the following warnings for Visual Studio

// This warning is:
//    "warning C4355: 'this' : used in base member initializer list"
// Which we get from this code but it is not an error so I'm turning this
// warning off and then turning it back on at the end of the file.
#pragma warning(disable : 4355)

#endif

namespace dlib
{

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

    // We want to return the compile time constant if our NR and NC dimensions
    // aren't zero but if they are then we want to call ref_.nx() and return
    // the correct values. 
39
    template < typename exp_type, long NR >
40
41
    struct get_nr_helper
    {
42
        static inline long get(const exp_type&) { return NR; }
43
44
    };

45
46
    template < typename exp_type >
    struct get_nr_helper<exp_type,0>
47
    {
48
        static inline long get(const exp_type& m) { return m.nr(); }
49
50
    };

51
    template < typename exp_type, long NC >
52
53
    struct get_nc_helper
    {
54
        static inline long get(const exp_type&) { return NC; }
55
56
    };

57
58
    template < typename exp_type >
    struct get_nc_helper<exp_type,0>
59
    {
60
        static inline long get(const exp_type& m) { return m.nc(); }
61
62
    };

63
64
65
66
67
    template <typename EXP>
    struct matrix_traits
    {
        typedef typename EXP::type type;
        typedef typename EXP::mem_manager_type mem_manager_type;
68
        typedef typename EXP::layout_type layout_type;
69
70
71
72
        const static long NR = EXP::NR;
        const static long NC = EXP::NC;
        const static long cost = EXP::cost;
    };
73
74
75
76

    template <
        typename EXP
        >
77
    class matrix_exp 
78
    {
79
80
81
82
83
84
        /*!
            REQUIREMENTS ON EXP
                EXP should be something convertible to a matrix_exp.  That is,
                it should inherit from matrix_exp
        !*/

85
    public:
86
87
        typedef typename matrix_traits<EXP>::type type;
        typedef typename matrix_traits<EXP>::mem_manager_type mem_manager_type;
88
        typedef typename matrix_traits<EXP>::layout_type layout_type;
89
90
91
        const static long NR = matrix_traits<EXP>::NR;
        const static long NC = matrix_traits<EXP>::NC;
        const static long cost = matrix_traits<EXP>::cost;
92

93
        typedef matrix<type,NR,NC,mem_manager_type,layout_type> matrix_type;
94
        typedef EXP exp_type;
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

        inline const type operator() (
            long r,
            long c
        ) const 
        { 
            DLIB_ASSERT(r < nr() && c < nc() && r >= 0 && c >= c, 
                "\tconst type matrix_exp::operator(r,c)"
                << "\n\tYou must give a valid row and column"
                << "\n\tr:    " << r 
                << "\n\tc:    " << c
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return ref_(r,c); 
        }

        const type operator() (
            long i
        ) const 
        {
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
            DLIB_ASSERT(nc() == 1 || nr() == 1, 
                "\tconst type matrix_exp::operator(i)"
                << "\n\tYou can only use this operator on column or row vectors"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, 
                "\tconst type matrix_exp::operator(i)"
                << "\n\tYou must give a valid row/column number"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            if (nc() == 1)
                return ref_(i,0);
            else
                return ref_(0,i);
        }

        long size (
        ) const { return nr()*nc(); }

        long nr (
144
        ) const { return get_nr_helper<exp_type,NR>::get(ref_); }
145
146

        long nc (
147
        ) const { return get_nc_helper<exp_type,NC>::get(ref_); }
148

149
        template <typename U, long iNR, long iNC, typename mm, typename l >
150
        bool aliases (
151
            const matrix<U,iNR,iNC,mm,l>& item
152
153
        ) const { return ref_.aliases(item); }

154
        template <typename U, long iNR, long iNC , typename mm, typename l>
155
        bool destructively_aliases (
156
            const matrix<U,iNR,iNC,mm,l>& item
157
158
        ) const { return ref_.destructively_aliases(item); }

159
        const exp_type& ref (
160
161
162
163
164
165
166
167
168
169
170
171
172
173
        ) const { return ref_; }

        inline operator const type (
        ) const 
        {
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0);
            COMPILE_TIME_ASSERT(NR == 1 || NR == 0);
            DLIB_ASSERT(nr() == 1 && nc() == 1, 
                "\tmatrix_exp::operator const type&() const"
                << "\n\tYou can only use this operator on a 1x1 matrix"
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
174
175
176
177
178
179
180

            // Put the expression contained in this matrix_exp into
            // a temporary 1x1 matrix so that the expression will encounter
            // all the overloads of matrix_assign() and have the chance to
            // go through any applicable optimizations.
            matrix<type,1,1> temp(ref_);
            return temp(0);
181
182
        }

183
184
185
186
187
    protected:
        explicit matrix_exp (
            const EXP& exp
        ) : ref_(exp) {}

188
189
190

    private:

191
192
193
194
        // you can't copy a matrix_exp at all.  Things that inherit from it must
        // define their own copy constructors that call the above protected 
        // constructor so that the reference below is maintained correctly.
        matrix_exp(const matrix_exp& item);
195
        matrix_exp& operator= (const matrix_exp&);
196

197
        const exp_type& ref_;
198
199
    };

200
201
// ----------------------------------------------------------------------------------------

202
// something is a matrix if it is convertible to a matrix_exp object
203
    template <typename T>
204
205
    struct is_matrix<T, typename enable_if<is_convertible<T, const matrix_exp<typename T::exp_type>& > >::type > 
    { static const bool value = true; }; 
206
207
208
209
    /*
        is_matrix<T>::value == 1 if T is a matrix type else 0
    */

210
// ----------------------------------------------------------------------------------------
211
212
213
214
215
216
217
218
219
220
221
222

    // This template will perform the needed loop for element multiplication using whichever
    // dimension is provided as a compile time constant (if one is at all).
    template <
        typename LHS,
        typename RHS,
        long lhs_nc = LHS::NC,
        long rhs_nr = RHS::NR
        >
    struct matrix_multiply_helper 
    {
        typedef typename LHS::type type;
223
        template <typename RHS_, typename LHS_>
224
        inline const static type  eval (
225
226
            const RHS_& rhs,
            const LHS_& lhs,
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
            long r, 
            long c
        )  
        { 
            type temp = type();
            for (long i = 0; i < rhs.nr(); ++i)
            {
                temp += lhs(r,i)*rhs(i,c);
            }
            return temp;
        }
    };

    template <
        typename LHS,
        typename RHS,
        long lhs_nc 
        >
    struct matrix_multiply_helper <LHS,RHS,lhs_nc,0>
    {
        typedef typename LHS::type type;
248
        template <typename RHS_, typename LHS_>
249
        inline const static type  eval (
250
251
            const RHS_& rhs,
            const LHS_& lhs,
252
253
254
255
256
257
258
259
260
261
262
263
264
            long r, 
            long c
        )  
        { 
            type temp = type();
            for (long i = 0; i < lhs.nc(); ++i)
            {
                temp += lhs(r,i)*rhs(i,c);
            }
            return temp;
        }
    };

265
266
267
268
269
270
271
272
    template <typename LHS, typename RHS>
    class matrix_multiply_exp;

    template <typename LHS, typename RHS>
    struct matrix_traits<matrix_multiply_exp<LHS,RHS> >
    {
        typedef typename LHS::type type;
        typedef typename LHS::mem_manager_type mem_manager_type;
273
        typedef typename LHS::layout_type layout_type;
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
        const static long NR = LHS::NR;
        const static long NC = RHS::NC;

        const static bool lhs_is_costly = (LHS::cost > 4)&&(RHS::NC != 1);
        const static bool rhs_is_costly = (RHS::cost > 4)&&(LHS::NR != 1);

        // Note that if we decide that one of the matrices is too costly we will evaluate it
        // into a temporary.  Doing this resets its cost back to 1.
        const static long lhs_cost = ((lhs_is_costly==true)? 1 : (LHS::cost));
        const static long rhs_cost = ((rhs_is_costly==true)? 1 : (RHS::cost));

        // The cost of evaluating an element of a matrix multiply is the cost of evaluating elements from
        // RHS and LHS times the number of rows/columns in the RHS/LHS matrix.  If we don't know the matrix
        // dimensions then just assume it is really large.
        const static long cost = (lhs_cost+rhs_cost)*
                                 ((tmax<LHS::NC,RHS::NR>::value!=0)? (tmax<LHS::NC,RHS::NR>::value):(10000));
    };

    template <typename T, bool is_ref> struct conditional_matrix_temp { typedef typename T::matrix_type type; };
    template <typename T> struct conditional_matrix_temp<T,true>      { typedef T& type; };

295
296
    template <
        typename LHS,
297
        typename RHS
298
        >
299
    class matrix_multiply_exp : public matrix_exp<matrix_multiply_exp<LHS,RHS> >
300
301
302
    {
        /*!
            REQUIREMENTS ON LHS AND RHS
303
                - must be matrix_exp objects.
304
305
306
        !*/
    public:

307
308
309
310
311
        typedef typename matrix_traits<matrix_multiply_exp>::type type;
        typedef typename matrix_traits<matrix_multiply_exp>::mem_manager_type mem_manager_type;
        const static long NR = matrix_traits<matrix_multiply_exp>::NR;
        const static long NC = matrix_traits<matrix_multiply_exp>::NC;
        const static long cost = matrix_traits<matrix_multiply_exp>::cost;
312
        typedef typename matrix_traits<matrix_multiply_exp>::layout_type layout_type;
313
314
315
316


        const static bool lhs_is_costly = matrix_traits<matrix_multiply_exp>::lhs_is_costly;
        const static bool rhs_is_costly = matrix_traits<matrix_multiply_exp>::rhs_is_costly;
317
318
        const static bool either_is_costly = lhs_is_costly || rhs_is_costly;
        const static bool both_are_costly = lhs_is_costly && rhs_is_costly;
319
320
321
322

        typedef typename conditional_matrix_temp<const LHS,lhs_is_costly == false>::type LHS_ref_type;
        typedef typename conditional_matrix_temp<const RHS,rhs_is_costly == false>::type RHS_ref_type;

323
324
325
326
327
328
329
330
        matrix_multiply_exp (
            const matrix_multiply_exp& item
        ) : matrix_exp<matrix_multiply_exp>(*this), lhs(item.lhs), rhs(item.rhs) {}

        // This constructor exists simply for the purpose of causing a compile time error if
        // someone tries to create an instance of this object with the wrong kind of objects.
        template <typename T1, typename T2>
        matrix_multiply_exp (T1,T2); 
331
332
333
334
335

        inline matrix_multiply_exp (
            const LHS& lhs_,
            const RHS& rhs_
        ) :
336
            matrix_exp<matrix_multiply_exp>(*this),
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
            lhs(lhs_),
            rhs(rhs_)
        {
            // You are trying to multiply two incompatible matrices together.  The number of columns 
            // in the matrix on the left must match the number of rows in the matrix on the right.
            COMPILE_TIME_ASSERT(LHS::NC == RHS::NR || LHS::NC*RHS::NR == 0);
            DLIB_ASSERT(lhs.nc() == rhs.nr(), 
                "\tconst matrix_exp operator*(const matrix_exp& lhs, const matrix_exp& rhs)"
                << "\n\tYou are trying to multiply two incompatible matrices together"
                << "\n\tlhs.nr(): " << lhs.nr()
                << "\n\tlhs.nc(): " << lhs.nc()
                << "\n\trhs.nr(): " << rhs.nr()
                << "\n\trhs.nc(): " << rhs.nc()
                << "\n\t&lhs: " << &lhs 
                << "\n\t&rhs: " << &rhs 
                );

            // You can't multiply matrices together if they don't both contain the same type of elements.
            COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
        }

        inline const type operator() (
            long r, 
            long c
        ) const 
        { 
            return matrix_multiply_helper<LHS,RHS>::eval(rhs,lhs,r,c);
        }

366
367
368
        inline const type operator() ( long i ) const 
        { return matrix_exp<matrix_multiply_exp>::operator()(i); }

369
370
371
372
373
374
        long nr (
        ) const { return lhs.nr(); }

        long nc (
        ) const { return rhs.nc(); }

375
        template <typename U, long iNR, long iNC, typename mm, typename l >
376
        bool aliases (
377
            const matrix<U,iNR,iNC,mm,l>& item
378
379
        ) const { return lhs.aliases(item) || rhs.aliases(item); }

380
        template <typename U, long iNR, long iNC , typename mm, typename l>
381
        bool destructively_aliases (
382
            const matrix<U,iNR,iNC,mm,l>& item
383
384
        ) const { return aliases(item); }

385
386
        LHS_ref_type lhs;
        RHS_ref_type rhs;
387
388
    };

389
390
    template < typename EXP1, typename EXP2 >
    inline const matrix_multiply_exp<EXP1, EXP2> operator* (
391
392
393
394
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
395
        return matrix_multiply_exp<EXP1, EXP2>(m1.ref(), m2.ref());
396
397
    }

398
399
400
401
402
403
404
405
406
407
408
409
410
411
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
    template <typename M, bool use_reference = true>
    class matrix_mul_scal_exp;

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

    // Now we declare some overloads that cause any scalar multiplications to percolate 
    // up and outside of any matrix multiplies.  Note that we are using the non-reference containing
    // mode of the matrix_mul_scal_exp object since we are passing in locally constructed matrix_multiply_exp 
    // objects.  So the matrix_mul_scal_exp object will contain copies of matrix_multiply_exp objects
    // rather than references to them.  This could result in extra matrix copies if the matrix_multiply_exp
    // decided it should evaluate any of its arguments.  So we also try to not apply this percolating operation 
    // if the matrix_multiply_exp would contain a fully evaluated copy of the original matrix_mul_scal_exp 
    // expression.
    // 
    // Also, the reason we want to apply this transformation in the first place is because it (1) makes
    // the expressions going into matrix multiply expressions simpler and (2) it makes it a lot more
    // straight forward to bind BLAS calls to matrix expressions involving scalar multiplies.
    template < typename EXP1, typename EXP2 >
    inline const typename disable_if_c< matrix_multiply_exp<matrix_mul_scal_exp<EXP1>, matrix_mul_scal_exp<EXP2> >::both_are_costly ,      
                                        matrix_mul_scal_exp<matrix_multiply_exp<EXP1, EXP2>,false> >::type operator* (
        const matrix_mul_scal_exp<EXP1>& m1,
        const matrix_mul_scal_exp<EXP2>& m2
    )
    {
        typedef matrix_multiply_exp<EXP1, EXP2> exp1;
        typedef matrix_mul_scal_exp<exp1,false> exp2;
        return exp2(exp1(m1.m, m2.m), m1.s*m2.s);
    }

    template < typename EXP1, typename EXP2 >
    inline const typename disable_if_c< matrix_multiply_exp<matrix_mul_scal_exp<EXP1>, EXP2 >::lhs_is_costly ,      
                                      matrix_mul_scal_exp<matrix_multiply_exp<EXP1, EXP2>,false> >::type operator* (
        const matrix_mul_scal_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
        typedef matrix_multiply_exp<EXP1, EXP2> exp1;
        typedef matrix_mul_scal_exp<exp1,false> exp2;
        return exp2(exp1(m1.m, m2.ref()), m1.s);
    }

    template < typename EXP1, typename EXP2 >
    inline const typename disable_if_c< matrix_multiply_exp<EXP1, matrix_mul_scal_exp<EXP2> >::rhs_is_costly ,      
                                      matrix_mul_scal_exp<matrix_multiply_exp<EXP1, EXP2>,false> >::type operator* (
        const matrix_exp<EXP1>& m1,
        const matrix_mul_scal_exp<EXP2>& m2
    )
    {
        typedef matrix_multiply_exp<EXP1, EXP2> exp1;
        typedef matrix_mul_scal_exp<exp1,false> exp2;
        return exp2(exp1(m1.ref(), m2.m), m2.s);
    }

451
// ----------------------------------------------------------------------------------------
452

453
454
    template <typename LHS, typename RHS>
    class matrix_add_exp;
455

456
457
    template <typename LHS, typename RHS>
    struct matrix_traits<matrix_add_exp<LHS,RHS> >
458
    {
459
460
        typedef typename LHS::type type;
        typedef typename LHS::mem_manager_type mem_manager_type;
461
        typedef typename LHS::layout_type layout_type;
462
463
464
465
        const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR;
        const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC;
        const static long cost = LHS::cost+RHS::cost;
    };
466
467
468
469
470

    template <
        typename LHS,
        typename RHS
        >
471
    class matrix_add_exp : public matrix_exp<matrix_add_exp<LHS,RHS> >
472
473
474
    {
        /*!
            REQUIREMENTS ON LHS AND RHS
475
                - must be matrix_exp objects. 
476
477
        !*/
    public:
478
479
480
481
482
        typedef typename matrix_traits<matrix_add_exp>::type type;
        typedef typename matrix_traits<matrix_add_exp>::mem_manager_type mem_manager_type;
        const static long NR = matrix_traits<matrix_add_exp>::NR;
        const static long NC = matrix_traits<matrix_add_exp>::NC;
        const static long cost = matrix_traits<matrix_add_exp>::cost;
483
        typedef typename matrix_traits<matrix_add_exp>::layout_type layout_type;
484
485
486
487

        matrix_add_exp (
            const matrix_add_exp& item
        ) : matrix_exp<matrix_add_exp>(*this), lhs(item.lhs), rhs(item.rhs) {}
488
489
490
491
492

        // This constructor exists simply for the purpose of causing a compile time error if
        // someone tries to create an instance of this object with the wrong kind of objects.
        template <typename T1, typename T2>
        matrix_add_exp (T1,T2); 
493
494

        matrix_add_exp (
495
496
497
            const LHS& lhs_,
            const RHS& rhs_
        ) :
498
            matrix_exp<matrix_add_exp>(*this),
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
            lhs(lhs_),
            rhs(rhs_)
        {
            // You can only add matrices together if they both have the same number of rows and columns.
            COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0);
            COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0);
            DLIB_ASSERT(lhs.nc() == rhs.nc() &&
                   lhs.nr() == rhs.nr(), 
                "\tconst matrix_exp operator+(const matrix_exp& lhs, const matrix_exp& rhs)"
                << "\n\tYou are trying to add two incompatible matrices together"
                << "\n\tlhs.nr(): " << lhs.nr()
                << "\n\tlhs.nc(): " << lhs.nc()
                << "\n\trhs.nr(): " << rhs.nr()
                << "\n\trhs.nc(): " << rhs.nc()
                << "\n\t&lhs: " << &lhs 
                << "\n\t&rhs: " << &rhs 
                );

            // You can only add matrices together if they both contain the same types of elements.
            COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
        }

        const type operator() (
            long r, 
            long c
        ) const { return lhs(r,c) + rhs(r,c); }

526
527
528
        inline const type operator() ( long i ) const 
        { return matrix_exp<matrix_add_exp>::operator()(i); }

529
        template <typename U, long iNR, long iNC , typename mm, typename l>
530
        bool aliases (
531
            const matrix<U,iNR,iNC,mm,l>& item
532
533
        ) const { return lhs.aliases(item) || rhs.aliases(item); }

534
        template <typename U, long iNR, long iNC, typename mm, typename l >
535
        bool destructively_aliases (
536
            const matrix<U,iNR,iNC,mm,l>& item
537
538
539
540
541
542
543
544
        ) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); }

        long nr (
        ) const { return lhs.nr(); }

        long nc (
        ) const { return lhs.nc(); }

545
546
        const LHS& lhs;
        const RHS& rhs;
547
548
549
550
551
552
    };

    template <
        typename EXP1,
        typename EXP2
        >
553
    inline const matrix_add_exp<EXP1, EXP2> operator+ (
554
555
556
557
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
558
        return matrix_add_exp<EXP1, EXP2>(m1.ref(),m2.ref());
559
560
561
562
    }

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

563
564
565
566
567
568
569
570
    template <typename LHS, typename RHS>
    class matrix_subtract_exp;

    template <typename LHS, typename RHS>
    struct matrix_traits<matrix_subtract_exp<LHS,RHS> >
    {
        typedef typename LHS::type type;
        typedef typename LHS::mem_manager_type mem_manager_type;
571
        typedef typename LHS::layout_type layout_type;
572
573
574
575
576
        const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR;
        const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC;
        const static long cost = LHS::cost+RHS::cost;
    };

577
578
579
580
    template <
        typename LHS,
        typename RHS
        >
581
    class matrix_subtract_exp : public matrix_exp<matrix_subtract_exp<LHS,RHS> >
582
583
584
    {
        /*!
            REQUIREMENTS ON LHS AND RHS
585
                - must be matrix_exp objects. 
586
587
        !*/
    public:
588
589
590
591
592
        typedef typename matrix_traits<matrix_subtract_exp>::type type;
        typedef typename matrix_traits<matrix_subtract_exp>::mem_manager_type mem_manager_type;
        const static long NR = matrix_traits<matrix_subtract_exp>::NR;
        const static long NC = matrix_traits<matrix_subtract_exp>::NC;
        const static long cost = matrix_traits<matrix_subtract_exp>::cost;
593
        typedef typename matrix_traits<matrix_subtract_exp>::layout_type layout_type;
594

595
596
597
598
599
600
601
602
        matrix_subtract_exp (
            const matrix_subtract_exp& item
        ) : matrix_exp<matrix_subtract_exp>(*this), lhs(item.lhs), rhs(item.rhs) {}

        // This constructor exists simply for the purpose of causing a compile time error if
        // someone tries to create an instance of this object with the wrong kind of objects.
        template <typename T1, typename T2>
        matrix_subtract_exp (T1,T2); 
603
604
605
606

        matrix_subtract_exp (
            const LHS& lhs_,
            const RHS& rhs_
607
        ) : matrix_exp<matrix_subtract_exp>(*this),
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
            lhs(lhs_),
            rhs(rhs_)
        {
            // You can only subtract one matrix from another if they both have the same number of rows and columns.
            COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0);
            COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0);
            DLIB_ASSERT(lhs.nc() == rhs.nc() &&
                   lhs.nr() == rhs.nr(), 
                "\tconst matrix_exp operator-(const matrix_exp& lhs, const matrix_exp& rhs)"
                << "\n\tYou are trying to add two incompatible matrices together"
                << "\n\tlhs.nr(): " << lhs.nr()
                << "\n\tlhs.nc(): " << lhs.nc()
                << "\n\trhs.nr(): " << rhs.nr()
                << "\n\trhs.nc(): " << rhs.nc()
                << "\n\t&lhs: " << &lhs 
                << "\n\t&rhs: " << &rhs 
                );

            // You can only subtract one matrix from another if they both contain elements of the same type.
            COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
        }

        const type operator() (
            long r, 
            long c
        ) const { return lhs(r,c) - rhs(r,c); }

635
636
637
        inline const type operator() ( long i ) const 
        { return matrix_exp<matrix_subtract_exp>::operator()(i); }

638
        template <typename U, long iNR, long iNC, typename mm, typename l >
639
        bool aliases (
640
            const matrix<U,iNR,iNC, mm,l>& item
641
642
        ) const { return lhs.aliases(item) || rhs.aliases(item); }

643
        template <typename U, long iNR, long iNC , typename mm, typename l>
644
        bool destructively_aliases (
645
            const matrix<U,iNR,iNC,mm,l>& item
646
647
648
649
650
651
652
653
        ) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); }

        long nr (
        ) const { return lhs.nr(); }

        long nc (
        ) const { return lhs.nc(); }

654
655
        const LHS& lhs;
        const RHS& rhs;
656
657
658
659
660
661
    };

    template <
        typename EXP1,
        typename EXP2
        >
662
    inline const matrix_subtract_exp<EXP1, EXP2> operator- (
663
664
665
666
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
667
        return matrix_subtract_exp<EXP1, EXP2>(m1.ref(),m2.ref());
668
669
670
671
    }

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

Davis King's avatar
Davis King committed
672
    template <typename M>
673
674
    class matrix_div_scal_exp;

Davis King's avatar
Davis King committed
675
676
    template <typename M>
    struct matrix_traits<matrix_div_scal_exp<M> >
677
678
679
    {
        typedef typename M::type type;
        typedef typename M::mem_manager_type mem_manager_type;
680
        typedef typename M::layout_type layout_type;
681
682
683
684
685
        const static long NR = M::NR;
        const static long NC = M::NC;
        const static long cost = M::cost+1;
    };

686
    template <
Davis King's avatar
Davis King committed
687
        typename M
688
        >
Davis King's avatar
Davis King committed
689
    class matrix_div_scal_exp : public matrix_exp<matrix_div_scal_exp<M> >
690
691
692
    {
        /*!
            REQUIREMENTS ON M 
693
                - must be a matrix_exp object.
694
695
        !*/
    public:
696
697
698
699
700
        typedef typename matrix_traits<matrix_div_scal_exp>::type type;
        typedef typename matrix_traits<matrix_div_scal_exp>::mem_manager_type mem_manager_type;
        const static long NR = matrix_traits<matrix_div_scal_exp>::NR;
        const static long NC = matrix_traits<matrix_div_scal_exp>::NC;
        const static long cost = matrix_traits<matrix_div_scal_exp>::cost;
701
        typedef typename matrix_traits<matrix_div_scal_exp>::layout_type layout_type;
702

703
704
705
706
707
708
709
        matrix_div_scal_exp (
            const matrix_div_scal_exp& item
        ) : matrix_exp<matrix_div_scal_exp>(*this), m(item.m), s(item.s) {}

        // This constructor exists simply for the purpose of causing a compile time error if
        // someone tries to create an instance of this object with the wrong kind of objects.
        template <typename T1>
Davis King's avatar
Davis King committed
710
        matrix_div_scal_exp (T1, const type&); 
711

712
        matrix_div_scal_exp (
713
            const M& m_,
Davis King's avatar
Davis King committed
714
            const type& s_
715
        ) :
716
            matrix_exp<matrix_div_scal_exp>(*this),
717
718
719
720
721
722
723
724
725
            m(m_),
            s(s_)
        {}

        const type operator() (
            long r, 
            long c
        ) const { return m(r,c)/s; }

726
727
728
        inline const type operator() ( long i ) const 
        { return matrix_exp<matrix_div_scal_exp>::operator()(i); }

729
        template <typename U, long iNR, long iNC, typename mm , typename l>
730
        bool aliases (
731
            const matrix<U,iNR,iNC,mm,l>& item
732
733
        ) const { return m.aliases(item); }

734
        template <typename U, long iNR, long iNC, typename mm, typename l >
735
        bool destructively_aliases (
736
            const matrix<U,iNR,iNC,mm,l>& item
737
738
739
740
741
742
743
744
        ) const { return m.destructively_aliases(item); }

        long nr (
        ) const { return m.nr(); }

        long nc (
        ) const { return m.nc(); }

745
        const M& m;
Davis King's avatar
Davis King committed
746
        const type s;
747
748
749
750
    };

    template <
        typename EXP,
751
        typename S
752
        >
753
    inline const typename enable_if_c<std::numeric_limits<typename EXP::type>::is_integer, matrix_div_scal_exp<EXP> >::type operator/ (
754
755
756
757
        const matrix_exp<EXP>& m,
        const S& s
    )
    {
758
        return matrix_div_scal_exp<EXP>(m.ref(),static_cast<typename EXP::type>(s));
759
760
761
762
    }

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

763
764
    template <typename M, bool use_reference >
    struct matrix_traits<matrix_mul_scal_exp<M,use_reference> >
765
766
767
    {
        typedef typename M::type type;
        typedef typename M::mem_manager_type mem_manager_type;
768
        typedef typename M::layout_type layout_type;
769
770
771
772
773
        const static long NR = M::NR;
        const static long NC = M::NC;
        const static long cost = M::cost+1;
    };

774
775
776
777
    template <typename T, bool is_ref> struct conditional_reference { typedef T type; };
    template <typename T> struct conditional_reference<T,true>      { typedef T& type; };


778
    template <
779
780
        typename M,
        bool use_reference
781
        >
782
    class matrix_mul_scal_exp : public matrix_exp<matrix_mul_scal_exp<M,use_reference> >
783
784
785
    {
        /*!
            REQUIREMENTS ON M 
786
                - must be a matrix_exp object.
787
788
789

        !*/
    public:
790
791
792
793
794
        typedef typename matrix_traits<matrix_mul_scal_exp>::type type;
        typedef typename matrix_traits<matrix_mul_scal_exp>::mem_manager_type mem_manager_type;
        const static long NR = matrix_traits<matrix_mul_scal_exp>::NR;
        const static long NC = matrix_traits<matrix_mul_scal_exp>::NC;
        const static long cost = matrix_traits<matrix_mul_scal_exp>::cost;
795
        typedef typename matrix_traits<matrix_mul_scal_exp>::layout_type layout_type;
796

797
        matrix_mul_scal_exp (
798
799
800
801
802
803
            const matrix_mul_scal_exp& item
        ) : matrix_exp<matrix_mul_scal_exp>(*this), m(item.m), s(item.s) {}

        // This constructor exists simply for the purpose of causing a compile time error if
        // someone tries to create an instance of this object with the wrong kind of objects.
        template <typename T1>
Davis King's avatar
Davis King committed
804
        matrix_mul_scal_exp (T1, const type&); 
805
806

        matrix_mul_scal_exp (
807
            const M& m_,
Davis King's avatar
Davis King committed
808
            const type& s_
809
        ) :
810
            matrix_exp<matrix_mul_scal_exp>(*this),
811
812
813
814
815
816
817
818
819
            m(m_),
            s(s_)
        {}

        const type operator() (
            long r, 
            long c
        ) const { return m(r,c)*s; }

820
821
822
        inline const type operator() ( long i ) const 
        { return matrix_exp<matrix_mul_scal_exp>::operator()(i); }

823
        template <typename U, long iNR, long iNC , typename mm, typename l>
824
        bool aliases (
825
            const matrix<U,iNR,iNC,mm,l>& item
826
827
        ) const { return m.aliases(item); }

828
        template <typename U, long iNR, long iNC, typename mm, typename l >
829
        bool destructively_aliases (
830
            const matrix<U,iNR,iNC,mm,l>& item
831
832
833
834
835
836
837
838
        ) const { return m.destructively_aliases(item); }

        long nr (
        ) const { return m.nr(); }

        long nc (
        ) const { return m.nc(); }

839
840
841
        typedef typename conditional_reference<const M,use_reference>::type M_ref_type;

        M_ref_type m;
Davis King's avatar
Davis King committed
842
        const type s;
843
844
845
846
847
848
    };

    template <
        typename EXP,
        typename S 
        >
Davis King's avatar
Davis King committed
849
    inline typename disable_if<is_matrix<S>, const matrix_mul_scal_exp<EXP> >::type operator* (
850
851
852
853
        const matrix_exp<EXP>& m,
        const S& s
    )
    {
Davis King's avatar
Davis King committed
854
        return matrix_mul_scal_exp<EXP>(m.ref(),s);
855
856
    }

857
858
859
860
861
862
863
864
865
866
867
868
869
    template <
        typename EXP,
        typename S,
        bool B
        >
    inline typename disable_if<is_matrix<S>, const matrix_mul_scal_exp<EXP> >::type operator* (
        const matrix_mul_scal_exp<EXP,B>& m,
        const S& s
    )
    {
        return matrix_mul_scal_exp<EXP>(m.m,s*m.s);
    }

870
871
872
873
    template <
        typename EXP,
        typename S 
        >
Davis King's avatar
Davis King committed
874
    inline typename disable_if<is_matrix<S>, const matrix_mul_scal_exp<EXP> >::type operator* (
875
876
877
878
        const S& s,
        const matrix_exp<EXP>& m
    )
    {
Davis King's avatar
Davis King committed
879
        return matrix_mul_scal_exp<EXP>(m.ref(),s);
880
881
882
    }

    template <
883
884
885
        typename EXP,
        typename S,
        bool B
886
        >
887
888
889
    inline typename disable_if<is_matrix<S>, const matrix_mul_scal_exp<EXP> >::type operator* (
        const S& s,
        const matrix_mul_scal_exp<EXP,B>& m
890
891
    )
    {
892
        return matrix_mul_scal_exp<EXP>(m.m,s*m.s);
893
894
895
    }

    template <
896
897
        typename EXP ,
        typename S
898
        >
899
    inline const typename disable_if_c<std::numeric_limits<typename EXP::type>::is_integer, matrix_mul_scal_exp<EXP> >::type operator/ (
900
        const matrix_exp<EXP>& m,
901
        const S& s
902
903
    )
    {
904
905
906
        typedef typename EXP::type type;
        const type one = 1;
        return matrix_mul_scal_exp<EXP>(m.ref(),one/static_cast<type>(s));
907
908
909
    }

    template <
910
911
912
        typename EXP,
        bool B,
        typename S
913
        >
914
915
916
    inline const typename disable_if_c<std::numeric_limits<typename EXP::type>::is_integer, matrix_mul_scal_exp<EXP> >::type operator/ (
        const matrix_mul_scal_exp<EXP,B>& m,
        const S& s
917
918
    )
    {
919
920
        typedef typename EXP::type type;
        return matrix_mul_scal_exp<EXP>(m.m,m.s/static_cast<type>(s));
921
922
923
924
925
    }

    template <
        typename EXP
        >
Davis King's avatar
Davis King committed
926
    inline const matrix_mul_scal_exp<EXP> operator- (
927
928
929
        const matrix_exp<EXP>& m
    )
    {
Davis King's avatar
Davis King committed
930
        return matrix_mul_scal_exp<EXP>(m.ref(),-1);
931
932
    }

933
934
935
936
937
938
939
940
941
942
943
    template <
        typename EXP,
        bool B
        >
    inline const matrix_mul_scal_exp<EXP> operator- (
        const matrix_mul_scal_exp<EXP,B>& m
    )
    {
        return matrix_mul_scal_exp<EXP>(m.m,-1*m.s);
    }

944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
// ----------------------------------------------------------------------------------------

    template <
        typename EXP1,
        typename EXP2
        >
    bool operator== (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
        if (m1.nr() == m2.nr() && m1.nc() == m2.nc())
        {
            for (long r = 0; r < m1.nr(); ++r)
            {
                for (long c = 0; c < m1.nc(); ++c)
                {
                    if (m1(r,c) != m2(r,c))
                        return false;
                }
            }
            return true;
        }
        return false;
    }

    template <
        typename EXP1,
        typename EXP2
        >
    inline bool operator!= (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    ) { return !(m1 == m2); }

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

    template <
        typename T,
        long num_rows,
        long num_cols,
987
988
        typename mem_manager,
        typename layout
989
        >
990
    struct matrix_traits<matrix<T,num_rows, num_cols, mem_manager, layout> >
991
992
993
    {
        typedef T type;
        typedef mem_manager mem_manager_type;
994
        typedef layout layout_type;
995
996
        const static long NR = num_rows;
        const static long NC = num_cols;
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
        const static long cost = 1;

    };

    template <
        typename T,
        long num_rows,
        long num_cols,
        typename mem_manager,
        typename layout
        >
    class matrix : public matrix_exp<matrix<T,num_rows,num_cols, mem_manager,layout> > 
    {

        COMPILE_TIME_ASSERT(num_rows >= 0 && num_cols >= 0); 
1012

1013
1014
1015
1016
1017
1018
1019
1020
1021
    public:
        typedef typename matrix_traits<matrix>::type type;
        typedef typename matrix_traits<matrix>::mem_manager_type mem_manager_type;
        typedef typename matrix_traits<matrix>::layout_type layout_type;
        const static long NR = matrix_traits<matrix>::NR;
        const static long NC = matrix_traits<matrix>::NC;
        const static long cost = matrix_traits<matrix>::cost;

        matrix () : matrix_exp<matrix>(*this) 
1022
1023
1024
1025
1026
        {
        }

        explicit matrix (
            long length 
1027
        ) : matrix_exp<matrix>(*this) 
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
        {
            // This object you are trying to call matrix(length) on is not a column or 
            // row vector.
            COMPILE_TIME_ASSERT(NR == 1 || NC == 1);
            DLIB_ASSERT( length >= 0, 
                "\tmatrix::matrix(length)"
                << "\n\tlength must be at least 0"
                << "\n\tlength: " << length 
                << "\n\tNR:     " << NR 
                << "\n\tNC:     " << NC 
                << "\n\tthis:   " << this
                );

            if (NR == 1)
            {
                DLIB_ASSERT(NC == 0 || NC == length,
                    "\tmatrix::matrix(length)"
Davis King's avatar
Davis King committed
1045
                    << "\n\tSince this is a statically sized matrix length must equal NC"
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                data.set_size(1,length);
            }
            else
            {
                DLIB_ASSERT(NR == 0 || NR == length,
                    "\tvoid matrix::set_size(length)"
Davis King's avatar
Davis King committed
1058
                    << "\n\tSince this is a statically sized matrix length must equal NR"
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                data.set_size(length,1);
            }
        }

        matrix (
            long rows,
            long cols 
1072
        ) : matrix_exp<matrix>(*this) 
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
        {
            DLIB_ASSERT( (NR == 0 || NR == rows) && ( NC == 0 || NC == cols) && 
                    rows >= 0 && cols >= 0, 
                "\tvoid matrix::matrix(rows, cols)"
                << "\n\tYou have supplied conflicting matrix dimensions"
                << "\n\trows: " << rows
                << "\n\tcols: " << cols
                << "\n\tNR:   " << NR 
                << "\n\tNC:   " << NC 
                );
            data.set_size(rows,cols);
        }

        template <typename EXP>
        matrix (
            const matrix_exp<EXP>& m
1089
        ): matrix_exp<matrix>(*this) 
1090
        {
1091
1092
1093
1094
1095
            // You get an error on this line if the matrix m contains a type that isn't
            // the same as the type contained in the target matrix.
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true) ||
                                (is_matrix<typename EXP::type>::value == true));

1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
            // The matrix you are trying to assign m to is a statically sized matrix and 
            // m's dimensions don't match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT((NR == 0 || NR == m.nr()) && (NC == 0 || NC == m.nc()), 
                "\tmatrix& matrix::matrix(const matrix_exp& m)"
                << "\n\tYou are trying to assign a dynamically sized matrix to a statically sized matrix with the wrong size"
                << "\n\tNR:     " << NR
                << "\n\tNC:     " << NC
                << "\n\tm.nr(): " << m.nr()
                << "\n\tm.nc(): " << m.nc()
                << "\n\tthis:   " << this
                );

            data.set_size(m.nr(),m.nc());

1112
            matrix_assign(*this, m);
1113
1114
1115
1116
        }

        matrix (
            const matrix& m
1117
        ): matrix_exp<matrix>(*this) 
1118
1119
        {
            data.set_size(m.nr(),m.nc());
1120
            matrix_assign(*this, m);
1121
1122
1123
1124
1125
        }

        template <typename U, size_t len>
        matrix (
            U (&array)[len]
1126
        ): matrix_exp<matrix>(*this) 
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
        {
            COMPILE_TIME_ASSERT(NR*NC == len && len > 0);
            size_t idx = 0;
            for (long r = 0; r < NR; ++r)
            {
                for (long c = 0; c < NC; ++c)
                {
                    data(r,c) = static_cast<T>(array[idx]);
                    ++idx;
                }
            }
        }

        T& operator() (
            long r, 
            long c
        ) 
        { 
            DLIB_ASSERT(r < nr() && c < nc() &&
                   r >= 0 && c >= 0, 
                "\tT& matrix::operator(r,c)"
                << "\n\tYou must give a valid row and column"
                << "\n\tr:    " << r 
                << "\n\tc:    " << c
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return data(r,c); 
        }

        const T& operator() (
            long r, 
            long c
        ) const 
        { 
            DLIB_ASSERT(r < nr() && c < nc() &&
                   r >= 0 && c >= 0, 
                "\tconst T& matrix::operator(r,c)"
                << "\n\tYou must give a valid row and column"
                << "\n\tr:    " << r 
                << "\n\tc:    " << c
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return data(r,c);
        }

        T& operator() (
            long i
        ) 
        {
            // You can only use this operator on column vectors.
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
            DLIB_ASSERT(nc() == 1 || nr() == 1, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou can only use this operator on column or row vectors"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou must give a valid row/column number"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            return data(i);
        }

        const T& operator() (
            long i
        ) const
        {
            // You can only use this operator on column vectors.
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
            DLIB_ASSERT(nc() == 1 || nr() == 1, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou can only use this operator on column or row vectors"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou must give a valid row/column number"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            return data(i);
        }

        inline operator const type (
        ) const 
        {
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0);
            COMPILE_TIME_ASSERT(NR == 1 || NR == 0);
            DLIB_ASSERT( nr() == 1 && nc() == 1 , 
                "\tmatrix::operator const type"
Davis King's avatar
Davis King committed
1233
                << "\n\tYou can only attempt to implicit convert a matrix to a scalar if"
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
                << "\n\tthe matrix is a 1x1 matrix"
                << "\n\tnr(): " << nr() 
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return data(0);
        }

        void set_size (
            long rows,
            long cols
        )
        {
            DLIB_ASSERT( (NR == 0 || NR == rows) && ( NC == 0 || NC == cols) &&
                    rows >= 0 && cols >= 0, 
                "\tvoid matrix::set_size(rows, cols)"
                << "\n\tYou have supplied conflicting matrix dimensions"
                << "\n\trows: " << rows
                << "\n\tcols: " << cols
                << "\n\tNR:   " << NR 
                << "\n\tNC:   " << NC 
                << "\n\tthis: " << this
                );
            if (nr() != rows || nc() != cols)
                data.set_size(rows,cols);
        }

        void set_size (
            long length
        )
        {
            // This object you are trying to call set_size(length) on is not a column or 
            // row vector.
            COMPILE_TIME_ASSERT(NR == 1 || NC == 1);
            DLIB_ASSERT( length >= 0, 
                "\tvoid matrix::set_size(length)"
                << "\n\tlength must be at least 0"
                << "\n\tlength: " << length 
                << "\n\tNR:     " << NR 
                << "\n\tNC:     " << NC 
                << "\n\tthis:   " << this
                );

            if (NR == 1)
            {
                DLIB_ASSERT(NC == 0 || NC == length,
                    "\tvoid matrix::set_size(length)"
Davis King's avatar
Davis King committed
1281
                    << "\n\tSince this is a statically sized matrix length must equal NC"
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                if (nc() != length)
                    data.set_size(1,length);
            }
            else
            {
                DLIB_ASSERT(NR == 0 || NR == length,
                    "\tvoid matrix::set_size(length)"
Davis King's avatar
Davis King committed
1295
                    << "\n\tSince this is a statically sized matrix length must equal NR"
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                if (nr() != length)
                    data.set_size(length,1);
            }
        }

        long nr (
        ) const { return data.nr(); }

        long nc (
        ) const { return data.nc(); }

        long size (
        ) const { return data.nr()*data.nc(); }

        template <typename U, size_t len>
        matrix& operator= (
            U (&array)[len]
        )
        {
            COMPILE_TIME_ASSERT(NR*NC == len && len > 0);
            size_t idx = 0;
            for (long r = 0; r < NR; ++r)
            {
                for (long c = 0; c < NC; ++c)
                {
                    data(r,c) = static_cast<T>(array[idx]);
                    ++idx;
                }
            }
            return *this;
        }

        template <typename EXP>
        matrix& operator= (
            const matrix_exp<EXP>& m
        )
        {
1339
1340
1341
            // You get an error on this line if the matrix you are trying to 
            // assign m to is a statically sized matrix and  m's dimensions don't 
            // match that of *this. 
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT((NR == 0 || nr() == m.nr()) && 
                   (NC == 0 || nc() == m.nc()), 
                "\tmatrix& matrix::operator=(const matrix_exp& m)"
                << "\n\tYou are trying to assign a dynamically sized matrix to a statically sized matrix with the wrong size"
                << "\n\tnr():   " << nr()
                << "\n\tnc():   " << nc()
                << "\n\tm.nr(): " << m.nr()
                << "\n\tm.nc(): " << m.nc()
                << "\n\tthis:   " << this
                );
1354
1355
1356
1357
1358

            // You get an error on this line if the matrix m contains a type that isn't
            // the same as the type contained in the target matrix.
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true) ||
                                (is_matrix<typename EXP::type>::value == true));
1359
1360
            if (m.destructively_aliases(*this) == false)
            {
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
                // This if statement is seemingly unnecessary since set_size() contains this
                // exact same if statement.  However, structuring the code this way causes
                // gcc to handle the way it inlines this function in a much more favorable way.
                if (data.nr() == m.nr() && data.nc() == m.nc())
                {
                    matrix_assign(*this, m);
                }
                else
                {
                    set_size(m.nr(),m.nc());
                    matrix_assign(*this, m);
                }
1373
1374
1375
            }
            else
            {
1376
1377
1378
                // we have to use a temporary matrix object here because
                // *this is aliased inside the matrix_exp m somewhere.
                matrix temp;
1379
                temp.set_size(m.nr(),m.nc());
1380
                matrix_assign(temp, m);
1381
                temp.swap(*this);
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
            }
            return *this;
        }

        template <typename EXP>
        matrix& operator += (
            const matrix_exp<EXP>& m
        )
        {
            // The matrix you are trying to assign m to is a statically sized matrix and 
            // m's dimensions don't match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT(this->nr() == m.nr() && this->nc() == m.nc(), 
                "\tmatrix& matrix::operator+=(const matrix_exp& m)"
                << "\n\tYou are trying to add a dynamically sized matrix to a statically sized matrix with the wrong size"
                << "\n\tthis->nr(): " << nr()
                << "\n\tthis->nc(): " << nc()
                << "\n\tm.nr():     " << m.nr()
                << "\n\tm.nc():     " << m.nc()
                << "\n\tthis:       " << this
                );
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true));
            if (m.destructively_aliases(*this) == false)
            {
1407
                matrix_assign(*this, *this + m);
1408
1409
1410
            }
            else
            {
1411
                // we have to use a temporary matrix object here because
1412
                // this->data is aliased inside the matrix_exp m somewhere.
1413
                matrix temp;
1414
                temp.set_size(m.nr(),m.nc());
1415
                matrix_assign(temp, *this + m);
1416
                temp.swap(*this);
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
            }
            return *this;
        }


        template <typename EXP>
        matrix& operator -= (
            const matrix_exp<EXP>& m
        )
        {
            // The matrix you are trying to assign m to is a statically sized matrix and 
            // m's dimensions don't match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT(this->nr() == m.nr() && this->nc() == m.nc(), 
                "\tmatrix& matrix::operator-=(const matrix_exp& m)"
                << "\n\tYou are trying to subtract a dynamically sized matrix from a statically sized matrix with the wrong size"
                << "\n\tthis->nr(): " << nr()
                << "\n\tthis->nc(): " << nc()
                << "\n\tm.nr():     " << m.nr()
                << "\n\tm.nc():     " << m.nc()
                << "\n\tthis:       " << this
                );
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true));
            if (m.destructively_aliases(*this) == false)
            {
1443
                matrix_assign(*this, *this - m);
1444
1445
1446
            }
            else
            {
1447
                // we have to use a temporary matrix object here because
1448
                // this->data is aliased inside the matrix_exp m somewhere.
1449
                matrix temp;
1450
                temp.set_size(m.nr(),m.nc());
1451
                matrix_assign(temp, *this - m);
1452
                temp.swap(*this);
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
            }
            return *this;
        }

        matrix& operator += (
            const matrix& m
        )
        {
            const long size = m.nr()*m.nc();
            for (long i = 0; i < size; ++i)
                data(i) += m.data(i);
            return *this;
        }

        matrix& operator -= (
            const matrix& m
        )
        {
            const long size = m.nr()*m.nc();
            for (long i = 0; i < size; ++i)
                data(i) -= m.data(i);
            return *this;
        }

        matrix& operator *= (
            const T& a
        )
        {
            const long size = data.nr()*data.nc();
            for (long i = 0; i < size; ++i)
                data(i) *= a;
            return *this;
        }

        matrix& operator /= (
            const T& a
        )
        {
            const long size = data.nr()*data.nc();
            for (long i = 0; i < size; ++i)
                data(i) /= a;
            return *this;
        }

        matrix& operator= (
            const matrix& m
        )
        {
            if (this != &m)
            {
                set_size(m.nr(),m.nc());
                const long size = m.nr()*m.nc();
                for (long i = 0; i < size; ++i)
                    data(i) = m.data(i);
            }
            return *this;
        }

        void swap (
            matrix& item
        )
        {
            data.swap(item.data);
        }

1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const  { return false; }

        template <typename U, long iNR, long iNC, typename mm, typename l>
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return false; }

        bool aliases (
            const matrix& item
        ) const { return (this == &item); }

1532
    private:
1533
1534
1535
1536
1537
1538
1539
1540
        struct literal_assign_helper
        {
            /*
                This struct is a helper struct returned by the operator<<() function below.  It is
                used primarily to enable us to put DLIB_CASSERT statements on the usage of the
                operator<< form of matrix assignment.
            */

1541
1542
            literal_assign_helper(const literal_assign_helper& item) : m(item.m), r(item.r), c(item.c), has_been_used(false) {}
            literal_assign_helper(matrix* m_): m(m_), r(0), c(0),has_been_used(false) {next();}
1543
1544
            ~literal_assign_helper()
            {
1545
1546
                DLIB_CASSERT(!has_been_used || r == m->nr(),
                             "You have used the matrix comma based assignment incorrectly by failing to\n"
1547
1548
1549
1550
                             "supply a full set of values for every element of a matrix object.\n");
            }

            const literal_assign_helper& operator, (
1551
                const T& val
1552
1553
1554
            ) const
            {
                DLIB_CASSERT(r < m->nr() && c < m->nc(),
1555
                             "You have used the matrix comma based assignment incorrectly by attempting to\n" <<
1556
1557
                             "supply more values than there are elements in the matrix object being assigned to.\n\n" <<
                             "Did you forget to call set_size()?\n");
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
                (*m)(r,c) = val;
                next();
                has_been_used = true;
                return *this;
            }

        private:

            void next (
            ) const
            {
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
                ++c;
                if (c == m->nc())
                {
                    c = 0;
                    ++r;
                }
            }

            matrix* m;
            mutable long r;
            mutable long c;
1580
            mutable bool has_been_used;
1581
1582
1583
1584
        };

    public:

1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
        const literal_assign_helper operator = (
            const T& val
        ) 
        {  
            // assign the given value to every spot in this matrix
            for (long r = 0; r < nr(); ++r)
            {
                for (long c = 0; c < nc(); ++c)
                {
                    data(r,c) = val;
                }
            }

            // Now return the literal_assign_helper so that the user
            // can use the overloaded comma notation to initialize 
            // the matrix if they want to.
            return literal_assign_helper(this); 
        }
1603
1604
1605
1606

    private:


1607
        typename layout::template layout<T,NR,NC,mem_manager> data;
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
    };

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

    template <
        typename T,
        long NR,
        long NC,
1618
1619
        typename mm,
        typename l
1620
1621
        >
    void swap(
1622
1623
        matrix<T,NR,NC,mm,l>& a,
        matrix<T,NR,NC,mm,l>& b
1624
1625
1626
1627
1628
1629
    ) { a.swap(b); }

    template <
        typename T,
        long NR,
        long NC,
1630
1631
        typename mm,
        typename l
1632
1633
        >
    void serialize (
1634
        const matrix<T,NR,NC,mm,l>& item, 
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
        std::ostream& out
    )
    {
        try
        {
            serialize(item.nr(),out);
            serialize(item.nc(),out);
            for (long r = 0; r < item.nr(); ++r)
            {
                for (long c = 0; c < item.nc(); ++c)
                {
                    serialize(item(r,c),out);
                }
            }
        }
        catch (serialization_error& e)
        {
            throw serialization_error(e.info + "\n   while serializing dlib::matrix");
        }
    }

    template <
        typename T,
        long NR,
        long NC,
1660
1661
        typename mm,
        typename l
1662
1663
        >
    void deserialize (
1664
        matrix<T,NR,NC,mm,l>& item, 
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
        std::istream& in
    )
    {
        try
        {
            long nr, nc;
            deserialize(nr,in); 
            deserialize(nc,in); 

            if (NR != 0 && nr != NR)
                throw serialization_error("Error while deserializing a dlib::matrix.  Invalid rows");
            if (NC != 0 && nc != NC)
                throw serialization_error("Error while deserializing a dlib::matrix.  Invalid columns");

            item.set_size(nr,nc);
            for (long r = 0; r < nr; ++r)
            {
                for (long c = 0; c < nc; ++c)
                {
                    deserialize(item(r,c),in);
                }
            }
        }
        catch (serialization_error& e)
        {
            throw serialization_error(e.info + "\n   while deserializing a dlib::matrix");
        }
    }

    template <
        typename EXP
        >
    std::ostream& operator<< (
        std::ostream& out,
        const matrix_exp<EXP>& m
    )
    {
        using namespace std;
        const streamsize old = out.width();

        // first figure out how wide we should make each field
        string::size_type w = 0;
        ostringstream sout;
        for (long r = 0; r < m.nr(); ++r)
        {
            for (long c = 0; c < m.nc(); ++c)
            {
                sout << m(r,c); 
                w = std::max(sout.str().size(),w);
                sout.str("");
            }
        }

        // now actually print it
        for (long r = 0; r < m.nr(); ++r)
        {
            for (long c = 0; c < m.nc(); ++c)
            {
                out.width(static_cast<streamsize>(w));
                out << m(r,c) << " ";
            }
            out << "\n";
        }
        out.width(old);
        return out;
    }

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

}

#ifdef _MSC_VER
// put that warning back to its default setting
#pragma warning(default : 4355)
#endif

#endif // DLIB_MATRIx_