lib/e_adp.cc

   1 /*
   2  * Copyright (C) 2010 Felix Salfelder
   3  * Author: felix@salfelder.org
   4  *
   5  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20  * 02110-1301, USA.
  21  */
  22
  23 #include "e_adp.h"
  24 #include "m_divdiff.h"
  25 #ifdef HAVE_GSL_FIT_H
  26 # include <gsl/gsl_fit.h>
  27 #endif
  28 #define DEBUG_ADP
  29 using namespace std;
  30
  31 /*
  32 ''=working; accessed by get.
  33 0=current time step.
  34 1=old time step. accepted value only
  35 2=history
  36 */
  37
  38 /*----------------------------------------------------------------------------*/
  39 uint_t ADP_NODE::order() const{
  40   return( min(_order, CKT_BASE::_sim->get_tt_order() ) );
  41 } // order used for extrapolation.
  42 /*----------------------------------------------------------------------------*/
  43 //ADP_NODE::ADP_NODE( const ADP_NODE& p ) :
  44 //  // NODE(p),
  45 //  _number(p._number)
  46 //{
  47 //  unreachable();
  48 //  std::cout << "copy?? (should not happen)\n";
  49 //  assert(false);
  50 //  tr_lo = inf;
  51 //  tr_hi = -inf;
  52 //}
  53 /*----------------------------------------------------------------------------*/
  54 ADP_NODE::ADP_NODE( const std::string n, const CARD_LIST* c):
  55   NODE_BASE(n, unsigned(-1), c),
  56   dbg(0),
  57   tr_value(0.),
  58   tr_noise(0.),
  59   tt_value(0.),
  60   _abs_tt_err(0.),
  61   _abs_tr_err(0.),
  62   _rel_tt_err(0.),
  63   _rel_tr_err(0.),
  64   trhack(0.)
  65 {
  66   trace1("ADP_NODE::ADP_NODE ", n);
  67   set_label( n );
  68   set_owner( 0 );
  69
  70   _delta_expect = 0.;
  71   _positive = true; //hack
  72   //_debug = 0;
  73   _order = 0;
  74   _wdT = 0;
  75   dbg++;
  76
  77   ADP_NODE_LIST::adp_node_list.push_back( this );
  78
  79   trace2("ADP_NODE::init", long_label(), _user_number );
  80
  81   // tt_value3=NAN; //FIXME
  82   tt_expect=NAN; //FIXME
  83   tr_value3=NAN; //FIXME
  84
  85   // _n[n_ic].new_model_node("." + long_label() + ".ic", this);
  86   assert(c);
  87   tr_lo = inf;
  88   tr_hi = -inf;
  89 }
  90 /*----------------------------------------------------------------------------*/
  91 //ADP_NODE::ADP_NODE( const COMPONENT* c, const std::string n ) :
  92 //  NODE(),
  93 //  _number(0),
  94 //  dbg(0),
  95 //  tr_noise(0)
  96 //{
  97 //  assert(false);
  98 //  trace0("ADP_NODE::ADP_NODE ");
  99 //  init(c,n);
 100 //  assert(c);
 101 //  tr_lo = inf;
 102 //  tr_hi = -inf;
 103 //}
 104 /*----------------------------------------------------------------------------*/
 105 // ADP_NODE_UDC::ADP_NODE_UDC( const COMPONENT* c ) : ADP_NODE(c, "udc") { }
 106 /*----------------------------------------------------------------------------*/
 107 // ADP_NODE_RCD::ADP_NODE_RCD( const COMPONENT* c ) : ADP_NODE(c, "rcd"), udc(.2) {}
 108 /*----------------------------------------------------------------------------*/
 109 // obsolete
 110 void ADP_NODE::init(const COMPONENT* c, const std::string name_in)
 111 {
 112   assert(0);
 113   set_label( c->long_label() + "." +  name_in );
 114   tr_value = (0.);
 115   tr_noise = 0;
 116   dbg=0;
 117
 118   tt_value = 0.;
 119   _delta_expect = 0.;
 120   _positive = true; //hack
 121   //_debug = 0;
 122   _order = 0;
 123   _wdT = 0;
 124   _rel_tr_err = 0;
 125   _rel_tt_err = 0;
 126   _abs_tr_err = 0;
 127   _abs_tt_err = 0;
 128   dbg++;
 129
 130   // _number=_sim->newnode_adp();
 131   ADP_NODE_LIST::adp_node_list.push_back( this );
 132
 133   trace2("ADP_NODE::init", long_label(), _user_number );
 134
 135   // tt_value3=NAN; //FIXME
 136   tt_expect=NAN; //FIXME
 137   tr_value3=NAN; //FIXME
 138
 139   // _n[n_ic].new_model_node("." + long_label() + ".ic", this);
 140 }
 141 /*----------------------------------------------------------------------------*/
 142 ADP_NODE::~ADP_NODE(){
 143 }
 144 /*----------------------------------------------------------------------------*/
 145 void ADP_NODE::set_tr( hp_float_t x ){
 146   if(x>tr())  _region= 2;
 147   if(x<tr())  _region=-2;
 148   if(x==tr()) _region= 0;
 149   tr() = x;
 150 }
 151 /*----------------------------------------------------------------------------*/
 152 void ADP_NODE::set_tt( hp_float_t x )
 153 {
 154   assert(is_number(x));
 155   if( x<0 && _positive){
 156     error( bDANGER, "%s, expecting postive tt value, got %E\n", long_label().c_str(), x);
 157   }
 158   assert(x >=-1e-19 || !_positive);
 159   assert(m_() != INVALID_NODE);
 160   assert(m_() < _sim->_adp_nodes);
 161   tt() = x;
 162   // _sim->_tt1[m_()] = x;
 163 }
 164 /*----------------------------------------------------------------------------*/
 165 double ADP_NODE::State() const
 166 {
 167   // _sim->_tt is state @ begin of frame.
 168   // how to avoid trhack?
 169   assert(m_() < _sim->_adp_nodes);
 170   assert (_sim->_tt);
 171   return _sim->_tt[m_()] + trhack;
 172 }
 173 /*----------------------------------------------------------------------------*/
 174 hp_float_t ADP_NODE::get_aft_1() const { return tt1() + _delta[1] ; }
 175 /*----------------------------------------------------------------------------*/
 176 double ADP_NODE::tr_duration()const{ return CKT_BASE::_sim->last_time(); }
 177 /*----------------------------------------------------------------------------*/
 178 // called right after tr
 179 TIME_PAIR ADP_NODE_RCD::tt_review( ) {
 180   unreachable(); // no longer used
 181   return TIME_PAIR(0,NEVER);
 182 }
 183 /*----------------------------------------------------------------------------*/
 184 TIME_PAIR ADP_NODE::tt_preview( ) {
 185   TIME_PAIR _time_pre;
 186     _time_pre.reset();
 187     return _time_pre;
 188 }
 189 /*----------------------------------------------------------------------------*/
 190 void ADP_NODE::tr_add( double x ) const {
 191   assert(is_number(tr()));
 192   tr() += x;
 193 }
 194 /*----------------------------------------------------------------------------*/
 195 // FIXME: remove
 196 void ADP_NODE::reset() {
 197   trace2( "ADP_NODE::reset()", CKT_BASE::_sim->_Time0, tr_value );
 198   std::cout << "rst " << label() << "\n";
 199   assert(false); //BUG
 200   tr_value = 0; // tt1() = 0; tt() =0;
 201 }
 202 /*----------------------------------------------------------------------------*/
 203 double* ADP_NODE::derivatives(double* c)
 204 {
 205
 206   double Time[order()+1];
 207   switch(order()){
 208     case 4: incomplete();
 209     case 3:
 210       c[3] = tr3();
 211       Time[3] = Time3();
 212     case 2:
 213       c[2] = tr2();
 214       Time[2] = Time2();
 215     case 1:
 216       c[1] = tr1();
 217       Time[1] = Time1();
 218     case 0:
 219       c[0] = tr();
 220       Time[0] = Time0();
 221       break;
 222     default:
 223       incomplete();
 224   }
 225   ::derivatives(c, order()+1, Time);
 226   return c;
 227 }
 228 /*----------------------------------------------------------------------------*/
 229 void ADP_NODE::tt_accept()
 230 {
 231
 232   // sets t{t,r}_value0 to the accepted values.
 233   // should compute new wanted_dT.
 234
 235   trace4(("ADP_NODE::tt_accept() " + long_label()).c_str(),
 236       tr_value, tt_value, CKT_BASE::_sim->tt_iteration_number(), CKT_BASE::_sim->get_tt_order());
 237   assert( is_number(tr_value) );
 238   assert( CKT_BASE::_sim->last_time() >0 );
 239
 240   uint_t order = CKT_BASE::_sim->get_tt_order();
 241
 242   // BUG? why +1?
 243 //  set_order(min(_order+1, order)); hmmm
 244
 245   return;
 246
 247   switch (order) { untested();
 248     case 0:
 249       // fixme. merge
 250       return tt_accept_first();
 251       break;
 252     case 1:
 253       tr_dd12 = ( tr() - tr1() ) / dT0();
 254       break;
 255     case 2:
 256       tr_dd23 = tr_dd12;
 257       tr_dd123 = (tr_dd12 - tr_dd23)/( dT1() + dT0() );
 258       break;
 259     case 3:
 260       break;
 261     default:
 262       assert(false);
 263   }
 264
 265
 266   if (_integrator && order){
 267     trace0("ADP_NODE::tt_accept reintegration");
 268 //    tr_value = _delta[0];
 269     assert(tr()<10);
 270 //    tt_value = (this->*_integrator)(_delta[0]);
 271     assert( !_positive || tt_value > E_min );
 272
 273     if( tt_value > 1 ){
 274       error(bDANGER, "ADP_NODE::tt_accept overshoot? %f\n", label().c_str(), tt_value );
 275       tt_value=1;
 276     }
 277   } else {
 278     untested();
 279   }
 280
 281   assert(tt_value <=1 );
 282
 283   // std::cout << "acc after int " << label() << tr_value << "\n";
 284
 285   assert( tt_value == tt_value );
 286
 287   tt() = tt_value;
 288
 289   //tt_value += 13; // (printed) output value
 290   if (tt_value <0 && _positive){
 291     untested();
 292     tt_value=0;
 293 //    _delta[0] -= tt();
 294   }
 295
 296   trace2(("ADP_NODE::tt_accept" + short_label()).c_str(), tt_value, tr_value);
 297   trace1(("ADP_NODE::tt_accept done " + short_label()).c_str(), tt());
 298   assert(tt() >=0 || !_positive);
 299
 300   if (( dT0() / tr_duration() ) < 1.000000001){
 301     trace1(("ADP_NODE::tt_accept" + short_label()).c_str(), dT0() / tr_duration() );
 302     trace6(("ADP_NODE::tt_accept" + short_label()).c_str(), tt() ,
 303         tt1(), tr(), tt() - tt1() - tr1(), tr1(),
 304          ( tt() - tt1() - tr1()) / (tr1() + 1e-10)  );
 305   }
 306
 307 }
 308 /*----------------------------------------------------------------------------*/
 309 void ADP_NODE::tt_last( ) {
 310   unreachable();
 311  // tt() = tt() + tr(); //FIXME. add tt_keep_values.
 312  // if( tt() < 0 && _positive ){    //FIXME. add tt_keep_values.
 313  //   tr()+=tt();
 314  //   tt() = 0;                       //FIXME. add tt_keep_values.
 315  // }
 316
 317 }
 318 /*----------------------------------------------------------------------------*/
 319 // fixme: merge
 320 void ADP_NODE::tt_accept_first( )
 321 {
 322
 323   return;
 324   tt_value = tt() + tr_value;
 325   assert(tt() >=0 || !_positive);
 326   assert (tt()<=1);
 327
 328   trace3(("ADP_NODE::tt_accept_first" + short_label()).c_str(),
 329       tr_value,CKT_BASE::_sim->tt_iteration_number(), tt_value );
 330   assert( tr_value == tr_value );
 331   assert(order()==0);
 332   assert( CKT_BASE::_sim->get_tt_order() == 0 );
 333   assert( CKT_BASE::_sim->tt_iteration_number() == 0 );
 334
 335
 336   // assert(0  == CKT_BASE::_sim->_dT0);
 337
 338 //  tr_value = 0;
 339
 340
 341  // FIXME tt_value??
 342
 343   assert(tt_value == tt_value);
 344   if (_positive && tt_value < 0 ){
 345     untested();
 346     // assert(tt() >= 0);
 347   }
 348   trace4(("ADP_NODE::tt_accept_first" + short_label()).c_str(),
 349       tr_value,CKT_BASE::_sim->tt_iteration_number(), tt_value, tt_value);
 350
 351 }
 352 /*----------------------------------------------------------------------------*/
 353 // stress_last, tt_review, tt_accept, outdata, tt_advance, apply
 354 /*----------------------------------------------------------------------------*/
 355 void ADP_NODE::tt_advance()
 356 {
 357   if (_sim->_dT0) {
 358     _order = min(_new_order, CKT_BASE::_sim->get_tt_order());
 359   }else{
 360     _order=CKT_BASE::_sim->get_tt_order();
 361     assert(_order==0);
 362   }
 363   _new_order = 99;
 364 }
 365 /*----------------------------------------------------------------------------*/
 366 void ADP_NODE::tt_regress()
 367 {
 368   incomplete();
 369 }
 370 /*----------------------------------------------------------------------------*/
 371 hp_float_t ADP_NODE::get_total() const{
 372   return( NAN );
 373 }
 374 /*----------------------------------------------------------------------------*/
 375 void ADP_NODE::tr_stress_last( ) { unreachable(); }
 376 /*----------------------------------------------------------------------------*/
 377 #if 0
 378 void ADP_NODE::tr_stress_last( double val ) {
 379
 380   tr_value = val - get();
 381   trace3(("ADP_NODE::tr_stress_last" + label()).c_str(), val, get(), CKT_BASE::_sim->tt_iteration_number());
 382   _delta[0] = tr_value; // good idea? no.
 383
 384   trace4(("ADP_NODE::tr_stress_last(double) " + short_label()).c_str(),
 385       tr_value,CKT_BASE::_sim->tt_iteration_number(), val, get());
 386
 387  // std::cout << "ADP_NODE::tr_stress_last " << label() << tr_value << "\n";
 388  //
 389  //
 390 }
 391 #endif
 392 /*----------------------------------------------------------------------------*/
 393 void ADP_NODE::tr_expect_2_avg(){
 394   assert(order()==2);
 395
 396   _delta_expect = (_delta[2]+_delta[1])/2.0;
 397   tt_expect = tt_integrate_2_linear(_delta_expect);
 398 }
 399 /*----------------------------------------------------------------------------*/
 400 double ADP_NODE::tt_integrate_2_linear(double tr_){
 401   incomplete(); // regression?
 402   hp_float_t delta = (_delta[1] + tr_)/2 * ( dT0() / tr_duration() -1 );
 403   double ret = get_aft_1() + delta; // RET
 404   assert (ret == ret);
 405   // _debug+=2;
 406
 407   trace6(( "ADP_NODE::tt_integrate_2_linear "+short_label()).c_str(), get_tr(),
 408       get1(), ret, delta, tr_, tr_duration() );
 409
 410
 411   if (ret < -1e-10 && _positive ){
 412     unsigned order = CKT_BASE::_sim->get_tt_order();
 413     error(bDANGER, "* ADP_NODE::tt_integrate_2_linear neg step %i, Time0=%f (%f,%f,%f), %s, tt_value = %g, ( %g, %g; %g) %i tt: %f, expecting %f\n", \
 414         CKT_BASE::_sim->tt_iteration_number(),
 415         CKT_BASE::_sim->_Time0,   dT0(), Time_delta(), dT1(),
 416         short_label().c_str(), tt_value,
 417         _delta[2], _delta[1], tr_,
 418         order,
 419          get1(),
 420         ret);
 421 //    assert(false);
 422   }
 423   if (ret < 0 && _positive ) ret =0;
 424   _integrator = &ADP_NODE::tt_integrate_2_linear;
 425   return ret;
 426 }
 427 /*---------------------------------*/
 428 void ADP_NODE::tr_expect_2_linear(){
 429   assert(order()==2);
 430   trace4(( "ADP_NODE::tt_expect_2_linear " + short_label()).c_str(),
 431       tr_value, tr(), tr1(), tr2() );
 432   trace3(( "ADP_NODE::tt_expect_2_linear " + short_label()).c_str(), Time_delta(), get1(), dT1() );
 433   // _delta_expect = fabs( (_delta[1]) + ( (_delta[1]) - (_delta[2])) * (hp_float_t) ((Time_delta() )/dT1()));
 434   //
 435   // expected tr for time0
 436   _delta_expect = ( tr1() + ( (tr1()) - (tr2())) * (hp_float_t) ((Time_delta() )/dT1()));
 437   // _debug+=3;
 438
 439   if(_delta_expect != _delta_expect) assert(false);
 440   if(tr1() != tr1()) assert(false);
 441
 442 //  tt_expect = tt_integrate_2_linear(_delta_expect);
 443
 444   if (tt_expect < -1e-10 && _positive ){
 445     unsigned order = CKT_BASE::_sim->get_tt_order();
 446     error(bDANGER, "* ADP_NODE::tt_expect_2_linear neg step %i, Time0=%f (%f,%f,%f), %s, tt_value = %g, ( %g, %g; %g) %i tt: %f, expecting %f\n", \
 447         CKT_BASE::_sim->tt_iteration_number(),
 448         CKT_BASE::_sim->_Time0,   dT0(), Time_delta(), dT1(),
 449         short_label().c_str(), tt_value,
 450         tr2(), tr1(), _delta_expect,
 451         order,
 452          get1(),
 453         tt_expect);
 454 //    assert(false);
 455   }
 456   if (tt_expect < 0 && _positive ) tt_expect =0;
 457
 458   _corrector = NULL; // &ADP_NODE::tr_correct_1_exp; // HACK. ouch. move to ADP_EXP
 459
 460   tr()=_delta_expect;
 461 //  _integrator = &ADP_NODE::tt_integrate_2_linear;
 462 }
 463 /*---------------------------------*/
 464 // 3 points denoising expect
 465 void ADP_NODE::tr_expect_3_linear(){
 466   assert(order()==3);
 467   trace4(( "ADP_NODE::tt_expect3_linear " + short_label()).c_str(), tr_value, _delta[0], _delta[1], _delta[2] );
 468   trace3(( "ADP_NODE::tt_expect3_linear " + short_label()).c_str(), Time_delta(), get1(), dT1() );
 469   // tr_expect = fabs( (_delta[1]) + ( (_delta[1]) - (_delta[2])) * (hp_float_t) ((Time_delta() )/dT1()));
 470   //
 471
 472   hp_float_t t1 = dT2();
 473   hp_float_t t2 = dT1()+dT2();
 474   hp_float_t t3 = dT0()+dT1()+dT2();
 475
 476   // expected tr for time0
 477
 478   const double x[3]={0, t1, t2};
 479   const double y[3]={tr_value3, _delta[2], _delta[1]};
 480
 481   double c0=0, c1=0;
 482   double cov00, cov01, cov11, sumsq;
 483 #ifdef HAVE_GSL_FIT_H
 484   gsl_fit_linear( x, 1, y, 1, 3, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
 485 #else
 486   USE(x); USE(y); USE(cov00); USE(cov01); USE(cov11); USE(sumsq);
 487   incomplete(); untested();
 488 #endif
 489
 490   _delta_expect = c0 + c1 * t3;
 491
 492   if(_delta_expect != _delta_expect) _delta_expect=0;
 493   if(_delta[1] != _delta[1]) assert(false);
 494
 495   hp_float_t diff = (_delta[1] + _delta_expect)/2 * Time_delta() / tr_duration();
 496   tt_expect = get1() + diff;
 497   assert (tt_expect == tt_expect);
 498
 499   trace6(( "ADP_NODE::tt_expect3_linear "+short_label()).c_str(), get_tr(),
 500       get1(), tt_expect, diff, _delta_expect, tr_duration() );
 501
 502   if (tt_expect < 0 ){
 503     //positive?
 504     unsigned order = CKT_BASE::_sim->get_tt_order();
 505     error(bDANGER, "* ADP_NODE::tt_expect3_linear neg error step %i, Time0=%f  \
 506        * (%f,%f,%f), %s, tt_value = %g, ( %g, %g, %g)-> %g --  %i tt: %f, \
 507         expecting %f\n", \
 508         CKT_BASE::_sim->tt_iteration_number(),
 509         CKT_BASE::_sim->_Time0,   dT0(), dT1(), dT2(),
 510         short_label().c_str(),
 511         tt_value, tr_value3, _delta[2], _delta[1], _delta_expect,
 512         order,
 513          get1(),
 514         tt_expect);
 515 //    assert(false);
 516   }
 517 }
 518 /*---------------------------------*/
 519 /*---------------------------------*/
 520
 521 void ADP_NODE::tr_expect_2_square(){
 522   assert(order()==2);
 523   _delta_expect = square( sqrt(_delta[1]) + ( sqrt(_delta[1]) - sqrt(_delta[2])) * (hp_float_t) ((Time_delta() )/dT1()));
 524 }
 525 /*---------------------------------*/
 526 double foo(double x){
 527   // FIxme: horner.
 528   if(fabs(x)<1e-4)
 529     return (1-x/2 + 11/24*x*x -7/16*x*x*x )*exp(1);
 530   return pow((1+x),(1/x));
 531
 532 }
 533 /*----------------------------------------------------------------------*/
 534 double my_log(double x){
 535   if(x>0) return log(x);
 536   return -log(x);
 537 }
 538 /*---*/
 539 double ADP_NODE::tt_integrate_2_linear2(double tr_) {
 540   trace2("ADP_NODE::tt_integrate_2_exp", _delta[1], tr_);
 541   trace2("ADP_NODE::tt_integrate_2_exp", _delta[1],  _sign);
 542   hp_float_t h = tr_duration();
 543
 544   double del0 = tr_;
 545   double del1 = _delta[1];
 546
 547   double x = (del0 - del1) / 2;
 548   double delta = (_delta[1]  + tr_ )/2 * dT0()/h   + x ;
 549
 550   if (delta != delta || delta==-inf || delta==inf){
 551     delta= (tr_ + _delta[1])/2 * (dT0()-h) / h ;
 552     incomplete(); // call another thing.
 553   }
 554
 555   assert(delta != -inf);
 556   assert(delta != inf);
 557
 558   double ret = get_aft_1() + _sign * delta; // RET
 559
 560   trace5(( "ADP_NODE::tt_integrate_2_exp "+short_label()).c_str(),
 561       _delta[2], _delta[1], tr_, tt_expect, delta );
 562
 563   assert(ret == ret);
 564   _integrator = &ADP_NODE::tt_integrate_2_linear2;
 565   return ret;
 566 }
 567 /*----------------------------------------------------------------------*/
 568 void ADP_NODE::tr_expect_3_exp(){
 569   assert(order()==3);
 570   _sign=1;
 571
 572   if ( ( _delta[2] > 0 && _delta[1] > 0 && _delta[1] < _delta[2] && tr_value3 < _delta[2] )
 573       || ( _delta[2] < 0 && _delta[1] < 0 && _delta[1] > _delta[2] && tr_value3 > _delta[2] ) ){
 574     _order--;
 575     return tr_expect_2_exp();
 576   }
 577
 578   if (_delta[1] < 0 && _delta[2] < 0 && tr_value3 <0){
 579         _sign=-1;
 580   } else if ( ! (_delta[1] > 0 && _delta[2] > 0 && tr_value3 >0) ){
 581      return tr_expect_3_linear();
 582   }
 583
 584   if ( tr_value3 > _delta[2] && _delta[2] > _delta[1] ){
 585   }
 586
 587
 588   double ln3= log(_sign * tr_value3/tr_duration());
 589   double ln2= log(_sign * _delta[2]/tr_duration());
 590   double ln1= log(_sign * _delta[1]/tr_duration());
 591
 592   hp_float_t tdh=tr_duration()/2;
 593
 594   hp_float_t t3 = tdh;
 595   hp_float_t t2 = t3+dT2();
 596   hp_float_t t1 = t2+dT1();
 597
 598   hp_float_t x =  dT2()+dT1()+dT0()+tdh;
 599
 600
 601  //   exponent= ax^2 + bx + c
 602  //
 603   // neville @ x.
 604   double ln0 = ln3;
 605   ln0 += (ln2-ln3) * ( x - t3) / (t2-t3);
 606   ln0 += (ln1 - ln3 - ( ln2 - ln3 ) * (t1-t3)/ (t2-t3) ) * (x-t3) * ( x-t2) / (t1-t3) / (t1-t3);
 607
 608
 609   // int_0^x exp( ln1 + t/dT0*(ln0-ln1) ) @ x=dT0
 610   _delta_expect =  _sign*exp(ln0) * tr_duration();
 611   _corrector = &ADP_NODE::tr_correct_3_exp;
 612   tt_expect = tt_integrate_2_exp(_delta_expect); // 3 doesnt make sense yet?
 613 }
 614 /*----------------------------------------------------------------------*/
 615 double ADP_NODE::tr_correct_3_exp()
 616 {
 617   trace0(("ADP_NODE::tr_correct_3_exp correction " + label()).c_str());
 618   return (tr_value+_delta_expect)/2;
 619 }
 620 /*----------------------------------------------------------------------*/
 621 double ADP_NODE::tr_correct_generic()
 622 {
 623   // return(18);
 624   trace0(("ADP_NODE::tr_correct_generic correction " + label()).c_str());
 625   return (tr_value+_delta_expect)/2;
 626 }
 627 /*----------------------------------------------------------------------*/
 628 double ADP_NODE::tt_integrate_3_exp(double tr_)
 629 {
 630   assert(order()==3);
 631  // hp_float_t sign=1;
 632   double h=tr_duration();
 633
 634   if ( ( _delta[2] > 0 && _delta[1] > 0 && _delta[1] < _delta[2] && tr_value3 < _delta[2] )
 635       || ( _delta[2] < 0 && _delta[1] < 0 && _delta[1] > _delta[2] && tr_value3 > _delta[2] ) ){
 636     _order--;
 637     assert(false);
 638     //return tr_expect_2_exp();
 639   }
 640
 641   if (_delta[1] < 0 && _delta[2] < 0 && tr_value3 <0){
 642 //        sign=-1;
 643   } else if ( ! (_delta[1] > 0 && _delta[2] > 0 && tr_value3 >0) ){
 644 //    fixme.
 645     assert(false);
 646 //     return tr_expect_3_linear();
 647   }
 648
 649   if ( tr_value3 > _delta[2] && _delta[2] > _delta[1] ){
 650   }
 651
 652   double ln1= log(_sign*_delta[1]/tr_duration());
 653   double ln0 = log(_sign*tr_/tr_duration());
 654
 655   hp_float_t delta = _sign* (exp( ln1 + (ln0-ln1) ) - (_sign*_delta[1]/tr_duration()) )/(ln0-ln1)* dT0();
 656
 657   if (delta != delta || delta==inf || delta==-inf){
 658
 659     // fixme: cll anothert integrator.
 660     hp_float_t delta_trapez = (_delta[1] + tr_)/2 * ( dT0()-h) / tr_duration() ;
 661     trace2("ADP_NODE::tt_expect3_exp fallback to delta_trapez", delta_trapez, delta);
 662     // _debug+=1000;
 663
 664     _integrator = &ADP_NODE::tt_integrate_1_const;
 665     return( get1() + delta_trapez);
 666   }
 667
 668   tt_value = get1() + delta * _sign;
 669   assert(tt_value == tt_value);
 670   assert(tt_value >=0 || !_positive );
 671
 672   _integrator = &ADP_NODE::tt_integrate_3_exp;
 673   return tt_value;
 674 }
 675 /*----------------------------------------------------------------------*/
 676 void ADP_NODE::tr_expect_3_exp_fit(){
 677   assert(order()==3);
 678   hp_float_t sign=1;
 679
 680   if ( ( _delta[2] > 0 && _delta[1] > 0 && _delta[1] < _delta[2] && tr_value3 < _delta[2] )
 681       || ( _delta[2] < 0 && _delta[1] < 0 && _delta[1] > _delta[2] && tr_value3 > _delta[2] ) ){
 682     _order--;
 683     return tr_expect_2_exp();
 684   }
 685
 686   if (_delta[1] < 0 && _delta[2] < 0 && tr_value3 <0){
 687         sign=-1;
 688   } else if ( ! (_delta[1] > 0 && _delta[2] > 0 && tr_value3 >0) ){
 689      return tr_expect_3_linear();
 690   }
 691
 692   if ( tr_value3 > _delta[2] && _delta[2] > _delta[1] ){
 693   }
 694
 695   hp_float_t f0 = sign*tr_value3;
 696   hp_float_t f1 = sign*_delta[2];
 697   hp_float_t f2 = sign*_delta[1];
 698
 699   hp_float_t t1 = dT2();
 700   hp_float_t t2 = dT1()+dT2();
 701   hp_float_t t3 = dT0()+dT1()+dT2();
 702
 703   hp_float_t a = ( (log(f0) + log(f1) + log(f2)) * ( t1*t1 + t2*t2 ) - (t1 + t2) * ( t1*log(f1) + t2*log(f2) ) ) /
 704                          (3 * ( t1*t1 + t2*t2 ) - (t1+t2)*(t1+t2) );
 705   hp_float_t alpha = ( 3*( t1*log(f1) + t2*log(f2)) - (t1 + t2) * ( log(f0)+log(f1) +log(f2) ) ) /
 706                          (3 * ( t1*t1 + t2*t2 ) - (t1+t2)*(t1+t2) );
 707
 708 //  hack...
 709   if (alpha != alpha) {
 710     _order--;
 711     return tr_expect_2_linear();
 712   }
 713
 714   hp_float_t c=exp(a);
 715
 716   // _delta_expect =
 717   _delta_expect = sign*c * exp(alpha*t3);
 718
 719   trace6(( "ADP_NODE::tt_expect3_exp() "+short_label()).c_str(),  Time0(), tr_value3, _delta[2], _delta[1], get1(), _delta_expect );
 720   assert( _delta_expect == _delta_expect );
 721
 722   hp_float_t E = - c/alpha;
 723   hp_float_t delta = E * ( - exp( t3*alpha ) + exp(t2*alpha) )/ CKT_BASE::_sim->last_time() ;
 724   hp_float_t delta_trapez = (_delta[1] + _delta_expect)/2 * Time_delta() / CKT_BASE::_sim->last_time();
 725
 726   trace4("ADP_NODE::tt_expect3_exp", delta, delta_trapez, tt1(), get() );
 727   tt_expect = get1() +  sign * delta;
 728
 729   if ( tt_expect != tt_expect){
 730      std::cerr << "tt_expect is broken:\n"
 731         << "c=" << c << "\n"
 732         << "E=" << E << "\n"
 733         << "alpha=" << alpha<< "\n"
 734         << "( 0.0, " <<                 f0 <<" )\n"
 735         << "( "      << t1  << ", " <<  f1 <<" )\n"
 736         << "( "      << t2  << ", " <<  f2 <<" )\n"
 737         << "( "      << t3  << ", ?    )\n";
 738      _order--;
 739      return tr_expect_2_linear();
 740   }
 741
 742   if (tt_expect < -1e-8 && _positive){
 743     error(bDANGER, "* ADP_NODE::tt_expect3_exp neg error step %i, Time0=%f, %s,"
 744         " tt_value = %g, ( %g, %g, %g; %g) tte: %g \n", \
 745         CKT_BASE::_sim->tt_iteration_number(),
 746          dT0(), short_label().c_str(), tt_value,
 747          tr_value3, _delta[2], _delta[1], _delta_expect,
 748         tt_expect
 749         );
 750 //    assert(false);
 751   }
 752
 753   //sanitycheck
 754   if ( f0<0 && f1<0 && f2 <0)
 755     assert (_delta_expect <= 0);
 756
 757 }
 758 /*---------------------------------*/
 759 void ADP_NODE::tr_expect_3_quadratic(){
 760   unreachable();
 761   assert(order()==3);
 762   trace3(( "ADP_NODE::tr_expect_3_quadratic() "+short_label()).c_str(), tr_value3, _delta[2],  _delta[1] );
 763   // uses 1 2 3 to set expect.
 764   assert( tr_dd12 == ( _delta[1] - _delta[2] ) / dT1());
 765   hp_float_t sum = dT0() + dT1() + dT2();
 766   hp_float_t tr_expect_kaputt;
 767
 768 //        doch nicht kaputt.
 769   tr_expect_kaputt = ( \
 770       sum * ( _delta[1]  + tr_dd12 * dT0() ) \
 771       - dT0() * (( - dT1()* tr_dd12 + _delta[1] ) + (dT0()+dT1()) * ( tr_dd12 + tr_dd123* ( -dT2()-dT1() )  ) ) \
 772       ) / (dT2()+dT1());
 773   USE(tr_expect_kaputt);
 774
 775
 776   _delta_expect = -(((_delta[2] - tr_value3)/dT2() - (_delta[1] - _delta[2])/dT1())*(dT0() + dT1())/(dT1() + dT2()) -
 777     (_delta[2] - tr_value3)/dT2())*(dT0() + dT1() + dT2()) + tr_value3;
 778
 779   hp_float_t diff = (_delta[1] + _delta_expect)/2 * Time_delta() / CKT_BASE::_sim->last_time();
 780   tt_expect = get() + diff;
 781   assert (tt_expect == tt_expect);
 782 }
 783 /*--------------------------------------------------------------------------*/
 784 // evaluate tt at Time0 using parabola through tt1, tr1 and tr0
 785 double ADP_NODE::tt_integrate_2_poly(double offset) const
 786 {
 787   // tt1
 788   trace5("integrate_2_poly", tr1(), tr2(), dT1(), dT0(), tt1());
 789   double a = .5*(tr1() - tr2()) / dT1();
 790   double b = tr1();
 791   double T = dT0() - offset;
 792     // hack. _dTmin is (fixed) period width.
 793     // tt1() is at end of period(!)
 794
 795   double ret = a;
 796   ret *= T;
 797   ret += b;
 798   ret *= T;
 799   ret += tt1();
 800
 801   return ret;
 802 }
 803 /*--------------------------------------------------------------------------*/
 804 double ADP_NODE::tt_integrate_3_poly() const
 805 { incomplete();
 806   assert(order()>=3);
 807
 808   trace3("integrate_3_poly", tr1(), tr2(), tr3());
 809   trace4("integrate_3_poly", tt1(), dT0(), dT1(), dT2());
 810   double tr_dd12 = ( tr1() - tr2() ) / dT1();
 811   double tr_dd23 = ( tr2() - tr3() ) / dT2();
 812   double tr_dd123 = (tr_dd12 - tr_dd23)/( dT1() + dT2() );
 813   double sum = dT0() + dT1() + dT2();
 814   double integral;
 815
 816   integral = tt1() + ( \
 817       sum * ( tr1()  + tr_dd12 * dT0() ) \
 818       - dT0() * (( - dT1()* tr_dd12 + tr1() ) + (dT0()+dT1()) * ( tr_dd12 + tr_dd123* ( -dT2()-dT1() )  ) ) \
 819       ) / (dT2()+dT1());
 820
 821
 822   assert(is_number(integral));
 823   return integral;
 824 }
 825 /*--------------------------------------------------------------------------*/
 826 double ADP_NODE::debug(){
 827   return _debug;
 828 }
 829 /*---------------------------------*/
 830 void ADP_NODE::tr_expect_2_something(){
 831   assert(order()==2);
 832   _delta_expect = _delta[1] + (_delta[1] - _delta[2]) * dT0()/dT1();
 833   tt_expect = get1() + ( Time_delta() ) / CKT_BASE::_sim->last_time();
 834   trace2(( "ADP_NODE::tt_expect2_something "+short_label()).c_str(), get1(), tt_expect );
 835 }
 836 /*---------------------------------*/
 837 void ADP_NODE::tr_expect_1_const(){
 838   assert(order()==1);
 839   trace2(( "ADP_NODE::tt_expect_1_const() "+short_label()).c_str(),
 840        tr(), tr1() );
 841   trace3(( "ADP_NODE::tt_expect_1_const() "+short_label()).c_str(),
 842       Time_delta(), CKT_BASE::_sim->_time0, CKT_BASE::_sim->last_time() );
 843
 844   tr() = tr1();
 845 //  _corrector = &ADP_NODE::tr_correct_1_const;
 846   // tr() = tt_integrate_1_const(_delta_expect);
 847
 848
 849   assert (!_positive || tt_expect >=0);
 850 }
 851 /*---------------------------------*/
 852 double ADP_NODE::tt_integrate_1_exp( double tr_ ){
 853   // FIXME implement quadratic model
 854   return( tt_integrate_1_const(tr_) );
 855 }
 856 /*---------------------------------*/
 857 double ADP_NODE::tt_integrate_1_const( double tr_){
 858   _sign = 1;
 859   // _debug=3;
 860   double delta = ( tr1() + tr_ )/2 *  (dT0() / tr_duration() - 1) ;
 861 //   double ret = tt1() + delta + tr_; // RET
 862
 863   long double ret = get_aft_1() + delta; // RET
 864
 865   trace5(("ADP_NODE::tt_integrate_1_const" + label()).c_str(), tt1(),
 866       delta , tr_, ret, dT0() / tr_duration());
 867   trace1(("ADP_NODE::tt_integrate_1_const" + label()).c_str(),tr_noise);
 868   // trace1("ADP_NODE::tt_integrate_1_const", tt());
 869
 870   if (ret < 0 && _positive){
 871     error(bDANGER, "* ADP_NODE::tt_integrate_1_const neg error step %i,"
 872         " Time0=%f (%f), %s, tt_value = %g, ( %g, %g) \n",
 873         CKT_BASE::_sim->tt_iteration_number(),
 874         CKT_BASE::_sim->_Time0, dT0(),
 875         short_label().c_str(), tt_expect,
 876         _delta[0], _delta[1]);
 877 //    assert(false);
 878      ret=0;
 879   }
 880
 881   if (_delta[1] < 0 && tr_ < 0 ){
 882     _sign = -1;
 883   }
 884
 885   _integrator = &ADP_NODE::tt_integrate_1_const;
 886   assert (!_positive || ret >=0);
 887   return( (double) ret);
 888 }
 889 /*---------------------------------*/
 890 double ADP_NODE::tr_correct_1_const(){
 891   trace0(("ADP_NODE::tr_correct_1_const correction " + label()).c_str());
 892   return  (_delta_expect + tr_value)/2;
 893 }
 894 /*---------------------------------*/
 895 double ADP_NODE::tt_integrate_( double ){
 896 // take order and tr_value
 897 // integrate things. write to tt_value.
 898   return 0;
 899 }
 900 /*---------------------------------*/
 901 /*----------------------------------------------------------------------------*/
 902 //  tt1() += diff;
 903 void ADP_NODE::tt_commit_first( )
 904 {
 905   trace2(("ADP_NODE::tt_commit_first() " + short_label()).c_str(), tt(), tr());
 906   assert (order()==0);
 907   _integrator = 0;
 908   _corrector = 0;
 909
 910   tt_first_time = CKT_BASE::_sim->_Time0;
 911   assert(_delta_expect != HUGE_VAL );
 912 }
 913 /*---------------------------------*/
 914 void ADP_NODE::tt_commit( )
 915 {
 916   assert(0);
 917  //  called before CARD_LIST apply...
 918  //  extrapolates tr values for guessed dT0
 919  //
 920  //   FIXME: cleanup.
 921
 922   if (order()==0)
 923     tt_commit_first();
 924
 925   trace3(( "ADP_NODE::tt_commit "+short_label()).c_str(), tt_value, tr(),
 926      CKT_BASE::_sim->tt_iteration_number() );
 927   trace6(("ADP_NODE::tt_commit " + short_label()).c_str(),
 928       dT0(), dT1(),
 929       CKT_BASE::_sim->get_tt_order(),
 930       CKT_BASE::_sim->_Time0,
 931       tt(), tt1()
 932       );
 933
 934   // sets _delta_expect and tt_expect
 935   if(!is_number( tr1()) && CKT_BASE::_sim->tt_iteration_number()>=2 )
 936   {
 937     error(bDANGER,"ADP_NODE::tt_commit history broken %s, %i, step %i\n", long_label().c_str(), m_(), _sim->tt_iteration_number());
 938     assert( is_number( tr1()) || CKT_BASE::_sim->tt_iteration_number()<2   );
 939   }
 940   tr_expect_();
 941
 942
 943   // this is before sweep....
 944   assert( CKT_BASE::_sim->last_time() >0 );
 945   assert( tt_value == tt_value );
 946   assert( tt_value >=0 || !_positive );
 947   // tr_value will be printed as "tr"
 948
 949   if (tt_expect < -1e-8 ){
 950     //positive?
 951     error(bDANGER, "* ADP_NODE::tt_commit neg error step %i, order %i, tt_expect=%e\n", \
 952         CKT_BASE::_sim->tt_iteration_number(), _order,tt_expect);
 953   }
 954 }
 955 /*---------------------------------*/
 956 // tr relative to Time1
 957 hp_float_t ADP_NODE::tr_rel( double dT ) const{
 958   long double now_rel = dT;
 959
 960   switch(order()){
 961     case 0: // transient sim
 962       if(_sim->analysis_is_tran() ){
 963         return tr();
 964       }else{incomplete();
 965         return 0;
 966       }
 967     case 1:
 968       return tr1();
 969     case 2:
 970       trace1("ADP_NODE::tr 2", now_rel/ dT1());
 971       return tr1() + ( tr1() - tr2()) * double ((now_rel )/(long double) dT1());
 972     case 3:
 973       return  -(((tr2() - tr3())/dT2() - (tr1() - tr2())/dT1())*(hp_float_t(now_rel) + dT1())/(dT1() + dT2()) -
 974         (tr2() - tr3())/dT2())*(hp_float_t(now_rel) + dT1() + dT2()) + tr3();
 975     default:
 976       assert(false);
 977   }
 978   return(0);
 979 }
 980 /*---------------------------------*/
 981 // evaluate parabola through (y1-dy1-dy2, x3), (y1-dy1,x2) , (y1,x1) at y
 982 static double interpoly2(double x3, double x2, double x1, double dy2, double dy1, double y1, double y)
 983 {
 984   double rel = y - y1;
 985
 986   return 0;
 987   return  -(((x2 - x3)/dy2 - (x1 - x2)/dy1)*((rel) + dy1)/(dy1 + dy2) -
 988       (x2 - x3)/dy2)*(rel + dy1 + dy2) + x3;
 989
 990 }
 991 // order:
 992 // 0 during 1st timeframe
 993 // 1 during 2nd timeframe.
 994 double ADP_NODE::tr(double Time, const double* guess)const
 995 {
 996   long double Time1 = Time0() - dT0();
 997   long double now_rel = Time - Time1; // dT0, in case Time=Time0;
 998
 999   if (guess) {
1000     switch (order()) {
1001       case 0: untested(); unreachable();
1002         assert( _sim->analysis_is_tran() );
1003         return *guess;
1004       case 1:
1005         return tr1() + ( *guess - tr1()) * double((now_rel )/(long double) dT0());
1006       case 2: incomplete();
1007               // return interpoly3
1008       default:
1009         return(NAN);
1010     }
1011   }
1012
1013   switch (order()) {
1014     case 0:
1015       // transient sim
1016       assert( _sim->analysis_is_tran() );
1017       return tr();
1018     case 1:
1019       assert(is_number(tr1()));
1020       return tr1();
1021     case 2:
1022       trace1("ADP_NODE::tr 2", now_rel/ dT1());
1023       return tr1() + ( tr1() - tr2()) * double((now_rel )/(long double) dT1());
1024     case 3:
1025       return  -(((tr2() - tr3())/dT2() - (tr1() - tr2())/dT1())*(double(now_rel) + dT1())/(dT1() + dT2()) -
1026         (tr2() - tr3())/dT2())*double((now_rel) + dT1() + dT2()) + tr3();
1027       return interpoly2(tr3(), tr2(), tr1(), dT2(), dT1(), (double) Time1, (double)Time);
1028     default: unreachable();
1029       return 0;
1030   }
1031 }
1032 /*---------------------------------*/
1033 void ADP_NODE::tr_expect_()
1034 {
1035   _integrator = 0;
1036   _corrector = 0;
1037   _order = min(_order+1,CKT_BASE::_sim->get_tt_order());
1038
1039   trace2("ADP_NODE::tr_expect_", _order, CKT_BASE::_sim->tt_iteration_number());
1040   trace4("ADP_NODE::tr_expect_ ", _sim->_Time0, tr(_sim->_Time0), tr1(), tr2() );
1041   assert(_order <= CKT_BASE::_sim->tt_iteration_number());
1042   assert ( isnan(tr()) || order()==0 );
1043
1044   switch(_order){
1045     case 0:
1046       trace0("ADP_NODE::tr_expect_, not doing anything");
1047       return;
1048     case 2:
1049       trace3(("ADP_NODE::tr_expect_ extradebug" + long_label()).c_str(), tr2(), tr1(), tr(Time0()) );
1050       if ( fabs(tr2()-tr_rel(-dT1())) > 1e-6   && !(is_almost(tr2(), tr(Time2())))){
1051         error( bDANGER, "ADP_NODE::tr_expect_ tt_iteration_number %i\n", tt_iteration_number());
1052         error( bDANGER, "ADP_NODE::tr_expect_ mismatch, T0: %.20E, T1 %.20E, Time2 %.20E\n", Time0(), Time1(), Time2());
1053         error( bDANGER, "ADP_NODE::tr_expect_ mismatch dT1() %E dT0() %E\n", dT1(), dT0() );
1054         error( bDANGER, "ADP_NODE::tr_expect_ mismatch, tr1:   %E, tr2: %E d %E  \n", tr1(), tr2(), tr1()-tr2() );
1055         error( bDANGER, "ADP_NODE::tr_expect_ mismatch, tr2():tr(Time2()=%E))= %E : %E\n", Time2(), tr2(), tr(Time2()) );
1056         error( bDANGER, "ADP_NODE::tr_expect_ delta, tr2()-tr(Time2()))= %E\n", tr2()- tr(Time2()) );
1057         error( bDANGER, "ADP_NODE::tr_expect_ mismatch in %s\n", long_label().c_str() );
1058
1059         throw(Exception(" mismatch " + long_label()));
1060         assert( false);
1061       }
1062     case 1:
1063       if (! is_almost(tr1(), tr(Time1()))) {
1064         error( bDANGER, "ADP_NODE::tr_sanity check failed %1.19g %1.19g \n",tr1(),tr(Time1()));
1065         untested();
1066       }
1067       // assert(is_almost(tr1(), tr(Time1())));
1068       //tr()=tr(_sim->_Time0);
1069       tr()=tr_rel(dT0());
1070
1071       break;
1072     case 3:
1073       // FIXME assert( tr_dd12 == ( _delta[1] - _delta[2] ) / dT1());
1074       tr_dd12 = ( tr1() - tr2() ) / dT1();
1075       tr()=tr(_sim->_Time0);
1076
1077       break;
1078     default:
1079       assert(false);
1080   }
1081   trace6( "ADP_NODE::tr_expect_ done ", tr1(), _delta_expect, _order,
1082       tt(), tt1(), _user_number);
1083   assert(tt1() >=0 || !_positive);
1084
1085   // assert(is_number(_tt[0]));
1086
1087 //  return _delta_expect;
1088  // tt_value  << " " << get1() << " order " << _order << " val1 " <<
1089  // tt1() <<
1090  //  "ttex " << tt_expect <<  "\n";
1091 }
1092 /*-----------------------------------------*/
1093 int ADP_NODE::region() const{
1094         return _region;
1095 }
1096 /*-----------------------------------------*/
1097 void ADP_NODE::tr_expect_1_exp(){
1098   // constant
1099   assert(order()==1);
1100   _sign=1;
1101   if(_delta[1] < 0) _sign=-1;
1102   trace3(( "ADP_NODE::tt_expect_1_exp() "+short_label()).c_str(),
1103      tr_value, _delta[0], _delta[1] );
1104   trace3(( "ADP_NODE::tt_expect_1_exp() "+short_label()).c_str(),
1105       Time_delta(), CKT_BASE::_sim->_time0, CKT_BASE::_sim->last_time() );
1106
1107   _delta_expect = _delta[1];
1108   assert (_delta_expect < 0.1);
1109
1110   _corrector = &ADP_NODE::tr_correct_1_exp;
1111   tt_expect = tt_integrate_1_const(_delta_expect);
1112
1113   _integrator = &ADP_NODE::tt_integrate_1_exp;
1114 }
1115 /*---------------------------------*/
1116 void ADP_NODE::tr_expect_2_exp(){
1117   trace2("ADP_NODE::tr_expect_2_exp()", _delta[1], _delta[2]);
1118   double h = tr_duration();
1119   // _debug+=1000;
1120   assert( order() == 2 );
1121   assert( is_number(_delta[1]) );
1122   assert( is_number(_delta[2]) );
1123   _sign = 1;
1124   if (_delta[1] <= 0 && _delta[2] < 0 ){
1125     _sign=-1;
1126   }
1127
1128   // need total ordering.
1129   if ( !( ( _delta[2] <= _delta[1] && _delta[1] <= 0)
1130        || ( _delta[2] >= _delta[1] && _delta[1] >= 0))){
1131 //    _order--;
1132     trace2(("ADP_NODE::tr_expect_2_exp no total ord linear... " +
1133           short_label()).c_str(),  _delta[2], _delta[1]);
1134     incomplete();
1135 //    _debug+=5000;
1136 //    tr_expect_2_linear();
1137     _order--;
1138     tr_expect_1_exp();
1139     _corrector = &ADP_NODE::tr_correct_1_exp;
1140     return;
1141   }
1142   if ( fabs(_delta[1] - _delta[2]) < tr_noise ){
1143     trace3(("ADP_NODE::tr_expect_2_exp just noise " +
1144           short_label()).c_str(),  _delta[2], _delta[1], tr_noise);
1145     incomplete();
1146     _order--;
1147     tr_expect_1_exp();
1148   //  _corrector = &ADP_NODE::tr_correct_1_exp;
1149     return;
1150   }
1151   //hp_float_t t2 = tr_duration()/2;
1152   //hp_float_t t1 = t2+dT1();
1153   //dashier geht nicht. kleine schritte zu schlecht.
1154   double d1 = _delta[1];
1155   double d2 = _delta[2];
1156
1157   d1/=h;
1158   d2/=h;
1159   // double alpha = - log ( d2/d1 ) / dT1();
1160
1161   //d1=_der_aft[1]; // dashier muesste man noch ausrechnen.
1162   //
1163   //
1164   // expect using deltas as derivatives.
1165   d1=_delta[1] / h;
1166   d2=_delta[2] / h;
1167   //  d2=_der_aft[2];
1168
1169   // alpha = - log ( d2/d1 ) / dT1();
1170   double q =  ( d1/d2 ) ;
1171
1172   if ( _sign*d1 <= _sign*d2 && q != q ) q=1;
1173
1174  //expect using deltas.
1175   d1=_delta[1] / h;
1176   d2=_delta[2] / h;
1177
1178  // _delta_expect = pow(d1,((dT0() + dT1())/dT1())) / pow(d2,(dT0()/dT1()))*h;
1179
1180   _delta_expect = pow(q, (dT0()/dT1()))*d1*h;
1181
1182   trace6("ADP_NODE::tr_expect_2_exp", _delta_expect - _delta[1], _der_aft[2], _der_aft[1],
1183       _delta[1], _delta[2], dT0());
1184   trace6("ADP_NODE::tr_expect_2_exp", _delta_expect, _der_aft[2], _der_aft[1],
1185       _delta[1] - _delta[2], _delta_expect, _sign);
1186   //trace6("ADP_NODE::tr_expect_2_exp", _delta[1], _delta[2], _delta_expect, p, r, q);
1187
1188   assert(_delta_expect == _delta_expect);
1189   assert(_delta_expect < 1e6);
1190   _corrector = &ADP_NODE::tr_correct_1_exp;
1191   assert (_sign*_delta[1] >= 0 );
1192
1193   assert( 1e-15 >= _sign*(_delta_expect - _delta[1]) );
1194
1195   tt_expect = tt_integrate_2_exp(_delta_expect);
1196 }
1197 /*---------------------------------*/
1198 double ADP_NODE::tr_correct_1_exp(){
1199
1200   assert(order()!=0);
1201   trace0(("ADP_NODE::tr_correct_1_exp correction " + label()).c_str());
1202
1203   double new_delta; // return value. new value for tr_value.
1204   double h = tr_duration();
1205
1206   //assert(tt1() == tt1());
1207
1208   long double a = _delta[1];
1209   //long double c = tr_value;
1210   long double b = tt_expect;
1211   long double d = tt1();
1212   long double dT = dT0();
1213
1214   trace4(("ADP_NODE::tr_correct_1_exp, _der_aft? " + label()).c_str(),
1215       _delta[1], tt_expect, tr_value, tt1() );
1216
1217   //model value as (does it make sense?)
1218   /// t \mapsto K* ( 1- exp (tau*t) )+ d
1219   // long double K = (b-d)/(c/a-1);
1220   // long double tau = log((c-a)/(b-d) + 1) / h;
1221
1222   //_delta_expect = K *(- exp(tau* (dT0()+h) )  + exp ((tau*dT0())) );
1223
1224   // a*(a - b - c + d)^(dT0/h)/(-b + d)^(dT0/h)
1225   //
1226
1227   long double d_b = d-b;
1228   long double a_c = _delta[1] - tr_value;
1229
1230
1231   if (  fabs ( a_c ) < tr_noise )  {
1232     trace2("ADP_NODE::tr_correct_1_exp ", a_c, tr_noise );
1233     _integrator=&ADP_NODE::tt_integrate_1_const;
1234     return(tr_value + _delta[1])/2;
1235   }
1236
1237   if (  fabs ( d_b )  < 1e-20 )  {
1238     trace6("ADP_NODE::tr_correct_1_exp ", _delta[1], tr_value, tt1(),
1239         tt_expect, 1e-13, fabs ( d_b / (fabs(d)+fabs(b))));
1240     _integrator=&ADP_NODE::tt_integrate_1_const;
1241     return(tr_value + _delta[1])/2;
1242   }
1243   if ( ( fabs ( d_b ) / (fabs(d)+fabs(b) ) < 1e-13 ) ) {
1244     trace6("ADP_NODE::tr_correct_1_exp ", _delta[1], tr_value, tt1(),
1245         tt_expect, 1e-13, fabs ( d_b / (fabs(d)+fabs(b))));
1246     _integrator=&ADP_NODE::tt_integrate_1_const;
1247     return(tr_value + _delta[1])/2;
1248   }
1249   long double B = (1 + (a_c)/(d_b));
1250   long double gain = ( dT/h );
1251   assert(is_number(gain));
1252
1253   if (B<0){
1254     B=0;
1255     untested();
1256   }
1257   long double P = pow(B,gain);
1258
1259   if ( !is_number((double)P) ) {
1260     trace0("ADP_NODE::tr_correct_1_exp P is not number. probably just noise");
1261     return(tr_value + _delta[1])/2;
1262   }
1263
1264 //  new_delta = (double)( a*pow((a - c + d-b)/(d-b),(dT/h)) );
1265   new_delta = (double)( a*P  );
1266
1267   trace6(("ADP_NODE::tr_correct_1_exp, " + label()).c_str(), \
1268       a, tr_value, P, dT0(), B, gain );
1269   trace6(("ADP_NODE::tr_correct_1_exp, " + label()).c_str(), tt_expect,
1270       tt1(), tt_value, _delta[1], tr_value, new_delta);
1271   trace4( "ADP_NODE::tr_correct_1_exp " , CKT_BASE::_sim->tt_iteration_number(), d_b, a_c,
1272       CKT_BASE::_sim->_dT0/CKT_BASE::_sim->_dTmin );
1273
1274   assert(is_number(B));
1275   assert(is_number(P));
1276   assert(is_number(new_delta));
1277   assert(new_delta * _delta[1] >= 0);
1278   assert(new_delta<1);
1279
1280   if( (_delta[1] >= new_delta && new_delta >= 0)
1281       ||  (_delta[1] <= new_delta && new_delta <= 0)){
1282   } else {
1283     _integrator=&ADP_NODE::tt_integrate_1_const;
1284   }
1285
1286   return new_delta;
1287 }
1288 /*----------------------------------------------------------------------*/
1289 double ADP_NODE::tt_integrate_2_exp(double tr_) {
1290   trace4("ADP_NODE::tt_integrate_2_exp", _delta[1], tr_ , _sign, tr_value);
1291   trace4("ADP_NODE::tt_integrate_2_exp", tr_duration(), dT0(), tt1() +
1292       _delta[1], tt1());
1293   hp_float_t h = tr_duration();
1294
1295   assert (_sign * tr_       >= 0 );
1296   assert (_sign * _delta[1] >= 0 );
1297   trace4( "ADP_NODE::tt_integrate_2_exp ", _delta[1],  tr_, _delta[1]- tr_, _sign);
1298
1299   if (_sign*_delta[1] < _sign*tr_) {
1300     untested();
1301   }
1302
1303   long double ln1 = log(_sign*_delta[1]/h);
1304   long double ln0 = log(_sign*tr_/h);
1305
1306   //int_{h}^{dT0} exp( ln1 + v/dT0 * (ln0 - ln1) ) dv
1307   long double delta = dT0()/(ln0-ln1) * ( _sign*tr_/h - exp( ln1 + (h/dT0())
1308         *(ln0-ln1) ) ) * _sign;
1309
1310
1311   //int_{h/2}^{dT0-h/2} exp( ln1 + v/dT0 * (ln0 - ln1) ) dv
1312   delta = dT0()/(ln0-ln1) * ( exp( ln1 + (1-h/2/dT0()) * (ln0-ln1) )
1313                             - exp( ln1 + (  h/2/dT0()) * (ln0-ln1) ) ) * _sign;
1314
1315   if (delta != delta || delta==-inf || delta==inf){
1316     return ( tt_integrate_1_const( tr_) );
1317     trace0( "FIXME" ); // call another thiongh.
1318   }
1319
1320   assert(delta != -inf);
1321   assert(delta != inf);
1322
1323   long double ret = get1() + _delta[1] + delta; // RET
1324
1325   if ( h == dT0() && fabs(delta)> 1e-14 ){
1326     trace2("ADP_NODE::tt_integrate_1_exp 0", _delta[1], delta);
1327     assert(false); // throw something?
1328   }
1329
1330   trace6(( "ADP_NODE::tt_integrate_1_exp "+short_label()).c_str(),
1331       _delta[2], _delta[1], tr_, tt_expect, delta, ret );
1332
1333   assert(ret == ret);
1334
1335   assert (_sign*tr_ >= 0 );
1336   _integrator = &ADP_NODE::tt_integrate_1_exp;
1337
1338   if( _positive && ret < 0){
1339         ret=0;
1340   }
1341
1342   assert(tt_expect >=0 || !_positive);
1343   assert(ret >=0 || !_positive);
1344   assert(ret<2); // RCD?
1345
1346   return double(ret);
1347 }
1348 /*----------------------------------------------------------------------*/
1349 // ADP_CARD::ADP_CARD() {unreachable();}
1350 /*----------------------------------------------------------------------*/
1351 ADP_CARD::ADP_CARD(const ADP_CARD& p) : COMPONENT(p) {unreachable();}
1352 /*----------------------------------------------------------------------*/
1353
1354 TIME_PAIR ADP_NODE_UDC::tt_review( ) {
1355   return ADP_NODE::tt_review();
1356 }
1357 /*----------------------------------------------------------------------*/
1358 TIME_PAIR ADP_NODE::tt_review( ) {
1359   hp_float_t myreltol = OPT::adpreltol;
1360   hp_float_t myabstol = OPT::adpabstol;
1361   double h = tr_duration();
1362   double delta_model;
1363   //fixme: where corrector?
1364   //if (_corrector){
1365   //  assert(order()>0);
1366   //  trace1(("ADP_NODE::tt_review: correction " + label()).c_str(), _delta_expect);
1367   //  delta_model = (this->*_corrector)(); // value predicted by model.
1368   //} else {
1369   delta_model = 0;
1370   //  trace1(("ADP_NODE::tt_review: no corrector " + label()).c_str(), delta_model );
1371   //}
1372
1373   trace3("ADP_NODE::tt_review", CKT_BASE::_sim->tt_iteration_number(), myabstol, myreltol);
1374   assert(delta_model == delta_model);
1375   if( myreltol == 0 ) {_wdT=0; return TIME_PAIR(0,0); }
1376   if( myabstol == 0 ) {_wdT=0; return TIME_PAIR(0,0); }
1377
1378   if ( ( tr_value * delta_model ) < 0 ) {
1379     trace2(("ADP_NODE::tt_review: oops, sign has changed "+ label()).c_str(), tr_value,delta_model);
1380   }
1381   // FIXME: implement abs err!
1382   double tr_high = tr_value + tr_noise;
1383   double tr_low = tr_value - tr_noise;
1384
1385   if ( delta_model > tr_high ) {
1386     _abs_tr_err = tr_high - delta_model;
1387   }  else if (delta_model < tr_low ){
1388     _abs_tr_err = delta_model - tr_low;
1389   }else {
1390     _abs_tr_err = fabs(delta_model - tr_value);
1391   }
1392
1393   trace2("ADP_NODE::tt_review", long_label(), tr_noise);
1394
1395   if (tr_noise<0 || !is_number(tr_noise)){
1396     error(bDANGER, "ADP_NODE::tt_review noise error: %s\n", long_label().c_str());
1397     assert (tr_noise >=0);
1398   }
1399
1400   _abs_tr_err = fabs (tr_value - delta_model);
1401   _abs_tr_err = max ( fabs (tr_value - delta_model) - tr_noise ,0.0 );
1402 //  _abs_tr_err = fabs (tr_value - delta_model) / tr_noise ;
1403
1404   if (_abs_tr_err < 1e-40 ) {
1405     _rel_tr_err = 0;
1406     _abs_tr_err = 0;
1407   } else {
1408     _rel_tr_err = _abs_tr_err; //  / max( fabs(tr_value) , fabs(delta_model));
1409
1410     // use noise for predictor only?
1411     // _rel_tr_err = fabs( tr_value-delta_model )/ fabs(tr_noise ) ;
1412   }
1413   _rel_tt_err = fabs (tt_value - tt_expect) / (fabs(tt_value) + fabs(tt_expect));
1414
1415   tr_value = delta_model;
1416
1417   trace3(("ADP_NODE::tt_review" + label() + "got tol").c_str(), myreltol, _abs_tr_err, _rel_tr_err);
1418   if( _abs_tr_err == 0 ) {
1419     trace1(("ADP_NODE::tt_review" + label() + "noabs").c_str(), TIME_PAIR(0, NEVER)._event );
1420     _wdT = inf; return TIME_PAIR(0, NEVER);
1421   }
1422   if( myreltol == inf && myabstol == inf ) { _wdT = inf; return TIME_PAIR();}
1423   if( myreltol == 0 && myabstol == 0 ) {_wdT=0; return TIME_PAIR(0,0);}
1424
1425 // FIXME: _order.
1426 //
1427   _wdT = dT0() * sqrt( myreltol / _rel_tr_err )  + h;
1428   _wdT = dT0() * log( 1+ myreltol / _rel_tr_err )  + h;
1429
1430
1431   if (_rel_tr_err >= 1 ){
1432     error( bDANGER, "stepdebug %i dev: %s reltr %f>1 model %E, tr_ %E fill %E\n",
1433         CKT_BASE::_sim->tt_iteration_number(),
1434         label().c_str(),   _rel_tr_err, delta_model, tr_value, tt_value );
1435   }
1436   if (_wdT < 1.1*h){
1437     error( bDANGER, "stepdebug  dev: _wdT %s wdT: %g %g %f rel_tr %f abserr: %g, corr %i ord %i\n",
1438         label().c_str(), _wdT, dT0(), myreltol,  _rel_tr_err, _abs_tr_err, (bool)_corrector, order() );
1439   }
1440
1441   trace4("ADP_NODE::tt_review adaptive _wdT", myreltol, inf, h, _rel_tr_err);
1442   if( _wdT < 0  || _wdT != _wdT) {
1443     error(bDANGER, " dev: _wdT %s %f %f %f\n", label().c_str(), _wdT, dT0(), _rel_tr_err );
1444   }
1445
1446   if( _rel_tr_err != _rel_tr_err ) _wdT = NEVER;
1447   if( !(_rel_tr_err == _rel_tr_err) ) _wdT = NEVER;
1448
1449   trace3( ("ADP_NODE::tt_review " + label()).c_str(), CKT_BASE::_sim->_Time0 , _wdT , myreltol );
1450
1451   assert(_wdT >= 0 );
1452
1453   assert(tr_value == tr_value);
1454   assert(_delta_expect == _delta_expect);
1455   tt_value += tr_value; // only printed if OPT::printrejected
1456
1457   return( TIME_PAIR(0, _wdT));
1458 }
1459
1460 void ADP_NODE_RCD::tr_expect_1(){ return tr_expect_1_exp();}
1461 void ADP_NODE_RCD::tr_expect_2(){ return tr_expect_2_exp();}
1462 void ADP_NODE_RCD::tr_expect_3(){ return tr_expect_3_exp();}
1463
1464 void ADP_NODE_UDC::tr_expect_1(){ return tr_expect_1_const();}
1465 void ADP_NODE_UDC::tr_expect_2(){ return tr_expect_2_linear();}
1466 void ADP_NODE_UDC::tr_expect_3(){ return tr_expect_3_quadratic();}
1467 /*----------------------------------------------------------------------------*/
1468 // vim:ts=8:sw=2:et: