1.1 --- a/Optimal_transportation.cpp Tue May 03 14:34:31 2011 +0200
1.2 +++ b/Optimal_transportation.cpp Tue May 03 14:38:27 2011 +0200
1.3 @@ -340,8 +340,15 @@
1.4 mass * avg_color * (QImage_area(image) / poly_to.area());
1.5 QPolygonF qpoly;
1.6 polygon_to_QPolygonF(poly_to, qpoly);
1.7 -
1.8 +
1.9 +#if 1
1.10 painter.setPen(QPen(double_to_QColor(color)));
1.11 +#else
1.12 + QPen pen(QColor(0,255,255));
1.13 + pen.setWidth(2);
1.14 +
1.15 + painter.setPen(pen);
1.16 +#endif
1.17 painter.setBrush(QBrush(double_to_QColor(color)));
1.18 painter.drawConvexPolygon(qpoly);
1.19 }
1.20 @@ -380,7 +387,7 @@
1.21
1.22 //_weights.resize(0);
1.23 decomposition.compute_weights(density_from, dt, _level, result, _weights,
1.24 - 0.05);
1.25 + 1e-6);
1.26 _weights = result;
1.27 _level ++;
1.28
2.1 --- a/include/OT/Density_base_2.hpp Tue May 03 14:34:31 2011 +0200
2.2 +++ b/include/OT/Density_base_2.hpp Tue May 03 14:38:27 2011 +0200
2.3 @@ -80,7 +80,7 @@
2.4 {
2.5 return convex_intersection (bound, poly, clipped_poly);
2.6 }
2.7 -
2.8 +#if 0
2.9 template<class K>
2.10 bool clip (const typename CGAL::Polygon_2<K> &bound,
2.11 const typename K::Segment_2 &segment,
2.12 @@ -100,7 +100,7 @@
2.13 CGAL::internal::intersection(ray, bound, K());
2.14 return CGAL::assign(clipped_segment, isect_object);
2.15 }
2.16 -
2.17 +#endif
2.18
2.19 template <class K, class Density>
2.20 double
3.1 --- a/include/OT/bfgs_descent.hpp Tue May 03 14:34:31 2011 +0200
3.2 +++ b/include/OT/bfgs_descent.hpp Tue May 03 14:38:27 2011 +0200
3.3 @@ -139,6 +139,7 @@
3.4 lbfgs_parameter_t param;
3.5 lbfgs_parameter_init(¶m);
3.6
3.7 + param.m = 10;
3.8 param.epsilon = 0;
3.9 param.linesearch = linesearch;
3.10 param.max_iterations = max_iterations;
3.11 @@ -148,8 +149,9 @@
3.12
3.13 size_t linesearches [] = {LBFGS_LINESEARCH_MORETHUENTE,
3.14 LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE,
3.15 - LBFGS_LINESEARCH_BACKTRACKING_WOLFE,
3.16 - LBFGS_LINESEARCH_BACKTRACKING_ARMIJO};
3.17 + //LBFGS_LINESEARCH_BACKTRACKING_WOLFE,
3.18 + //LBFGS_LINESEARCH_BACKTRACKING_ARMIJO
3.19 + };
3.20 size_t num_linesearches = sizeof(linesearches)/sizeof(linesearches[0]);
3.21
3.22 int current_line_search = 0;
3.23 @@ -189,6 +191,28 @@
3.24 } while (current_line_search < num_linesearches);
3.25 }
3.26
3.27 + template <class Function, class StoppingCriterion, class Progress>
3.28 + void
3.29 + bfgs_descent_zm(const Function &function,
3.30 + const StoppingCriterion &stop,
3.31 + Progress &prog,
3.32 + uvector &weights,
3.33 + int linesearch = LBFGS_LINESEARCH_MORETHUENTE,
3.34 + int max_iterations = 0)
3.35 + {
3.36 + size_t N = function.number_of_vertices();
3.37 + OT_ASSERT(weights.size() == N);
3.38 +
3.39 + for (size_t i = 0; i < weights.size(); ++i)
3.40 + weights[i] -= weights[N-1];
3.41 + weights.resize(N-1);
3.42 +
3.43 + bfgs_descent (function, stop, prog, weights, linesearch, max_iterations);
3.44 +
3.45 + weights.resize(N);
3.46 + weights[N-1] = 0;
3.47 + }
3.48 +
3.49 template <class Function>
3.50 void
3.51 bfgs_descent(const Function &function,
3.52 @@ -202,6 +226,20 @@
3.53 weights,
3.54 linesearch, max_iterations);
3.55 }
3.56 +
3.57 + template <class Function>
3.58 + void
3.59 + bfgs_descent_zm(const Function &function,
3.60 + uvector &weights,
3.61 + double eps,
3.62 + int linesearch = LBFGS_LINESEARCH_MORETHUENTE,
3.63 + int max_iterations = 0)
3.64 + {
3.65 + Display_value_gradient dvg;
3.66 + bfgs_descent_zm (function, Stop_norm_inf_gradient(eps), dvg,
3.67 + weights,
3.68 + linesearch, max_iterations);
3.69 + }
3.70 }
3.71
3.72 #endif
4.1 --- a/include/OT/objective_function.hpp Tue May 03 14:34:31 2011 +0200
4.2 +++ b/include/OT/objective_function.hpp Tue May 03 14:38:27 2011 +0200
4.3 @@ -191,12 +191,13 @@
4.4 typedef typename RT::Geom_traits::Kernel K;
4.5 typedef typename RT::Vertex_handle Vertex_handle;
4.6 typedef typename RT::Bare_point Bare_point;
4.7 + typedef typename RT::Point Point;
4.8
4.9 protected:
4.10 RT _rt;
4.11 const Function *_f;
4.12 std::vector<Vertex_handle> _vertices;
4.13 - std::vector<double> _diff_masses;
4.14 + std::vector<double> _diff_masses, _gradient;
4.15 size_t _num_active;
4.16
4.17 mutable bool _has_wasserstein;
4.18 @@ -212,14 +213,9 @@
4.19 template <class InputIterator>
4.20 void prepare(const Function *f, InputIterator begin, InputIterator end);
4.21
4.22 - size_t dimension()
4.23 - {
4.24 - return _diff_masses.size();
4.25 - }
4.26 -
4.27 const Gradient_type &gradient() const
4.28 {
4.29 - return _diff_masses;
4.30 + return _gradient;
4.31 }
4.32
4.33 double value() const
4.34 @@ -400,7 +396,7 @@
4.35 wsp.prepare (this, begin, end);
4.36 }
4.37
4.38 - size_t size() const
4.39 + size_t number_of_vertices () const
4.40 {
4.41 return _positions.size();
4.42 }
4.43 @@ -432,23 +428,31 @@
4.44 InputIterator begin,
4.45 InputIterator end)
4.46 {
4.47 - OT_ASSERT (end - begin == f->size());
4.48 -
4.49 _f = f;
4.50 _num_active = 0;
4.51 _has_wasserstein = false;
4.52 _wasserstein = 0;
4.53 +
4.54 + size_t N = _f->number_of_vertices();
4.55 + size_t Nw = end - begin;
4.56 +
4.57 + OT_ASSERT (Nw <= N);
4.58
4.59 // Recompute regular triangulation with new weights.
4.60 - rt_build (_rt, _f->_positions.begin(), _f->_positions.end(),
4.61 + rt_build (_rt,
4.62 + _f->_positions.begin(), _f->_positions.begin() + Nw,
4.63 begin, end);
4.64
4.65 + for (size_t i = Nw; i < N; ++i)
4.66 + _rt.insert (Point (_f->_positions[i], 0.0));
4.67 +
4.68 // Extract vertices and reorder
4.69 OT::rt_copy_vertices_in_order(_rt, _f->_indices, _vertices);
4.70
4.71 // Pre-compute gradient
4.72 typename std::vector<Vertex_handle>::iterator jt = _vertices.begin();
4.73 - _diff_masses.resize(_f->size());
4.74 +
4.75 + _diff_masses.resize (N);
4.76
4.77 for (size_t i = 0; jt != _vertices.end(); ++jt, ++i)
4.78 {
4.79 @@ -461,6 +465,11 @@
4.80 if (current_mass != 0.0)
4.81 _num_active++;
4.82 }
4.83 +
4.84 + _gradient.resize(Nw);
4.85 + for (size_t i = 0; i < Nw; ++i)
4.86 + _gradient[i] = _diff_masses[i]; // - _diff_masses[N-1];
4.87 +
4.88 //std::cerr << n_active << " active sites\n";
4.89 }
4.90 }
5.1 --- a/tests/CMakeLists.txt Tue May 03 14:34:31 2011 +0200
5.2 +++ b/tests/CMakeLists.txt Tue May 03 14:38:27 2011 +0200
5.3 @@ -1,4 +1,4 @@
5.4 add_executable(test_rasterization test_rasterization.cpp)
5.5 -add_executable(test_wasserstein test_wasserstein.cpp)
5.6 +#add_executable(test_wasserstein test_wasserstein.cpp)
5.7 # ADD_TEST(test_rasterization ${opttransport_BINARY_DIR}/tests/test_rasterization)
5.8 #message(${opttransport_BINARY_DIR}/tests/test_rasterization)
5.9 \ No newline at end of file