PageRenderTime 514ms CodeModel.GetById 60ms app.highlight 385ms RepoModel.GetById 39ms app.codeStats 1ms

/@QMITL_Formula/private/src/robusthom/robustness.cpp

https://bitbucket.org/donze/breach
C++ | 797 lines | 586 code | 140 blank | 71 comment | 145 complexity | beebe878d863a4436b239163904a61fe MD5 | raw file
  1#include "robustness.h"
  2#include <list>
  3
  4using namespace std;
  5
  6/*---------------------------------------------------------------------------*
  7 *                           MAIN FUNCTIONS                                  *
  8 *---------------------------------------------------------------------------*/
  9
 10Signal * computeNot(Signal * y) {
 11  Signal * z = new Signal();
 12  Signal::iterator i;
 13	
 14  z->beginTime=y->beginTime;
 15  z->endTime=y->endTime;
 16
 17  for(i = y->begin(); i != y->end(); i++) {
 18    z->push_back(Sample(i->time, -(i->value), -(i->derivative)));	        
 19  }
 20  return z;
 21}
 22
 23
 24Signal * computeAnd(Signal * x, Signal * y) {
 25  Signal::reverse_iterator i=x->rbegin();
 26  Signal::reverse_iterator j=y->rbegin();
 27	
 28  Signal * z=new Signal();	
 29
 30  z->beginTime=fmax(x->beginTime, y->beginTime);
 31  z->endTime=fmin(x->endTime, y->endTime);
 32
 33  while(i->time >= z->endTime) i++;
 34  while(j->time >= z->endTime) j++;
 35
 36  computePartialAnd(z, i, j, z->beginTime, z->endTime);
 37  z->simplify();
 38	
 39  return z;
 40
 41}
 42
 43
 44Signal * computeOr(Signal * x, Signal * y) {
 45  Signal::reverse_iterator i=x->rbegin();
 46  Signal::reverse_iterator j=y->rbegin();
 47	
 48  Signal * z=new Signal();	
 49
 50  z->beginTime=fmax(x->beginTime, y->beginTime);
 51  z->endTime=fmin(x->endTime, y->endTime);
 52
 53  while(i->time >= z->endTime) i++;
 54  while(j->time >= z->endTime) j++;
 55
 56  computePartialOr(z, i, j, z->beginTime, z->endTime);
 57  z->simplify();
 58	
 59  return z;
 60
 61}
 62
 63
 64Signal * computeEventually(Signal * x) {
 65#ifdef DEBUG__
 66  cout << "Entering computeEventually" << endl;
 67#endif
 68  
 69  Signal * z=new Signal();
 70
 71  z->beginTime=x->beginTime;
 72  z->endTime=x->endTime;
 73
 74  Signal::reverse_iterator i=x->rbegin();	
 75
 76  computePartialEventually(z,i,z->beginTime, z->endTime);
 77
 78  return z;
 79}
 80
 81
 82Signal * computeUntil(Signal * x, Signal * y) {
 83#ifdef DEBUG__
 84  cout << "Entering computeUntil" << endl;
 85#endif
 86
 87  double s, t;
 88  double z_max=BOTTOM;	
 89  Signal * z=new Signal();
 90	
 91  Signal::reverse_iterator i=x->rbegin();
 92  Signal::reverse_iterator j=y->rbegin();
 93
 94  z->beginTime=fmax(x->beginTime, y->beginTime);
 95  z->endTime=fmin(x->endTime, y->endTime);
 96
 97  while(i->time >= z->endTime) i++;
 98  while(j->time >= z->endTime) j++;
 99
100  s=z->beginTime;
101  t=z->endTime;		
102
103  while(i->time > s) {
104
105    computeSegmentUntil(z, *i, t, j, z_max);
106    z_max=z->front().value;
107
108    if(j->time == i->time) j++; 
109    t=i->time;
110    i++;
111  }
112
113  if(i->time == s) computeSegmentUntil(z, *i, t, j, z_max);
114  else computeSegmentUntil(z, Sample(s, i->valueAt(s), i->derivative), t, j, z_max);
115
116  z->simplify();
117	
118  return z;
119}
120
121
122Signal * computeBoundedEventually(Signal * x, double a) {
123#ifdef DEBUG__
124  cout << "Entering computeBoundedEventually" << endl;
125#endif
126
127  Signal *z, *z1, *z2, *z3;
128
129  z1=plateauMax(x, a);
130  z2=new Signal(*x);
131  z2->resize(x->beginTime + a, x->endTime + a, BOTTOM);
132  z2->shift(-a);
133  z3=computeOr(z2, z1);
134    
135  delete z1;
136  delete z2;
137
138  z=new Signal();
139  z=computeOr(x, z3);
140    
141  delete z3;
142  z->simplify();
143	
144  return z;
145}
146
147
148Signal * computeBoundedGlobally(Signal * x, double a) {
149  Signal *z, *z1, *z2, *z3;
150
151  z1=plateauMin(x, a);
152
153  z2=new Signal(*x);
154  z2->resize(x->beginTime + a, x->endTime + a, TOP);
155  z2->shift(-a);
156
157  z3=computeAnd(z2, z1);
158  delete z1;
159  delete z2;
160
161  z=new Signal();
162  z=computeAnd(x, z3);
163  delete z3;
164		
165  z->simplify();
166	
167  return z;
168}
169
170
171Signal * computeTimedUntil(Signal * x, Signal * y, double a, double b) {
172  Signal *z, *z1, *z2, *z3, *z4;
173
174  z=new Signal();
175  if (a>0)
176    z1=computeBoundedGlobally(x,a);
177	
178  z2=computeBoundedEventually(y,b-a);
179  z3=computeUntil(x,y);
180  z4=computeAnd(z2,z3);
181	
182  if (a>0) {
183    z4->shift(-a);
184    z=computeAnd(z1,z4);
185  }
186  else
187    z=computeAnd(x,z4);
188
189  return z;
190}
191
192/*---------------------------------------------------------------------------*
193 *                         CONJUNCTION SUBROUTINES                           *
194 *---------------------------------------------------------------------------*/
195
196//PRECONDITIONS: j->time < t, i.time < t.
197//POSTCONDITIONS: j->time <= i.time.
198void computeSegmentAnd(Signal * z, const Sample & i, double t, Signal::reverse_iterator & j) {
199  bool continued=false;
200  double s=j->time;
201
202  // for every sample *j in (i.time, t) 
203  while(s > i.time) {
204    if(i.valueAt(t) < j->valueAt(t)) {
205      if (i.valueAt(s) > j->value) {
206	t=i.timeIntersect(*j);
207	z->push_front(Sample(t,i.valueAt(t),i.derivative));
208	z->push_front(Sample(s,j->value,j->derivative));
209	continued=false;
210      }
211      else continued=true;
212    }
213    else if (i.valueAt(t) == j->valueAt(t)) {
214      if (i.valueAt(s) > j->value) {
215	if(continued) {
216	  z->push_front(Sample(t,i.valueAt(t), i.derivative));
217	  continued=false;
218	}
219	z->push_front(Sample(s,j->value,j->derivative));
220      }
221      else continued=true;
222    }
223    else {
224      if (i.valueAt(s) < j->value) {
225	if(continued) {
226	  z->push_front(Sample(t,i.valueAt(t),i.derivative));
227	}
228	t=i.timeIntersect(*j);
229	z->push_front(Sample(t,j->valueAt(t),j->derivative));
230	continued=true;
231      }
232      else {
233	if(continued) {
234	  z->push_front(Sample(t,i.valueAt(t),i.derivative));
235	  continued=false;
236	}
237	z->push_front(Sample(s,j->value,j->derivative));	
238      }
239    }
240		
241    //increment reverse iterator j
242    t = s;
243    j++;
244    s = j->time;
245  }
246
247
248  //here we may have j->time < i.time
249  // "i" values of z are no longer "continued"
250  s=i.time;
251  if(i.valueAt(t) < j->valueAt(t)) {
252    if (i.value > j->valueAt(s)) {
253      t=i.timeIntersect(*j);
254      z->push_front(Sample(t,i.valueAt(t),i.derivative));
255      z->push_front(Sample(s,j->valueAt(s),j->derivative));
256    }
257    else {
258      z->push_front(i);
259    }
260  }
261  else if (i.valueAt(t) == j->valueAt(t)) {
262    if (i.value > j->valueAt(s)) {
263      if(continued) {
264	z->push_front(Sample(t,i.valueAt(t), i.derivative));
265      }
266      z->push_front(Sample(s,j->valueAt(s),j->derivative));
267    }
268    else {
269      z->push_front(i);
270    }
271  }
272  else {
273    if (i.value < j->valueAt(s)) {
274      if(continued) {
275          z->push_front(Sample(t,i.valueAt(t),i.derivative));
276      }
277      t=i.timeIntersect(*j);
278      z->push_front(Sample(t,j->valueAt(t),j->derivative));
279      z->push_front(i);
280    }
281    else {
282      if(continued) {
283          z->push_front(Sample(t,i.valueAt(t),i.derivative));
284      }
285      z->push_front(Sample(s,j->valueAt(s),j->derivative));	
286    }
287  }
288
289}
290
291
292//void computeConstantSegmentAnd(Signal *, double, Signal::reverse_iterator &, double, double) { }
293
294void computePartialAnd(Signal * z, Signal::reverse_iterator & i, Signal::reverse_iterator & j, double s, double t) {
295
296  while(i->time > s) {
297    computeSegmentAnd(z,*i,t,j);
298    if (j->time == i->time) j++; 
299    t=i->time;
300    i++;
301  }
302
303  if(i->time == s) computeSegmentAnd(z, *i, t, j);
304  else computeSegmentAnd(z, Sample(s,i->valueAt(s),i->derivative), t, j);	
305
306}
307
308
309/*---------------------------------------------------------------------------*
310 *                         DISJUNCTION SUBROUTINES                           *
311 *---------------------------------------------------------------------------*/
312
313// copy of computeSegmentAnd, operator "<" switched with operator ">"
314void computeSegmentOr(Signal * z, const Sample & i, double t, Signal::reverse_iterator & j) {
315  bool continued=false;
316  double s=j->time;
317
318  // for every sample *j in (i.time, t) 
319  while(s > i.time) {
320    if(i.valueAt(t) > j->valueAt(t)) {
321      if (i.valueAt(s) < j->value) {
322          t=i.timeIntersect(*j);
323          z->push_front(Sample(t,i.valueAt(t),i.derivative));
324          z->push_front(Sample(s,j->value,j->derivative));
325          continued=false;
326      }
327      else continued=true;
328    }
329    else if (i.valueAt(t) == j->valueAt(t)) {
330      if (i.valueAt(s) < j->value) {
331	if(continued) {
332	  z->push_front(Sample(t,i.valueAt(t), i.derivative));
333	  continued=false;
334	}
335	z->push_front(Sample(s,j->value,j->derivative));
336      }
337      else continued=true;
338    }
339    else {
340      if (i.valueAt(s) > j->value) {
341	if(continued) {
342	  z->push_front(Sample(t,i.valueAt(t),i.derivative));
343	}
344	t=i.timeIntersect(*j);
345	z->push_front(Sample(t,j->valueAt(t),j->derivative));
346	continued=true;
347      }
348      else {
349	if(continued) {
350	  z->push_front(Sample(t,i.valueAt(t),i.derivative));
351	  continued=false;
352	}
353	z->push_front(Sample(s,j->value,j->derivative));	
354      }
355    }
356		
357    //increment reverse iterator j
358    t = s;
359    j++;
360    s = j->time;
361  }
362
363
364  //here we may have j->time < i.time
365  // "i" values of z are no longer "continued"
366  s=i.time;
367  if(i.valueAt(t) > j->valueAt(t)) {
368    if (i.value < j->valueAt(s)) {
369      t=i.timeIntersect(*j);
370      z->push_front(Sample(t,i.valueAt(t),i.derivative));
371      z->push_front(Sample(s,j->valueAt(s),j->derivative));
372    }
373    else {
374      z->push_front(i);
375    }
376  }
377  else if (i.valueAt(t) == j->valueAt(t)) {
378    if (i.value < j->valueAt(s)) {
379      if(continued) {
380	z->push_front(Sample(t,i.valueAt(t), i.derivative));
381      }
382      z->push_front(Sample(s,j->valueAt(s),j->derivative));
383    }
384    else {
385      z->push_front(i);
386    }
387  }
388  else {
389    if (i.value > j->valueAt(s)) {
390      if(continued) {
391	z->push_front(Sample(t,i.valueAt(t),i.derivative));
392      }
393      t=i.timeIntersect(*j);
394      z->push_front(Sample(t,j->valueAt(t),j->derivative));
395      z->push_front(i);
396    }
397    else {
398      if(continued) {
399	z->push_front(Sample(t,i.valueAt(t),i.derivative));
400      }
401      z->push_front(Sample(s,j->valueAt(s),j->derivative));	
402    }
403  }
404
405}
406
407
408void computePartialOr(Signal * z, Signal::reverse_iterator & i, Signal::reverse_iterator & j, double s, double t) {
409
410  while(i->time > s) {
411    computeSegmentOr(z,*i,t,j);
412    if (j->time == i->time) j++; 
413    t=i->time;
414    i++;
415  }
416
417  if(i->time == s) computeSegmentOr(z, *i, t, j);
418  else computeSegmentOr(z, Sample(s,i->valueAt(s),i->derivative), t, j);	
419
420}
421
422
423
424/*---------------------------------------------------------------------------*
425 *                              UNTIL SUBROUTINES                            *
426 *---------------------------------------------------------------------------*/
427
428void computePartialEventually(Signal* z, Signal::reverse_iterator & i, double s, double t) {
429  bool continued=false;
430  double z_max=BOTTOM;
431  while(i->time > s) {
432    if(i->derivative >= 0) {	
433      if(z_max < i->valueAt(t)) {
434	if(continued) {
435	  z->push_front(Sample(t,z_max,0));
436	}
437	z_max=i->valueAt(t);
438      }
439      continued=true;
440      //z->push_front(Sample(i->time, z_max, 0));
441    }
442    else if(i->valueAt(t) >= z_max) {
443      if(continued) {
444	z->push_front(Sample(t,z_max,0));
445	continued=false;
446      }
447      z_max=i->value;
448      z->push_front(*i);
449    }
450    else if(z_max >= i->value) {
451      continued=true;
452      //z->push_front(Sample(i->time, z_max, 0));
453    }
454    else {
455      z->push_front(Sample(i->time + (z_max-i->value)/i->derivative, z_max, 0)); //time at which y reaches value next_z
456      z->push_front(*i);
457      z_max=i->value;			
458      continued=false;
459    }
460
461    t=i->time;
462    i++;		
463  }
464
465  //leftmost sample *i may not be on s
466  //"z_max" values of z are not longer "continued".
467  if(i->derivative >= 0) {	
468    if(z_max < i->valueAt(t)) {
469      if(continued) {
470	z->push_front(Sample(t,z_max,0));
471      }
472      z_max=i->valueAt(t);
473    }
474    z->push_front(Sample(s, z_max, 0));
475  }
476  else if(i->valueAt(t) >= z_max) {
477    if(continued) {
478      z->push_front(Sample(t,z_max,0));
479    }
480    z->push_front(Sample(s, i->valueAt(s), i->derivative));
481  }
482  else if(z_max >= i->valueAt(s)) {
483    z->push_front(Sample(s, z_max, 0));
484  }
485  else {
486    z->push_front(Sample(s + (z_max-i->value)/i->derivative, z_max, 0)); //time at which y reaches value next_z
487    z->push_front(Sample(s, i->valueAt(s), i->derivative));
488  }
489	
490}
491
492void computeSegmentUntil(Signal * z, const Sample & i, double t, Signal::reverse_iterator & j, double z_max) {
493  Signal *z1, *z2, *z3;
494  Signal::reverse_iterator k, l;
495  double s=i.time;
496
497  if(i.derivative <= 0) {
498    z1=new Signal();
499    computeSegmentAnd(z1, i, t, j);
500
501    z2=new Signal();
502    k=z1->rbegin();
503    computePartialEventually(z2, k, s, t);
504    delete z1;
505
506    l=z2->rbegin();
507    computeSegmentOr(z, Sample(s, fmin(z_max, i.valueAt(t)), 0), t, l);
508    delete z2;
509  }
510  else {
511    z1=new Signal();
512    computePartialEventually(z1, j, s, t);
513
514    z2=new Signal();
515    k=z1->rbegin();
516    computeSegmentAnd(z2, i, t, k);
517    delete z1;
518
519    z1=new Signal();
520    z3=new Signal();
521    z3->push_front(Sample(s, z_max, 0));
522    k=z3->rbegin();
523    computeSegmentAnd(z1, i, t, k);
524    delete z3;
525
526    k=z1->rbegin();
527    l=z2->rbegin();
528    computePartialOr(z, k, l, s, t);
529    delete z1;
530    delete z2;
531  }
532}
533
534/*---------------------------------------------------------------------------*
535 *                           TIMED UNTIL SUBROUTINES                         *
536 *---------------------------------------------------------------------------*/
537
538//computation of the max induced by discontinuity points of x in t+(0,a] 
539//a is assumed to be smaller than length of x.
540//based on Lemire algorithm for "streaming min-max filter"
541Signal * plateauMax(Signal * x, double a) {
542#ifdef DEBUG__
543  cout << "Entering plateauMax" << endl;
544#endif
545  bool new_candidate, end_candidate;
546  double t, v;
547  Signal::reverse_iterator j;
548  Sequence M; //sorted in ascending times from front to back
549  Sequence y; //maximum of x(t) and x(t-) at discontinuity points of x
550  Sequence::iterator i;
551  
552  Signal * z=new Signal();
553  z->beginTime=x->beginTime;
554  z->endTime=x->endTime;
555
556
557  //PRECOMPUTATION: list the local maximums of x	
558  t=x->endTime;	
559  v=BOTTOM;
560  for(j=x->rbegin(); j!= x->rend(); j++) {
561    y.push_front(Point(t, fmax(v,j->valueAt(t))));
562    //    cout << "y.time:" << y.front().time << " y->value:" << y.front().value << endl;
563    t=j->time;
564    v=j->value;
565  }
566  y.push_front(Point(x->front().time, x->front().value));
567  
568  //INIT: read values in [0, a)
569  i=y.begin();
570  
571  while(i->time < x->beginTime + a) {
572    while(!M.empty() && i->value >= M.back().value) {
573      M.pop_back();
574    }
575    //    cout << "i->time:" << i->time << " i->value:" << i->value << endl;
576
577    M.push_back(*i);
578    i++;
579  }
580
581
582  if(i->time == x->beginTime + a) new_candidate=true;
583  else new_candidate=false;
584  
585  end_candidate=false;
586  t=x->beginTime;
587  
588  //  cout << "M.front "  << M.front().value << endl;; 
589
590  //STEP
591  bool cont=true;
592  while(cont) {
593    
594    //UPDATE OF CANDIDATE LIST
595    //candidate crosses t: remove it from list
596    if(end_candidate) {
597      //      cout << "end_candidate: M.pop front" << endl;
598      M.pop_front();
599    }
600
601    //sample crosses t+a: add it to candidate list
602    if(new_candidate) {
603      //cout << "new_candidate: doing stuff" << endl;
604      while(!M.empty() && i->value >= M.back().value) {
605	M.pop_back();
606      }
607      //detect if new maximum is found
608      if(!M.empty()) {
609	//if M non empty then t + a does not generate new maximum
610	new_candidate=false;
611      }
612      //add candidate
613      M.push_back(*i);
614
615      //increment iterator
616      i++;
617    }
618
619    //OUTPUT OF NEW MAXIMUM
620      //      cout << "Output new maximum" << endl;
621      //no candidate
622      if(M.empty()) {
623	//cout << "No candidate" << endl;
624	z->push_back(Sample(t, BOTTOM, 0));
625      }
626      //next best candidate
627      else {
628	z->push_back(Sample(t, M.front().value, 0));
629	//	cout << "Pushed " << M.front().value << endl;
630      }
631
632    //NEXT EVENT DETECTION
633    if(! M.empty()) {
634      if(i != y.end()) {
635	if(i->time - a == M.front().time) {
636	  t=M.front().time;
637	  new_candidate=true;
638	  end_candidate=true;
639	}
640	else if(i->time - a < M.front().time) {
641	  t=i->time - a;
642	  new_candidate=true;
643	  end_candidate=false;
644	}
645	else { //M.back().time < i->time - a
646	  t=M.front().time;
647	  new_candidate=false;
648	  end_candidate=true;
649					
650	}
651      }
652      else {
653	t=M.front().time;
654	new_candidate=false;
655	end_candidate=true;
656      }
657    }
658    else {
659      if(i != y.end()) {
660	t=i->time - a;
661	new_candidate=true;
662	end_candidate=false;
663      }
664      else {
665	new_candidate=false;
666	end_candidate=false;
667      }
668    }
669    cont = (new_candidate||end_candidate);
670 }
671
672  if(z->back().time==z->endTime) z->pop_back();
673  
674  return z;
675}
676
677
678//copy of plateauMax, operator < switched with operator >, TOP replaces BOTTOM
679Signal * plateauMin(Signal * x, double a) {
680  bool new_candidate, end_candidate;
681  double t, v;
682  Signal::reverse_iterator j;
683  Sequence M; //sorted in ascending times from front to back
684  Sequence y; //minimum of x(t) and x(t-) at discontinuity points of x
685  Sequence::iterator i;
686
687  Signal * z=new Signal();
688  z->beginTime=x->beginTime;
689  z->endTime=x->endTime;
690
691
692  //PRECOMPUTATION: list the local minimums of x	
693  t=x->endTime;	
694  v=TOP;
695  for(j=x->rbegin(); j!= x->rend(); j++) {
696    y.push_front(Point(t, fmin(v,j->valueAt(t))));
697    t=j->time;
698    v=j->value;
699  }
700  y.push_front(Point(x->front().time, x->front().value));
701
702  //INIT: read values in [0, a)
703  i=y.begin();
704  while(i->time < x->beginTime + a) {
705    while(!M.empty() && i->value <= M.back().value) {
706      M.pop_back();
707    }
708    M.push_back(*i);
709    i++;
710  }
711  if(i->time == x->beginTime + a) new_candidate=true;
712  else new_candidate = false;
713  end_candidate=false;
714  t=x->beginTime;
715
716  //STEP
717  bool cont=true;
718  while(cont) {
719    //UPDATE OF CANDIDATE LIST
720    //candidate crosses t: remove it from list
721    if(end_candidate) {
722      M.pop_front();
723    }
724
725    //sample crosses t+a: add it to candidate list
726    if(new_candidate) {
727      while(!M.empty() && i->value <= M.back().value) {
728	M.pop_back();
729      }
730      //detect if new minimum is found
731      if(!M.empty()) {
732	//if M non empty then t + a does not generate new minimum
733	new_candidate=false;
734      }
735      //add candidate
736      M.push_back(*i);
737
738      //increment iterator
739      i++;
740    }
741
742    //OUTPUT OF NEW MINIMUM
743      //no candidate
744      if(M.empty()) {
745          z->push_back(Sample(t, TOP, 0));
746      }
747      //next best candidate
748      else {
749          z->push_back(Sample(t, M.front().value, 0));
750      }
751
752    //NEXT EVENT DETECTION
753    if(! M.empty()) {
754      if(i != y.end()) {
755	if(i->time - a == M.front().time) {
756	  t=M.front().time;
757	  new_candidate=true;
758	  end_candidate=true;
759	}
760	else if(i->time - a < M.front().time) {
761	  t=i->time - a;
762	  new_candidate=true;
763	  end_candidate=false;
764	}
765	else { //M.back().time < i->time - a
766	  t=M.front().time;
767	  new_candidate=false;
768	  end_candidate=true;
769					
770	}
771      }
772      else {
773	t=M.front().time;
774	new_candidate=false;
775	end_candidate=true;
776      }
777    }
778    else {
779      if(i != y.end()) {
780	t=i->time - a;
781	new_candidate=true;
782	end_candidate=false;
783      }
784      else {
785	new_candidate=false;
786	end_candidate=false;
787      }
788    }
789    cont = (new_candidate||end_candidate);
790  }
791
792  if(z->back().time==z->endTime) z->pop_back();
793
794  return z;
795}
796
797