include/continewz/analytic_heuristic.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : analytic_heuristic.hpp
* DESCRIPTION: heuristic analytic contination
* 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__ANALYTIC_HEURISTIC__HPP
#define __MMX__ANALYTIC_HEURISTIC__HPP
#include <continewz/analytic_meta.hpp>

namespace mmx {
#define TMPL template<typename C,typename V>
#define R Abs_type(C)
#define Analytic analytic<C,V>
#define Analytic_rep analytic_rep<C,V>
#define Heuristic_rep heuristic_analytic_rep<C,V>
#define Series series<C>
#define Polynomial polynomial<C>
#define Pol polynomial<Fast_type(C) >
#define Assume assume_bound<typename unvectorize<R >::val>

// TODO:
// 1. Method for post-composition in analytic_rep
// 2. Routine conformal (f: Analytic, phi: Analytic) -> Analytic
//    for drawing f in the region phi^{inv} (disk_convergence (f o phi))
// 3. Routine for gluing several analytic functions together
//    (pick first one for which radial_eval returns a non-Nan)
// 4. Library of common phi-transforms and inverses

TMPL class heuristic_analytic_rep;
TMPL Analytic heuristic (const Polynomial& f);
TMPL Analytic heuristic (const Heuristic_rep* prev, const Polynomial& p,
                         const Analytic& phi, const C& t, const C& u);

template<typename C> polynomial<C>
compose (const polynomial<C>& p, const analytic<C>& post) {
  nat prec = precision (promote (0, CF(p)));
  nat order= (nat) (mmx_order_ratio * ((double) prec));
  nat n= max (N(p), order);
  polynomial<C> q= truncate (post, n);
  polynomial<C> r= promote (0, CF(p));
  for (int i= n-1; i>=0; i--) {
    polynomial<C> s= q * r;
    r= range (s, 0, min (n, N(s))) + polynomial<C> (p[i]);
  }
  return r;
}

template<typename C> void
heuristic_change (analytic<C>& phi, analytic<C>& psi,
                  analytic<C>& ret, const analytic<C>& src,
                  const C& v1, const C& v2, const R& r)
{
  C u= v2 - v1;
  u= (u / abs (u)) * r;
  analytic<C> z (polynomial<C> (promote (1, v1), (nat) 1));
  phi= radial (u * z / (analytic<C> (u) + z));
  psi= radial (u * z / (analytic<C> (u) - z));
  ret= radial (u * src / (analytic<C> (u) + src));
  /*
  phi= radial (u * (1 - exp (- z / u)));
  psi= radial (- u * log (1 - z / u));
  ret= radial (u * (1 - exp (- src / u)));
  */
}

TMPL
class heuristic_analytic_rep: public Analytic_rep {
protected:
  const Heuristic_rep* previous;
  Analytic   prev;
  Polynomial p;
  Pol        fp;
  R          r;
  Analytic   phi;
  C          t;
  C          u;
  vector<C>  u_path;
  Analytic   phi_at_u;

  void initialize () {
    nat order= (nat) (mmx_order_ratio * ((double) precision ((double) 0.0)));
    fp= fast (p);
    r = heuristic_radius (p, order);
    if (previous == NULL) {
      u_path= vec<C> (u);
      phi_at_u= radial_move (phi, u);
    }
    else {
      prev= Analytic (previous, true); // keep reference counter up to date
      u_path= copy (previous->u_path);
      u_path[N(u_path) - 1]= t - previous->t;
      u_path << (u - t);
      phi_at_u= move (phi, u_path);
    }
  }

  nat depth () const {
    if (previous == NULL) return 0;
    else return previous->depth () + 1;
  }

  bool check (const C& z, const R& alpha) const {
    C v2= radial_eval (phi_at_u, z-u);
    return abs (v2) < alpha * r;
  }

public:
  heuristic_analytic_rep (const Polynomial& p2):
    Analytic_rep (CF(p2)),
    previous (NULL),
    p (p2),
    phi (radial (Analytic ((Polynomial (promote (1, CF(p2)), (nat) 1))))),
    t (promote (0, CF(p2))),
    u (promote (0, CF(p2))) { initialize (); }
  heuristic_analytic_rep (const Heuristic_rep* previous2, const Polynomial& p2,
                          const Analytic& phi2, const C& t2, const C& u2):
    Analytic_rep (CF(p2)),
    previous (previous2),
    p (p2),
    phi (phi2),
    t (t2),
    u (u2) { initialize (); }
  syntactic expression (const syntactic& z) const {
    return apply ("heuristic", flatten (Expand (), z)); }
  Series Expand () const {
    return compose (p, phi_at_u); }
  R Radius_bound (nat order) const {
    ERROR ("not yet implemented");
    return R (0); }
  R Tail_bound (const R& r, nat order, Assume& a) const {
    ERROR ("not yet implemented");
    return R (0); }
  Analytic Default_move (const C& z) const {
    return heuristic<C,V> (previous, p, phi, t, u + z); }
  Analytic Move (const C& z) const {
    C u2= u + z;
    //mmout << u_path << ", " << z << "\n";
    R rr= shrink (promote (1, as<R> (abs (u2))));
    if (previous != NULL && previous->check (u2, rr))
      return previous->Move (u2 - previous->u);
    if (depth () >= 3 || check (u2, rr))
      return Default_move (z);
    C v = phi_at_u[0];
    C v2= radial_eval (phi_at_u, z);
    Analytic phi_post, psi_post;
    Analytic phi2;
    heuristic_change (phi_post, psi_post, phi2, phi, v, v2, r);
    Polynomial p2= compose (p, psi_post);
    //Analytic phi2= compose (phi, phi_post);
    Analytic phi2_at_u = move (phi2, u_path);
    Analytic phi2_at_u2= radial_move (phi2_at_u, u2 - u);
    nat order= (nat) (mmx_order_ratio * ((double) precision ((double) 0.0)));
    C w2= phi2_at_u2[0];
    R r2= heuristic_radius (p2, order);
    if (abs (w2) >= r2) return Default_move (z);
    //R d1= radial_distance (phi_at_u , u2 - u, r);
    //R d2= radial_distance (phi2_at_u, u2 - u, r2);
    //if (d1 >= d2) return Default_move (z);
    return heuristic<C,V> (this, p2, phi2, u2, u2); }
  Analytic Derive () const {
    ERROR ("not yet implemented");
    return Analytic (0, this->tfm ()); }
  C Radial_eval (const C& z) {
    if (!check (u+z, promote (1, as<R> (abs(z))))) return Nan (C);
    return slow<C> (eval (fp, fast (radial_eval (phi_at_u, z)))); }
};

TMPL Analytic
heuristic (const Polynomial& p) {
  return (Analytic_rep*) new heuristic_analytic_rep<C,V> (p);
}

TMPL Analytic
heuristic (const Heuristic_rep* prev, const Polynomial& p,
           const Analytic& phi, const C& t, const C& u) {
  return (Analytic_rep*) new heuristic_analytic_rep<C,V> (prev, p, phi, t, u);
}

template<typename C> analytic<C>
std_heuristic (const polynomial<C>& p) {
  return heuristic<C,std_analytic> (p);
}

#undef TMPL
#undef R
#undef Analytic
#undef Analytic_rep
#undef Heuristic_rep
#undef Series
#undef Polynomial
#undef Pol
#undef Assume
} // namespace mmx

#endif //__MMX__ANALYTIC_HEURISTIC__HPP

---