codac 1.5.6
Loading...
Searching...
No Matches
codac2_pow.h
Go to the documentation of this file.
1
9
10#pragma once
11
12#include "codac2_Interval.h"
13#include "codac2_AnalyticType.h"
15#include "codac2_sqr.h"
16
17namespace codac2
18{
19 struct PowOp
20 {
21 template<typename X1,typename X2>
22 static std::string str(const X1& x1, const X2& x2)
23 {
24 return x1->str(!x1->is_str_leaf()) + "^" + x2->str(!x2->is_str_leaf());
25 }
26
27 template<typename X1, typename X2>
28 static std::pair<Index,Index> output_shape([[maybe_unused]] const X1& s1, [[maybe_unused]] const X2& s2)
29 {
30 return {1,1};
31 }
32
33 static Interval fwd(const Interval& x1, const Interval& x2);
34 static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2);
35 static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2);
36 static void bwd(const Interval& y, Interval& x1, int x2);
37 static void bwd(const Interval& y, Interval& x1, Interval& x2);
38 };
39
40 // Analytic operator
41 // The following function can be used to build analytic expressions.
42
43 inline ScalarExpr
44 pow(const ScalarExpr& x1, const ScalarExpr& x2)
45 {
46 auto x2_as_const_value = std::dynamic_pointer_cast<ConstValueExpr<ScalarType>>(x2->copy());
47 if(x2_as_const_value && x2_as_const_value->value() == 2.)
48 {
49 // SqrOp may be better than the generic PowOp with a simple power of 2
50 return { std::make_shared<AnalyticOperationExpr<SqrOp,ScalarType,ScalarType>>(x1) };
51 }
52
53 return { std::make_shared<AnalyticOperationExpr<PowOp,ScalarType,ScalarType,ScalarType>>(x1,x2) };
54 }
55
56 inline ScalarExpr
57 operator^(const ScalarExpr& x1, const ScalarExpr& x2)
58 {
59 return pow(x1,x2);
60 }
61
62 // Inline functions
63
64 inline Interval PowOp::fwd(const Interval& x1, const Interval& x2)
65 {
66 return pow(x1,x2);
67 }
68
69 inline ScalarType PowOp::fwd_natural(const ScalarType& x1, const ScalarType& x2)
70 {
71 bool x2isint = x2.a.is_integer();
72 bool x2positive = x2.a.lb()>=0.0;
73 return {
74 fwd(x1.a, x2.a),
75 x1.def_domain && x2.def_domain
76 && (x2isint || x1.a.lb()>=0.0)
77 && (x2positive || !x1.a.contains(0.0))
78 };
79 }
80
81 inline ScalarType PowOp::fwd_centered(const ScalarType& x1, const ScalarType& x2)
82 {
83 if(centered_form_not_available_for_args(x1,x2))
84 return fwd_natural(x1,x2);
85
86 bool x2isint = x2.a.is_integer();
87 bool x2positive = x2.a.lb()>=0.0;
88 IntervalMatrix d(1,x1.da.size());
89
90 if (x2.a.is_degenerated())
91 {
92 // To avoid calling log(x1.a), as it would return emptyset if x1.a<=0
93 for(Index i = 0 ; i < d.size() ; i++)
94 d(0,i) = x2.a*x1.da(0,i)*pow(x1.a,x2.a-1.);
95 }
96
97 else
98 {
99 for(Index i = 0 ; i < d.size() ; i++)
100 {
101 d(0,i) = x2.a*x1.da(0,i)*pow(x1.a,x2.a-1.)
102 + x2.da(0,i)*log(x1.a)*pow(x1.a,x2.a);
103 // Not good when x1 close to 0
104 }
105 }
106
107 return {
108 fwd(x1.m, x2.m),
109 fwd(x1.a, x2.a),
110 d,
111 x1.def_domain && x2.def_domain
112 && (x2isint || x1.a.lb()>=0.0)
113 && (x2positive || !x1.a.contains(0.0))
114 };
115 }
116
117 inline void PowOp::bwd(const Interval& y, Interval& x1, int x2)
118 {
119 // The content of this function comes from the IBEX library.
120 // See ibex::Interval (IBEX lib, main author: Gilles Chabert)
121 // https://ibex-lib.readthedocs.io
122
123 if(x2 % 2 == 0)
124 {
125 Interval proj = root(y, x2);
126 Interval pos_proj = proj & x1;
127 Interval neg_proj = (-proj) & x1;
128 x1 = pos_proj | neg_proj;
129 }
130
131 else
132 x1 &= root(y, x2);
133 }
134
135 inline void PowOp::bwd(const Interval& y, Interval& x1, Interval& x2)
136 {
137 assert_release(x2.is_degenerated() && "PowOp::bwd(y,x1,x2) (with x1 and x2 intervals) not implemented yet with Gaol");
138 assert_release(x2 == (int)(x2.mid()) && "PowOp::bwd(y,x1,x2) (x2 not integer) not implemented yet with Gaol");
139 PowOp::bwd(y, x1, x2.mid());
140 }
141}
Interval class, for representing closed and connected subsets of .
Definition codac2_Interval.h:62
Interval pow(const Interval &x, int n)
Returns , .
Definition codac2_Interval_operations_impl.h:33
Interval log(const Interval &x)
Returns .
Definition codac2_Interval_operations_impl.h:85
Interval root(const Interval &x, int p)
Returns the p-th root: .
Definition codac2_Interval_operations_impl.h:60