codac 1.5.6
Loading...
Searching...
No Matches
codac2_det.h
Go to the documentation of this file.
1
9
10#pragma once
11
12#include "codac2_Interval.h"
13#include "codac2_IntervalVector.h"
14#include "codac2_IntervalMatrix.h"
15#include "codac2_AnalyticType.h"
16#include "codac2_AnalyticExprWrapper.h"
17#include "codac2_arith_sub.h"
18#include "codac2_arith_mul.h"
19
20namespace codac2
21{
22 struct DetOp
23 {
24 // For matrices
25 static Interval fwd(const IntervalMatrix& x);
26 static ScalarType fwd_natural(const MatrixType& x);
27 static ScalarType fwd_centered(const MatrixType& x);
28 static void bwd(const Interval& y, IntervalMatrix& x);
29
30 // For two vectors (merged into a 2×2 matrix)
31 static Interval fwd(const IntervalVector& x1, const IntervalVector& x2);
32 static ScalarType fwd_natural(const VectorType& x1, const VectorType& x2);
33 static ScalarType fwd_centered(const VectorType& x1, const VectorType& x2);
34 static void bwd(const Interval& y, IntervalVector& x1, IntervalVector& x2);
35
36 // For three vectors (merged into a 3×3 matrix)
37 static Interval fwd(const IntervalVector& x1, const IntervalVector& x2, const IntervalVector& x3);
38 static ScalarType fwd_natural(const VectorType& x1, const VectorType& x2, const VectorType& x3);
39 static ScalarType fwd_centered(const VectorType& x1, const VectorType& x2, const VectorType& x3);
40 static void bwd(const Interval& y, IntervalVector& x1, IntervalVector& x2, IntervalVector& x3);
41 };
42
43 // Analytic operator
44 // The following functions can be used to build analytic expressions.
45
46 inline ScalarExpr
47 det(const MatrixExpr& x1)
48 {
49 return { std::make_shared<AnalyticOperationExpr<DetOp,ScalarType,MatrixType>>(x1) };
50 }
51
52 inline ScalarExpr
53 det(const VectorExpr& x1, const VectorExpr& x2)
54 {
55 return { std::make_shared<AnalyticOperationExpr<DetOp,ScalarType,VectorType,VectorType>>(x1,x2) };
56 }
57
58 inline ScalarExpr
59 det(const VectorExpr& x1, const VectorExpr& x2, const VectorExpr& x3)
60 {
61 return { std::make_shared<AnalyticOperationExpr<DetOp,ScalarType,VectorType,VectorType,VectorType>>(x1,x2,x3) };
62 }
63
64 // Inline functions
65
66 inline Interval DetOp::fwd(const IntervalMatrix& x)
67 {
68 assert_release(x.is_squared() && "can only compute determinants for a square matrix");
69 assert_release((x.rows() == 1 || x.rows() == 2) && "determinant not yet computable for n×n matrices, n>2");
70
71 if(x.rows() == 1) // 1×1 matrix
72 return x(0,0);
73
74 else if(x.rows() == 2) // 2×2 matrix
75 return x(0,0)*x(1,1)-x(0,1)*x(1,0);
76
77 else
78 return Interval::empty(); // unhandled case
79 }
80
81 inline ScalarType DetOp::fwd_natural(const MatrixType& x)
82 {
83 return {
84 fwd(x.a),
85 x.def_domain
86 };
87 }
88
89 inline ScalarType DetOp::fwd_centered(const MatrixType& x)
90 {
91 return {
92 fwd(x.m),
93 fwd(x.a),
94 IntervalMatrix(0,0), // not supported yet for auto diff
95 x.def_domain
96 };
97 }
98
99 inline void DetOp::bwd(const Interval& y, IntervalMatrix& x)
100 {
101 assert_release(x.is_squared() && "can only compute determinants for a square matrix");
102 assert_release((x.rows() == 1 || x.rows() == 2) && "determinant not yet computable for n×n matrices, n>2");
103
104 if(x.rows() == 1) // 1×1 matrix
105 x(0,0) &= y;
106
107 else if(x.rows() == 2) // 2×2 matrix
108 {
109 Interval z1 = x(0,0)*x(1,1), z2 = x(1,0)*x(0,1);
110 SubOp::bwd(y, z1, z2);
111 MulOp::bwd(z1, x(0,0), x(1,1));
112 MulOp::bwd(z2, x(1,0), x(0,1));
113 }
114
115 else
116 {
117 // unhandled case
118 }
119 }
120
121 inline Interval DetOp::fwd(const IntervalVector& x1, const IntervalVector& x2)
122 {
123 assert_release(x1.size() == 2 && x2.size() == 2 && "determinant only computable for pairs of 2d vectors");
124 IntervalMatrix m(2,2);
125 m.col(0) = x1; m.col(1) = x2;
126 return DetOp::fwd(m);
127 }
128
129 inline ScalarType DetOp::fwd_natural(const VectorType& x1, const VectorType& x2)
130 {
131 IntervalMatrix a(2,2);
132 a.col(0) = x1.a; a.col(1) = x2.a;
133
134 return {
135 fwd(a),
136 x1.def_domain && x2.def_domain
137 };
138 }
139
140 inline ScalarType DetOp::fwd_centered(const VectorType& x1, const VectorType& x2)
141 {
142 if(centered_form_not_available_for_args(x1,x2))
143 return fwd_natural(x1,x2);
144
145 IntervalMatrix m(2,2);
146 m.col(0) = x1.m; m.col(1) = x2.m;
147 IntervalMatrix a(2,2);
148 a.col(0) = x1.a; a.col(1) = x2.a;
149
150 return {
151 fwd(m),
152 fwd(a),
153 IntervalMatrix(0,0), // not supported yet for auto diff
154 x1.def_domain && x2.def_domain
155 };
156 }
157
158 inline void DetOp::bwd(const Interval& y, IntervalVector& x1, IntervalVector& x2)
159 {
160 assert_release(x1.size() == 2 && x2.size() == 2 && "determinant only computable for pairs of 2d vectors");
161 IntervalMatrix m(2,2);
162 m.col(0) = x1; m.col(1) = x2;
163 DetOp::bwd(y,m);
164 x1 &= m.col(0);
165 x2 &= m.col(1);
166 }
167
168 inline Interval DetOp::fwd(const IntervalVector& x1, const IntervalVector& x2, const IntervalVector& x3)
169 {
170 assert_release(x1.size() == 3 && x2.size() == 3 && x3.size() == 3 && "determinant only computable for triplet of 3d vectors");
171 IntervalMatrix m(3,3);
172 m.col(0) = x1; m.col(1) = x2; m.col(2) = x3;
173 return DetOp::fwd(m);
174 }
175
176 inline ScalarType DetOp::fwd_natural(const VectorType& x1, const VectorType& x2, const VectorType& x3)
177 {
178 IntervalMatrix a(3,3);
179 a.col(0) = x1.a; a.col(1) = x2.a; a.col(2) = x3.a;
180
181 return {
182 fwd(a),
183 x1.def_domain && x2.def_domain && x3.def_domain
184 };
185 }
186
187 inline ScalarType DetOp::fwd_centered(const VectorType& x1, const VectorType& x2, const VectorType& x3)
188 {
189 if(centered_form_not_available_for_args(x1,x2,x3))
190 return fwd_natural(x1,x2,x3);
191
192 IntervalMatrix m(3,3);
193 m.col(0) = x1.m; m.col(1) = x2.m; m.col(2) = x3.m;
194 IntervalMatrix a(3,3);
195 a.col(0) = x1.a; a.col(1) = x2.a; a.col(2) = x3.a;
196
197 return {
198 fwd(m),
199 fwd(a),
200 IntervalMatrix(0,0), // not supported yet for auto diff
201 x1.def_domain && x2.def_domain && x3.def_domain
202 };
203 }
204
205 inline void DetOp::bwd(const Interval& y, IntervalVector& x1, IntervalVector& x2, IntervalVector& x3)
206 {
207 assert_release(x1.size() == 3 && x2.size() == 3 && x3.size() == 3 && "determinant only computable for triplet of 3d vectors");
208
209 IntervalMatrix m(3,3);
210 m.col(0) = x1; m.col(1) = x2; m.col(2) = x3;
211
212 DetOp::bwd(y,m);
213 x1 &= m.col(0);
214 x2 &= m.col(1);
215 x3 &= m.col(2);
216 }
217}
codac2::Interval
Interval class, for representing closed and connected subsets of .
Definition codac2_Interval.h:62
codac2::Interval::empty
static Interval empty()
Provides an empty interval.
Definition codac2_Interval_impl.h:535
codac2_AnalyticExprWrapper.h
codac2_AnalyticType.h
codac2_Interval.h
codac2_IntervalMatrix.h
codac2_IntervalVector.h
codac2_arith_mul.h
codac2_arith_sub.h
Generated by doxygen 1.13.0