[Olena-patches] r3022: [Markov] Speed up the algorithm

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2008-12-10 Matthieu Garrigues &lt;garrigues(a)lrde.epita.fr&gt; [Markov] Speed up the algorithm. * markov/approx_exp.hh: New, exp using a lookup table. * markov/markov.hh: compute delta u directly. * markov/random.hh, * markov/random.hxx: speed up random using a lookup table. --- approx_exp.hh | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ markov.hh | 76 +++++++++++++++++++++++++-------------------------------- random.hh | 6 +++- random.hxx | 22 ++++++++++++++-- 4 files changed, 135 insertions(+), 46 deletions(-) Index: trunk/milena/sandbox/markov/approx_exp.hh =================================================================== --- trunk/milena/sandbox/markov/approx_exp.hh (revision 0) +++ trunk/milena/sandbox/markov/approx_exp.hh (revision 3022) @@ -0,0 +1,77 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + + + +#ifndef MLN_APPROX_EXP_HH +# define MLN_APPROX_EXP_HH + +class approx_exp +{ +public: + + approx_exp(float min, float max, unsigned nvalues) + : min_(min), + max_(max), + nvalues_(nvalues), + values_(nvalues), + step_((max - min) / float(nvalues)) + { + + for (unsigned i = 0; i < nvalues_; ++i) + values_[i] = exp(min_ + i * step_); + } + + inline + float + get(float x) + { + if (x < min_) + return values_[0]; + if (x > max_) + return values_[nvalues_ - 1]; + + return values_[unsigned(float(x - min_) / step_)]; + } + + inline + float + operator() (float x) + { + return get(x); + } + +private: + const float min_; + const float max_; + const unsigned nvalues_; + const float step_; + + std::vector<float> values_; +}; + +#endif // ! MLN_APPROX_EXP_HH Index: trunk/milena/sandbox/markov/random.hxx =================================================================== --- trunk/milena/sandbox/markov/random.hxx (revision 3021) +++ trunk/milena/sandbox/markov/random.hxx (revision 3022) @@ -11,26 +11,42 @@ template <typename T> Random<T>::Random (T inf, T sup) - : inf_ (inf), sup_ (sup) + : inf_ (inf), + sup_ (sup), + i_ (0) { assert (sup >= inf); srand (time (0)); + for (int i = 0; i < size_; i++) + values_[i] = gen(); } template <typename T> +inline T -Random<T>::get () const +Random<T>::gen () const { double res = (sup_ - inf_) * (((double) rand ()) / RAND_MAX) + inf_; return res; } template <> +inline bool -Random<bool>::get () const +Random<bool>::gen () const { bool res = rand () > (RAND_MAX / 2); + return res; +} + +template <typename T> +inline +T +Random<T>::get () +{ + T res = values_[i_ % size_]; + ++i_; return res; } Index: trunk/milena/sandbox/markov/markov.hh =================================================================== --- trunk/milena/sandbox/markov/markov.hh (revision 3021) +++ trunk/milena/sandbox/markov/markov.hh (revision 3022) @@ -4,6 +4,7 @@ # include <cmath> # include <iomanip> # include <random.hh> +# include <approx_exp.hh> # include <T_gen.hh> # include <mln/binarization/threshold.hh> # include <mln/core/routine/clone.hh> @@ -11,36 +12,36 @@ namespace mln { - template <typename I, typename O, typename N> - double compute_energy(const I& ima, const O& out, const N& nbh, bool xi, const mln_site(I) &p) + template <typename I, typename N> + inline + float compute_du(const I& ima, const I& out, const mln_site(I)& p, const N& nbh) { + // Compute du : energy of clique with the new value minus with the + // old value (p_out.val()) + + bool old_val = out(p); // Compute u(x,y) - double u; - if (xi == ima(p)) - u = 0; + float du = 0; + if (old_val == ima(p)) + du += 1; else - u = 1; - + du -= 1; // u(x) is cst so we don't care - double diff_sum = 0; - double coeff = 0; + // sum the differences between new_val and the neighboors. + int diff_sum = 0; mln_niter(N) n(nbh, p); for_all(n) - if (ima.domain().has(n)) - { - diff_sum += abs(xi - out(n)); - coeff ++; - } - - diff_sum /= coeff; + if (old_val != out(n)) + ++diff_sum; -// std::cout << "energy : " << (u + diff_sum) << std::endl; + du -= float(5 * (2 * diff_sum - int(nbh.size()))) / nbh.size(); - return (u + diff_sum * 5); + return du; } + template <typename I> void dump(const Image<I>& ima) { @@ -58,50 +59,42 @@ const I &ima = exact(ima_); const N &nbh = exact(nbh_); - mln_ch_value(I, bool) bin = binarization::threshold(ima, 255 / 2); // FIXME : max - mln_ch_value(I, bool) out(bin.domain()); - - io::pbm::save(out, "threshold.pbm"); + typedef mln_ch_value(I, bool) O; + O bin = binarization::threshold(ima, 255 / 2); // FIXME : max + O out(bin.domain()); temperature_generator gtemp(start_temp, 0.8); - double temp = start_temp; + approx_exp my_exp(-1000, 5, 100000); + float temp = start_temp; Random<bool> v_random(0, 1); // mettre max et min ? - Random<double> p_random(0., 1.); // idem + Random<float> p_random(0., 1.); // idem unsigned modifications = 42; unsigned turn = 1; bool gradient = false; - int diffneg = 0; + while (!gradient || modifications) { // Trace. // dump(out); - mln_piter(I) p(bin.domain()); modifications = 0; - - for_all(p) + mln_piter(O) p_out(out.domain()); + for_all(p_out) if (v_random.get()) { - bool v = v_random.get(); + float d_u = compute_du(bin, out, p_out, nbh); - double u = compute_energy(bin, out, nbh, out(p), p); - double up = compute_energy(bin, out, nbh, v, p); - - double d_u = up - u; - double proba = exp(-d_u / temp); - - if ((d_u < 0 || !gradient && (p_random.get() < proba)) && out(p) != v) + if ((d_u < 0 || !gradient && (p_random.get() < my_exp(-d_u / temp)))) { - if (d_u < 0) - diffneg ++; - out(p) = v; - modifications ++; + out(p_out) = !out(p_out); + ++modifications; } } temp = gtemp; - std::cout << "Turn : " << turn << " Modifs : " << modifications << " DiffNeg : " << diffneg << " Temp : " << temp << std::endl; + + std::cout << "Turn : " << turn << " Temp : " << temp << " Modifications : " << modifications << std::endl; turn ++; if (!gradient && !modifications) { @@ -109,7 +102,6 @@ modifications = 1; gradient = true; } - diffneg = 0; } return out; Index: trunk/milena/sandbox/markov/random.hh =================================================================== --- trunk/milena/sandbox/markov/random.hh (revision 3021) +++ trunk/milena/sandbox/markov/random.hh (revision 3022) @@ -9,11 +9,15 @@ public: Random (T inf, T sup); - T get () const; + T get (); + T gen () const; public: + static const unsigned size_ = 10000; T inf_; T sup_; + int i_; + T values_[size_]; };