continewz: include/continewz/homotopy

00001 
00002 /******************************************************************************
00003 * MODULE     : homotopy_ball.hpp
00004 * DESCRIPTION: Certified homotopies
00005 * COPYRIGHT  : (C) 2009  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX_HOMOTOPY_BALL_HPP
00014 #define __MMX_HOMOTOPY_BALL_HPP
00015 #include <continewz/homotopy_euler.hpp>
00016 
00017 namespace mmx {
00018 #define Function slp_tangent<ball<C> >
00019 #define Function_C slp_tangent<C>
00020 
00021 struct homotopy_ball {};
00022 
00023 /*
00024 template<typename C> struct homotopy_variant_helper<ball<C> > {
00025   typedef homotopy_ball HV; };
00026 */
00027 
00028 /******************************************************************************
00029 * Subroutines
00030 ******************************************************************************/
00031 
00032 template<typename C> Abs_type(C)
00033 l2_norm (const vector<C>& v) {
00034   typedef Abs_type(C) R;
00035   R r= 0;
00036   for (nat i=0; i<N(v); i++)
00037     r += square (abs (v[i]));
00038   return sqrt (r);
00039 }
00040 
00041 /******************************************************************************
00042 * Certified homotopies
00043 ******************************************************************************/
00044 
00045 template<typename C>
00046 struct homotopy_stepper_rep<ball<C>,homotopy_ball>:
00047   public homotopy_stepper_rep<ball<C>,homotopy_abstract>
00048 {
00049 public:
00050   typedef ball<C> Ball;
00051   typedef Radius_type (Ball) R;
00052   typedef typename rounding_helper<R>::UV Up;
00053   typedef typename rounding_helper<R>::DV Down;
00054 
00055 public:
00056   Function_C   fc;
00057   vector<Ball> init;
00058   nat          k;
00059   Ball         goal;
00060   C            t0;
00061   C            t1;
00062   C            dt;
00063   vector<C>    x;
00064   vector<C>    fx;
00065   matrix<C>    Jx;
00066   matrix<C>    Ix;
00067   vector<Ball> tube;
00068   vector<C>    x2;
00069   xnat         coh;
00070 
00071 public:
00072   inline homotopy_stepper_rep (const Function& f):
00073     homotopy_stepper_rep<ball<C>,homotopy_abstract> (f),
00074     fc (center (f)) {}
00075 
00076   void
00077   prepare_step () {
00078     //mmerr << "-----------------------\n";
00079     //mmerr << "x = " << x << "\n";
00080     fx= eval (fc, x);
00081     //mmerr << "fx= " << fx << "\n";
00082     Jx= jacobian (fc, x);
00083     //mmerr << "Jx= " << Jx << "\n";
00084     Ix= invert (remove_column (Jx, k));
00085     //mmerr << "Ix= " << Ix << "\n";
00086   }
00087 
00088   bool
00089   newton_step () {
00090     vector<C> dd= -Ix * fx;
00091     vector<C> dx= append (range (dd, 0, k), cons (C(0), range (dd, k, N(dd))));
00092     vector<C> fy= eval (fc, x+dx);
00093     if (11 * l2_norm (fy) >= 10 * l2_norm (fx)) return false;
00094     x += dx;
00095     prepare_step ();
00096     return true;
00097   }
00098 
00099   void
00100   shrink_step () {
00101     newton_step ();
00102     dt= dt / 2;
00103   }
00104 
00105   void
00106   confirm_step () {
00107     x= x2;
00108     dt= dt * sqrt (C (2));
00109     if (abs (t1 - x[k]) < abs (dt)) dt= t1 - x[k];
00110     prepare_step ();
00111   }
00112 
00113   void
00114   propose_tube () {
00115     vector<C> fy= fx + column (Jx, k) * dt;
00116     vector<C> dd= -(Ix * fy);
00117     vector<C> dx= append (range (dd, 0, k), cons (dt, range (dd, k, N(dd))));
00118     x2= x + dx;
00119     R r= as<R> (l2_norm (dd));
00120     tube= fill<Ball> (N(x));
00121     for (nat i=0; i<N(tube); i++)
00122       if (i == k)
00123         tube[i]= Ball (x[k] + (dt / 2), as<R> (abs (dt)) / R(1.999999));
00124       else
00125         tube[i]= Ball (x[i] + (dx[i] / 2), r);
00126     //mmerr << "dx= " << dx << "\n";
00127     //mmerr << "r = " << r  << "\n";
00128     //mmerr << "ct= " << center (tube) << "\n";
00129     //mmerr << "rt= " << radius (tube) << "\n";
00130     if (dt != 0) coh= min (coh, coherence (fc, x, fx, Jx, k, dd));
00131   }
00132 
00133   void
00134   final_tube (const Ball& target, const R& scale) {
00135     R dt= Up::add (abs_up (target - x[k]), radius (target));
00136     vector<C> df= column (Jx, k) * as<C> (dt);
00137     vector<C> d1= -(Ix * fx);
00138     vector<C> d2= -(Ix * df);
00139     R r= scale * as<R> (l2_norm (d1) + l2_norm (d2));
00140     tube= fill<Ball> (N(x));
00141     for (nat i=0; i<N(tube); i++)
00142       if (i == k) tube[i]= Ball (x[k], dt);
00143       else tube[i]= Ball (x[i], r);
00144     //mmerr << "dt= " << dt << "\n";
00145     //mmerr << "ct= " << center (tube) << "\n";
00146     //mmerr << "rt= " << radius (tube) << "\n";
00147   }
00148 
00149   bool
00150   certify_tube () {
00151     //mmerr << "U = " << tube << "\n";
00152     matrix<Ball> JU = jacobian (this->f, tube);
00153     //mmerr << "JU= " << JU << "\n";
00154     matrix<Ball> J  = remove_column (JU, k);
00155     //mmerr << "J = " << J << "\n";
00156     //mmerr << "cJ= " << center (J) << "\n";
00157     vector<Ball> V  = column (JU, k);
00158     //mmerr << "V = " << V << "\n";
00159     matrix<Ball> T  = invert (as<matrix<Ball> > (center (J)));
00160     //mmerr << "T = " << T << "\n";
00161     matrix<Ball> Eps= T*J - 1;
00162     //mmerr << "E = " << Eps << "\n";
00163     R bnd1= Down::sub (R(1), l2_norm_up (Eps));
00164     //mmerr << "b1= " << bnd1 << "\n";
00165     if (bnd1 <= as<R> (0)) return false;
00166     vector<Ball> X= as<vector<Ball> > (center (tube));
00167     vector<Ball> Y= T * eval (this->f, X);
00168     vector<Ball> W= T * V;
00169     R bnd2= l2_norm_up (Y);
00170     //mmerr << "X = " << X << "\n";
00171     //mmerr << "Y = " << Y << "\n";
00172     //mmerr << "b2= " << bnd2 << "\n";
00173     R bnd3= l2_norm_up (W) * radius (tube[k]);
00174     //mmerr << "b3= " << bnd3 << "\n";
00175     R bnd4= big_min (radius (remove_entry (tube, k)));
00176     //mmerr << "b4= " << bnd4 << "\n";
00177     return Up::add (bnd2, bnd3) < Down::mul (bnd4, bnd1);
00178   }
00179 
00180   bool
00181   finish (const Ball& target) {
00182     // FIXME: it is not clear until which magnitude we should go
00183     // for the scale. Please investigate closer, so that we may
00184     // do with less than 10 iterations
00185     R scale= R (2.0);
00186     for (nat i=0; i<10; i++) {
00187       //mmerr << "scale= " << scale << "\n";
00188       final_tube (target, scale);
00189       if (certify_tube ()) {
00190         tube[k]= target;
00191         return true;
00192       }
00193       scale *= R (2.0);
00194       if (i >= 5) scale *= incexp2 (R (2.0), i-5);
00195     }
00196     return false;
00197   }
00198 
00199   vector<Ball>
00200   hom (const vector<Ball>& init2, nat k2, const Ball& t1b) {
00201     //mmerr << "From x = " << init2 << " to x[" << k2 << "] = " << t1b
00202     //<< " [ball stepper]" << "\n";
00203     ASSERT (N(init2) == N(this->f) + 1, "homotopy stepper not initialized");
00204     ASSERT (k2 < N(init2), "coordinate out of range");
00205     init= init2;
00206     k= k2;
00207     goal= t1b;
00208     t0= center (init[k]);
00209     t1= center (t1b);
00210     dt= t1-t0;
00211     x=center (init);
00212     add_multiplicative_error (goal);
00213     coh= precision (C(1));
00214 
00215     int credit= 1000;
00216     prepare_step (); 
00217     propose_tube ();
00218     while (abs (t1 - x[k]) > abs (sqrt (Accuracy (C)))) {
00219       //mmerr << "x= " << x << "\n";
00220       credit--;
00221       if (credit < 0) mmerr << "Failed at " << x << "\n";
00222       ASSERT (credit >= 0, "continuation failed");
00223       propose_tube ();
00224       if (coh <= 10) break;
00225       if (!certify_tube ()) shrink_step ();
00226       else confirm_step ();
00227     }
00228     //mmerr << "x= " << x << "\n";
00229     while (newton_step ());
00230     //mmerr << "x= " << x << "\n";
00231     prepare_step ();
00232 
00233     if (!finish (goal)) {
00234       Ball target= Ball (x[k]);
00235       add_multiplicative_error (target);
00236       ASSERT (finish (target), "continuation failed");
00237     }
00238     //mmerr << "Steps= " << 1000-credit << "\n";
00239     //mmerr << "Coherence= " << coh << "\n";
00240     return tube;
00241   }
00242 };
00243 
00244 #undef Function
00245 #undef Function_C
00246 } // namespace mmx
00247 #endif // __MMX_HOMOTOPY_BALL_HPP
include/continewz/homotopy_ball.hpp
