00001 /***************************************************************************
00002 rootscan.h
00003 -----------
00004
00005 begin : Fri Sep 15 12:22:11 CEST 2000
00006 author : (C) 2000 by Michael Peeters
00007 email : Michael.Peeters@vub.ac.be
00008 ***************************************************************************/
00009
00010 /***************************************************************************
00011 * *
00012 * This program is free software; you can redistribute it and/or modify *
00013 * it under the terms of the GNU General Public License as published by *
00014 * the Free Software Foundation; either version 2 of the License, or *
00015 * (at your option) any later version. *
00016 * *
00017 ***************************************************************************/
00018
00019 #ifndef ROOTSCAN_H
00020 #define ROOTSCAN_H
00021
00022 #include "numvector.h"
00023 #include "numvectorprint.h"
00024 #include "numerictraits.h"
00025 #include "vectorfunction.h"
00026 #include <string>
00027 #include "parameter.h"
00028 #include <list>
00029 #include <iostream>
00030 #include <map>
00031 #include "vfwithbump.h"
00032 #include "eigenvalues.h"
00033 #include "newtonroot.h"
00034 #include "jacobian.h"
00035
00036 namespace MODEL
00037 {
00041 template<integer dims, typename nelem=number, class NT =
00042 NumericTraits<nelem,dims> >
00043 class ScanList
00044 {
00045 public:
00046 typedef typename NT::vect vect;
00047
00053 typedef vector< NumVector<dims,nelem,NT> > parpoint;
00054
00059 typedef list<parpoint> parpointlist;
00060 typedef map<number,parpointlist> datalist;
00061
00063 typedef bool criteria(const parpoint&);
00064
00068 void clear() {dat.clear();}
00069
00071 parpointlist get_solution(number param)
00072 {
00073 datalist::iterator here=dat.upper_bound(param);
00074 return here->second;
00075 }
00076
00078 void add_point(const vect& v)
00079 {
00080 p.push_back(v);
00081 }
00082
00085 void add_param(const number& param)
00086 {
00087 dat[param].push_back(p);
00088 p.clear();
00089 }
00090
00094 ScanList select(criteria ifthis) const
00095 {
00096 // make a copy
00097 ScanList selected(*this);
00098
00099 // Run thru it, removing all bad elements
00100 for(datalist::iterator p=selected.dat.begin();
00101 p!=selected.dat.end();++p)
00102 {
00103 parpointlist::iterator pp=p->second.begin();
00104 while(pp!=p->second.end())
00105 if(!ifthis(*pp)) pp=p->second.erase(pp);
00106 else ++pp;
00107 }
00108
00109 return selected;
00110 }
00111
00115 void print_raw(ostream& out)
00116 {
00117 out << "# Raw data -----------------------" << endl;
00118
00119 for(datalist::iterator parr=dat.begin();
00120 parr!=dat.end();
00121 ++parr)
00122 {
00123
00124 for(parpointlist::iterator pointlist=parr->second.begin();
00125 pointlist!=parr->second.end();
00126 ++pointlist)
00127 {
00128 out << parr->first; // Write parameter value
00129 for(parpoint::iterator point=pointlist->begin();
00130 point!=pointlist->end();
00131 ++point)
00132 out << "\t" << *point;
00133 out << endl;
00134 }
00135 }
00136 }
00137
00140 ScanList operator+=(const ScanList<dims,nelem>& extra)
00141 {
00142 const datalist& other=extra.dat;
00143
00144 for(datalist::const_iterator parr=other.begin();
00145 parr!=other.end();
00146 ++parr)
00147 {
00148
00149 for(parpointlist::const_iterator pointlist=parr->second.begin();
00150 pointlist!=parr->second.end();
00151 ++pointlist)
00152 {
00153 // out << parr->first; // Write parameter value
00154 for(parpoint::const_iterator point=pointlist->begin();
00155 point!=pointlist->end();
00156 ++point)
00157 add_point(*point); // out << "\t" << *point;
00158 add_param(parr->first); // out << endl;
00159 }
00160 }
00161 return *this;
00162 }
00163
00164
00174 void print_list(ostream& out, const number& accuracy=1.)
00175 {
00176 datalist toprint(dat);
00177
00178 out << "# Scanlist output --------------------" << endl;
00179
00180 counter pieces(0);
00181
00182 // 1) Find a starting parameter value
00183
00184 datalist::iterator parrunner(toprint.begin());
00185 while(parrunner!=toprint.end()) { // while there are parameters left
00186
00187 // 2) Find a starting point
00188 parpointlist::iterator varrunner(parrunner->second.begin());
00189
00190 while (varrunner!=parrunner->second.end()) { // if there are any
00191 datalist::iterator nextrunner(parrunner);
00192 number param=nextrunner->first; // get first = KEY
00193
00194 // write comment
00195 out << "# -------------------------------"
00196 << " Start of list starting at " << param <<endl;
00197
00198 // We have a starting point, remove it from the list
00199 parpoint previous(*varrunner);
00200 number previousparam(param);
00201
00202 parpointlist::iterator diepiggy(varrunner);
00203
00204 parrunner->second.erase(diepiggy); // first erase,
00205 varrunner=parrunner->second.begin(); // then find another
00206
00207 // Add it to our list:
00208 out << param;
00209 for(parpoint::iterator i=previous.begin();
00210 i!=previous.end();++i)
00211 out << "\t" << *i;
00212 out << endl;
00213
00214 // Find a next point (a higher par value)
00215 ++nextrunner;
00216 while (nextrunner!=toprint.end()){ // loopy - we have points
00217 // find someone to love
00218 bool loved=false;
00219 number nextparam=nextrunner->first;
00220
00221 parpointlist::iterator inrunner=nextrunner->second.begin();
00222
00223 // calculate match requirements
00224 number maxchange=accuracy*nextparam/previousparam;
00225
00226 while (inrunner!=nextrunner->second.end()) { // find match
00227
00228 // Ok, we have one
00229 parpoint match(previous);
00230 parpoint gotcha(*inrunner);
00231 bool matched=true;
00232
00233 // Is it any good: calculate all differences
00234 for(parpoint::iterator
00235 i=match.begin(),
00236 j=gotcha.begin(),
00237 k=previous.begin();
00238 i!=match.end();
00239 ++i,++j,++k)
00240 {
00241 *i -= *j;
00242 *i /=length(*k); // normalized difference
00243
00244 // the abs is to make sure it works for complex
00245 // types
00248 matched= matched && abs(length(*i)) <maxchange;
00249 }
00250
00253 if(matched) { // we have a match
00254 previous=gotcha;
00255
00256 // Remove it from the list
00257 parrunner->second.erase(inrunner);
00258
00259 // Add it to our list:
00260 out << nextparam;
00261 for(parpoint::iterator i=previous.begin();
00262 i!=previous.end();++i)
00263 out << "\t" << *i;
00264 out << endl;
00265
00266 loved=true;
00267 break; // We are happy, go to next parval
00268 }
00269 else ++inrunner;
00270
00271 } // while find match
00272
00273 if(loved) {
00274 // try it with the next parval
00275 ++nextrunner;
00276 }
00277 else {
00278 // All exhausted, close our list
00279 // cout << "#----------------------------" << endl;
00280 ++nextrunner;
00281
00282 // break;
00283 }
00284
00285 previousparam=nextparam;
00286
00287 } // loopy
00288 // close the list
00289 out << "# -------------------------------"
00290 << " End of list starting at " << param <<endl<<endl<<endl;
00291 ++pieces;
00292
00293 // three endl's for gnuplot
00294
00295 } // if we had any
00296
00297
00298 // Go to next parameter value, as we are exhausted here
00299 ++parrunner;
00300
00301 } // while parameters left
00302 out << "# Number of distict pieces: " << pieces << endl;
00303 }
00304
00305 datalist& get_data(void)
00306 {
00307 return dat;
00308 }
00309
00310 private:
00311 datalist dat;
00312 parpoint p;
00313 };
00314
00315 //------------------------------------------------------------
00316
00325 template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
00326 class RootScan
00327 {
00328 public:
00329 typedef typename NT::vf vf;
00330 typedef typename NT::vect vect;
00331
00333 RootScan(vf& tobescanned, const string& param, vect starter=1.) : f(&tobescanned)
00334 {
00335 Parameter p=f->get_parameter(param);
00336 _p=&p;
00337
00338 // Generate some plausible starting points
00339 def_starters.push_back(starter);
00340 }
00341
00342
00346 RootScan(vf& tobescanned) : f(&tobescanned)
00347 {
00348 _p=&dummy;
00349 }
00350
00352 void add_start(const vect start)
00353 {
00354 def_starters.push_back(start);
00355 }
00356
00358 const ScanList<dims,nelem,NT>& scan(const number& from,
00359 const number& to,
00360 const counter& n)
00361 {
00362 // Clear out previous scan
00363 roots.clear();
00364
00365 number delta=(to-from)/number(n-1);
00366
00367 // A starting value list
00368 list< NumVector<dims,nelem,NT> > starters(def_starters);
00369
00370 for(*_p=from;*_p<=to;*_p+=delta) {
00371 #ifdef ROOTSCAN_DEBUG
00372 cerr << "PARAMETER SCAN: " << *_p << endl;
00373 #endif
00374 // Add a bumpy layer, to find other zeroes
00375 VFwithBump<dims,nelem,NT> bumpy(*f);
00376 // Find stationary solution (the zeroes)
00377 NewtonRoot<dims,nelem,NT> stat(bumpy);
00378 Jacobian<dims> J(*f,1E-6);
00379
00380 NumVector<dims,nelem,NT> solution;
00381
00382 // Get the list of starting points, defined from the
00383 // previous parameter value.
00384 list< NumVector<dims,nelem,NT> > oldstarters(starters);
00385
00386 starters.clear();
00387 starters=def_starters;
00388
00389 // Loop until no more points
00390
00391 list< NumVector<dims,nelem,NT> >::iterator
00392 start=oldstarters.begin();
00393 while(start!=oldstarters.end()){
00394 while(true){
00395 try {
00396 solution=stat(*start);
00397 }
00398 catch (logic_error le) {++start; break;}
00399
00400 //debug :
00401 // cerr << *_p << ":\t" << solution << endl;
00402
00403 if(stat.wrong_min()) {++start; break;}
00404 if(stat.no_root()) {++start; break;}
00405
00406 // Oh goody, a root !
00407
00408 // Find stability
00409 NumVector<dims,nelem,NT> reals=Eigenvalues<dims,nelem,NT>(J.calculate(solution)).real();
00410 NumVector<dims,nelem,NT> imags=Eigenvalues<dims,nelem,NT>(J.calculate(solution)).imag();
00411
00412 // Define the data
00413 #ifdef ROOTSCAN_DEBUG
00414 cerr << "STATPOINT:" << solution << endl;
00415 #endif
00416 roots.add_point(solution);
00417 roots.add_point(reals);
00418 roots.add_point(imags);
00419 // Save the data
00420 roots.add_param(*_p);
00421
00422 // add the solution to the bumplist
00423 bumpy.AddBump(solution);
00424
00425 // add to next try for next point
00426 starters.push_back( solution );
00427
00432 // // define new starting point
00433 // for(integer i=0;i<dims;i++)
00434 // {
00435 // NumVector<dims,nelem,NT> nsplus(solution);
00436 // NumVector<dims,nelem,NT> nsmin(solution);
00437 // nsplus[i]+= (delta / *_p)*solution[i]; // Linear
00438 // // Interpol
00439 // nsmin[i]-= (delta / *_p)*solution[i]; // Linear Interpol
00440 // starters.push_back(nsplus);
00441 // starters.push_back(nsmin);
00442
00443 // }
00444
00445 }
00446
00447 }
00448 }
00449 return roots;
00450 }
00451
00452 // A few examples of criteria
00454 static bool all(const ScanList<dims,nelem,NT>::parpoint& p) {return true;}
00458 static bool stable(const ScanList<dims,nelem,NT>::parpoint& p)
00459 {
00460 bool stab=true;
00461 for(counter i=0;i<dims;++i) stab=stab&&(p[1][i]<0);
00462
00463 return stab;
00464 }
00466 static bool unstable(const ScanList<dims,nelem,NT>::parpoint& p)
00467 {
00468 return !stable(p);
00469 }
00472 static bool positive(const ScanList<dims,nelem,NT>::parpoint &p)
00473 {
00474 bool pos=true;
00475 for(counter i=0;i<dims;++i) pos=pos&&(p[0][i]>-1E-4);
00476 // Agree, this is still negative, but still... Just to remove
00477 // rounding errors which play hell on the routines
00478
00479 return pos;
00480 }
00481
00482 static bool negative(const ScanList<dims,nelem,NT>::parpoint &p)
00483 {
00484 return !positive(p);
00485 }
00486
00487
00488 private:
00489 vf* f;
00490 number* _p;
00491 ScanList<dims,nelem,NT> roots;
00492
00493 // if we don't have a parameter to scan
00494 number dummy;
00495
00496 public:
00498 list< NumVector<dims,nelem,NT> > def_starters;
00499 };
00500
00501
00502
00503
00504
00505 } // end namespace
00506
00507 #endif
00508
00509 /*********************************************************************
00510 $Id: rootscan.h,v 1.3.2.2 2002/06/19 13:19:07 mpeeters Exp $
00511 **********************************************************************
00512
00513 $Log: rootscan.h,v $
00514 Revision 1.3.2.2 2002/06/19 13:19:07 mpeeters
00515 Made the select member const, so one can make a subselection of a const list.
00516
00517 Revision 1.3.2.1 2001/08/30 14:46:04 mpeeters
00518 Added #ifdef ROOTSCAN_DEBUG to surpress outputting ridiculous amounts of info to cerr.
00519
00520 Revision 1.3 2001/08/23 20:24:00 mpeeters
00521 Modified docs to provide correct version info.
00522
00523 Revision 1.2 2001/07/26 12:01:22 mpeeters
00524 Make start point selection easier (add_start). Removed some bugs.
00525
00526 Revision 1.1 2001/05/22 10:54:55 mpeeters
00527 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter".
00528
00529 Revision 1.5 2001/05/21 11:53:16 mpeeters
00530 Removed Makefile.in, which is automatically generated anyway.
00531
00532 Revision 1.4 2000/09/29 08:56:28 mpeeters
00533 Inserted into print_raw: write the parameter on EVERY line, otherwise gnuplot flips.
00534
00535 Revision 1.3 2000/09/22 08:50:03 mpeeters
00536 Changed dependencies causing trouble.
00537
00538 Revision 1.2 2000/09/15 10:26:31 mpeeters
00539 Cleaned out header and added CVS tails to files, removed superfuous
00540 @author comments, inserted dates
00541
00542 *********************************************************************/
More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be
Last update: June 2002.
Looking for Open Source? Check out
!