include/continewz/homotopy_post_certify.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : homotopy_post_certify.hpp
* DESCRIPTION: Homotopies with certified end-points (but not the entire paths)
*              This is OK, as long as we preserve the total number of
*              solutions of a system
* COPYRIGHT  : (C) 2011  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_HOMOTOPY_POST_CERTIFY_HPP
#define __MMX_HOMOTOPY_POST_CERTIFY_HPP
#include <continewz/homotopy_ball.hpp>
#include <continewz/homotopy_floating.hpp>
#include <basix/vector_sort.hpp>

namespace mmx {
#define Function slp_tangent<ball<C> >

template<typename V= homotopy_euler, typename W= homotopy_ball>
struct homotopy_post_certify {};

template<typename C>
struct homotopy_variant_helper<ball<C> > {
  typedef Homotopy_variant(C) V;
  typedef homotopy_post_certify<V> HV; };

/*
template<typename FV>
struct homotopy_variant_helper<ball<floating<FV> > > {
  typedef homotopy_floating<homotopy_euler,homotopy_euler> V;
  typedef homotopy_floating<homotopy_ball,homotopy_ball> W;
  typedef homotopy_post_certify<V,W> HV; };

template<typename FV>
struct homotopy_variant_helper<ball<complex<floating<FV> > > > {
  typedef homotopy_floating<homotopy_euler,homotopy_euler> V;
  typedef homotopy_floating<homotopy_ball,homotopy_ball> W;
  typedef homotopy_post_certify<V,W> HV; };
*/

/******************************************************************************
* Separating vectors of balls into disjoint groups
******************************************************************************/

template<typename R> vector<vector<nat> >
separate_balls (const vector<ball<R> >& v2) {
  typedef ball<R> B;
  nat n= N(v2);
  if (n == 0) return vector<vector<nat> > ();
  vector<B> v= copy (v2);
  vector<nat> sigma;
  sort (v, sigma);
  //mmerr << "v= " << v << "\n";
  //mmerr << "s= " << sigma << "\n";
  vector<R> lo= lower (v);
  vector<R> up= upper (v);
  //mmerr << "l= " << lo << "\n";
  //mmerr << "u= " << up << "\n";
  for (nat i=1  ; i<n; i++) up[i]  = max (up[i-1], up[i]);
  for (nat i=n-1; i>0; i--) lo[i-1]= min (lo[i-1], lo[i]);
  //mmerr << "l= " << lo << "\n";
  //mmerr << "u= " << up << "\n";
  vector<vector<nat> > r;
  vector<nat> current;
  current << sigma[0];
  for (nat i=1; i<n; i++) {
    if (up[i-1] < lo[i]) {
      r << current;
      current= vector<nat> ();
    }
    current << sigma[i];
  }
  r << current;
  //mmerr << "r= " << r << "\n";
  return r;
}

template<typename C> vector<vector<nat> >
separate_balls (const vector<ball<C> >& v, const vector<vector<nat> >& in) {
  vector<vector<nat> > out;
  for (nat i=0; i<N(in); i++) {
    nat n= N(in[i]);
    if (n <= 1) { out << in[i]; continue; }
    vector<ball<C> > w= extract (v, in[i]);
    vector<vector<nat> > sub= separate_balls (w);
    for (nat j=0; j<N(sub); j++) {
      //mmerr << "Extract " << sub[j] << " from " << in[i] << "\n";
      nat l= N(sub[j]);
      vector<nat> add= fill<nat> (l);
      for (nat k=0; k<l; k++) add[k]= in[i][sub[j][k]];
      out << add;
    }
  }
  return out;
}

template<typename R> vector<vector<nat> >
separate_balls (const vector<ball<complex<R> > >& v) {
  vector<vector<nat> > sub= separate_balls (Re (v));
  return separate_balls (Im (v), sub);
}

template<typename C> vector<vector<nat> >
separate_balls (const vector<vector<ball<C> > >& v) {
  nat n= N(v);
  if (n == 0) return vector<vector<nat> > ();
  vector<nat> cnt= fill<nat> (n);
  for (nat i=0; i<n; i++) cnt[i]= i;
  vector<vector<nat> > sub= vec<vector<nat> > (cnt);
  nat d= N(v[0]);
  if (d == 0) return vec<vector<nat> > (cnt);
  for (nat i=0; i<d; i++) {
    vector<ball<C> > w= fill<ball<C> > (ball<C> (0), n);
    for (nat j=0; j<n; j++) {
      ASSERT (N(v[j]) == d, "vectors of same length expected");
      w[j]= v[j][i];
    }
    sub= separate_balls (w, sub);
  }
  return sub;
}

/******************************************************************************
* Certified homotopies
******************************************************************************/

template<typename C, typename V, typename W>
struct homotopy_stepper_rep<ball<C>,homotopy_post_certify<V,W> >:
  public homotopy_stepper_rep<ball<C>,homotopy_abstract>
{
public:
  typedef homotopy_stepper<C,V> Numeric_stepper;
  typedef homotopy_stepper<ball<C>,W> Safe_stepper;
  typedef ball<C> Ball;

public:
  Numeric_stepper  numeric_stepper;
  Safe_stepper     safe_stepper;
  vector<Ball>     init;

public:
  inline homotopy_stepper_rep (const Function& f):
    homotopy_stepper_rep<ball<C>,homotopy_abstract> (f),
    numeric_stepper (center (f)),
    safe_stepper (f) {}

  vector<Ball>
  hom (const vector<Ball>& init, nat k, const Ball& t1) {
    return homx (safe_stepper, init, k, t1);
  }

  vector<vector<Ball> >
  homvp (const vector<vector<Ball> >& init, nat k, const vector<Ball>& p) {
    //mmerr << "From x = " << init << " via x[" << k << "] = " << p
    //      << " [post stepper]" << "\n";
    //double time_1= mmx_time ();
    vector<vector<C> > mid= center (init);
    mid= homx (numeric_stepper, mid, k, center (p));
    //mmerr << "mid= " << row_matrix (mid) << "\n";
    //double time_2= mmx_time ();
    //mmerr << "Initial approximation in " << time_2 - time_1 << " ms\n";
    vector<vector<Ball> > fin= as<vector<vector<Ball> > > (mid);
    fin= homx (safe_stepper, fin, k, p);
    //double time_3= mmx_time ();
    //mmerr << "Certification in " << time_3 - time_2 << " ms\n";

    nat n= N(init);
    bool busy= true;
    nat multiple= n;
    vector<bool> done= fill<bool> (false, n);
    while (busy) {
      busy= false;
      multiple= 0;
      vector<vector<nat> > parts= separate_balls (fin);
      for (nat i=0; i<N(parts); i++)
        if (N(parts[i]) > 1) {
          multiple += N(parts[i]);
          for (nat j=0; j<N(parts[i]); j++) {
            nat nr= parts[i][j];
            if (!done[nr]) {
              //mmerr << "Redo " << nr << "\n";
              //mmerr << "  Was: " << fin[nr] << "\n";
              fin[nr]= homx (safe_stepper, init[nr], k, p);
              //mmerr << "  Is : " << fin[nr] << "\n";
              done[nr]= true;
              busy= true;
            }
          }
        }
    }
    //if (multiple > 0) mmerr << "Multiple solutions: " << multiple << "\n";
    //double time_4= mmx_time ();
    //mmerr << "Error elimination in " << time_4 - time_3 << " ms\n";
    return fin;
  }
};

#undef Function
} // namespace mmx
#endif // __MMX_HOMOTOPY_POST_CERTIFY_HPP

---