Skip to content
This repository was archived by the owner on Aug 6, 2023. It is now read-only.

Commit 024808c

Browse files
luk036luk036
committed
remove xtensor dependency
1 parent 283335b commit 024808c

File tree

14 files changed

+559
-244
lines changed

14 files changed

+559
-244
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ option(INSTALL_ONLY "Enable for installation only" OFF)
1212
# Note: update this to your new project's name and version
1313
project(
1414
LmiSolver
15-
VERSION 1.3.5
15+
VERSION 1.3.6
1616
LANGUAGES CXX
1717
)
1818

‎bench/BM_lmi.cpp‎

Lines changed: 149 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,12 @@
11
#include <ellalgo/cutting_plane.hpp> // for cutting_plane_optim
22
#include <ellalgo/ell.hpp> // for Ell
3+
#include <ellalgo/ell_matrix.hpp> // for Matrix
34
#include <gsl/span> // for span
45
#include <lmisolver/lmi_old_oracle.hpp> // for LmiOldOracle
56
#include <lmisolver/lmi_oracle.hpp> // for LmiOracle
67
#include <tuple> // for tuple
78
#include <type_traits> // for move
89
#include <vector> // for vector
9-
#include <xtensor/xarray.hpp> // for xarray_container
10-
#include <xtensor/xlayout.hpp> // for layout_type, layout_type::row...
11-
#include <xtensor/xmath.hpp> // for sum
12-
#include <xtensor/xoperation.hpp> // for xfunction_type_t, operator*
13-
#include <xtensor/xreducer.hpp> // for xreducer
14-
#include <xtensor/xtensor_forward.hpp> // for xarray
1510

1611
#include "benchmark/benchmark.h" // for BENCHMARK, State, BENCHMARK_...
1712

@@ -21,14 +16,13 @@
2116
* @tparam Oracle
2217
*/
2318
template <typename Oracle> class my_oracle {
24-
using Arr = xt::xarray<double, xt::layout_type::row_major>;
2519
using Vec = std::valarray<double>;
26-
using Cut = std::pair<Arr, double>;
20+
using Cut = std::pair<Vec, double>;
2721

2822
private:
2923
Oracle lmi1;
3024
Oracle lmi2;
31-
const Arr c;
25+
const Vec c;
3226

3327
public:
3428
/**
@@ -40,9 +34,9 @@ template <typename Oracle> class my_oracle {
4034
* @param[in] B2
4135
* @param[in] c
4236
*/
43-
my_oracle(gsl::span<const Arr> F1, const Arr &B1, gsl::span<const Arr> F2,
44-
const Arr &B2, Arr c)
45-
: lmi1{F1, B1}, lmi2{F2, B2}, c{std::move(c)} {}
37+
my_oracle(size_t m1, gsl::span<const Matrix> F1, const Matrix &B1, size_t m2,
38+
gsl::span<const Matrix> F2, const Matrix &B2, Vec c)
39+
: lmi1{m1, F1, B1}, lmi2{m2, F2, B2}, c{std::move(c)} {}
4640

4741
/**
4842
* @brief
@@ -51,10 +45,10 @@ template <typename Oracle> class my_oracle {
5145
* @param[in,out] t
5246
* @return std::tuple<Cut, double>
5347
*/
54-
std::tuple<Cut, bool> assess_optim(const Arr &x, double &t) {
55-
const auto f0 = xt::sum(this->c * x)();
48+
std::tuple<Cut, bool> assess_optim(const Vec &x, double &t) {
49+
const auto f0 = (this->c * x).sum();
5650
const auto f1 = f0 - t;
57-
if (f1 > 0) {
51+
if (f1 > 0.0) {
5852
return {{this->c, f1}, false};
5953
}
6054
const auto cut1 = this->lmi1(x);
@@ -76,7 +70,7 @@ template <typename Oracle> class my_oracle {
7670
* @param[in,out] t
7771
* @return std::tuple<Cut, double>
7872
*/
79-
std::tuple<Cut, bool> operator()(const Arr &x, double &t) {
73+
std::tuple<Cut, bool> operator()(const Vec &x, double &t) {
8074
return this->assess_optim(x, t);
8175
}
8276
};
@@ -87,22 +81,55 @@ template <typename Oracle> class my_oracle {
8781
* @param[in,out] state
8882
*/
8983
static void LMI_Lazy(benchmark::State &state) {
90-
using Arr = xt::xarray<double, xt::layout_type::row_major>;
91-
92-
// auto c = Arr {1.0, -1.0, 1.0};
93-
const auto F1 = std::vector<Arr>{{{-7.0, -11.0}, {-11.0, 3.0}},
94-
{{7.0, -18.0}, {-18.0, 8.0}},
95-
{{-2.0, -8.0}, {-8.0, 1.0}}};
96-
const auto B1 = Arr{{33.0, -9.0}, {-9.0, 26.0}};
97-
const auto F2 = std::vector<Arr>{
98-
{{-21.0, -11.0, 0.0}, {-11.0, 10.0, 8.0}, {0.0, 8.0, 5.0}},
99-
{{0.0, 10.0, 16.0}, {10.0, -10.0, -10.0}, {16.0, -10.0, 3.0}},
100-
{{-5.0, 2.0, -17.0}, {2.0, -6.0, 8.0}, {-17.0, 8.0, 6.0}}};
101-
const auto B2 = Arr{{14.0, 9.0, 40.0}, {9.0, 91.0, 10.0}, {40.0, 10.0, 15.0}};
84+
using Vec = std::valarray<double>;
85+
using M_t = std::vector<Matrix>;
86+
87+
auto c = Vec{1.0, -1.0, 1.0};
88+
89+
auto m0F1 = Matrix(2);
90+
m0F1.row(0) = Vec{-7.0, -11.0};
91+
m0F1.row(1) = Vec{-11.0, 3.0};
92+
93+
auto m1F1 = Matrix(2);
94+
m1F1.row(0) = Vec{7.0, -18.0};
95+
m1F1.row(1) = Vec{-18.0, 8.0};
96+
97+
auto m2F1 = Matrix(2);
98+
m2F1.row(0) = Vec{-2.0, -8.0};
99+
m2F1.row(1) = Vec{-8.0, 1.0};
100+
101+
auto F1 = M_t{m0F1, m1F1, m2F1};
102+
103+
auto B1 = Matrix(2);
104+
B1.row(0) = Vec{33.0, -9.0};
105+
B1.row(1) = Vec{-9.0, 26.0};
106+
107+
auto m0F2 = Matrix(3);
108+
m0F2.row(0) = Vec{-21.0, -11.0, 0.0};
109+
m0F2.row(1) = Vec{-11.0, 10.0, 8.0};
110+
m0F2.row(2) = Vec{0.0, 8.0, 5.0};
111+
112+
auto m1F2 = Matrix(3);
113+
m1F2.row(0) = Vec{0.0, 10.0, 16.0};
114+
m1F2.row(1) = Vec{10.0, -10.0, -10.0};
115+
m1F2.row(2) = Vec{16.0, -10.0, 3.0};
116+
117+
auto m2F2 = Matrix(3);
118+
m2F2.row(0) = Vec{-5.0, 2.0, -17.0};
119+
m2F2.row(1) = Vec{2.0, -6.0, 8.0};
120+
m2F2.row(2) = Vec{-17.0, 8.0, 6.0};
121+
122+
auto F2 = M_t{m0F2, m1F2, m2F2};
123+
124+
auto B2 = Matrix(3);
125+
B2.row(0) = Vec{14.0, 9.0, 40.0};
126+
B2.row(1) = Vec{9.0, 91.0, 10.0};
127+
B2.row(2) = Vec{40.0, 10.0, 15.0};
102128

103129
while (state.KeepRunning()) {
104-
auto P = my_oracle<LmiOracle<Arr>>(F1, B1, F2, B2, Arr{1.0, -1.0, 1.0});
105-
auto E = Ell<Arr>(10.0, Arr{0.0, 0.0, 0.0});
130+
auto P =
131+
my_oracle<LmiOracle<Matrix>>(2, F1, B1, 3, F2, B2, Vec{1.0, -1.0, 1.0});
132+
auto E = Ell<Vec>(10.0, Vec{0.0, 0.0, 0.0});
106133
auto t = 1e100; // std::numeric_limits<double>::max()
107134
[[maybe_unused]] const auto rslt = cutting_plane_optim(P, E, t);
108135
}
@@ -119,22 +146,55 @@ BENCHMARK(LMI_Lazy);
119146
* @param[in,out] state
120147
*/
121148
static void LMI_old(benchmark::State &state) {
122-
using Arr = xt::xarray<double, xt::layout_type::row_major>;
123-
124-
// auto c = Arr {1.0, -1.0, 1.0};
125-
const auto F1 = std::vector<Arr>{{{-7.0, -11.0}, {-11.0, 3.0}},
126-
{{7.0, -18.0}, {-18.0, 8.0}},
127-
{{-2.0, -8.0}, {-8.0, 1.0}}};
128-
const auto B1 = Arr{{33.0, -9.0}, {-9.0, 26.0}};
129-
const auto F2 = std::vector<Arr>{
130-
{{-21.0, -11.0, 0.0}, {-11.0, 10.0, 8.0}, {0.0, 8.0, 5.0}},
131-
{{0.0, 10.0, 16.0}, {10.0, -10.0, -10.0}, {16.0, -10.0, 3.0}},
132-
{{-5.0, 2.0, -17.0}, {2.0, -6.0, 8.0}, {-17.0, 8.0, 6.0}}};
133-
const auto B2 = Arr{{14.0, 9.0, 40.0}, {9.0, 91.0, 10.0}, {40.0, 10.0, 15.0}};
149+
using Vec = std::valarray<double>;
150+
using M_t = std::vector<Matrix>;
151+
152+
auto c = Vec{1.0, -1.0, 1.0};
153+
154+
auto m0F1 = Matrix(2);
155+
m0F1.row(0) = Vec{-7.0, -11.0};
156+
m0F1.row(1) = Vec{-11.0, 3.0};
157+
158+
auto m1F1 = Matrix(2);
159+
m1F1.row(0) = Vec{7.0, -18.0};
160+
m1F1.row(1) = Vec{-18.0, 8.0};
161+
162+
auto m2F1 = Matrix(2);
163+
m2F1.row(0) = Vec{-2.0, -8.0};
164+
m2F1.row(1) = Vec{-8.0, 1.0};
165+
166+
auto F1 = M_t{m0F1, m1F1, m2F1};
167+
168+
auto B1 = Matrix(2);
169+
B1.row(0) = Vec{33.0, -9.0};
170+
B1.row(1) = Vec{-9.0, 26.0};
171+
172+
auto m0F2 = Matrix(3);
173+
m0F2.row(0) = Vec{-21.0, -11.0, 0.0};
174+
m0F2.row(1) = Vec{-11.0, 10.0, 8.0};
175+
m0F2.row(2) = Vec{0.0, 8.0, 5.0};
176+
177+
auto m1F2 = Matrix(3);
178+
m1F2.row(0) = Vec{0.0, 10.0, 16.0};
179+
m1F2.row(1) = Vec{10.0, -10.0, -10.0};
180+
m1F2.row(2) = Vec{16.0, -10.0, 3.0};
181+
182+
auto m2F2 = Matrix(3);
183+
m2F2.row(0) = Vec{-5.0, 2.0, -17.0};
184+
m2F2.row(1) = Vec{2.0, -6.0, 8.0};
185+
m2F2.row(2) = Vec{-17.0, 8.0, 6.0};
186+
187+
auto F2 = M_t{m0F2, m1F2, m2F2};
188+
189+
auto B2 = Matrix(3);
190+
B2.row(0) = Vec{14.0, 9.0, 40.0};
191+
B2.row(1) = Vec{9.0, 91.0, 10.0};
192+
B2.row(2) = Vec{40.0, 10.0, 15.0};
134193

135194
while (state.KeepRunning()) {
136-
auto P = my_oracle<LmiOldOracle<Arr>>(F1, B1, F2, B2, Arr{1.0, -1.0, 1.0});
137-
auto E = Ell<Arr>(10.0, Arr{0.0, 0.0, 0.0});
195+
auto P = my_oracle<LmiOldOracle<Matrix>>(2, F1, B1, 3, F2, B2,
196+
Vec{1.0, -1.0, 1.0});
197+
auto E = Ell<Vec>(10.0, Vec{0.0, 0.0, 0.0});
138198
auto t = 1e100; // std::numeric_limits<double>::max()
139199
[[maybe_unused]] const auto rslt = cutting_plane_optim(P, E, t);
140200
}
@@ -147,23 +207,54 @@ BENCHMARK(LMI_old);
147207
// * @param[in,out] state
148208
// */
149209
// static void LMI_No_Trick(benchmark::State &state) {
150-
// using Arr = xt::xarray<double, xt::layout_type::row_major>;
210+
// using Vec = std::valarray<double>;
211+
// using M_t = std::vector<Matrix>;
212+
//
213+
// auto c = Vec{1.0, -1.0, 1.0};
214+
//
215+
// auto m0F1 = Matrix(2);
216+
// m0F1.row(0) = Vec{-7.0, -11.0};
217+
// m0F1.row(1) = Vec{-11.0, 3.0};
218+
//
219+
// auto m1F1 = Matrix(2);
220+
// m1F1.row(0) = Vec{7.0, -18.0};
221+
// m1F1.row(1) = Vec{-18.0, 8.0};
222+
//
223+
// auto m2F1 = Matrix(2);
224+
// m2F1.row(0) = Vec{-2.0, -8.0};
225+
// m2F1.row(1) = Vec{-8.0, 1.0};
226+
//
227+
// auto F1 = M_t{m0F1, m1F1, m2F1};
228+
//
229+
// auto B1 = Matrix(2);
230+
// B1.row(0) = Vec{33.0, -9.0};
231+
// B1.row(1) = Vec{-9.0, 26.0};
232+
//
233+
// auto m0F2 = Matrix(3);
234+
// m0F2.row(0) = Vec{-21.0, -11.0, 0.0};
235+
// m0F2.row(1) = Vec{-11.0, 10.0, 8.0};
236+
// m0F2.row(2) = Vec{0.0, 8.0, 5.0};
237+
//
238+
// auto m1F2 = Matrix(3);
239+
// m1F2.row(0) = Vec{0.0, 10.0, 16.0};
240+
// m1F2.row(1) = Vec{10.0, -10.0, -10.0};
241+
// m1F2.row(2) = Vec{16.0, -10.0, 3.0};
242+
//
243+
// auto m2F2 = Matrix(3);
244+
// m2F2.row(0) = Vec{-5.0, 2.0, -17.0};
245+
// m2F2.row(1) = Vec{2.0, -6.0, 8.0};
246+
// m2F2.row(2) = Vec{-17.0, 8.0, 6.0};
247+
//
248+
// auto F2 = M_t{m0F2, m1F2, m2F2};
151249
//
152-
// // const auto c = Arr {1.0, -1.0, 1.0};
153-
// const auto F1 = std::vector<Arr>{{{-7.0, -11.0}, {-11.0, 3.0}},
154-
// {{7.0, -18.0}, {-18.0, 8.0}},
155-
// {{-2.0, -8.0}, {-8.0, 1.0}}};
156-
// const auto B1 = Arr{{33.0, -9.0}, {-9.0, 26.0}};
157-
// const auto F2 = std::vector<Arr>{
158-
// {{-21.0, -11.0, 0.0}, {-11.0, 10.0, 8.0}, {0.0, 8.0, 5.0}},
159-
// {{0.0, 10.0, 16.0}, {10.0, -10.0, -10.0}, {16.0, -10.0, 3.0}},
160-
// {{-5.0, 2.0, -17.0}, {2.0, -6.0, 8.0}, {-17.0, 8.0, 6.0}}};
161-
// const auto B2 = Arr{{14.0, 9.0, 40.0}, {9.0, 91.0, 10.0},
162-
// {40.0, 10.0, 15.0}};
250+
// auto B2 = Matrix(3);
251+
// B2.row(0) = Vec{14.0, 9.0, 40.0};
252+
// B2.row(1) = Vec{9.0, 91.0, 10.0};
253+
// B2.row(2) = Vec{40.0, 10.0, 15.0};
163254
//
164255
// while (state.KeepRunning()) {
165-
// auto P = my_oracle<LmiOracle<Arr>>(F1, B1, F2, B2, Arr{1.0, -1.0, 1.0});
166-
// auto E = Ell<Arr>(10.0, Arr{0.0, 0.0, 0.0});
256+
// auto P = my_oracle<LmiOracle<Matrix>>(2, F1, B1, 3, F2, B2, Vec{1.0,
257+
// -1.0, 1.0}); auto E = Ell<Vec>(10.0, Vec{0.0, 0.0, 0.0});
167258
// E.no_defer_trick = true;
168259
// auto t = 1e100; // std::numeric_limits<double>::max()
169260
// [[maybe_unused]] const auto rslt = cutting_plane_optim(P, E, t);

‎include/lmisolver/ldlt_ext.hpp‎

Lines changed: 22 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
* for all v in R^n.
1616
* - O(p^2) per iteration, independent of N
1717
*/
18-
template <typename Arr036> class ldlt_ext {
18+
class ldlt_ext {
1919
using Vec = std::valarray<double>;
2020
using Rng = std::pair<size_t, size_t>;
2121

@@ -132,9 +132,24 @@ template <typename Arr036> class ldlt_ext {
132132
*
133133
* @return auto
134134
*/
135-
auto witness() -> double;
135+
auto witness() -> double {
136+
assert(!this->is_spd());
136137

137-
auto set_witness_vec(Arr036 &v) const -> void {
138+
// const auto& [start, n] = this->p;
139+
const auto &start = this->p.first;
140+
const auto &n = this->p.second;
141+
auto m = n - 1; // assume stop > 0
142+
this->witness_vec[m] = 1.0;
143+
for (auto i = m; i > start; --i) {
144+
this->witness_vec[i - 1] = 0.0;
145+
for (auto k = i; k != n; ++k) {
146+
this->witness_vec[i - 1] -= this->T(k, i - 1) * this->witness_vec[k];
147+
}
148+
}
149+
return -this->T(m, m);
150+
}
151+
152+
template <typename Arr036> auto set_witness_vec(Arr036 &v) const -> void {
138153
for (auto i = 0U; i != this->_n; ++i) {
139154
v[i] = this->witness_vec[i];
140155
}
@@ -153,18 +168,19 @@ template <typename Arr036> class ldlt_ext {
153168
const auto &start = this->p.first;
154169
const auto &stop = this->p.second;
155170
for (auto i = start; i != stop; ++i) {
156-
auto s = double{};
171+
auto s = 0.0;
157172
for (auto j = i + 1; j != stop; ++j) {
158173
s += A(i, j) * v[j];
159174
}
160-
res += v[i] * (A(i, i) * v[i] + 2 * s);
175+
res += v[i] * (A(i, i) * v[i] + 2.0 * s);
161176
}
162177
return res;
163178
}
164179

165180
/**
166181
* @brief Return upper triangular matrix $R$ where $A = R^T R$
167182
*
183+
* Note: must input a zero matrix
168184
* @return typename ldlt_ext<Arr036>::Mat
169185
*/
170186
template <typename Mat> auto sqrt(Mat &M) -> void {
@@ -174,6 +190,7 @@ template <typename Arr036> class ldlt_ext {
174190
M(i, i) = std::sqrt(this->T(i, i));
175191
for (auto j = i + 1; j != this->_n; ++j) {
176192
M(i, j) = this->T(j, i) * M(i, i);
193+
M(j, i) = 0.0;
177194
}
178195
}
179196
}

0 commit comments

Comments
 (0)