modules/s_reach.cc

   1 /* Copyright (C) 2014 Felix Salfelder
   2  * Author: Felix Salfelder
   3  *
   4  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3, or (at your option)
   9  * any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  19  * 02110-1301, USA.
  20  *------------------------------------------------------------------
  21  * reachability stuff
  22  */
  23
  24 #include "config.h"
  25
  26 #ifdef HAVE_CLAPACK_H
  27 // #include "u_atlas.h"
  28 #endif
  29
  30
  31 #ifndef HAVE_CLAPACK_H
  32 #warning "untested"
  33 #endif
  34
  35 #include "u_status.h"
  36 #include "u_ssp.h"
  37 #include "s_ev.h"
  38 #include "u_gsl.h"
  39 #include "u_prblst.h"
  40 #include "u_cardst.h"
  41 #include "e_elemnt.h"
  42 #include "e_storag.h"
  43 #include "s_ev.h"
  44 #include "io_matrix.h"
  45 #include "m_matrix_extra.h"
  46 #include "e_aux.h"
  47 #include "s_tr.h"
  48 using namespace std;
  49 /*--------------------------------------------------------------------------*/
  50 /*--------------------------------------------------------------------------*/
  51 /*--------------------------------------------------------------------------*/
  52 namespace { //
  53 /*--------------------------------------------------------------------------*/
  54 class REACH : public EV_BASE, protected TRANSIENT
  55         { //
  56         public:
  57                 explicit REACH(): EV_BASE() {}
  58                 ~REACH() {}
  59                 void    do_it(CS&, CARD_LIST*);
  60         private:
  61                 void    setup(CS&);
  62                 explicit REACH(const REACH&): EV_BASE(), TRANSIENT() {unreachable(); incomplete();}
  63
  64                 COMMON_COMPONENT* _zap_bm[DCNEST]; /* pointer to thing to sweep, dc command */
  65                 void options_(CS&);
  66
  67         protected:
  68                 void sweep_recursive(int);
  69
  70                 void align_inf_eigenvectors(gsl_matrix_complex_view& EV_inf, size_t inputnumber=0);
  71                 unsigned sens_inf_vectors(gsl_matrix_complex_view& REACH);
  72
  73
  74                 void find_successors(index_t* swp, SSP_SPL& t, unsigned Nest, unsigned depth);
  75                 void find_successors(const SSP_VECTOR& inp, SSP_SPL& t, unsigned=0, unsigned depth=0){
  76                         index_t* swp = inp.clone_data();
  77                         find_successors(swp, t, _n_inputs, depth);
  78                         delete[] swp;
  79                 }
  80                 void find_successors(const gsl_vector* swp, SSP_SPL& t, unsigned Nest, unsigned depth){
  81                         SSP_VECTOR v(swp);
  82                         return find_successors(v, t, Nest, depth);
  83                 }
  84
  85                 void tran_step(double dt);
  86         private:
  87                 SSP_SPL* root_spl;
  88                 unsigned _depth;
  89                 unsigned _maxwidth;
  90
  91         private:
  92                 class SWEEP{
  93                         public:
  94                                 SWEEP(REACH& p, SSP_TRANS& t, unsigned depth) :
  95                                     _parent(p), _spl(t), _depth(depth) { itested();
  96                                         swp = t.input().clone_data();
  97
  98                                 }
  99                                 ~SWEEP(){}
 100
 101                                 const SSP_CHART& chart()const{return _spl.chart();}
 102                                 unsigned n_inputs()const{return _spl.chart().n_inputs();}
 103                                 void find(unsigned Nest);
 104                         private:
 105                                 double first(const SSP_SPL& c, index_t* swp, int Nest);
 106                                 bool next(const SSP_CHART& c, index_t* swp, int Nest) const;
 107                                 void do_it(int Nest, double dt, double bak);
 108                                 REACH &_parent;
 109                                 SSP_SPL& _spl;
 110                                 double* grid;
 111                                 index_t* swp;
 112                                 //index_t const* swp; ??
 113                                 unsigned _depth;
 114                                 index_t _last;
 115                 };
 116
 117 };
 118 /*--------------------------------------------------------------------------*/
 119 void REACH::SWEEP::do_it(int Nest, double dt, double bak)
 120 {
 121         size_t n = _sim->_total_nodes;
 122         gsl_vector_view v0__ = gsl_vector_view_array(_sim->_v0+1,n);
 123         _sim->tr_reset(); // required?
 124         assert(_sim->_time0 == 0.);
 125
 126         _sim->restore_voltages();
 127         CARD_LIST::card_list.tr_restore(); // done by TRANSIENT::sweep?
 128         trace7("transition", _spl.id(), Nest, Nest, _parent._zap_bm[Nest]->param_name(10), _parent._zap_bm[Nest]->param_value(10), swp[Nest], dt);
 129
 130         if(_parent.EV_BASE::_trace >= tVERBOSE){
 131                 _parent.EV_BASE::outdata(-double(_spl.id()));
 132         }
 133
 134         _parent.tran_step(dt);
 135         _parent.get_op_();
 136         double v0save[n];
 137         memcpy(v0save, _sim->_v0+1, n*sizeof(double));
 138
 139         if(_parent.EV_BASE::_trace >= tITERATION){
 140                 _parent.EV_BASE::outdata(1./0.);
 141         }
 142
 143         trace1("presolver", v0__);
 144         SSP_CHART* chart = NULL;
 145         SSP_STATE& state = _parent.ev_solver(chart);
 146         trace1("solved", v0__);
 147         gsl_vector_view v0bak__ = gsl_vector_view_array(v0save,n);
 148
 149         SSP_SPL& spl = state.insert_spl(v0__, chart);
 150         trace1("discretized", v0__);
 151
 152         gsl_vector_sub(v0bak__, v0__);
 153         double error = gsl_blas_dnrm2(&v0bak__.vector);
 154
 155         SSP_TRANS transition = _spl.insert_trans(swp, &spl);
 156
 157         if (spl.is_new()) { itested();
 158                 //                              _parent.EV_BASE::_out << double(spl.id()) << double(depth) << double(chart->id());
 159                 _parent.EV_BASE::_out << chart->id();
 160                 _parent.EV_BASE::outdata(double(spl.id()));
 161         }
 162
 163         _parent.EV_BASE::_out << _spl.chart().id() << "," << _spl.id() << " -> " << chart->id() << "," << spl.id();
 164         _parent.EV_BASE::_out << dt;
 165         for(unsigned i=0; i < n_inputs(); ++i){
 166                 _parent.EV_BASE::_out << double(swp[i]) * chart->grid(i);
 167         }
 168         assert(fabs(error)<10);
 169         _parent.EV_BASE::_out << " err " << error << "\n";
 170
 171         if (!spl.is_new()) { itested();
 172                 trace2("notnew", spl.id(), _depth);
 173         }else if(!_depth){ untested();
 174                 // border.
 175                 // hmm push to stack?
 176         }else{ itested();
 177                 trace2("isnew", spl.id(), swp);
 178                 _sim->keep_voltages(true);
 179                 SWEEP sweeper(_parent, transition, _depth-1);
 180                 sweeper.find(spl.n_inputs());
 181                 _sim->pop_voltages();
 182         }
 183
 184         trace2("resetting iv", Nest, bak);
 185         _parent._zap_bm[Nest]->set_param_by_name("iv", ::to_string(bak));
 186 }
 187 /*--------------------------------------------------------------------------*/
 188 // TODO: breadth-first search...?
 189 void REACH::SWEEP::find(unsigned Nest)
 190 { itested();
 191         --Nest;
 192
 193         trace3("find_successors", _spl.id(), Nest, swp[Nest]);
 194         double bak = first(_spl, swp, Nest);
 195         do { itested();
 196                 if (Nest == 0) {
 197
 198 //                      double dt = 1e-8;
 199 //                      do_it(Nest, dt, bak);
 200                         do_it(Nest, 1e-2, bak);
 201
 202
 203                 }else{ untested();
 204                         find(Nest);
 205                 }
 206         } while (next(_spl.chart(), swp, Nest));
 207 }
 208 /*--------------------------------------------------------------------------*/
 209 /*--------------------------------------------------------------------------*/
 210 void REACH::do_it(CS& Cmd, CARD_LIST* Scope)
 211 {
 212
 213         trace0("REACH::do_it(CS&, CARD_LIST*)");
 214
 215         trace0("doing ddc");
 216         EV_BASE::_scope = Scope;
 217         _sim->_time0 = 0.;
 218         _sim->set_command_ddc();
 219         _sim->_phase = p_INIT_DC;
 220         ::status.ddc.reset().start();
 221         _sim->_temp_c = temp_c_in;
 222         _do_tran_step = 0;
 223         EV_BASE::_dump_matrix = 0;
 224         EV_BASE::command_base(Cmd);
 225         ::status.ddc.stop();
 226 }
 227 /*--------------------------------------------------------------------------*/
 228 void REACH::setup(CS& Cmd)
 229 {
 230         EV_BASE::_cont = false;
 231         EV_BASE::_trace = tNONE;
 232         EV_BASE::_out = IO::mstdout;
 233         EV_BASE::_out.reset(); //BUG// don't know why this is needed */
 234         bool ploton = IO::plotset  &&  EV_BASE::plotlist().size() > 0;
 235         _n_inputs = 0;
 236
 237         if (Cmd.more()) {
 238                 for (_n_sweeps = 0; Cmd.more() && _n_sweeps < DCNEST; ++_n_sweeps) {
 239                         if (Cmd.umatch("to")) {
 240                                 break;
 241                         }
 242                         CARD_LIST::fat_iterator ci = findbranch(Cmd, &CARD_LIST::card_list);
 243                         if (!ci.is_end()) {
 244                                 if (ELEMENT* c = dynamic_cast<ELEMENT*>(*ci)) {
 245                                         _zap[_n_sweeps] = c;
 246                                         c->tr_begin();
 247                                         trace2("REACH::setup", c->long_label(), c->value());
 248                                 } else { untested();
 249                                         throw Exception("dc/op: can't sweep " + (**ci).long_label() + '\n');
 250                                 }
 251                                 ++_n_inputs;
 252                         }else if (Cmd.is_float()) {             // sweep the generator
 253                                 _zap[_n_sweeps] = NULL;
 254                         }else if (Cmd.is_alpha()) {
 255                                 string pname;
 256                                 unsigned here = Cmd.cursor();
 257                                 Cmd >> pname;
 258                                 try {
 259                                         _param[_n_sweeps] = &(EV_BASE::_scope->params()->find(pname));
 260                                 } catch(Exception) {
 261                                         Cmd.reset(here);
 262                                 }
 263                                 throw Exception("reach: can't sweep parameter " + pname + '\n');
 264                                 _zap[_n_sweeps] = NULL;
 265                         }else{ untested();
 266                                 // leave as it was .. repeat Cmd with no args
 267                         }
 268                         if (Cmd.match1("'\"({") || Cmd.is_float()) {
 269                                 _start[_n_sweeps] = "NA";
 270                                 _stop[_n_sweeps] = "NA";
 271                                 Cmd >> _start[_n_sweeps] >> _stop[_n_sweeps];
 272                                 _step[_n_sweeps] = 0.;
 273                         }else{
 274                                 // leave it as it was .. repeat Cmd with no args
 275                         }
 276                         _sim->_genout = 0.;
 277                         temp_c_in = OPT::temp_c;
 278                         _sim->_temp_c = temp_c_in;
 279                         EV_BASE::options(Cmd,_n_sweeps);
 280
 281                 }
 282                 unsigned here = Cmd.cursor();
 283                 string output;
 284                 bool newprobes = true;
 285                 do{ itested();
 286                         // FIXME: can we use the probe parser?!
 287                         unsigned here2 = Cmd.cursor();
 288                         if( Cmd >> "tran" ){ itested();
 289                                 Cmd.reset(here2);
 290                                 break;
 291                         }else if( Cmd >> "v" ){
 292                                 output_t t;
 293                                 t.brh[1] = 0;
 294                                 trace1("SENS::setup, have v", Cmd.tail());
 295
 296                                 int paren = Cmd.skip1b('(');
 297                                 Cmd >> output;
 298                                 t.label = "v("+output;
 299
 300                                 CKT_NODE* node = dynamic_cast<CKT_NODE*>((*EV_BASE::_scope).node(output));
 301                                 if(node)
 302                                         t.brh[0] = node;
 303                                 else{
 304                                         continue;
 305                                 }
 306
 307                                 trace2("SENS::setup, have", output, Cmd.tail());
 308                                 if(!Cmd.match1(')')){
 309                                         output = Cmd.ctos(")");
 310                                         trace1("SENS:skip1b:setup, have2", output);
 311                                         t.label += "," + output;
 312                                         node = dynamic_cast<CKT_NODE*>((*EV_BASE::_scope).node(output));
 313                                         if(node)
 314                                                 t.brh[1] = node;
 315                                         else{
 316                                                 Cmd.warn(bWARNING, "probelist: what's this?");
 317                                         }
 318
 319                                 }else{
 320                                         trace0("SENS::setup no comma");
 321                                 }
 322
 323                                 paren -= Cmd.skip1b(')');
 324                                 if (paren != 0) {untested();
 325                                         Cmd.warn(bWARNING, "need )");
 326                                 }else if (output.empty()) {untested();
 327                                         Cmd.warn(bWARNING, "probelist: what's this?");
 328                                 }else{
 329                                 }
 330
 331                                 t.label+=")";
 332                                 if (newprobes){
 333                                         _output.clear();
 334                                         newprobes = false;
 335                                 }
 336                                 _output.push_back(t);
 337                         }
 338                         //try{
 339                         //  _output.add_list(Cmd);
 340                         //}
 341                         //catch(Exception_Cant_Find)
 342                         //{}
 343                         ONE_OF
 344                                 || Get(Cmd, "dm",                 &(EV_BASE::_dump_matrix))
 345                                 || EV_BASE::_out.outset(Cmd);
 346                         ;
 347                 }while (Cmd.more() && !Cmd.stuck(&here));
 348         }else{ untested();
 349         }
 350
 351         TRANSIENT::_scope = EV_BASE::_scope;
 352         options_(Cmd);
 353         Cmd.check(bWARNING, "what's this?");
 354
 355         IO::plotout = (ploton) ? IO::mstdout : OMSTREAM();
 356         initio(EV_BASE::_out);
 357
 358         assert(_n_sweeps > 0);
 359         for (int ii = 0;  ii < _n_sweeps;  ++ii) { untested();
 360                 trace1("REACH setup", ii);
 361                 _start[ii].e_val(0., EV_BASE::_scope);
 362                 fix_args(ii);
 363
 364                 if (_zap[ii]) {
 365                         trace2("zap inc_probes" + _zap[ii]->long_label(), ii, _zap[ii]->value());
 366                         _stash[ii] = _zap[ii];                  // stash the std value
 367                         _zap[ii]->inc_probes();                 // we need to keep track of it
 368
 369                         // urghs. hack
 370                         STORAGE* s = dynamic_cast<STORAGE*>(_zap[ii]);
 371                         if (s) { unreachable(); incomplete();
 372                                 // trace2("REACH::setup ", _zap[ii]->long_label(), _zap[ii]->is_constant());
 373                                 // _zap[ii]->set_constant(false);                  // so it will be updated
 374                                 // _pushel[ii] = _zap[ii];                 // point to value to patch
 375                                 // _uic_caplist.push_back(s);
 376                                 // _sweepval[ii] = s->set__ic();
 377                         }else{
 378                                 COMMON_COMPONENT* pulse = bm_dispatcher.clone("pulse");
 379                                 trace2("REACH::setup", ii, hp(pulse));
 380                                 pulse->set_param_by_name("iv", "0"); // BUG?.
 381                                 pulse->set_param_by_name("pv", "0"); // BUG?.
 382                                 _zap_bm[ii] = pulse;
 383                                 _zap[ii]->set_value(_zap[ii]->value(), pulse);
 384                                 assert(_zap[ii]->common() == pulse);
 385                                 _zap[ii]->set_constant(false);                 // so it will be updated
 386                                 _sweepval[ii] = _zap[ii]->set__value();
 387                         }
 388                 } else if (_param[ii]) {
 389                         _sweepval[ii] = _param[ii]->pointer_hack();
 390                 } else {
 391                         trace1("REACH::setup does not exist", ii);
 392                         //_sweepval[ii] = 0;
 393                         _sweepval[ii] = &_sim->_genout;
 394                         // throw(Exception("something went wrong\n"));
 395                 }
 396         }
 397         _sim->_freq = 0;
 398
 399         untested();
 400 }
 401 /*--------------------------------------------------------------------------*/
 402 // unnested options
 403 void REACH::options_(CS& Cmd)
 404 {
 405         trace1("REACH::options(CS&, int)", Cmd.tail());
 406         _quantize_states = false;
 407         _depth = -1;
 408         _maxwidth = 10;
 409
 410         _sim->_uic = false;
 411         unsigned here = Cmd.cursor();
 412         bool tr=false;
 413         do{ itested();
 414                 ONE_OF
 415                         ||(Get(Cmd, "uic",            &_uic))
 416                         || Get(Cmd, "sr",             &_slew_rate)
 417                         || Get(Cmd, "dt{emp}",        &temp_c_in,   mOFFSET, OPT::temp_c)
 418                         || Get(Cmd, "te{mperature}",  &temp_c_in)
 419                         || Get(Cmd, "de{pth}",        &_depth)
 420                         || Get(Cmd, "max{width}",     &_maxwidth)
 421                         || (Get(Cmd, "qu{antize}",    &_quantize_states))
 422                         || (Cmd.umatch("tr{ace} {=}") &&
 423                                         (ONE_OF
 424                                          || Set(Cmd, "n{one}",      &(EV_BASE::_trace), tNONE)
 425                                          || Set(Cmd, "o{ff}",       &(EV_BASE::_trace), tNONE)
 426                                          || Set(Cmd, "w{arnings}",  &(EV_BASE::_trace), tUNDER)
 427                                          || Set(Cmd, "i{terations}",&(EV_BASE::_trace), tITERATION)
 428                                          || Set(Cmd, "v{erbose}",   &(EV_BASE::_trace), tVERBOSE)
 429                                          || Cmd.warn(bWARNING,
 430                                                  "need none, off, warnings, iterations, verbose")
 431                                         )
 432                                 )
 433                         // || ( Get(Cmd , "tran",      &tr) &&
 434                         || ( Cmd >> "tran" &&
 435                                         (trace1("TR", Cmd.tail()), true) &&
 436                                         (TRANSIENT::options(Cmd), tr=true)
 437                                 )
 438
 439                         || (trace1("options_", Cmd.tail()), EV_BASE::_out.outset(Cmd));
 440         }while (Cmd.more() && !Cmd.stuck(&here) && !tr);
 441
 442         if (!tr){untested();
 443                 CS c(CS::_STRING,"");
 444                 TRANSIENT::options(c);
 445         }
 446 }
 447 /*--------------------------------------------------------------------------*/
 448 static REACH p2;
 449 EV_BASE* p3 = dynamic_cast<EV_BASE*>(&p2);
 450 static DISPATCHER<CMD>::INSTALL d2(&command_dispatcher, "reach", p3);
 451 /*--------------------------------------------------------------------------*/
 452 void REACH::sweep_recursive(int Nest)
 453 { untested();
 454         size_t n = _sim->_total_nodes;
 455         trace2("REACH::sweep_recursive(int)", Nest, _sim->_uic);
 456         //  unsigned n = _sim->_total_nodes;
 457
 458         --Nest;
 459         assert(Nest >= 0);
 460         assert(Nest < DCNEST);
 461
 462         // double iddc[d];
 463
 464         OPT::ITL itl = OPT::DCBIAS;
 465
 466         _sim->_uic = _uic;
 467
 468         CARD_LIST::card_list.precalc_last();
 469         double inputvals[_n_inputs+1];
 470
 471         for(unsigned i=0; i < _n_inputs; ++i) { untested();
 472                 inputvals[i+1] = _zap[i]->value();
 473                 trace2("setup pulse input", _zap[i]->long_label(), _zap[i]->value());
 474                 _zap_bm[i]->set_param_by_name("iv", ::to_string(inputvals[i+1]));
 475                 _zap_bm[i]->set_param_by_name("pv", ::to_string(inputvals[i+1]));
 476                 _zap[i]->precalc_first();
 477                 _zap[i]->precalc_last();
 478                 _zap[i]->tr_begin();
 479
 480                 if (ELEMENT* c = dynamic_cast<ELEMENT*>(_zap[i])) { untested();
 481                         c->set_constant(false);
 482                 }
 483         }
 484
 485         get_op(itl); // FIXME: only one root op.
 486                      // need to untangle input range from dc-sweep init points
 487         _sim->_uic = false;
 488
 489         if(EV_BASE::_trace >= tVERBOSE){
 490                 EV_BASE::outdata(1./0.);
 491         }
 492
 493         SSP_CHART* chart = NULL;
 494         SSP_STATE& state = ev_solver(chart);
 495         assert(chart);
 496         assert(chart->id()<10);
 497         gsl_vector_view v0__ = gsl_vector_view_array(_sim->_v0+1,n);
 498
 499         double probes[EV_BASE::printlist().size()];
 500
 501         size_t k = 0;
 502         for (PROBELIST::const_iterator
 503                         p=EV_BASE::printlist().begin();  p!=EV_BASE::printlist().end();  ++p) {
 504                 double v= (*p)->value();
 505                 probes[k++] = v;
 506         }
 507         USE(probes);
 508
 509         trace1("to", v0__);
 510         untested();
 511         SSP_CHART* ch = chart; // new SSP_CHART(_n_inputs);
 512         SSP_SPL& spl = state.insert_spl(v0__, ch);
 513         trace1("ins", v0__);
 514
 515         inputvals[_n_inputs] = 1./0.;
 516         for(unsigned i=0; i < _n_inputs; ++i){
 517                 inputvals[i] = *(_sweepval[i]);
 518                 ch->_inputgrid[i] = _step[i];
 519         }
 520         gsl_vector_view inp__ = gsl_vector_view_array(inputvals,_n_inputs+1);
 521         root_spl = &spl;
 522         SSP_TRANS transition = root_spl->insert_trans(&inp__.vector, &spl);
 523
 524         SSP_VECTOR inp(inputvals,_n_inputs+1);
 525
 526         if (spl.is_new()){
 527                 EV_BASE::_out << chart->id();
 528                 EV_BASE::outdata(spl.id());
 529         }
 530
 531         if (!_depth){ untested();
 532         }else   if (spl.is_new()){ untested();
 533                 trace1("isnew", v0__.vector);
 534                 _sim->keep_voltages(true);
 535                 SSP_VECTOR swp(transition.input());
 536                 SWEEP sweeper(*this, transition, _depth-1);
 537                 sweeper.find(spl.n_inputs());
 538                 _sim->pop_voltages();
 539         }else{ unreachable(); // for now.
 540                 trace1("notnew", v0__.vector);
 541         }
 542 }
 543 /*--------------------------------------------------------------------------*/
 544 //double REACH::first(SSP_VECTOR& t, int Nest) const
 545 double REACH::SWEEP::first(const SSP_SPL& s, index_t* swp, int Nest)
 546 { itested();
 547         assert(Nest >= 0);
 548         assert(Nest < DCNEST);
 549         const SSP_CHART& ch = s.chart();
 550         assert (unsigned(Nest) < ch.n_inputs());
 551
 552         double start = -.2; // FIXME!
 553         double stop = -start;
 554
 555         int width = _parent._maxwidth; // FIXME: must depend on SR
 556
 557         double from = ch.grid(Nest) * double(swp[Nest]);
 558
 559         index_t min_index = index_t(floor( (start+ch.grid(Nest)*.5) /(ch.grid(Nest))));
 560         index_t max_index = index_t(floor( (stop+ch.grid(Nest)*.5) /(ch.grid(Nest))));
 561         index_t here = swp[Nest];
 562
 563         assert(min_index<=max_index);
 564
 565         _last = min(max_index, here+width);
 566         swp[Nest] = max(min_index, here-width);
 567
 568         double to = ch.grid(Nest) * double(swp[Nest]);
 569         trace5("REACH::first, pulse", ch.id(), Nest, ch.grid(Nest), from, to);
 570         trace3("REACH::first, pulse", ch.id(), swp[Nest], _last);
 571         trace2("REACH::first, pulse", min_index, max_index);
 572
 573         _parent._zap_bm[Nest]->set_param_by_name("iv", ::to_string(from));
 574         _parent._zap_bm[Nest]->set_param_by_name("pv", ::to_string(to));
 575
 576         *(_parent._sweepval[Nest]) = to; // needed by inputsens scaling
 577
 578         return from;
 579 }
 580 /*--------------------------------------------------------------------------*/
 581 // FIXME: individual precisions
 582 // FIXME: individual inputno?
 583 // FIXME: break upon leaving chart.
 584 bool REACH::SWEEP::next(const SSP_CHART& c, index_t* swp, int Nest) const
 585 {
 586         ++swp[Nest];
 587         double p(c.grid(Nest));
 588         double q = double(swp[Nest]);
 589         assert(c.grid(Nest));
 590         double to = q*p;
 591         _parent._zap_bm[Nest]->set_param_by_name("pv", ::to_string(to));
 592         *(_parent._sweepval[Nest]) = to; // BUG? used by inputsens scaling
 593
 594         trace5("EV_BASE::next", Nest, swp[Nest], p, to, _last);
 595         return swp[Nest] <= _last;
 596 }
 597 /*--------------------------------------------------------------------------*/
 598 // take some eigenvectors (alpha!=0, beta==0)
 599 // rearrange considering input sensitivities.
 600 void REACH::align_inf_eigenvectors(gsl_matrix_complex_view& EV_inf, size_t i)
 601 {
 602         if (!_zap[i]){ untested();
 603                 return;
 604         }
 605         size_t n = _sim->_total_nodes;
 606         trace3("fixing\n", EV_inf.matrix, EV_inf.matrix.size1, EV_inf.matrix.size2);
 607
 608
 609
 610         trace1("lu", _sim->_acx);
 611
 612         assert(_zap[i]); // incomplete
 613         clear_sens();
 614         _zap[i]->sens_load("dummy");
 615         complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
 616         trace1("2 sensstamp", *sens_);
 617
 618         ::status.back.start();
 619         _sim->_acx.fbsub(_sim->_sens);
 620         ::status.back.stop();
 621
 622         complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
 623         trace2("sens", _zap[i]->long_label(), *sens_);
 624
 625         gsl_vector_complex_view first_EV = firstcol(EV_inf);
 626         gsl_vector_complex_memcpy(&first_EV.vector, sens_);
 627
 628
 629         if(EV_inf.matrix.size2 > 1) { untested();
 630                 EV_inf = gsl_matrix_complex_submatrix (&EV_inf.matrix, 0, 1, n, EV_inf.matrix.size2-1);
 631                 align_inf_eigenvectors(EV_inf, i+1);
 632         }
 633
 634 #if 0
 635         relevant_EV = gsl_matrix_complex_submatrix (EV_, 0, 0, n, fin_dim);
 636         gsl_vector_complex_view relevant_v0prime_ = gsl_vector_complex_subvector(v0prime_,0,fin_dim);
 637         trace1("preimg cut ", relevant_v0prime_.vector);
 638         gsl_matrix_complex_view relevant_EV = gsl_matrix_complex_submatrix (EV_, 0, 0, n, fin_dim);
 639         trace1("relevant EV\n", (relevant_EV.matrix));
 640         err+= gsl_blas_zgemv( CblasNoTrans, one, &relevant_EV.matrix, &relevant_v0prime_.vector, zero, ztmp_ );
 641         trace1("proj\n", *ztmp_);
 642
 643         //extra vector: normalize(v0 - ztmp)
 644         err+= gsl_vector_complex_sub(v0_, ztmp_);
 645         double len = gsl_blas_dznrm2(v0_);
 646         trace1("res\n", *v0_);
 647         gsl_vector_complex_view v = gsl_matrix_complex_column (EV_, fin_dim);
 648         gsl_vector_complex_memcpy (&v.vector, v0_);
 649
 650         //    gsl_vector_complex_set(relevant_v0prime.vector,0,len);
 651         gsl_complex c; GSL_REAL(c)=1./len; GSL_IMAG(c)=0;
 652         gsl_vector_complex_scale (&v.vector, c);
 653         fin_dim+= 1;
 654 #endif
 655
 656 }
 657 /*--------------------------------------------------------------------------*/
 658 // FIXME: this is ad-hoc.
 659 void REACH::tran_step(double dt)
 660 { itested();
 661
 662         for(unsigned i=0; i < _n_inputs; ++i) {
 663                 assert(_zap[i]);
 664                 assert(_zap_bm[i]);
 665                 _zap_bm[i]->set_param_by_name("rise", ::to_string(dt));
 666                 _zap[i]->precalc_last();
 667                 _zap[i]->tr_begin(); // hmmm?
 668         }
 669
 670         size_t n = _sim->_total_nodes;
 671         _sim->set_command_tran();
 672         _tstart = 0.;
 673         TRANSIENT::_time1 = 0.;
 674         TRANSIENT::_cont_dc = true;
 675         TRANSIENT::_trace = tALLTIME;
 676         _sim->_time0 = 0.;
 677         _sim->_dtmin = dt/10.;
 678         _dtmax = dt/2.;
 679         _tstop = dt;
 680         _tstrobe = dt;
 681         _sim->_freezetime = true;
 682         TRANSIENT::sweep();
 683         _sim->_freezetime = false;
 684         gsl_vector_view v0__ = gsl_vector_view_array(_sim->_v0+1,n);
 685         trace1("swept", v0__);
 686 }
 687 /*--------------------------------------------------------------------------*/
 688 }// namespace
 689 /*--------------------------------------------------------------------------*/
 690 /*--------------------------------------------------------------------------*/