src/m_wave.h

   1 /*                          -*- C++ -*-
   2  * Copyright (C) 2001 Albert Davis
   3  * Author: Albert Davis <aldavis@gnu.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  * "wave" class, for transmission lines and delays
  23  */
  24 //testing=script 2006.07.13
  25 #include "l_denoise.h"
  26 #include "m_interp.h"
  27 #include "io_misc.h"
  28 #include <boost/range/iterator_range.hpp>
  29 /*--------------------------------------------------------------------------*/
  30 class WAVE {
  31 private:
  32   std::deque<DPAIR> _w;
  33   double _delay;
  34 public:
  35   typedef std::deque<DPAIR>::iterator iterator;
  36   typedef std::deque<DPAIR>::const_iterator const_iterator;
  37
  38   explicit WAVE(double d=0);
  39   explicit WAVE(const WAVE&);
  40           ~WAVE() {}
  41   WAVE&    set_delay(double d);
  42   WAVE&    initialize();
  43   WAVE&    push(double t, double v);
  44   FPOLY1   v_out(double t)const;
  45   FPOLY1   operator()(double xx)const {return v_out(xx);}
  46   double   v_reflect(double t, double v_total)const;
  47   WAVE&    operator+=(const WAVE& x);
  48   WAVE&    operator+=(double x);
  49   WAVE&    operator*=(const WAVE& x);
  50   WAVE&    operator*=(double x);
  51   WAVE&    warp(double x);
  52   double   dhd_linear(const WAVE& x, std::pair<DPAIR, DPAIR>* where=NULL) const;
  53   double   dhd_discrete(const WAVE& x, std::pair<DPAIR, DPAIR>* where=NULL) const;
  54   const_iterator begin()const {return _w.begin();}
  55   const_iterator end()const {return _w.end();}
  56   size_t size()const {return _w.size();}
  57
  58 private:
  59   double segment2point(const_iterator l, const DPAIR& p, DPAIR* where=NULL) const;
  60   double norm2(const DPAIR&, const DPAIR&) const;
  61   double norm2(const DPAIR&) const;
  62   double point2chain(DPAIR p, const_iterator B, DPAIR* w=NULL) const; // FIXME: pass right end
  63   double edge2chain(DPAIR l, DPAIR r, const_iterator& L, const_iterator& R, std::pair<DPAIR, DPAIR>* where=NULL) const;
  64   FPOLY1 bisector_dual(const const_iterator& a, const const_iterator& b) const;
  65   double bisector_distance(const const_iterator& a, const const_iterator& b, DPAIR l, DPAIR r) const;
  66   double endpoints2chain(DPAIR l, DPAIR r, const_iterator& L, const_iterator& R, std::pair<DPAIR, DPAIR>* where=NULL) const;
  67 private:
  68   class RANGE {
  69   private:
  70     const_iterator _a;
  71     const_iterator _b;
  72   public:
  73     RANGE(const const_iterator& a, const const_iterator& b) : _a(a), _b(b){ untested();}
  74     const_iterator begin() const {return _a;}
  75     const_iterator end() const {return _b;}
  76   };
  77 };
  78 /*--------------------------------------------------------------------------*/
  79 /*--------------------------------------------------------------------------*/
  80 // push: insert a signal on the "input" end.
  81 // args: t = the time now
  82 //       v = the value to push
  83 //
  84 inline WAVE& WAVE::push(double t, double v)
  85 {
  86   _w.push_back(DPAIR(t+_delay, v));
  87   return *this;
  88 }
  89 /*--------------------------------------------------------------------------*/
  90 // initialize: remove all info, fill it with all 0.
  91 //
  92 inline WAVE& WAVE::initialize()
  93 {
  94   _w.clear();
  95   return *this;
  96 }
  97 /*--------------------------------------------------------------------------*/
  98 inline WAVE::WAVE(const WAVE& w)
  99   :_w(w._w),
 100    _delay(w._delay)
 101 {
 102 }
 103 /*--------------------------------------------------------------------------*/
 104 // constructor -- argument is the delay
 105 //
 106 inline WAVE::WAVE(double d)
 107   :_w(),
 108    _delay(d)
 109 {
 110   initialize();
 111 }
 112 /*--------------------------------------------------------------------------*/
 113 inline WAVE& WAVE::set_delay(double d)
 114 {
 115   _delay = d;
 116   return *this;
 117 }
 118 /*--------------------------------------------------------------------------*/
 119 // v_out: return the value at the "output" end
 120 // args: t = the time now
 121 //
 122 inline FPOLY1 WAVE::v_out(double t)const
 123 {
 124   return interpolate(_w.begin(), _w.end(), t, 0., 0.);
 125 }
 126 /*--------------------------------------------------------------------------*/
 127 // reflect: calculate a reflection
 128 // args: t = the time now
 129 //       v_total = actual voltage across the termination
 130 // returns: the value (voltage) to send back as the reflection
 131 //
 132 inline double WAVE::v_reflect(double t, double v_total)const
 133 {
 134   // return (v_total*2 - v_out(t)); // de-noised
 135   return dn_diff(v_total*2, v_out(t).f0);
 136 }
 137 /*--------------------------------------------------------------------------*/
 138 inline WAVE& WAVE::operator+=(const WAVE& x)
 139 {
 140   untested();
 141   for (std::deque<DPAIR>::iterator
 142          i = _w.begin(); i != _w.end(); ++i) {
 143     untested();
 144     (*i).second += x.v_out((*i).first).f0;
 145   }
 146   return *this;
 147 }
 148 /*--------------------------------------------------------------------------*/
 149 inline WAVE& WAVE::operator+=(double x)
 150 {
 151   untested();
 152   for (std::deque<DPAIR>::iterator
 153          i = _w.begin(); i != _w.end(); ++i) {
 154     untested();
 155     (*i).second += x;
 156   }
 157   return *this;
 158 }
 159 /*--------------------------------------------------------------------------*/
 160 inline WAVE& WAVE::operator*=(const WAVE& x)
 161 {
 162   untested();
 163   for (std::deque<DPAIR>::iterator
 164          i = _w.begin(); i != _w.end(); ++i) {
 165     untested();
 166     (*i).second *= x.v_out((*i).first).f0;
 167   }
 168   return *this;
 169 }
 170 /*--------------------------------------------------------------------------*/
 171 inline WAVE& WAVE::operator*=(double x)
 172 {
 173   for (std::deque<DPAIR>::iterator
 174          i = _w.begin(); i != _w.end(); ++i) {
 175     (*i).second *= x;
 176   }
 177   return *this;
 178 }
 179 /*--------------------------------------------------------------------------*/
 180 inline WAVE& WAVE::warp(double t)
 181 {
 182   assert(t>0);
 183   for (std::deque<DPAIR>::iterator
 184          i = _w.begin(); i != _w.end(); ++i) {
 185     (*i).first *= t;
 186   }
 187   return *this;
 188 }
 189 /*--------------------------------------------------------------------------*/
 190 inline double WAVE::norm2(const DPAIR& a) const
 191 {
 192   double sum = (a.first) * (a.first);
 193   assert(is_number(sum));
 194   sum += (a.second) * (a.second);
 195   assert(is_number(sum));
 196   return sqrt(sum);
 197 }
 198 /*--------------------------------------------------------------------------*/
 199 inline double WAVE::norm2(const DPAIR& a, const DPAIR& b) const
 200 {
 201   double sum = (a.first - b.first) * (a.first - b.first);
 202   assert(is_number(sum));
 203   sum += (a.second - b.second) * (a.second - b.second);
 204   assert(is_number(sum));
 205   return sqrt(sum);
 206 }
 207 /*--------------------------------------------------------------------------*/
 208 inline DPAIR operator+(const DPAIR& a, const DPAIR& b)
 209 {
 210   DPAIR sum(a);
 211   sum.first += b.first;
 212   sum.second += b.second;
 213   return sum;
 214 }
 215 /*--------------------------------------------------------------------------*/
 216 inline DPAIR operator-(const DPAIR& a, const DPAIR& b)
 217 {
 218   DPAIR sub(a);
 219   sub.first -= b.first;
 220   sub.second -= b.second;
 221   return sub;
 222 }
 223 /*--------------------------------------------------------------------------*/
 224 inline DPAIR operator*(const double& t, const DPAIR& a)
 225 {
 226   DPAIR s(a);
 227   s.first *= t;
 228   s.second *= t;
 229   return s;
 230 }
 231 #if __cplusplus > 199711L
 232 // compute segment to point distance
 233 // take left endpoint of segment
 234 inline double WAVE::segment2point(const_iterator li, const DPAIR& p, DPAIR* where) const
 235 {
 236   WAVE::const_iterator ri = next(li);
 237   const DPAIR l = *li;
 238   const DPAIR r = *ri;
 239 //  trace3("norm2", l, r, p);
 240   double ret = 0;
 241
 242   if(li == end()){ unreachable();
 243   }else if(ri == end()){
 244     ret = norm2(l,p);
 245     if(where) {
 246       *where = *li;
 247     }
 248   }else{
 249     assert(l != r);
 250     if(l.first >= r.first){
 251       error(bDANGER, "something wrong %E %E\n", l.first, r.first);
 252     }
 253     assert(l.first < r.first);
 254     // solve l + t*(r-l) closest to p
 255     double t;
 256     double p1 = p.first;
 257     double p2 = p.second;
 258     double l1 = l.first;
 259     double l2 = l.second;
 260     double r1 = r.first;
 261     double r2 = r.second;
 262     t = - (p1 - l1) * (l1 - r1) - (p2 - l2) * (l2 - r2);
 263     assert(is_number(t));
 264     t /= (l1 - r1) * (l1 - r1) + (l2 - r2) * (l2 - r2);
 265     assert(is_number(t));
 266
 267     DPAIR op;
 268     if(t<=0) {
 269       op = l;
 270     } else if(t>=1) {
 271       op = r;
 272     }else{
 273       op = l+(t*(r-l));
 274     }
 275     ret = norm2(op,p);
 276     if(where){
 277       *where = op;
 278     }
 279 //    trace2("norm2 found", t, ret);
 280   }
 281   return ret;
 282 }
 283 /*--------------------------------------------------------------------------*/
 284 // compute directed distance from point to chain
 285 // distance is min_{x \in chain} |b-x|
 286 // chain is iterator of this starting at B
 287 // edge has vertices l and r
 288 // FIXME: pass right end.
 289 inline double WAVE::point2chain(DPAIR p, const_iterator L, DPAIR* w) const
 290 {
 291   assert(L!=end());
 292
 293   double ret = norm2(*L,p);
 294   if(w){
 295     *w = *L;
 296     trace2("point2chain", ret, *w);
 297   }
 298
 299   while(next(L) != end() && L->first < p.first+ret){
 300     trace2("point2chain", *L, *next(L));
 301     assert(next(L)!=end());
 302     DPAIR w2;
 303     DPAIR* pw2 = (w)? &w2 : NULL;
 304     double distance = segment2point(L,p,pw2);
 305     trace4("point2chain", p, distance, *L, w2);
 306     if(distance<ret){
 307       ret = distance;
 308       if (w){
 309         *w = w2;
 310       }
 311     }
 312     ++L;
 313   }
 314
 315   trace4("point2chain", p, *L, *next(L), ret);
 316   if(w) trace1("point2chain", *w);
 317   return(ret);
 318 }
 319 /*--------------------------------------------------------------------------*/
 320 // compute bisector of 2 edges
 321 inline FPOLY1 WAVE::bisector_dual(const const_iterator& a, const const_iterator& b) const
 322 {
 323
 324   FPOLY1 fa(a->first, a->second, (next(a)->second - a->second)/(next(a)->first - a->first));
 325   FPOLY1 fb(b->first, b->second, (next(b)->second - b->second)/(next(b)->first - b->first));
 326
 327 //  trace2("bisector_dual", *a, *next(a));
 328 //  trace2("bisector_dual", *b, *next(b));
 329 //  trace2("bisector_dual", fa, fb);
 330 //  trace2("bisector_dual", CPOLY1(fa), CPOLY1(fb));
 331
 332   DPAIR t = CPOLY1(fa).intersect(CPOLY1(fb));
 333 //   trace1("bisector_dual", t);
 334
 335   double norma = 1./norm2(*next(a) - *a);
 336   DPAIR A = norma*(*next(a) - *a);
 337   double normb = 1./norm2(*next(b) - *b);
 338   DPAIR B = normb*(*next(b) - *b);
 339   DPAIR sum = A+B;
 340
 341   return FPOLY1(t.first, t.second, sum.second/sum.first );
 342 }
 343 /*--------------------------------------------------------------------------*/
 344 // find distance from endpoints to chain
 345 // update chain iterators to relevant region for inner points
 346 inline double WAVE::endpoints2chain(DPAIR l, DPAIR r, const_iterator& L, const_iterator& R, std::pair<DPAIR, DPAIR>* where) const
 347 {
 348   trace4("endpoints2chain start....", l, r, *L, *R);
 349   if (L == R){
 350     // empty range
 351   }
 352   assert(R != end()); // bad right end.
 353   auto Lsave = L; // mark start, so we don't sweep twice
 354   auto Rsave = R;
 355   auto Lmin = L; // found minimal distance here (points to left vertex of edge)
 356   DPAIR hmmmL, hmmmR; // FIXME: merge. later.
 357   DPAIR* hmmmLp = NULL;
 358   DPAIR* hmmmRp = NULL;
 359   if(where){
 360     hmmmL = *L;
 361     hmmmR = *R;
 362     hmmmLp = &hmmmL;
 363     hmmmRp = &hmmmR;
 364   }
 365   double distL = norm2(*L, l);
 366   hmmmL = *L;
 367
 368   // shift L to left if necessary
 369   while(L!=begin()){
 370     double tmp = segment2point(prev(L), l, hmmmLp);
 371     if(tmp < distL){ untested();
 372       distL = tmp;
 373       if(where){ untested();
 374         where->second = hmmmL;
 375       }
 376       Lmin = L;
 377       --L;
 378     } else if( prev(L)->first > l.first - distL) {
 379       --L;
 380     }else{
 381       break;
 382     }
 383   }
 384
 385   trace4("endpoints2chain first L pass", *Lsave, *Lmin, distL, hmmmL);
 386   // shift L to the right as much as possible
 387   while(next(Lsave)!=end()){
 388     trace3("L ->", *Lsave, *next(Lsave), *R);
 389     double tmp = segment2point(Lsave, l, hmmmLp);
 390     if(tmp < distL){
 391       distL = tmp;
 392       if(where){
 393         where->second = hmmmL;
 394       }
 395       Lmin = Lsave;
 396       ++Lsave;
 397     } else if( Lsave->first < l.first + distL) {
 398       ++Lsave;
 399     }else{
 400       break;
 401     }
 402   }
 403
 404   trace4("endpoints2chain second L pass", *Lsave, *Lmin, distL, hmmmL);
 405   L = Lmin;
 406
 407   double distR = norm2(*R, r);
 408   if (distR>distL) {
 409     if(where){
 410       where->second = *R;
 411     }
 412   }
 413   auto Rmin = R; // found minimal distance here (points to right vertex of edge)
 414
 415   trace2("guess", distR, *R);
 416   // shift R to right if neccessary
 417   while(next(R)!=end()){
 418     trace3("R ->", *R, r, *prev(R));
 419     double tmp = segment2point(R, r, hmmmRp); // check edge right to range
 420     if(tmp < distR){
 421       distR = tmp;
 422       if(where && (distR > distL)){
 423         trace1("endpoints2chain R->", hmmmR);
 424         where->second = hmmmR;
 425       }else if(where){
 426         trace3("endpoints2chain R-> nosecond", hmmmR, distL, distR);
 427       }
 428       ++R;
 429       Rmin = R;
 430     } else if( R->first < r.first + distR) {
 431       ++R;
 432     }else{
 433       break;
 434     }
 435   }
 436
 437   trace6("endpoints2chain done first R pass", *Rsave, *Rmin, *next(Rmin), distR, *R, hmmmR);
 438
 439
 440   // shift R to the left as much as possible
 441   while(Rsave != Lmin){
 442     assert(Rsave!=begin());
 443     double tmp = segment2point(prev(Rsave), r, hmmmRp);
 444     if(tmp < distR){
 445       Rmin = Rsave;
 446       --Rsave;
 447       distR = tmp;
 448       if(where && (distR > distL)){
 449         trace1("endpoints2chain <-R", hmmmR);
 450         where->second = hmmmR;
 451       }else if(where){
 452         trace3("endpoints2chain <-R nosecond", hmmmR, distL, distR);
 453       }
 454     } else if( prev(Rsave)->first < r.first + distR) {
 455       --Rsave;
 456     }else{ untested();
 457       break;
 458     }
 459   }
 460   trace6("endpoints2chain done second R pass", *Rsave, *Rmin, *next(Rmin), distR, *R, hmmmR);
 461
 462   R = Rmin;
 463
 464   trace4("endpoints2chain check", *L, *R, distL, distR);
 465   assert(L->first <= R->first);
 466   assert(L != end());
 467
 468 //  auto x = L;
 469 //  for(; x!=R; ++x){
 470 //    assert(x!=end());
 471 //  }
 472 //  assert(x==R);
 473
 474   trace8("endpoints2chain done", l, r, *Lmin, *Rmin, distL, distR, *L, *R);
 475   if (distL > distR){
 476     if(where){
 477       where->first = l;
 478       //where->second = hmmmL;
 479       trace1("endpoints2chain done l", where->second);
 480     }
 481     return distL;
 482   }else{
 483     if(where){
 484       where->first = r;
 485       trace1("endpoints2chain done r", where->second);
 486       //where->second = hmmmR;
 487     }
 488     return distR;
 489   }
 490 }
 491 /*--------------------------------------------------------------------------*/
 492 // compute worst case distance from edge to chain, leaving out right (?) endpoint
 493 // chain is subsegment of this incluging L .. R
 494 // edge has vertices l and r
 495 // B will be incremented as needed
 496 inline double WAVE::edge2chain(DPAIR l, DPAIR r, const_iterator& L, const_iterator& R, std::pair<DPAIR, DPAIR>* where) const
 497 {
 498   double ret = 0;
 499   trace4("edge2chain...", l, r, *L, *R);
 500
 501   // need to check (left) endpoint of edge against range
 502   if(l == r){ incomplete();
 503     return point2chain(l,L);
 504   }
 505
 506   ret = endpoints2chain(l, r, L, R, where);
 507   trace5("edge2chain new endpoints", l, r, *L, *R, ret);
 508   if (where) trace2("", where->first, where->second);
 509
 510   std::set<DPAIR> intersections;
 511
 512   for(auto a = L; a != R; ++a){
 513     trace2("iterating left...", *a, *next(a));
 514     if (next(a)==end()){
 515     }else {
 516       for(auto b = next(a); b != R; ++b){
 517         assert (next(b)!=end());
 518         trace2("iterating right...", *b, *next(b));
 519
 520         if (b==end()){ incomplete();
 521
 522         } else if (next(b)==end()){ incomplete();
 523           // we have only one edge left
 524
 525         } else {
 526           // need to check bisector
 527
 528           FPOLY1 bis = bisector_dual(a,b); // cache?
 529           trace5("bisector...", *a, *next(a), *b, *next(b), bis);
 530           FPOLY1 edge(l,r);
 531
 532           double t = CPOLY1(edge).zero(bis);
 533
 534           if(t>=l.first && t<=r.first){
 535             double ft = edge(t);
 536             trace6("found intersection...", t, ft, *a, *b, l, r);
 537             DPAIR ww;
 538             DPAIR* w2 = (where)? &ww : NULL;
 539             double d = point2chain(DPAIR(t,ft), a, w2); // FIXME. pass right end
 540             if(w2) trace1("point2chain done", *w2);
 541             if(d<ret){
 542               trace1("intersection is closer", d);
 543             }else{
 544               if(where){
 545                 trace2("further intersection...", *w2, d);
 546                 where->first = DPAIR(t,ft);
 547                 where->second = ww;
 548               }
 549               ret = d;
 550             }
 551           }else{
 552             trace3("no intersection", t, bis, edge);
 553           }
 554         }
 555       }
 556     }
 557   }
 558   return ret;
 559 }
 560 /*--------------------------------------------------------------------------*/
 561 // compute directed hausdorff distance from this to that using linear
 562 // interpolations
 563 // == max_{i \in lin(this)} min_{j \in lin(that)} d(i,j)
 564 inline double WAVE::dhd_linear(const WAVE& that, std::pair<DPAIR,DPAIR>* wp) const
 565 {
 566   auto L = that.begin();
 567   auto R = L;
 568   assert(L!=that.end()); // for now.
 569
 570   double ret = 0.; // fabs(that.v_out(begin()->first).f0 - begin()->second);
 571
 572   for(auto ii = begin(); ; ++ii){
 573     auto other = next(ii);
 574
 575 //    ret = std::max(ret, fabs(other->second - that.v_out(other->first).f0)); // hmmm
 576     trace3("dhd outer loop", ret, *ii, *other);
 577
 578     if (other==end()){
 579       // necessary?!
 580       trace3("dhd, end", *ii, ret, *L);
 581       DPAIR ww;
 582       auto wp2 = (wp)? &ww : NULL;
 583       double r = that.point2chain(*ii, L, wp2);
 584       if(ret<r){
 585         ret = r;
 586         if(wp){
 587           wp->first = *ii;
 588           wp->second = ww;
 589         }
 590       }else{
 591       }
 592       trace3("dhd end", ii->first, ret, *L);
 593       if(wp) trace3("dhd end", ii->second, wp->first, wp->second);
 594       return ret;
 595     } else if (R==that.end()) { untested();
 596       DPAIR ww;
 597       auto wp2 = (wp)? &ww : NULL;
 598       double r = that.segment2point(ii, *L, wp2);
 599       if(ret<r){untested();
 600         ret = r;
 601         if(wp){untested();
 602           wp->first = *ii;
 603           wp->second = ww;
 604         }
 605       }
 606       trace3("dhd singleton chain", ret, *L, *R);
 607     } else {
 608       // perform windowed sweep.
 609       std::pair<DPAIR,DPAIR> ww;
 610       auto wp2 = (wp)? &ww : NULL;
 611       double r = that.edge2chain(*ii, *other, L, R, wp2);
 612
 613       if(r>ret){
 614         ret = r;
 615         if(wp){
 616           *wp = ww;
 617         }
 618         trace5("dhd windowed sweep, found", *ii, *other, ret, *L, *R);
 619         if(wp) trace3("dhd windowed sweep", ret, wp->first, wp->second);
 620       }else{
 621       }
 622     }
 623
 624   }
 625
 626   return ret;
 627 }
 628 /*--------------------------------------------------------------------------*/
 629 // compute directed hausdorff distance from points of this to points of that
 630 // this is useless, so lets not try to be efficient.
 631 inline double WAVE::dhd_discrete(const WAVE& that, std::pair<DPAIR, DPAIR>* where) const
 632 {
 633   double ret = 0;
 634   DPAIR here;
 635   for(auto ii : _w) {
 636     double dist = inf;
 637     for(auto j : that._w) {
 638       double norm = norm2(ii,j);
 639       if(norm < dist){
 640         dist = norm;
 641         if(where) {
 642           here = j;
 643         }
 644       }
 645     }
 646     if(ret<dist){
 647       ret = dist;
 648       if(where){
 649         where->first = ii;
 650         where->second = here;
 651       }
 652     }
 653   }
 654   return ret;
 655 }
 656 #endif
 657 /*--------------------------------------------------------------------------*/
 658 // vim:ts=8:sw=2:noet: