Lattice.h
Go to the documentation of this file.
1//===========================================================================
2/*!
3 *
4 * \brief Various functions for generating n-dimensional grids
5 * (simplex-lattices).
6 *
7 *
8 * \author Bjørn Bugge Grathwohl
9 * \date February 2017
10 *
11 * \par Copyright 1995-2017 Shark Development Team
12 *
13 * <BR><HR>
14 * This file is part of Shark.
15 * <https://shark-ml.github.io/Shark/>
16 *
17 * Shark is free software: you can redistribute it and/or modify
18 * it under the terms of the GNU Lesser General Public License as published
19 * by the Free Software Foundation, either version 3 of the License, or
20 * (at your option) any later version.
21 *
22 * Shark is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU Lesser General Public License for more details.
26 *
27 * You should have received a copy of the GNU Lesser General Public License
28 * along with Shark. If not, see <http://www.gnu.org/licenses/>.
29 *
30 */
31//===========================================================================
32
33#ifndef SHARK_ALGORITHMS_DIRECT_SEARCH_OPERATORS_LATTICE
34#define SHARK_ALGORITHMS_DIRECT_SEARCH_OPERATORS_LATTICE
35
36#include <shark/LinAlg/Base.h>
38#include <shark/Core/Random.h>
39
40#include <numeric>
41
42namespace shark {
43
44namespace detail {
45
46/*
47 * An n-dimensional point sums to 's' if the sum of the parts equal 's',
48 * e.g. the point (x_0,x_1,x_2) sums to x_0+x_1+x+2 etc. The number of
49 * n-dimensional points that sum to 's' is given by the formula "N over K" where
50 * N is n - 2 + s + 1 and K is s.
51 */
52std::size_t sumlength(std::size_t const n, std::size_t const sum);
53
54void pointLattice_helper(
55 UIntMatrix & pointMatrix,
56 std::size_t const rowidx,
57 std::size_t const colidx,
58 std::size_t const sum_rest);
59
60/// \brief A corner is a point where exactly one dimension is non-zero.
61template <typename Iterator>
62bool isLatticeCorner(Iterator begin, Iterator end){
63 std::size_t nonzero = 0;
64 for(auto iter = begin; iter != end; ++iter)
65 {
66 if(nonzero > 1)
67 {
68 return false;
69 }
70 if(*iter > 0)
71 {
72 ++nonzero;
73 }
74 }
75 return nonzero == 1;
76}
77
78} // namespace detail
79
80
81
82/*
83 * Sample points uniformly and uniquely from the simplex given in the matrix.
84 * Corners are always included in the sampled point set (unless excplicitly
85 * turned off with keep_corners set to false). The returned matrix will always
86 * have n points or the same number of points as the original matrix if n is
87 * smaller.
88 */
89template <typename Matrix, typename randomType = shark::random::rng_type>
91 randomType & rng, Matrix const & matrix,
92 std::size_t const n,
93 bool const keep_corners = true
94){
95 // No need to do all the below stuff if we're gonna grab it all anyway.
96 if(matrix.size1() <= n){
97 return matrix;
98 }
99 Matrix sampledMatrix(n, matrix.size2());
100 std::set<std::size_t> added_rows;
101 // First find all the corners and add them to our set of sampled points.
102 if(keep_corners){
103 for(std::size_t row = 0; row < matrix.size1(); ++row){
104 if(detail::isLatticeCorner(matrix.major_begin(row), matrix.major_end(row))){
105 added_rows.insert(row);
106 }
107 }
108 }
109 while(added_rows.size() < n){
110 // If we sample an existing index it doesn't alter the set.
111 added_rows.insert(random::discrete(rng, std::size_t(0), matrix.size1() - 1));
112 }
113 std::size_t i = 0;
114 for(std::size_t major_idx : added_rows)
115 {
116 std::copy(
117 matrix.major_begin(major_idx), matrix.major_end(major_idx),
118 sampledMatrix.major_begin(i)
119 );
120 ++i;
121 }
122 return sampledMatrix;
123}
124
125/// \brief A preferred region in a lattice-sampled unit sphere.
126///
127/// A preferred region is a pair with a radius and a vector. The intersection
128/// of the vector and the unit sphere denotes the center of the preferred
129/// region, and the radius denotes the radius of a sphere constructed such that
130/// the center point (intersection of vector and unit sphere) is on its
131/// periphery. Points are then sampled on this smaller sphere and projected up
132/// on the unit sphere, thus restricting the covered
133/// segment/area/volume/hypervolume of the unit sphere. See Figure 1 in
134/// "Evolutionary Many-objective Optimization of Hybrid Electric Vehicle
135/// Control: From General Optimization to Preference Articulation"
136typedef std::pair<double, RealVector> Preference;
137
138/// \brief Return a set of evenly spaced n-dimensional points on the unit sphere
139/// clustered around the specified preference points.
141 std::size_t const n,
142 std::size_t const sum,
143 std::vector<Preference> const & preferences);
144
145/// \brief Return a set of of evenly spaced n-dimensional points on the "unit
146/// simplex" clustered around the specified preference points.
148 std::size_t const n,
149 std::size_t const sum,
150 std::vector<Preference> const & preferences);
151
152///\brief Computes the number of Ticks for a grid of a certain size
153///
154///
155/// Computes the least number of ticks in each dimension required to make an
156/// n-dimensional simplex grid at least as many points as a target number points.
157/// For example, the points in a two-dimensional
158/// grid -- a line -- with size n are the points (0,n-1), (1,n-2), ... (n-1,0).
160 std::size_t const n, std::size_t const target_count
161);
162
163/// \brief Returns a set of evenly spaced n-dimensional points on the "unit simplex".
164RealMatrix weightLattice(std::size_t const n, std::size_t const sum);
165
166/// \brief Return a set of evenly spaced n-dimensional points on the unit sphere.
167RealMatrix unitVectorsOnLattice(std::size_t const n, std::size_t const sum);
168
169/*
170 * Computes the pairwise euclidean distance between all row vectors in the
171 * matrix and returns a matrix containing, for each row vector, the indices of
172 * the 'n' closest row vectors.
173 */
174template <typename Matrix>
176 Matrix const & m, std::size_t const n){
177 const RealMatrix distances = distanceSqr(m, m);
178 UIntMatrix neighbourIndices(m.size1(), n);
179 // For each vector we are interested in indices of the t closest vectors.
180 for(std::size_t i = 0; i < m.size1(); ++i)
181 {
182 // Make some indices we can sort.
183 std::vector<std::size_t> indices(distances.size2());
184 std::iota(indices.begin(), indices.end(), 0);
185 // Sort indices by the distances.
186 std::sort(indices.begin(), indices.end(),
187 [&](std::size_t a, std::size_t b)
188 {
189 return distances(i, a) < distances(i, b);
190 });
191 // Copy the T closest indices into B.
192 std::copy_n(indices.begin(), n, neighbourIndices.major_begin(i));
193 }
194 return neighbourIndices;
195}
196
197
198} // namespace shark
199
200#endif // SHARK_ALGORITHMS_DIRECT_SEARCH_OPERATORS_LATTICE