1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
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
|
#pragma once
#include <array>
#include <functional>
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
template <class Scalar>
using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
template <class Scalar>
using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
template <class Scalar, int p>
using Tensor = Eigen::Tensor<Scalar, p>;
template <class Scalar, int p, std::size_t... Indices>
Tensor<Scalar, p> initializeJHelper(unsigned N, std::index_sequence<Indices...>) {
std::array<unsigned, p> Ns;
std::fill_n(Ns.begin(), p, N);
return Tensor<Scalar, p>(std::get<Indices>(Ns)...);
}
template <class Scalar, int p>
Tensor<Scalar, p> initializeJ(unsigned N) {
return initializeJHelper<Scalar, p>(N, std::make_index_sequence<p>());
}
template <class Scalar, int p, std::size_t... Indices>
void setJHelper(Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) {
J(std::get<Indices>(ind)...) = val;
}
template <class Scalar, int p>
void setJ(Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) {
do {
setJHelper<Scalar, p>(J, ind, val, std::make_index_sequence<p>());
} while (std::next_permutation(ind.begin(), ind.end()));
}
template <class Scalar, int p, std::size_t... Indices>
Scalar getJHelper(const Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
return J(std::get<Indices>(ind)...);
}
template <class Scalar, int p>
Scalar getJ(const Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) {
return getJHelper<Scalar, p>(J, ind, std::make_index_sequence<p>());
}
template <class Scalar, int p>
void iterateOverHelper(Tensor<Scalar, p>& J,
std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)>& f,
unsigned l, std::array<unsigned, p> is) {
if (l == 0) {
f(J, is);
} else {
unsigned iMin;
if (l == p) {
iMin = 0;
} else {
iMin = is[p - l - 1];
}
for (unsigned i = iMin; i < J.dimension(p - 1); i++) {
std::array<unsigned, p> js = is;
js[p - l] = i;
iterateOverHelper<Scalar, p>(J, f, l - 1, js);
}
}
}
template <class Scalar, int p>
void iterateOver(Tensor<Scalar, p>& J, std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) {
std::array<unsigned, p> is;
iterateOverHelper<Scalar, p>(J, f, p, is);
}
template <class Scalar, int p, class Distribution, class Generator>
Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
Tensor<Scalar, p> J = initializeJ<Scalar, p>(N);
std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom =
[&d, &r] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
setJ<Scalar, p>(JJ, ind, d(r));
};
iterateOver<Scalar, p>(J, setRandom);
return J;
}
template <class Scalar, int p, class Distribution, class Generator>
Tensor<Scalar, p> plantState(const Tensor<Scalar, p>& J, const Vector<Scalar>& z, double β) {
Tensor<Scalar, p> JPlanted = J;
std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> plant =
[&z, β] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
Scalar Ji = getJ<Scalar, p>(JJ, ind);
Scalar zzz = 1;
for (unsigned i : ind) {
zzz *= z(i);
}
setJ<Scalar, p>(JJ, ind, Ji - β * zzz / pow(zzz.size(), p));
};
iterateOver<Scalar, p>(JPlanted, plant);
return JPlanted;
}
template <class Scalar>
Matrix<Scalar> contractDown(const Tensor<Scalar, 2>& J, const Vector<Scalar>& z) {
return Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
}
const std::array<Eigen::IndexPair<int>, 1> ip00 = {Eigen::IndexPair<int>(0, 0)};
template <class Scalar, int r>
Matrix<Scalar> contractDown(const Tensor<Scalar, r>& J, const Vector<Scalar>& z) {
Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
Tensor<Scalar, r - 1> Jz = J.contract(zT, ip00);
return contractDown(Jz, z);
}
template <int f, class Scalar>
Tensor<Scalar, f> contractDownTo(const Tensor<Scalar, f>& J, const Vector<Scalar>& z) {
return J;
}
template <int f, class Scalar, int r>
Tensor<Scalar, f> contractDownTo(const Tensor<Scalar, r>& J, const Vector<Scalar>& z) {
Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
Tensor<Scalar, r - 1> Jz = J.contract(zT, ip00);
return contractDownTo<f>(Jz, z);
}
|
