ViSP  2.8.0
vpHinkley.cpp
1 /****************************************************************************
2  *
3  * $Id: vpHinkley.cpp 4056 2013-01-05 13:04:42Z fspindle $
4  *
5  * This file is part of the ViSP software.
6  * Copyright (C) 2005 - 2013 by INRIA. All rights reserved.
7  *
8  * This software is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * ("GPL") version 2 as published by the Free Software Foundation.
11  * See the file LICENSE.txt at the root directory of this source
12  * distribution for additional information about the GNU GPL.
13  *
14  * For using ViSP with software that can not be combined with the GNU
15  * GPL, please contact INRIA about acquiring a ViSP Professional
16  * Edition License.
17  *
18  * See http://www.irisa.fr/lagadic/visp/visp.html for more information.
19  *
20  * This software was developed at:
21  * INRIA Rennes - Bretagne Atlantique
22  * Campus Universitaire de Beaulieu
23  * 35042 Rennes Cedex
24  * France
25  * http://www.irisa.fr/lagadic
26  *
27  * If you have questions regarding the use of this file, please contact
28  * INRIA at visp@inria.fr
29  *
30  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
31  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
32  *
33  *
34  * Description:
35  * Hinkley's cumulative sum test implementation.
36  *
37  * Authors:
38  * Fabien Spindler
39  *
40  *****************************************************************************/
41 
51 #include <visp/vpHinkley.h>
52 #include <visp/vpDebug.h>
53 #include <visp/vpIoTools.h>
54 #include <visp/vpMath.h>
55 
56 #include <stdio.h>
57 #include <stdlib.h>
58 #include <iostream>
59 #include <cmath> // std::fabs
60 #include <limits> // numeric_limits
61 
62 /* VP_DEBUG_MODE fixed by configure:
63  1:
64  2:
65  3: Print data
66 */
67 
68 
81 {
82  init();
83 
84  setAlpha(0.2);
85  setDelta(0.2);
86 }
87 
104 vpHinkley::vpHinkley(double alpha, double delta)
105 {
106  init();
107 
108  setAlpha(alpha);
109  setDelta(delta);
110 }
111 
125 void
126 vpHinkley::init(double alpha, double delta)
127 {
128  init();
129 
130  setAlpha(alpha);
131  setDelta(delta);
132 }
133 
140 {
141 }
142 
150 {
151  nsignal = 0;
152  mean = 0.0;
153 
154  Sk = 0;
155  Mk = 0;
156 
157  Tk = 0;
158  Nk = 0;
159 }
160 
169 void vpHinkley::setDelta(double delta)
170 {
171  dmin2 = delta / 2;
172 }
173 
181 void vpHinkley::setAlpha(double alpha)
182 {
183  this->alpha = alpha;
184 }
185 
197 {
198 
199  vpHinkleyJumpType jump = noJump;
200 
201  nsignal ++; // Signal length
202 
203  if (nsignal == 1) mean = signal;
204 
205  // Calcul des variables cumulées
206  computeSk(signal);
207 
208  computeMk();
209 
210  vpCDEBUG(2) << "alpha: " << alpha << " dmin2: " << dmin2
211  << " signal: " << signal << " Sk: " << Sk << " Mk: " << Mk;
212 
213  // teste si les variables cumulées excèdent le seuil
214  if ((Mk - Sk) > alpha)
215  jump = downwardJump;
216 
217 #ifdef VP_DEBUG
218  if (VP_DEBUG_MODE >=2) {
219  switch(jump) {
220  case noJump:
221  std::cout << "noJump " << std::endl;
222  break;
223  case downwardJump:
224  std::cout << "downWardJump " << std::endl;
225  break;
226  case upwardJump:
227  std::cout << "upwardJump " << std::endl;
228  break;
229  }
230  }
231 #endif
232 
233  computeMean(signal);
234 
235  if (jump == downwardJump) {
236  vpCDEBUG(2) << "\n*** Reset the Hinkley test ***\n";
237 
238  Sk = 0; Mk = 0; nsignal = 0;
239  }
240 
241  return (jump);
242 }
243 
255 {
256 
257  vpHinkleyJumpType jump = noJump;
258 
259  nsignal ++; // Signal length
260 
261  if (nsignal == 1) mean = signal;
262 
263  // Calcul des variables cumulées
264  computeTk(signal);
265 
266  computeNk();
267 
268  vpCDEBUG(2) << "alpha: " << alpha << " dmin2: " << dmin2
269  << " signal: " << signal << " Tk: " << Tk << " Nk: " << Nk;
270 
271  // teste si les variables cumulées excèdent le seuil
272  if ((Tk - Nk) > alpha)
273  jump = upwardJump;
274 
275 #ifdef VP_DEBUG
276  if (VP_DEBUG_MODE >= 2) {
277  switch(jump) {
278  case noJump:
279  std::cout << "noJump " << std::endl;
280  break;
281  case downwardJump:
282  std::cout << "downWardJump " << std::endl;
283  break;
284  case upwardJump:
285  std::cout << "upWardJump " << std::endl;
286  break;
287  }
288  }
289 #endif
290  computeMean(signal);
291 
292  if (jump == upwardJump) {
293  vpCDEBUG(2) << "\n*** Reset the Hinkley test ***\n";
294 
295  Tk = 0; Nk = 0; nsignal = 0;
296  }
297 
298  return (jump);
299 }
300 
312 {
313 
314  vpHinkleyJumpType jump = noJump;
315 
316  nsignal ++; // Signal length
317 
318  if (nsignal == 1) mean = signal;
319 
320  // Calcul des variables cumulées
321  computeSk(signal);
322  computeTk(signal);
323 
324  computeMk();
325  computeNk();
326 
327  vpCDEBUG(2) << "alpha: " << alpha << " dmin2: " << dmin2
328  << " signal: " << signal
329  << " Sk: " << Sk << " Mk: " << Mk
330  << " Tk: " << Tk << " Nk: " << Nk << std::endl;
331 
332  // teste si les variables cumulées excèdent le seuil
333  if ((Mk - Sk) > alpha)
334  jump = downwardJump;
335  else if ((Tk - Nk) > alpha)
336  jump = upwardJump;
337 
338 #ifdef VP_DEBUG
339  if (VP_DEBUG_MODE >= 2) {
340  switch(jump) {
341  case noJump:
342  std::cout << "noJump " << std::endl;
343  break;
344  case downwardJump:
345  std::cout << "downWardJump " << std::endl;
346  break;
347  case upwardJump:
348  std::cout << "upwardJump " << std::endl;
349  break;
350  }
351  }
352 #endif
353 
354  computeMean(signal);
355 
356  if ((jump == upwardJump) || (jump == downwardJump)) {
357  vpCDEBUG(2) << "\n*** Reset the Hinkley test ***\n";
358 
359  Sk = 0; Mk = 0; Tk = 0; Nk = 0; nsignal = 0;
360  // Debut modif FS le 03/09/2003
361  mean = signal;
362  // Fin modif FS le 03/09/2003
363  }
364 
365  return (jump);
366 }
367 
376 void vpHinkley::computeMean(double signal)
377 {
378  // Debut modif FS le 03/09/2003
379  // Lors d'une chute ou d'une remontée lente du signal, pariculièrement
380  // après un saut, la moyenne a tendance à "dériver". Pour réduire ces
381  // dérives de la moyenne, elle n'est remise à jour avec la valeur
382  // courante du signal que si un début de saut potentiel n'est pas détecté.
383  //if ( ((Mk-Sk) == 0) && ((Tk-Nk) == 0) )
384  if ( ( std::fabs(Mk-Sk) <= std::fabs(vpMath::maximum(Mk,Sk))*std::numeric_limits<double>::epsilon() )
385  &&
386  ( std::fabs(Tk-Nk) <= std::fabs(vpMath::maximum(Tk,Nk))*std::numeric_limits<double>::epsilon() ) )
387  // Fin modif FS le 03/09/2003
388 
389  // Mise a jour de la moyenne.
390  mean = (mean * (nsignal - 1) + signal) / (nsignal);
391 
392 }
400 void vpHinkley::computeSk(double signal)
401 {
402 
403  // Calcul des variables cumulées
404  Sk += signal - mean + dmin2;
405 }
411 void vpHinkley::computeMk()
412 {
413  if (Sk > Mk) Mk = Sk;
414 }
421 void vpHinkley::computeTk(double signal)
422 {
423 
424  // Calcul des variables cumulées
425  Tk += signal - mean - dmin2;
426 }
432 void vpHinkley::computeNk()
433 {
434  if (Tk < Nk) Nk = Tk;
435 }
436 
438 {
439  switch(jump)
440  {
441  case noJump :
442  std::cout << " No jump detected " << std::endl ;
443  break ;
444  case downwardJump :
445  std::cout << " Jump downward detected " << std::endl ;
446  break ;
447  case upwardJump :
448  std::cout << " Jump upward detected " << std::endl ;
449  break ;
450  default:
451  std::cout << " Jump detected " << std::endl ;
452  break ;
453 
454  }
455 }
456 
vpHinkleyJumpType testDownUpwardJump(double signal)
Definition: vpHinkley.cpp:311
vpHinkleyJumpType
Definition: vpHinkley.h:106
vpHinkleyJumpType testUpwardJump(double signal)
Definition: vpHinkley.cpp:254
static Type maximum(const Type &a, const Type &b)
Definition: vpMath.h:137
static void print(vpHinkleyJumpType jump)
Definition: vpHinkley.cpp:437
void init()
Definition: vpHinkley.cpp:149
#define vpCDEBUG(niv)
Definition: vpDebug.h:478
vpHinkleyJumpType testDownwardJump(double signal)
Definition: vpHinkley.cpp:196
void setAlpha(double alpha)
Definition: vpHinkley.cpp:181
void setDelta(double delta)
Definition: vpHinkley.cpp:169