include/continewz/riemann_bundle.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : riemann_bundle.hpp
* DESCRIPTION: digital Riemann bundles
* COPYRIGHT  : (C) 2008  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_RIEMANN_BUNDLE_HPP
#define __MMX_RIEMANN_BUNDLE_HPP
#include <continewz/riemann_surface.hpp>

namespace mmx {
#define TMPL_DEF template<typename Germ,typename R,typename V= std_riemann>
#define TMPL template<typename Germ,typename R,typename V>
#define C complex<R>
#define Square riemann_square<R,V>
#define Point riemann_point<R,V>
#define Squares table<bool,Square >
#define Iterator iterator<Square >
#define Surface riemann_surface<R,V>
#define Bundle riemann_bundle<Germ,R,V>
#define Bundle_rep riemann_bundle_rep<Germ,R,V>
#define Table table<vector<Germ>,Square >

/******************************************************************************
* Required functionality from Germ class
******************************************************************************/

template<typename P, typename Germ> bool inside (const P& p, const Germ& g);

template<typename P, typename Germ> inline bool
better (const P& p, const Germ& g1, const Germ& g2) {
  return inside (p, g1);
}

TMPL bool
inside (const Square& s, const Germ& g) {
  // Is the germ defined on the whole square s?
  R x= Re (center (s));
  R y= Im (center (s));
  R r= radius (s);
  return
    inside (C (x-r, y-r), g) &&
    inside (C (x-r, y+r), g) &&
    inside (C (x+r, y-r), g) &&
    inside (C (x+r, y+r), g);
}

/******************************************************************************
* Digital Riemann bundles associate to each square in the underlying
* digital Riemann surface a vector of "germs", which, glued together,
* constitute a bundle
******************************************************************************/

TMPL_DEF
class riemann_bundle_rep REP_STRUCT {
public:
  Surface surface;
  Table bundle;
public:
  inline riemann_bundle_rep () {}
  inline riemann_bundle_rep (const format<R>& fm):
    surface (fm), bundle (format<vector<Germ > > (format<Germ > (fm)),
                          format<Square > (fm)) {}
  inline riemann_bundle_rep (const Surface& rs):
    surface (rs), bundle (vector<Germ > (format<Germ > (CF(rs))),
                          format<Square > (CF(rs))) {}
  inline riemann_bundle_rep (const Surface& rs, const Table& bun):
    surface (rs), bundle (bun) {}
  void upgrade (const Squares& ss);
};

TMPL_DEF
class riemann_bundle {
INDIRECT_PROTO_3 (riemann_bundle, riemann_bundle_rep, Germ, R, V)
public:
  inline riemann_bundle ():
    rep (new Bundle_rep ()) {}
  inline riemann_bundle (const format<R>& fm):
    rep (new Bundle_rep (fm)) {}
  inline riemann_bundle (const Surface& rs):
    rep (new Bundle_rep (rs)) {}
  riemann_bundle (const Surface& rs, const Germ& g):
    rep (new Bundle_rep (rs))
  {
    for (Iterator it= entries (rs->squares); busy (it); ++it)
      rep->bundle[*it]= vec<Germ> (g);
  }
  riemann_bundle (const Surface& rs, const Table& bun):
    rep (new Bundle_rep (rs, bun)) {}
};
INDIRECT_IMPL_3 (riemann_bundle, riemann_bundle_rep,
                 typename Germ, Germ, typename R, R, typename V, V)
DEFINE_UNARY_FORMAT_3_bis (riemann_bundle)

TMPL inline format<R> CF (const Bundle& b) {
  return CF (b->surface); }

/******************************************************************************
* Required routines
******************************************************************************/

template<typename Op, typename Germ, typename R, typename V> nat
unary_hash (const Bundle& b) {
  return hard_hash (b);
}

template<typename Op, typename Germ, typename R, typename V> bool
binary_test (const Bundle& b1, const Bundle& b2) {
  return hard_eq (b1, b2);
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Bundle)
EXACT_IDENTITY_OP_SUGAR(TMPL,Bundle)

TMPL syntactic
flatten (const Bundle& b) {
  return apply ("bundle", flatten (b->surface), flatten (b->bundle));
}

TMPL
struct binary_helper<Bundle >: public void_binary_helper<Bundle > {
  static inline string short_type_name () {
    return "Rbun" * Short_type_name (Germ) * Short_type_name (R); }
  static inline generic full_type_name () {
    return gen ("Riemann_bundle", Full_type_name (R), Full_type_name (Germ)); }
  static inline generic disassemble (const Bundle& b) {
    return gen_vec (as<generic> (b->surface), as<generic> (b->bundle)); }
  static inline Bundle assemble (const generic& v) {
    return Bundle (as<Surface > (vector_access (v, 0)),
                   as<Table > (vector_access (v, 1))); }
  static inline void write (const port& out, const Bundle& b) {
    binary_write<Surface > (out, b->surface);
    binary_write<Table > (out, b->bundle); }
  static inline Bundle read (const port& in) {
    Surface rs= binary_read<Surface > (in);
    Table bun= binary_read<Table > (in);
    return Bundle (rs, bun); }
};

/******************************************************************************
* Upgrade after gluing
******************************************************************************/

TMPL vector<Germ >
merge (const Square& s, const vector<Germ >& v1, const vector<Germ >& v2) {
  if (N(v1) == 0) return v2;
  if (N(v2) == 0) return v1;
  if (N(v1) == 1 && inside (s, v1[0])) return v1;
  if (N(v2) == 1 && inside (s, v2[0])) return v2;
  vector<Germ > r= v1;
  nat i, j;
  for (i=0; i<N(v2); i++) {
    for (j=0; j<N(r); j++)
      if (v2[i] == r[j]) break;
    if (j == N(r)) r << v2[i];
  }
  return r;
}

TMPL void
Bundle_rep::upgrade (const Squares& ss) {
  for (Iterator it= entries (ss); busy (it); ++it) {
    Square s1= *it;
    Square s2= inside (surface)->upgrade (s1);
    vector<Germ> v1= bundle[s1];
    vector<Germ> v2= bundle[s2];
    reset (bundle, s1);
    bundle[s2]= merge (s2, v1, v2);
  }
}

/******************************************************************************
* Public interface
******************************************************************************/

TMPL Squares
glue (Bundle& b1, const Point& p1, const Bundle& b2, const Point& p2) {
  Squares replaced= glue (inside (b1)->surface, p1, b2->surface, p2);
  for (Iterator it= entries (b2->bundle); busy (it); ++it)
    inside (b1)->bundle[*it]= merge (*it, b1->bundle[*it], b2->bundle[*it]);
  inside (b1)->upgrade (replaced);
  return replaced;
}

TMPL void
glue (Bundle& b1, const Point& p1, const Bundle& b2) {
  Point p2= lift (b2->surface, project (p1));
  glue (b1, p1, b2, p2);
}

TMPL inline vector<Germ >
germs (const Bundle& b, const Point& p_old) {
  Point p= upgrade (b->surface, p_old);
  return b->bundle[owner (p)];
}

TMPL Germ
germ (const Bundle& b, const Point& p_old) {
  Point p= upgrade (b->surface, p_old);
  vector<Germ > v= germs (b, p);
  nat best= N(v);
  for (nat i=0; i<N(v); i++)
    if (inside (project (p), v[i]))
      if (best == N(v) || better (project (p), v[i], v[best]))
        best= i;
  if (best != N(v)) return v[best];
  mmout << "p= " << p << "\n";
  mmout << "v= " << v << "\n";
  ERROR ("germ not found");
  return v[0];
}

TMPL inline Point lift (const Bundle& b, const C& z) {
  return lift (b->surface, z); }
TMPL inline vector<Point > lifts (const Bundle& b, const C& z) {
  return lifts (b->surface, z); }
TMPL inline Point move (const Bundle& b, const Point& p, const C& delta) {
  return move (b->surface, p, delta); }

TMPL inline Surface bundle_surface (const Bundle& b) { // for glue
  return b->surface; }

#undef TMPL_DEF
#undef TMPL
#undef C
#undef Square
#undef Point
#undef Squares
#undef Iterator
#undef Surface
#undef Bundle
#undef Bundle_rep
#undef Table
} // namespace mmx
#endif // __MMX_RIEMANN_HPP

---