1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup glue_polyfit
20 //! @{
21
22
23
24 template<typename eT>
25 inline
26 bool
apply_noalias(Mat<eT> & out,const Col<eT> & X,const Col<eT> & Y,const uword N)27 glue_polyfit::apply_noalias(Mat<eT>& out, const Col<eT>& X, const Col<eT>& Y, const uword N)
28 {
29 arma_extra_debug_sigprint();
30
31 // create Vandermonde matrix
32
33 Mat<eT> V(X.n_elem, N+1, arma_nozeros_indicator());
34
35 V.tail_cols(1).ones();
36
37 for(uword i=1; i <= N; ++i)
38 {
39 const uword j = N-i;
40
41 Col<eT> V_col_j (V.colptr(j ), V.n_rows, false, false);
42 Col<eT> V_col_jp1(V.colptr(j+1), V.n_rows, false, false);
43
44 V_col_j = V_col_jp1 % X;
45 }
46
47 Mat<eT> Q;
48 Mat<eT> R;
49
50 const bool status1 = auxlib::qr_econ(Q, R, V);
51
52 if(status1 == false) { return false; }
53
54 const bool status2 = auxlib::solve_trimat_fast(out, R, (Q.t() * Y), uword(0));
55
56 if(status2 == false) { return false; }
57
58 return true;
59 }
60
61
62
63 template<typename T1, typename T2>
64 inline
65 bool
apply_direct(Mat<typename T1::elem_type> & out,const Base<typename T1::elem_type,T1> & X_expr,const Base<typename T1::elem_type,T2> & Y_expr,const uword N)66 glue_polyfit::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& X_expr, const Base<typename T1::elem_type, T2>& Y_expr, const uword N)
67 {
68 arma_extra_debug_sigprint();
69
70 typedef typename T1::elem_type eT;
71
72 const quasi_unwrap<T1> UX(X_expr.get_ref());
73 const quasi_unwrap<T2> UY(Y_expr.get_ref());
74
75 const Mat<eT>& X = UX.M;
76 const Mat<eT>& Y = UY.M;
77
78 arma_debug_check
79 (
80 ( ((X.is_vec() == false) && (X.is_empty() == false)) || ((Y.is_vec() == false) && (Y.is_empty() == false)) ),
81 "polyfit(): given object must be a vector"
82 );
83
84 arma_debug_check( (X.n_elem != Y.n_elem), "polyfit(): given vectors must have the same number of elements" );
85
86 if(X.n_elem == 0)
87 {
88 out.reset();
89 return true;
90 }
91
92 arma_debug_check( (N >= X.n_elem), "polyfit(): N must be less than the number of elements in X" );
93
94 const Col<eT> X_as_colvec( const_cast<eT*>(X.memptr()), X.n_elem, false, false);
95 const Col<eT> Y_as_colvec( const_cast<eT*>(Y.memptr()), Y.n_elem, false, false);
96
97 bool status = false;
98
99 if(UX.is_alias(out) || UY.is_alias(out))
100 {
101 Mat<eT> tmp;
102 status = glue_polyfit::apply_noalias(tmp, X_as_colvec, Y_as_colvec, N);
103 out.steal_mem(tmp);
104 }
105 else
106 {
107 status = glue_polyfit::apply_noalias(out, X_as_colvec, Y_as_colvec, N);
108 }
109
110 return status;
111 }
112
113
114
115 template<typename T1, typename T2>
116 inline
117 void
apply(Mat<typename T1::elem_type> & out,const Glue<T1,T2,glue_polyfit> & expr)118 glue_polyfit::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyfit>& expr)
119 {
120 arma_extra_debug_sigprint();
121
122 const bool status = glue_polyfit::apply_direct(out, expr.A, expr.B, expr.aux_uword);
123
124 if(status == false)
125 {
126 out.soft_reset();
127 arma_stop_runtime_error("polyfit(): failed");
128 }
129 }
130
131
132
133 //! @}
134