diff --git a/.github/workflows/ubuntu-sanitizers.yml b/.github/workflows/ubuntu-sanitizers.yml index 9894441bac..e215a425c0 100644 --- a/.github/workflows/ubuntu-sanitizers.yml +++ b/.github/workflows/ubuntu-sanitizers.yml @@ -82,4 +82,4 @@ jobs: # SUMMARY: AddressSanitizer: odr-violation: global 'ALIGNMENT' at /home/runner/work/visp/visp/3rdparty/simdlib/Simd/SimdLib.cpp:82:18 ASAN_OPTIONS: detect_odr_violation=0 working-directory: build - run: ctest -j$(nproc) --output-on-failure + run: ctest -j$(nproc) --output-on-failure -V diff --git a/doc/image/tutorial/misc/img-tutorial-spc-run.jpg b/doc/image/tutorial/misc/img-tutorial-spc-run.jpg new file mode 100644 index 0000000000..ccbcce61ea Binary files /dev/null and b/doc/image/tutorial/misc/img-tutorial-spc-run.jpg differ diff --git a/doc/tutorial/misc/tutorial-spc.dox b/doc/tutorial/misc/tutorial-spc.dox new file mode 100644 index 0000000000..d4df27c12c --- /dev/null +++ b/doc/tutorial/misc/tutorial-spc.dox @@ -0,0 +1,89 @@ +/** + \page tutorial-spc Tutorial: Using Statistical Process Control to monitor your signal + \tableofcontents + +\section tuto-spc-intro Introduction + +Statistical Process Control (SPC) is defined as the use of statistical methods to monitor +if a signal is "in control". + +In this tutorial, we will use a Statistical Process Control method to monitor if a +random signal following a normal distribution is "in control". + +\subsection tuto-spc-intro-methods Available methods + +The different methods available in ViSP aim at detecting if the mean of a signal is +changing, either due to an abrupt jump or due to a slow drift. + +The different methods that are available are the following: +- *Exponentially Weighted Moving Average* (EWMA), implementent in the vpStatisticalTestEWMA class. +- *Hinkley's test*, implemented in the vpStatisticalTestHinkley class. +- *Mean Adjusted Cumulative Sum* (mean adjusted CUSUM), implemented in the vpStatisticalTestMeanAdjustedCUSUM class. +- *Shewhart's test*, implemented in the vpStatisticalTestShewhart class. +- *Sigma test*, implemented in the vpStatisticalTestSigma class. + +We refer the reader to the documentation of each class to have more detailed information on +each method. + +\section tuto-spc-tutorial Explanations about the tutorial + +\subsection tuto-spc-tutorial-howtorun How to run the tutorial + +To see the different options of the tutorial, please run the following commands: + +``` +$ cd $VISP_WS/visp-build/tutorial/mean-drift +$ ./tutorial-meandrift -h +``` + +If you run the program without argument, you should see something similar to the following image: + +\image html img-tutorial-spc-run.jpg + +A Gaussian signal of mean equal to 6 and of standard deviation equal to 2 is generated, +without any mean drift. The program tells you which method has been chosen in the +console, and which are its parameters. A monitoring loop stops once an alarm +is raised. When the alarm is raised, some information about the alarm and the +test signal(s) + limits of the SPC method are given. Press `Return` to leave the program. + +\subsection tuto-spc-tutorial-explained Detailed explanations about the SPC tutorial + +For this tutorial, we use the main program tutorial-meandrift.cpp . + +It uses the following enumeration to let the user choose which SPC method to use to monitor +the signal: + +\snippet tutorial-meandrift.cpp Enum_For_Test_Choice + +The program arguments are parsed and fill the following structure that stores the SPC methods +parameters: + +\snippet tutorial-meandrift.cpp Structure_Parameters + +It is possible to choose to monitor only upward mean drifts, only downward mean drifts or both. +To do so, use the `--alarms` option with the name of the alarm(s) you want to monitor. + +First, the plot that will show the data is created in the following section of code: + +\snippet tutorial-meandrift.cpp Plot_Init + +Then, the desired method is created in the following section of code, with the desired parameters: + +\snippet tutorial-meandrift.cpp Test_Creat + +Then, the method is filled with "in control" signals to initialize the expected mean and the standard deviation: + +\snippet tutorial-meandrift.cpp Test_Init + +Then, the monitoring loop is run, where the signal is randomly generated with a potential mean drift +(if desired). This random signal is then tested, and the loop is stopped if an alarm we +want to monitor is raised: + +\snippet tutorial-meandrift.cpp Loop_Monitor + +Finally, some information about why the loop was stopped is displayed in the console: + +\snippet tutorial-meandrift.cpp Failure_Debrief + +The program stops once the `Return` key is pressed. +*/ diff --git a/doc/tutorial/tutorial-users.dox b/doc/tutorial/tutorial-users.dox index 43e6e80996..9dc2c0359f 100644 --- a/doc/tutorial/tutorial-users.dox +++ b/doc/tutorial/tutorial-users.dox @@ -169,6 +169,7 @@ This page introduces the user to other tools that may be useful. - \subpage tutorial-pcl-viewer
This tutorial explains how to use a threaded PCL-based point cloud visualizer. - \subpage tutorial-json
This tutorial explains how to read and save data in the portable JSON format. It focuses on saving the data generated by a visual servoing experiment and exporting it to Python in order to generate plots. - \subpage tutorial-synthetic-blenderproc
This tutorial shows you how to easily generate synthetic data from the 3D model of an object and obtain various modalities. This data can then be used to train a neural network for your own task. +- \subpage tutorial-spc
This tutorial shows you how to monitor if a signal is "in control" using Statistical Process Control methods. */ /*! \page tutorial_munkres Munkres Assignment Algorithm diff --git a/modules/core/include/visp3/core/vpHinkley.h b/modules/core/include/visp3/core/vpHinkley.h index d24760c2d6..6d9ee47e07 100644 --- a/modules/core/include/visp3/core/vpHinkley.h +++ b/modules/core/include/visp3/core/vpHinkley.h @@ -40,6 +40,7 @@ */ #include +#if defined(VISP_BUILD_DEPRECATED_FUNCTIONS) /*! \class vpHinkley @@ -88,13 +89,14 @@ N_{k^{'}} = 0 \f$. */ -class VISP_EXPORT vpHinkley +class vp_deprecated vpHinkley { public: /*! \enum vpHinkleyJumpType Indicates if a jump is detected by the Hinkley test. */ - typedef enum { + typedef enum + { noJump, /*!< No jump is detected by the Hinkley test. */ downwardJump, /*!< A downward jump is detected by the Hinkley test. */ upwardJump /*!< An upward jump is detected by the Hinkley test. */ @@ -162,5 +164,5 @@ class VISP_EXPORT vpHinkley double Tk; double Nk; }; - +#endif #endif diff --git a/modules/core/include/visp3/core/vpStatisticalTestAbstract.h b/modules/core/include/visp3/core/vpStatisticalTestAbstract.h new file mode 100644 index 0000000000..dbf9604a2d --- /dev/null +++ b/modules/core/include/visp3/core/vpStatisticalTestAbstract.h @@ -0,0 +1,266 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ +/*! + * \file vpStatisticalTestAbstract.h + * \brief Base class for Statistical Process Control methods implementation. + */ + +#ifndef _vpStatisticalTestAbstract_h_ +#define _vpStatisticalTestAbstract_h_ + +#include +#include +#include +#include + +#include + +/** + * \ingroup group_core_math_tools + * \brief Base class for methods detecting the drift of the mean of a process. + * + * To detect only downward drifts of the input signal \f$ s(t) \f$ use + * testDownwardMeanDrift().To detect only upward drifts in \f$ s(t) \f$ use + * testUpwardMeanDrift(). To detect both, downward and upward drifts use + * testDownUpwardMeanDrift(). + */ +class VISP_EXPORT vpStatisticalTestAbstract +{ +public: + /** + * \brief Enum that indicates if a drift of the mean occurred. + */ + typedef enum vpMeanDriftType + { + MEAN_DRIFT_NONE = 0, /*!< No mean drift occurred*/ + MEAN_DRIFT_DOWNWARD = 1, /*!< A downward drift of the mean occurred.*/ + MEAN_DRIFT_UPWARD = 2, /*!< An upward drift of the mean occurred.*/ + MEAN_DRIFT_BOTH = 3, /*!< Both an aupward and a downward drifts occurred.*/ + MEAN_DRIFT_COUNT = 4, + MEAN_DRIFT_UNKNOWN = MEAN_DRIFT_COUNT + } vpMeanDriftType; + + /** + * \brief Cast a \b vpMeanDriftType into a string. + * + * \param[in] type The type of mean drift. + * \return std::string The corresponding message. + */ + static std::string vpMeanDriftTypeToString(const vpMeanDriftType &type); + + /** + * \brief Cast a string into a \b vpMeanDriftType. + * + * \param[in] name The name of the mean drift. + * \return vpMeanDriftType The corresponding \b vpMeanDriftType. + */ + static vpMeanDriftType vpMeanDriftTypeFromString(const std::string &name); + + /** + * \brief Get the list of available \b vpMeanDriftType objects that are handled. + * + * \param[in] prefix The prefix that should be placed before the list. + * \param[in] sep The separator between each element of the list. + * \param[in] suffix The suffix that should terminate the list. + * \return std::string The list of handled type of process tests, presented as a string. + */ + static std::string getAvailableMeanDriftType(const std::string &prefix = "<", const std::string &sep = " , ", + const std::string &suffix = ">"); + + /** + * \brief Print the message corresponding to the type of mean drift. + * + * \param[in] type The type of mean drift. + */ + static void print(const vpMeanDriftType &type); + +protected: + bool m_areStatisticsComputed; /*!< Set to true once the mean and the standard deviation are available.*/ + float m_count; /*!< Current number of data used to compute the mean and the standard deviation.*/ + float m_limitDown; /*!< Upper limit for the test signal m_wt.*/ + float m_limitUp; /*!< Lower limit for the test signal m_wt*/ + float m_mean; /*!< Mean of the monitored signal.*/ + unsigned int m_nbSamplesForStatistics; /*!< Number of samples to use to compute the mean and the standard deviation.*/ + float *m_s; /*!< Array that keeps the samples used to compute the mean and standard deviation.*/ + float m_stdev; /*!< Standard deviation of the monitored signal.*/ + float m_stdevmin; /*!< Minimum allowed standard deviation of the monitored signal.*/ + float m_sumForMean; /*!< Sum of the samples used to compute the mean and standard deviation.*/ + + /** + * \brief Detects if a downward mean drift occurred. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_DOWNWARD if a downward mean drift occurred, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectUpwardMeanDrift() + */ + virtual vpMeanDriftType detectDownwardMeanDrift() = 0; + + /** + * \brief Detects if a upward mean drift occurred. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_UPWARD if an upward mean drift occurred, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectDownwardMeanDrift() + */ + virtual vpMeanDriftType detectUpwardMeanDrift() = 0; + + /** + * \brief Update \b m_s and if enough values are available, compute the mean, the standard + * deviation and the limits. + * + * \param[in] signal The new value of the signal to monitor. + * + * \return true if the statistics have been computed, false if data are missing. + */ + virtual bool updateStatistics(const float &signal); + + /** + * \brief Update the test signals. + * + * \param[in] signal The new value of the signal to monitor. + */ + virtual void updateTestSignals(const float &signal) = 0; +public: + /** + * \brief Construct a new vpStatisticalTestAbstract object. + */ + vpStatisticalTestAbstract(); + + /** + * \brief Construct by copy a new vpStatisticalTestAbstract object. + */ + vpStatisticalTestAbstract(const vpStatisticalTestAbstract &other); + + /** + * \brief Destroy the vpStatisticalTestAbstract object. + */ + virtual ~vpStatisticalTestAbstract(); + + /** + * \brief Get the upper and lower limits of the test signal. + * + * \param[out] limitDown The lower limit. + * \param[out] limitUp The upper limit. + */ + inline void getLimits(float &limitDown, float &limitUp) const + { + limitDown = m_limitDown; + limitUp = m_limitUp; + } + + /** + * \brief Get the mean used as reference. + * + * \return float The mean. + */ + inline float getMean() const + { + return m_mean; + } + + /** + * \brief Get the standard deviation used as reference. + * + * \return float The standard deviation. + */ + inline float getStdev() const + { + return m_stdev; + } + + /** + * \brief (Re)Initialize the algorithm. + */ + void init(); + + /** + * \brief Copy operator of a vpStatisticalTestAbstract. + * + * \param[in] other The vpStatisticalTestAbstract to copy. + * \return vpStatisticalTestAbstract& *this after copy. + */ + vpStatisticalTestAbstract &operator=(const vpStatisticalTestAbstract &other); + + /** + * \brief Set the minimum value of the standard deviation that is expected. + * The computed standard deviation cannot be lower this value if set. + * + * \param[in] stdevmin The minimum value of the standard deviation that is expected. + */ + void setMinStdev(const float &stdevmin) + { + m_stdevmin = stdevmin; + } + + /** + * \brief Set the number of samples required to compute the mean and standard deviation + * of the signal and allocate the memory accordingly. + * + * \param[in] nbSamples The number of samples we want to use. + */ + void setNbSamplesForStat(const unsigned int &nbSamples); + + /** + * \brief Test if a downward or an upward mean drift occurred + * according to the new value of the signal. + * + * \param[in] signal The new value of the signal. + * \return vpMeanDriftType The type of mean drift that occurred. + * + * \sa testDownwardMeanDrift() testUpwardMeanDrift() + */ + vpMeanDriftType testDownUpwardMeanDrift(const float &signal); + + /** + * \brief Test if a downward mean drift occurred + * according to the new value of the signal. + * + * \param[in] signal The new value of the signal. + * \return vpMeanDriftType The type of mean drift that occurred. + * + * \sa testUpwardMeanDrift() + */ + vpMeanDriftType testDownwardMeanDrift(const float &signal); + + /** + * \brief Test if an upward mean drift occurred + * according to the new value of the signal. + * + * \param[in] signal The new value of the signal. + * \return vpMeanDriftType The type of mean drift that occurred. + * + * \sa testDownwardMeanDrift() + */ + vpMeanDriftType testUpwardMeanDrift(const float &signal); +}; + +#endif diff --git a/modules/core/include/visp3/core/vpStatisticalTestEWMA.h b/modules/core/include/visp3/core/vpStatisticalTestEWMA.h new file mode 100644 index 0000000000..134debf418 --- /dev/null +++ b/modules/core/include/visp3/core/vpStatisticalTestEWMA.h @@ -0,0 +1,182 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ +/*! + * \file vpStatisticalTestEWMA.h + * \brief Statistical Process Control Exponentially Weighted Moving Average implementation. + */ + +#ifndef _vpStatisticalTestEWMA_h_ +#define _vpStatisticalTestEWMA_h_ + +#include + +#include + +/** + * \ingroup group_core_math_tools + * \brief Class that permits to perform Exponentially Weighted Moving Average mean drft tests. + * + * The EWMA test is designed to detect drift in the mean \f$ \mu \f$ + * of an observed signal \f$ s(t) \f$. + * + * The test signal \f$ w(t) \f$ is computed as follow: + * + * \f$ w(0) = \mu \f$ + * + * \f$ w(t) = \alpha s(t) + ( 1 - \alpha ) * w(t-1) \f$ + * + * Be \f$ \sigma \f$ the standard deviation of the input signal \f$ s(t) \f$. + * + * A downward alarm is raised if: + * \f$ w(t) <= \mu - 3 * \sigma * \sqrt{ \frac{\alpha}{2 - \alpha}}\f$ + * + * An upward alarm is raised if: + * \f$ w(t) >= \mu + 3 * \sigma * \sqrt{ \frac{\alpha}{2 - \alpha}}\f$ + * + * To detect only downward drifts of the input signal \f$ s(t) \f$ use + * testDownwardMeanDrift().To detect only upward drifts in \f$ s(t) \f$ use + * testUpwardMeanDrift(). To detect both, downward and upward drifts use + * testDownUpwardMeanDrift(). + */ +class VISP_EXPORT vpStatisticalTestEWMA : public vpStatisticalTestAbstract +{ +protected: + float m_alpha; /*!< Forgetting factor: the higher, the more weight the current signal value has.*/ + float m_wt; /*!< Test signal that permits to raise an alarm.*/ + float m_wtprev; /*!< Previous value of the test signal.*/ + + /** + * \brief Compute the upper and lower limits of the test signal. + */ + virtual void computeDeltaAndLimits(); + + /** + * \brief Detects if a downward mean drift occured. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_DOWNWARD if a downward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectUpwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectDownwardMeanDrift() override; +#else + virtual vpMeanDriftType detectDownwardMeanDrift(); +#endif + + /** + * \brief Detects if an upward mean drift occured on the mean. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_UPWARD if an upward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectDownwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectUpwardMeanDrift() override; +#else + virtual vpMeanDriftType detectUpwardMeanDrift(); +#endif + + /** + * \brief Update m_s and if enough values are available, compute the mean, the standard + * deviation and the limits. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual bool updateStatistics(const float &signal) override; +#else + virtual bool updateStatistics(const float &signal); +#endif + + /** + * \brief Update the test signals. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual void updateTestSignals(const float &signal) override; +#else + virtual void updateTestSignals(const float &signal); +#endif + +public: + /** + * \brief Construct a new vpStatisticalTestEWMA object. + * + * \param[in] alpha The forgetting factor. + */ + vpStatisticalTestEWMA(const float &alpha = 0.1f); + + /** + * \brief Get the forgetting factor of the algorithm. + * + * \return float The forgetting factor. + */ + inline float getAlpha() const + { + return m_alpha; + } + + /** + * \brief Get the current value of the test signal. + * + * \return float The current value of the test signal. + */ + inline float getWt() const + { + return m_wt; + } + + /** + * \brief Initialize the EWMA algorithm. + * + * \param[in] alpha The forgetting factor. + */ + void init(const float &alpha); + + /** + * \brief Initialize the EWMA algorithm. + * + * \param[in] alpha The forgetting factor. + * \param[in] mean The expected mean of the signal to monitor. + * \param[in] stdev The expected standard deviation of the signal to monitor. + */ + void init(const float &alpha, const float &mean, const float &stdev); + + /** + * \brief Set the forgetting factor. + * + * \param[in] alpha The forgetting factor. + */ + void setAlpha(const float &alpha); +}; +#endif diff --git a/modules/core/include/visp3/core/vpStatisticalTestHinkley.h b/modules/core/include/visp3/core/vpStatisticalTestHinkley.h new file mode 100644 index 0000000000..e07d4b1591 --- /dev/null +++ b/modules/core/include/visp3/core/vpStatisticalTestHinkley.h @@ -0,0 +1,333 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ +/*! + * \file vpStatisticalTestHinkley.h + * \brief Statistical Process Control Hinkley's test implementation. + */ + +#ifndef _vpStatisticalTestHinkley_h_ +#define _vpStatisticalTestHinkley_h_ + +#include + +#include + +/** + * \ingroup group_core_math_tools + * \brief This class implements the Hinkley's cumulative sum test. + * + * The Hinkley's cumulative sum test is designed to detect drift in the mean + of an observed signal \f$ s(t) \f$. It is known to be robust (by + taking into account all the past of the observed quantity), + efficient, and inducing a very low computational load. The other + attractive features of this test are two-fold. First, it can + straightforwardly and accurately provide the drift instant. Secondly, + due to its formulation (cumulative sum test), it can simultaneously + handle both very abrupt and important changes, and gradual smaller + ones without adapting the involved thresholds. + + Two tests are performed in parallel to look for downwards or upwards + drifts in \f$ s(t) \f$, respectively defined by: + + \f[ S_k = \sum_{t=0}^{k} (s(t) - m_0 + \frac{\delta}{2}) \f] + \f[ M_k = \max_{0 \leq i \leq k} S_i\f] + \f[ T_k = \sum_{t=0}^{k} (s(t) - m_0 - \frac{\delta}{2}) \f] + \f[ N_k = \min_{0 \leq i \leq k} T_i\f] + + In which \f$m_o\f$ is computed on-line and corresponds to the mean + of the signal \f$ s(t) \f$ we want to detect a drift. \f$m_o\f$ is + re-initialized at zero after each drift detection. \f$\delta\f$ + denotes the drift minimal magnitude that we want to detect and + \f$\alpha\f$ is a predefined threshold. These values are set by + default to 0.2 in the default constructor vpHinkley(). To modify the + default values use setAlpha() and setDelta() or the + vpHinkley(double alpha, double delta) constructor. + + A downward drift is detected if \f$ M_k - S_k > \alpha \f$. + A upward drift is detected if \f$ T_k - N_k > \alpha \f$. + + To detect only downward drifts in \f$ s(t) \f$ use + testDownwardMeanDrift().To detect only upward drifts in \f$ s(t) \f$ use + testUpwardMeanDrift(). To detect both, downward and upward drifts use + testDownUpwardMeanDrift(). + + If a drift is detected, the drift location is given by the last instant + \f$k^{'}\f$ when \f$ M_{k^{'}} - S_{k^{'}} = 0 \f$, or \f$ T_{k^{'}} - + N_{k^{'}} = 0 \f$. + */ +class VISP_EXPORT vpStatisticalTestHinkley : public vpStatisticalTestAbstract +{ +protected: + float m_dmin2; /*!< Half of \f$\delta\f$, the drift minimal magnitude that we want to detect.*/ + float m_alpha; /*!< The \f$\alpha\f$ threshold indicating that a mean drift occurs. */ + float m_Sk; /*!< Test signal for downward mean drift.*/ + float m_Mk; /*!< Maximum of the test signal for downward mean drift \f$S_k\f$ .*/ + float m_Tk; /*!< Test signal for upward mean drift.*/ + float m_Nk; /*!< Minimum of the test signal for upward mean drift \f$T_k\f$*/ + bool m_computeDeltaAndAlpha; /*!< If true, compute \f$\delta\f$ and \f$\alpha\f$ from the standard deviation, + the alarm factor and the detection factor.*/ + float m_h; /*!< The alarm factor, that permits to compute \f$\alpha\f$ from the standard deviation of the signal.*/ + float m_k; /*!< The detection factor, that permits to compute \f$\delta\f$ from the standard deviation of the signal.*/ + + /** + * \brief Compute \f$\delta\f$ and \f$\alpha\f$ from the standard deviation of the signal. + */ + virtual void computeAlphaDelta(); + + /** + * \brief Compute the mean value \f$m_0\f$ of the signal. The mean value must be + * computed before the mean drift is estimated on-line. + * + * \param[in] signal The new value of the signal to monitor. + */ + void computeMean(double signal); + + /** + * \brief Compute \f$S_k = \sum_{t=0}^{k} (s(t) - m_0 + \frac{\delta}{2})\f$ + * + * \param[in] signal The new value of the signal to monitor. + */ + void computeSk(double signal); + + /** + * \brief Compute \f$M_k\f$, the maximum value of \f$S_k\f$. + */ + void computeMk(); + + /** + * \brief Compute \f$T_k = \sum_{t=0}^{k} (s(t) - m_0 - \frac{\delta}{2})\f$ + * + * \param[in] signal The new value of the signal to monitor. + */ + void computeTk(double signal); + + /** + * \brief Compute \f$N_k\f$, the minimum value of \f$T_k\f$. + */ + void computeNk(); + + /** + * \brief Detects if a downward mean drift occurred. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_DOWNWARD if a downward mean drift occurred, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectDownwardMeanDrift() override; +#else + virtual vpMeanDriftType detectDownwardMeanDrift(); +#endif + + /** + * \brief Detects if an upward mean drift occured on the mean. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_UPWARD if an upward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectDownwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectUpwardMeanDrift() override; +#else + virtual vpMeanDriftType detectUpwardMeanDrift(); +#endif + + /** + * \brief Update m_s and if enough values are available, compute the mean, the standard + * deviation and the limits. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual bool updateStatistics(const float &signal) override; +#else + virtual bool updateStatistics(const float &signal); +#endif + +/** + * \brief Update the test signals. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual void updateTestSignals(const float &signal) override; +#else + virtual void updateTestSignals(const float &signal); +#endif +public: + /** + * @brief Construct a new vpStatisticalTestHinkley object. + * Call init() to initialise the Hinkley's test and set \f$\alpha\f$ + * and \f$\delta\f$ to default values. + * + * By default \f$ \delta = 0.2 \f$ and \f$ \alpha = 0.2\f$. Use + * setDelta() and setAlpha() to modify these values. + */ + vpStatisticalTestHinkley(); + + /** + * \brief Call init() to initialise the Hinkley's test and set \f$\alpha\f$ + * and \f$\delta\f$ thresholds. + * \param[in] alpha : \f$\alpha\f$ threshold indicating that a mean drift occurs. + * \param[in] delta : \f$\delta\f$ denotes the drift minimal magnitude that + * we want to detect. + * \param[in] nbSamplesForInit : number of signal samples to initialize the mean of the signal. + * + * \sa setAlpha(), setDelta() + */ + vpStatisticalTestHinkley(const float &alpha, const float &delta, const unsigned int &nbSamplesForInit = 30); + + /** + * \brief Construct a new vpStatisticalTestHinkley object. \f$\alpha\f$ and \f$\delta\f$ will be computed + * from the standard deviation of the signal. + * + * \param[in] h : the alarm factor that permits to compute \f$\alpha\f$ from the standard deviation. + * \param[in] k : the detection factor that permits to compute \f$\delta\f$ from the standard deviation. + * \param[in] computeAlphaDeltaFromStdev : must be equal to true, otherwise throw a vpException. + * \param[in] nbSamplesForInit : number of signal samples to initialize the mean of the signal. + */ + vpStatisticalTestHinkley(const float &h, const float &k, const bool &computeAlphaDeltaFromStdev, const unsigned int &nbSamplesForInit = 30); + + /** + * \brief Construct a new vpStatisticalTestHinkley object. \f$\alpha\f$ and \f$\delta\f$ will be computed + * from the standard deviation of the signal. + * + * \param[in] h : the alarm factor that permits to compute \f$\alpha\f$ from the standard deviation. + * \param[in] k : the detection factor that permits to compute \f$\delta\f$ from the standard deviation. + * \param[in] mean : the expected mean of the signal. + * \param[in] stdev : the expected standard deviation of the signal. + */ + vpStatisticalTestHinkley(const float &h, const float &k, const float &mean, const float &stdev); + + /** + * \brief Destroy the vpStatisticalTestHinkley object. + */ + virtual ~vpStatisticalTestHinkley(); + + /** + * \brief Get the \f$\alpha\f$ threshold indicating that a mean drift occurs. + * + * \return The \f$\alpha\f$ threshold. + */ + inline float getAlpha() const { return m_alpha; } + + /*! + * \brief Get the test signal for downward mean drift. + * + * \return The value of \f$S_k = \sum_{t=0}^{k} (s(t) - m_0 + \frac{\delta}{2})\f$ . + */ + inline float getSk() const { return m_Sk; } + + /*! + * \brief Get the maximum of the test signal for downward mean drift \f$S_k\f$ . + * + * \return The value of \f$M_k\f$, the maximum value of \f$S_k\f$. + */ + inline float getMk() const { return m_Mk; } + + /*! + * \brief Get the test signal for upward mean drift.. + * + * \return The value of \f$T_k = \sum_{t=0}^{k} (s(t) - m_0 - \frac{\delta}{2})\f$ . + + */ + inline float getTk() const { return m_Tk; } + + /*! + * \brief Get the minimum of the test signal for upward mean drift \f$T_k\f$ + * + * \return The value of \f$N_k\f$, the minimum value of \f$T_k\f$. + */ + inline float getNk() const { return m_Nk; } + + /** + * \brief Initialise the Hinkley's test by setting the mean signal value + * \f$m_0\f$ to zero as well as \f$S_k, M_k, T_k, N_k\f$. + */ + void init(); + + /** + * \brief Call init() to initialise the Hinkley's test and set \f$\alpha\f$ + * and \f$\delta\f$ thresholds. + * + * \param[in] alpha The threshold indicating that a mean drift occurs. + * \param[in] delta The drift minimal magnitude that we want to detect. + * \param[in] nbSamplesForInit : number of signal samples to initialize the mean of the signal. + */ + void init(const float &alpha, const float &delta, const unsigned int &nbSamplesForInit); + + /** + * \brief (Re)Initialize a new vpStatisticalTestHinkley object. \f$\alpha\f$ and \f$\delta\f$ will be computed + * from the standard deviation of the signal. + * + * \param[in] h : the alarm factor that permits to compute \f$\alpha\f$ from the standard deviation. + * \param[in] k : the detection factor that permits to compute \f$\delta\f$ from the standard deviation. + * \param[in] computeAlphaDeltaFromStdev : must be equal to true, otherwise throw a vpException. + * \param[in] nbSamplesForInit : number of signal samples to initialize the mean of the signal. + */ + void init(const float &h, const float &k, const bool &computeAlphaDeltaFromStdev, const unsigned int &nbSamplesForInit); + + /** + * \brief Call init() to initialise the Hinkley's test, set \f$\alpha\f$ + * and \f$\delta\f$ thresholds, and the mean of the signal \f$m_0\f$. + * + * \param[in] alpha The threshold indicating that a mean drift occurs. + * \param[in] delta The drift minimal magnitude that we want to detect. + * \param[in] mean The expected value of the mean. + */ + void init(const float &alpha, const float &delta, const float &mean); + + /** + * \brief (Re)Initialize a new vpStatisticalTestHinkley object. \f$\alpha\f$ and \f$\delta\f$ will be computed + * from the standard deviation of the signal. + * + * \param[in] h : the alarm factor that permits to compute \f$\alpha\f$ from the standard deviation. + * \param[in] k : the detection factor that permits to compute \f$\delta\f$ from the standard deviation. + * \param[in] mean : the expected mean of the signal. + * \param[in] stdev : the expected standard deviation of the signal. + */ + void init(const float &h, const float &k, const float &mean, const float &stdev); + + /** + * \brief Set the drift minimal magnitude that we want to detect. + * + * \param[in] delta The drift magnitude. + */ + void setDelta(const float &delta); + + /** + * \brief The threshold indicating that a mean drift occurs. + * + * \param[in] alpha The threshold. + */ + void setAlpha(const float &alpha); +}; + +#endif diff --git a/modules/core/include/visp3/core/vpStatisticalTestMeanAdjustedCUSUM.h b/modules/core/include/visp3/core/vpStatisticalTestMeanAdjustedCUSUM.h new file mode 100644 index 0000000000..fdb9f09ebf --- /dev/null +++ b/modules/core/include/visp3/core/vpStatisticalTestMeanAdjustedCUSUM.h @@ -0,0 +1,274 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ +/*! + * \file vpStatisticalTestMeanAdjustedCUSUM.h + * \brief Statistical Process Control mean adjusted CUSUM implementation. + */ + +#ifndef _vpStatisticalTestMeanAdjustedCUSUM_h_ +#define _vpStatisticalTestMeanAdjustedCUSUM_h_ + +#include + +#include + +/** + * \ingroup group_core_math_tools + * \brief Class that permits to perform a mean adjusted Cumulative Sum test. + * + * The mean adjusted CUSUM test is designed to detect drift in the mean \f$ \mu \f$ + * of an observed signal \f$ s(t) \f$. + * + * Be \f$ \delta \f$ the amplitude of the mean drift we want to detect. + * Two test signals are computed at each iteration: + * + * \f$ S_-(t) = max\{0, S_-(t-1) - (s(t) - \mu) - \frac{\delta}{2}\} \f$ + * + * \f$ S_+(t) = max\{0, S_+(t-1) + (s(t) - \mu) - \frac{\delta}{2}\} \f$ + * + * A downward alarm is raised if: + * \f$ S_-(t) >= thresh\f$ + * + * An upward alarm is raised if: + * \f$ S_+(t) >= thresh\f$ + * + * To ease the understanding of the detection threshold \f$ \delta \f$ and the + * alarm threshold \f$ thresh \f$, ViSP implemented these two thresholds as + * a multiple of the standard deviation of the signal \f$ \sigma \f$: + * + * \f$ \delta = k \sigma , k \in R^{+*} \f$ + * + * \f$ thresh = h \sigma , h \in R^{+*} \f$ + * + * To have an Average Run Lenght of ~374 samples for a detection threshold \f$ \delta \f$ + * of 1 standard deviation \f$ \sigma \f$, set \f$ h \f$ to 4.76 . + * + * To detect only downward drifts of the input signal \f$ s(t) \f$ use + * testDownwardMeanDrift().To detect only upward drifts in \f$ s(t) \f$ use + * testUpwardMeanDrift(). To detect both, downward and upward drifts use + * testDownUpwardMeanDrift(). + */ +class VISP_EXPORT vpStatisticalTestMeanAdjustedCUSUM : public vpStatisticalTestAbstract +{ +protected: + float m_delta; /*!< Slack of the CUSUM test, i.e. amplitude of mean shift we want to be able to detect.*/ + float m_h; /*!< Alarm factor that permits to determine the limit telling when a mean shift occurs: \f$thresh = h * \sigma \f$ . + To have an Average Run Lenght of ~374 samples for a detection of 1 stdev, set it to 4.76f*/ + float m_half_delta; /*!< Half of the amplitude we want to detect.*/ + float m_k; /*!< Detection factor that permits to determine the slack: \f$\delta = k * \sigma\f$ .*/ + float m_sminus; /*!< Test signal for downward mean shift: \f$ S_-(t) = max\{0, S_-(t-1) - (s(t) - \mu) - \frac{\delta}{2}\} \f$.*/ + float m_splus; /*!< Test signal for upward mean shift: \f$ S_+(t) = max\{0, S_+(t-1) + (s(t) - \mu) - \frac{\delta}{2}\} \f$.*/ + + /** + * \brief Compute the upper and lower limits of the test signal. + */ + virtual void computeDeltaAndLimits(); + + /** + * \brief Detects if a downward mean drift occured. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_DOWNWARD if a downward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectUpwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectDownwardMeanDrift() override; +#else + virtual vpMeanDriftType detectDownwardMeanDrift(); +#endif + +/** + * \brief Detects if a upward jump occurred on the mean. + * + * \return upwardJump if a upward jump occurred, noJump otherwise. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectUpwardMeanDrift() override; +#else + virtual vpMeanDriftType detectUpwardMeanDrift(); +#endif + + /** + * \brief Update m_s and if enough values are available, compute the mean, the standard + * deviation and the limits. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual bool updateStatistics(const float &signal) override; +#else + virtual bool updateStatistics(const float &signal); +#endif + + /** + * \brief Update the test signals. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual void updateTestSignals(const float &signal) override; +#else + virtual void updateTestSignals(const float &signal); +#endif +public: + /** + * \brief Construct a new vpStatisticalTestMeanAdjustedCUSUM object. + * + * \param[in] h The alarm factor that permits to determine when the process is out of control from the standard + * deviation of the signal. + * \param[in] k The detection factor that permits to determine the slack of the CUSUM test, i.e. the + * minimum value of the jumps we want to detect, from the standard deviation of the signal. + * \param[in] nbPtsForStats The number of samples to use to compute the mean and the standard deviation of the signal + * to monitor. + */ + vpStatisticalTestMeanAdjustedCUSUM(const float &h = 4.f, const float &k = 1.f, const unsigned int &nbPtsForStats = 30); + + /** + * \brief Get the slack of the CUSUM test, + * i.e. amplitude of mean shift we want to be able to detect. + * + * \return float The slack of the CUSUM test. + */ + inline float getDelta() const + { + return m_delta; + } + + /** + * \brief Get the alarm factor. + * + * \return float The alarm factor. + */ + inline float getH() const + { + return m_h; + } + + /** + * \brief Get the detection factor. + * + * \return float The detection factor. + */ + inline float getK() const + { + return m_k; + } + + /** + * \brief Get the latest value of the test signal for downward jumps of the mean. + * + * \return float Its latest value. + */ + inline float getTestSignalMinus() const + { + return m_sminus; + } + + /** + * \brief Get the latest value of the test signal for upward jumps of the mean. + * + * \return float Its latest value. + */ + inline float getTestSignalPlus() const + { + return m_splus; + } + + /** + * \brief (Re)Initialize the mean adjusted CUSUM test. + * + * \param[in] h The alarm factor that permits to determine when the process is out of control from the standard + * deviation of the signal. + * \param[in] k The detection factor that permits to determine the slack of the CUSUM test, i.e. the + * minimum value of the jumps we want to detect, from the standard deviation of the signal. + * \param[in] nbPtsForStats The number of points to use to compute the mean and the standard deviation of the signal + */ + void init(const float &h, const float &k, const unsigned int &nbPtsForStats); + + /** + * \brief Initialize the mean adjusted CUSUM test. + * + * \param[in] delta The slack of the CUSUM test, i.e. the minimum value of the jumps we want to detect. + * \param[in] limitDown The lower limit of the CUSUM test, for the downward jumps. + * \param[in] limitUp The upper limit of the CUSUM test, for the upward jumps. + * \param[in] nbPtsForStats The number of points to use to compute the mean and the standard deviation of the signal + * to monitor. + */ + void init(const float &delta, const float &limitDown, const float &limitUp, const unsigned int &nbPtsForStats); + + /** + * \brief Initialize the mean adjusted CUSUM test. + * + * \param[in] h The alarm factor that permits to determine when the process is out of control from the standard + * deviation of the signal. + * \param[in] k The detection factor that permits to determine the slack of the CUSUM test, i.e. the + * minimum value of the jumps we want to detect, from the standard deviation of the signal. + * \param[in] mean The expected mean of the signal to monitor. + * \param[in] stdev The expected standard deviation of the signal to monitor. + */ + void init(const float &h, const float &k, const float &mean, const float &stdev); + + /** + * \brief Initialize the mean adjusted CUSUM test. + * + * \param[in] delta The slack of the CUSUM test, i.e. the minimum value of the jumps we want to detect. + * \param[in] limitDown The lower limit of the CUSUM test, for the downward jumps. + * \param[in] limitUp The upper limit of the CUSUM test, for the upward jumps. + * \param[in] mean The expected mean of the signal to monitor. + * \param[in] stdev The expected standard deviation of the signal to monitor. + */ + void init(const float &delta, const float &limitDown, const float &limitUp, const float &mean, const float &stdev); + + /** + * \brief Set the slack of the CUSUM test, i.e. the minimum value of the jumps we want to detect. + * + * \param[in] delta The slack of the CUSUM test. + */ + inline void setDelta(const float &delta) + { + m_delta = delta; + m_half_delta = 0.5f * delta; + } + + /** + * \brief Set the upward and downward jump limits. + * + * \param[in] limitDown The limit for the downward jumps. + * \param[in] limitUp The limit for the upward jumps. + */ + inline void setLimits(const float &limitDown, const float &limitUp) + { + m_limitDown = limitDown; + m_limitUp = limitUp; + } +}; +#endif diff --git a/modules/core/include/visp3/core/vpStatisticalTestShewhart.h b/modules/core/include/visp3/core/vpStatisticalTestShewhart.h new file mode 100644 index 0000000000..32c543086a --- /dev/null +++ b/modules/core/include/visp3/core/vpStatisticalTestShewhart.h @@ -0,0 +1,242 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ +/*! + * \file vpStatisticalTestShewhart.h + * \brief Statistical Process Control Shewhart's test implementation. + */ + +#ifndef _vpStatisticalTestShewhartTest_h_ +#define _vpStatisticalTestShewhartTest_h_ + +#include + +#include +#include + +/** + * \ingroup group_core_math_tools + * \brief Class that permits a Shewhart's test. + * + * Be \f$ s(t) \f$ the signal to monitor, \f$ \mu \f$ and \f$ \sigma \f$ the mean and standard deviation + * of this signal when it is "in control". + * + * A downward alarm is raised if: + * \f$ s(t) >= \mu - 3 \sigma \f$ + * + * An upward alarm is raised if: + * \f$ s(t) >= \mu + 3 \sigma \f$ + * + * Additionnally, we can activate the WECO's rules that have been + * proposed by the Western Electric Company to add additionnal verifications: + * - An alarm is raised if two out of three consecutive points fall beyond the \f$2\sigma\f$-limit, on the same side of the mean \f$ \mu \f$ + * - An alarm is raised if four out of five consecutive points fall beyond the \f$1\sigma\f$-limit, on the same side of the mean \f$ \mu \f$ + * - An alarm is raised if eight consecutive points fall on the same side of the mean \f$ \mu \f$. + * + * The user can decide to use or not the WECO's rules. Additionnally, the user can choose which WECO's + * rule(s) to activate. + * + * To detect only downward drifts of the input signal \f$ s(t) \f$ use + * testDownwardMeanDrift().To detect only upward drifts in \f$ s(t) \f$ use + * testUpwardMeanDrift(). To detect both, downward and upward drifts use + * testDownUpwardMeanDrift(). + */ + +class VISP_EXPORT vpStatisticalTestShewhart : public vpStatisticalTestSigma +{ +public: + typedef enum vpWecoRulesAlarm + { + THREE_SIGMA_WECO = 0, /*!< When a \f$ 3\sigma \f$ alarm was raised.*/ + TWO_SIGMA_WECO = 1, /*!< When a \f$ 2\sigma \f$ alarm was raised.*/ + ONE_SIGMA_WECO = 2, /*!< When a \f$ 1\sigma \f$ alarm was raised*/ + SAME_SIDE_WECO = 3, /*!< When a alarm raised when 8 consecutive points lie on the same side of the mean \f$ \mu \f$ was raised.*/ + NONE_WECO = 4, /*!< When no WECO's rule alarm was raised.*/ + COUNT_WECO = 5 /*!< Number of WECO's rules that are implemented.*/ + } vpWecoRulesAlarm; + + static std::string vpWecoRulesAlarmToString(const vpWecoRulesAlarm &alarm); + + static const bool CONST_ALL_WECO_ACTIVATED[COUNT_WECO - 1]; + static const int NB_DATA_SIGNAL = 8; + +protected: + unsigned int m_nbDataInBuffer; /*!< Indicate how many data are available in the circular buffer.*/ + float m_signal[NB_DATA_SIGNAL]; /*!< The last values of the signal.*/ + bool m_activateWECOrules; /*!< If true, activate the WECO's rules (NB: it increases the sensitivity of the Shewhart + control chart but the false alarm frequency is also increased.)*/ + bool m_activatedWECOrules[COUNT_WECO - 1]; /*!< The WECO's rules that are activated. The more are activated, the higher the + sensitivity of the Shewhart control chart is but the higher the false + alarm frequency is.*/ + int m_idCurrentData; /*!< The index of the current data in m_signal.*/ + vpWecoRulesAlarm m_alarm; /*!< The type of alarm raised due to WECO's rules.*/ + float m_oneSigmaNegLim; /*!< The \f$ \mu - \sigma \f$ threshold.*/ + float m_oneSigmaPosLim; /*!< The \f$ \mu + \sigma \f$ threshold.*/ + float m_twoSigmaNegLim; /*!< The \f$ \mu - 2 \sigma \f$ threshold.*/ + float m_twoSigmaPosLim; /*!< The \f$ \mu + 2 \sigma \f$ threshold.*/ + + /** + * \brief Compute the upper and lower limits of the test signal. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual void computeLimits() override; +#else + virtual void computeLimits(); +#endif + +/** + * \brief Detects if a downward mean drift occured. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_DOWNWARD if a downward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectUpwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectDownwardMeanDrift() override; +#else + virtual vpMeanDriftType detectDownwardMeanDrift(); +#endif + + /** + * \brief Detects if an upward mean drift occured on the mean. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_UPWARD if an upward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectDownwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectUpwardMeanDrift() override; +#else + virtual vpMeanDriftType detectUpwardMeanDrift(); +#endif + + /** + * \brief Update m_s and if enough values are available, compute the mean, the standard + * deviation and the limits. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual bool updateStatistics(const float &signal) override; +#else + virtual bool updateStatistics(const float &signal); +#endif + + /** + * \brief Update the test signals. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual void updateTestSignals(const float &signal) override; +#else + virtual void updateTestSignals(const float &signal); +#endif +public: + /** + * \brief Construct a new vpStatisticalTestShewhart object. + * + * \param[in] activateWECOrules If true, activate the WECO's rules (NB: it increases the sensitivity of the Shewhart + * control chart but the false alarm frequency is also increased.) + * \param[in] activatedRules An array where true means that the corresponding WECO's rule is activated and false means + * that it is not. + * \param[in] nbSamplesForStats The number of samples to compute the statistics of the signal. + */ + vpStatisticalTestShewhart(const bool &activateWECOrules = true, const bool activatedRules[COUNT_WECO - 1] = CONST_ALL_WECO_ACTIVATED, const unsigned int &nbSamplesForStats = 30); + + /** + * \brief Construct a new vpStatisticalTestShewhart object. + * + * \param[in] activateWECOrules If true, activate the WECO's rules (NB: it increases the sensitivity of the Shewhart + * control chart but the false alarm frequency is also increased.) + * \param[in] activatedRules An array where true means that the corresponding WECO's rule is activated and false means + * that it is not. + * \param[in] mean The expected mean of the signal. + * \param[in] stdev The expected standard deviation of the signal. + */ + vpStatisticalTestShewhart(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1], const float &mean, const float &stdev); + + /** + * \brief Get the alarm raised by the last test due to the WECO's rules. + * + * \return vpWecoRulesAlarm The type of raised alarm. + */ + vpWecoRulesAlarm getAlarm() const + { + return m_alarm; + } + + /** + * \brief Get the last value of the signal. + * + * \return float The signal. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + inline virtual float getSignal() const override +#else + inline virtual float getSignal() const +#endif + { + return m_signal[m_idCurrentData]; + } + + /** + * \brief Get the NB_DATA_SIGNAL last signal values, sorted from the latest [0] to the newest [NB_DATA_SIGNAL - 1]. + * + * \return std::vector The last NB_DATA_SIGNAL values. + */ + std::vector getSignals() const; + + /** + * \brief (Re)Initialize the test. + * + * \param[in] activateWECOrules If true, activate the WECO's rules (NB: it increases the sensitivity of the Shewhart + * control chart but the false alarm frequency is also increased.) + * \param[in] activatedRules An array where true means that the corresponding WECO's rule is activated and false means + * that it is not. + * \param[in] nbSamplesForStats The number of samples to compute the statistics of the signal. + */ + void init(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1] = CONST_ALL_WECO_ACTIVATED, const unsigned int &nbSamplesForStats = 30); + + /** + * \brief (Re)Initialize the test. + * + * \param[in] activateWECOrules If true, activate the WECO's rules (NB: it increases the sensitivity of the Shewhart + * control chart but the false alarm frequency is also increased.) + * \param[in] activatedRules An array where true means that the corresponding WECO's rule is activated and false means + * that it is not. + * \param[in] mean The expected mean of the signal. + * \param[in] stdev The expected standard deviation of the signal. + */ + void init(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1], const float &mean, const float &stdev); +}; + +#endif diff --git a/modules/core/include/visp3/core/vpStatisticalTestSigma.h b/modules/core/include/visp3/core/vpStatisticalTestSigma.h new file mode 100644 index 0000000000..3c0f10facd --- /dev/null +++ b/modules/core/include/visp3/core/vpStatisticalTestSigma.h @@ -0,0 +1,173 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ +/*! + * \file vpStatisticalTestSigma.h + * \brief Statistical Process Control sigma test implementation. + */ + +#ifndef _vpStatisticalTestSigma_h_ +#define _vpStatisticalTestSigma_h_ + +#include + +#include + +/** + * \ingroup group_core_math_tools + * \brief Class that permits a simple test comparing the current value to the + * standard deviation of the signal. + * + * Be \f$ s(t) \f$ the signal to monitor, \f$ \mu \f$ and \f$ \sigma \f$ the mean and standard deviation + * of this signal when it is "in control". + * + * Be \f$ h \f$ a user-defined alarm factor. + * + * A downward alarm is raised if: + * \f$ s(t) >= \mu - h \sigma \f$ + * + * An upward alarm is raised if: + * \f$ s(t) >= \mu - h \sigma \f$ + * + * \f$ h \f$ is often set to 3 if we assume the \f$ s(t) \f$ follows a normal distribution. + * + * To detect only downward drifts of the input signal \f$ s(t) \f$ use + * testDownwardMeanDrift().To detect only upward drifts in \f$ s(t) \f$ use + * testUpwardMeanDrift(). To detect both, downward and upward drifts use + * testDownUpwardMeanDrift(). + */ + +class VISP_EXPORT vpStatisticalTestSigma : public vpStatisticalTestAbstract +{ +protected: + float m_h; /*!< The alarm factor applied to the standard deviation to compute the limits.*/ + float m_s; /*!< The last value of the signal.*/ + + /** + * \brief Compute the upper and lower limits of the test signal. + */ + virtual void computeLimits(); + + /** + * \brief Detects if a downward mean drift occured. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_DOWNWARD if a downward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectUpwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectDownwardMeanDrift() override; +#else + virtual vpMeanDriftType detectDownwardMeanDrift(); +#endif + + /** + * \brief Detects if an upward mean drift occured on the mean. + * + * \return \b vpMeanDriftType::MEAN_DRIFT_UPWARD if an upward mean drift occured, \b vpMeanDriftType::MEAN_DRIFT_NONE otherwise. + * + * \sa detectDownwardMeanDrift() + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual vpMeanDriftType detectUpwardMeanDrift() override; +#else + virtual vpMeanDriftType detectUpwardMeanDrift(); +#endif + + /** + * \brief Update m_s and if enough values are available, compute the mean, the standard + * deviation and the limits. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual bool updateStatistics(const float &signal) override; +#else + virtual bool updateStatistics(const float &signal); +#endif + + /** + * \brief Update the test signals. + * + * \param[in] signal The new value of the signal to monitor. + */ +#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11) + virtual void updateTestSignals(const float &signal) override; +#else + virtual void updateTestSignals(const float &signal); +#endif +public: + /** + * \brief Construct a new vpStatisticalTestSigma object. + * + * \param[in] h The alarm factor applied to the standard deviation to compute the limits. + * \param[in] nbSamplesForStats The number of samples to compute the statistics of the signal. + */ + vpStatisticalTestSigma(const float &h = 3.f, const unsigned int &nbSamplesForStats = 30); + + /** + * \brief Construct a new vpStatisticalTestSigma object. + * + * \param[in] h The alarm factor applied to the standard deviation to compute the limits. + * \param[in] mean The expected mean of the signal. + * \param[in] stdev The expected standard deviation of the signal. + */ + vpStatisticalTestSigma(const float &h, const float &mean, const float &stdev); + + /** + * \brief Get the last value of the signal. + * + * \return float The signal. + */ + inline virtual float getSignal() const + { + return m_s; + } + + /** + * \brief (Re)Initialize the test. + * + * \param[in] h The alarm factor applied to the standard deviation to compute the limits. + * \param[in] nbSamplesForStats The number of samples to compute the statistics of the signal. + */ + void init(const float &h = 3.f, const unsigned int &nbSamplesForStats = 30); + + /** + * \brief (Re)Initialize the test. + * + * \param[in] h The alarm factor applied to the standard deviation to compute the limits. + * \param[in] mean The expected mean of the signal. + * \param[in] stdev The expected standard deviation of the signal. + */ + void init(const float &h, const float &mean, const float &stdev); +}; + +#endif diff --git a/modules/core/src/math/misc/vpHinkley.cpp b/modules/core/src/math/misc/vpHinkley.cpp index 40b06647ef..fc534c3e7e 100644 --- a/modules/core/src/math/misc/vpHinkley.cpp +++ b/modules/core/src/math/misc/vpHinkley.cpp @@ -42,9 +42,9 @@ */ -#include #include -//#include +#if defined(VISP_BUILD_DEPRECATED_FUNCTIONS) +#include #include #include // std::fabs @@ -70,7 +70,7 @@ setDelta() and setAlpha() to modify these values. */ -vpHinkley::vpHinkley() : dmin2(0.1), alpha(0.2), nsignal(0), mean(0), Sk(0), Mk(0), Tk(0), Nk(0) {} +vpHinkley::vpHinkley() : dmin2(0.1), alpha(0.2), nsignal(0), mean(0), Sk(0), Mk(0), Tk(0), Nk(0) { } /*! @@ -90,8 +90,7 @@ vpHinkley::vpHinkley() : dmin2(0.1), alpha(0.2), nsignal(0), mean(0), Sk(0), Mk( vpHinkley::vpHinkley(double alpha_val, double delta_val) : dmin2(delta_val / 2.), alpha(alpha_val), nsignal(0), mean(0), Sk(0), Mk(0), Tk(0), Nk(0) -{ -} +{ } /*! @@ -119,7 +118,7 @@ void vpHinkley::init(double alpha_val, double delta_val) Destructor. */ -vpHinkley::~vpHinkley() {} +vpHinkley::~vpHinkley() { } /*! @@ -305,9 +304,9 @@ vpHinkley::vpHinkleyJumpType vpHinkley::testDownUpwardJump(double signal) computeNk(); vpCDEBUG(2) << "alpha: " << alpha << " dmin2: " << dmin2 << " signal: " << signal << " Sk: " << Sk << " Mk: " << Mk - << " Tk: " << Tk << " Nk: " << Nk << std::endl; + << " Tk: " << Tk << " Nk: " << Nk << std::endl; - // teste si les variables cumulees excedent le seuil +// teste si les variables cumulees excedent le seuil if ((Mk - Sk) > alpha) jump = downwardJump; else if ((Tk - Nk) > alpha) @@ -433,3 +432,7 @@ void vpHinkley::print(vpHinkley::vpHinkleyJumpType jump) break; } } +#elif !defined(VISP_BUILD_SHARED_LIBS) + // Work around to avoid warning: libvisp_core.a(vpHinkley.cpp.o) has no symbols +void dummy_vpRHinkley() { }; +#endif diff --git a/modules/core/src/math/misc/vpStatisticalTestAbstract.cpp b/modules/core/src/math/misc/vpStatisticalTestAbstract.cpp new file mode 100644 index 0000000000..399b4b3ddd --- /dev/null +++ b/modules/core/src/math/misc/vpStatisticalTestAbstract.cpp @@ -0,0 +1,247 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +/** +* +* \file vpStatisticalTestAbstract.cpp +* +* \brief Definition of the vpStatisticalTestAbstract base class. +*/ + +#include + +std::string vpStatisticalTestAbstract::vpMeanDriftTypeToString(const vpStatisticalTestAbstract::vpMeanDriftType &type) +{ + std::string name; + switch (type) { + case MEAN_DRIFT_NONE: + name = "no_drift"; + break; + case MEAN_DRIFT_DOWNWARD: + name = "downward_drift"; + break; + case MEAN_DRIFT_UPWARD: + name = "upward_drift"; + break; + case MEAN_DRIFT_BOTH: + name = "both_drift"; + break; + case MEAN_DRIFT_UNKNOWN: + default: + name = "undefined_drift"; + break; + } + return name; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestAbstract::vpMeanDriftTypeFromString(const std::string &name) +{ + vpMeanDriftType result = MEAN_DRIFT_UNKNOWN; + unsigned int count = static_cast(MEAN_DRIFT_COUNT); + unsigned int id = 0; + bool hasNotFound = true; + while ((id < count) && hasNotFound) { + vpMeanDriftType temp = static_cast(id); + if (vpMeanDriftTypeToString(temp) == name) { + result = temp; + hasNotFound = false; + } + ++id; + } + return result; +} + +std::string vpStatisticalTestAbstract::getAvailableMeanDriftType(const std::string &prefix, const std::string &sep, + const std::string &suffix) +{ + std::string msg(prefix); + unsigned int count = static_cast(MEAN_DRIFT_COUNT); + unsigned int lastId = count - 1; + for (unsigned int i = 0; i < lastId; i++) { + msg += vpMeanDriftTypeToString(static_cast(i)) + sep; + } + msg += vpMeanDriftTypeToString(static_cast(lastId)) + suffix; + return msg; +} + +void vpStatisticalTestAbstract::print(const vpStatisticalTestAbstract::vpMeanDriftType &type) +{ + std::cout << vpMeanDriftTypeToString(type) << " detected" << std::endl; +} + +bool vpStatisticalTestAbstract::updateStatistics(const float &signal) +{ + m_s[static_cast(m_count)] = signal; + m_count += 1.f; + m_sumForMean += signal; + if (static_cast (m_count) >= m_nbSamplesForStatistics) { + // Computation of the mean + m_mean = m_sumForMean / m_count; + + // Computation of the stdev + float sumSquaredDiff = 0.f; + unsigned int count = static_cast(m_nbSamplesForStatistics); + for (unsigned int i = 0; i < count; ++i) { + sumSquaredDiff += (m_s[i] - m_mean) * (m_s[i] - m_mean); + } + float stdev = std::sqrt(sumSquaredDiff / m_count); + if (m_stdevmin < 0) { + m_stdev = stdev; + } + else { + m_stdev = std::max(m_stdev, m_stdevmin); + } + + m_areStatisticsComputed = true; + } + return m_areStatisticsComputed; +} + +vpStatisticalTestAbstract::vpStatisticalTestAbstract() + : m_areStatisticsComputed(false) + , m_count(0.f) + , m_limitDown(0.f) + , m_limitUp(0.f) + , m_mean(0.f) + , m_nbSamplesForStatistics(0) + , m_s(nullptr) + , m_stdev(0.f) + , m_stdevmin(-1.f) + , m_sumForMean(0.f) +{ } + +vpStatisticalTestAbstract::vpStatisticalTestAbstract(const vpStatisticalTestAbstract &other) +{ + *this = other; +} + +vpStatisticalTestAbstract::~vpStatisticalTestAbstract() +{ + if (m_s != nullptr) { + delete[] m_s; + m_s = nullptr; + } +} + +void vpStatisticalTestAbstract::init() +{ + m_areStatisticsComputed = false; + m_count = 0.f; + m_limitDown = 0.f; + m_limitUp = 0.f; + m_mean = 0.f; + m_nbSamplesForStatistics = 0; + if (m_s != nullptr) { + delete[] m_s; + m_s = nullptr; + } + m_stdev = 0.f; + m_sumForMean = 0.f; +} + +vpStatisticalTestAbstract &vpStatisticalTestAbstract::operator=(const vpStatisticalTestAbstract &other) +{ + m_areStatisticsComputed = other.m_areStatisticsComputed; + m_count = other.m_count; + m_limitDown = other.m_limitDown; + m_limitUp = other.m_limitUp; + m_mean = other.m_mean; + if (other.m_nbSamplesForStatistics > 0) { + setNbSamplesForStat(other.m_nbSamplesForStatistics); + } + else if (m_s != nullptr) { + m_nbSamplesForStatistics = 0; + delete[] m_s; + m_s = nullptr; + } + m_stdev = 0.f; + m_sumForMean = 0.f; + return *this; +} + +void vpStatisticalTestAbstract::setNbSamplesForStat(const unsigned int &nbSamples) +{ + m_nbSamplesForStatistics = nbSamples; + if (m_s != nullptr) { + delete[] m_s; + } + m_s = new float[nbSamples]; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestAbstract::testDownUpwardMeanDrift(const float &signal) +{ + if (m_areStatisticsComputed) { + updateTestSignals(signal); + vpMeanDriftType jumpDown = detectDownwardMeanDrift(); + vpMeanDriftType jumpUp = detectUpwardMeanDrift(); + if ((jumpDown != MEAN_DRIFT_NONE) && (jumpUp != MEAN_DRIFT_NONE)) { + return MEAN_DRIFT_BOTH; + } + else if (jumpDown != MEAN_DRIFT_NONE) { + return jumpDown; + } + else if (jumpUp != MEAN_DRIFT_NONE) { + return jumpUp; + } + else { + return MEAN_DRIFT_NONE; + } + } + else { + updateStatistics(signal); + return MEAN_DRIFT_NONE; + } +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestAbstract::testDownwardMeanDrift(const float &signal) +{ + if (m_areStatisticsComputed) { + updateTestSignals(signal); + return detectDownwardMeanDrift(); + } + else { + updateStatistics(signal); + return MEAN_DRIFT_NONE; + } +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestAbstract::testUpwardMeanDrift(const float &signal) +{ + if (m_areStatisticsComputed) { + updateTestSignals(signal); + return detectUpwardMeanDrift(); + } + else { + updateStatistics(signal); + return MEAN_DRIFT_NONE; + } +} diff --git a/modules/core/src/math/misc/vpStatisticalTestEWMA.cpp b/modules/core/src/math/misc/vpStatisticalTestEWMA.cpp new file mode 100644 index 0000000000..6cb06fee31 --- /dev/null +++ b/modules/core/src/math/misc/vpStatisticalTestEWMA.cpp @@ -0,0 +1,129 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +/** +* +* \file vpStatisticalTestEWMA.cpp +* +* \brief Definition of the vpStatisticalTestEWMA class that implements Exponentially Weighted Moving Average +* mean drift test. +*/ + +#include + +void vpStatisticalTestEWMA::computeDeltaAndLimits() +{ + float delta = 3.f * m_stdev * std::sqrt(m_alpha / (2.f - m_alpha)); + m_limitDown = m_mean - delta; + m_limitUp = m_mean + delta; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestEWMA::detectDownwardMeanDrift() +{ + if (m_wt <= m_limitDown) { + return MEAN_DRIFT_DOWNWARD; + } + else { + return MEAN_DRIFT_NONE; + } +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestEWMA::detectUpwardMeanDrift() +{ + if (m_wt >= m_limitUp) { + return MEAN_DRIFT_UPWARD; + } + else { + return MEAN_DRIFT_NONE; + } +} + +bool vpStatisticalTestEWMA::updateStatistics(const float &signal) +{ + bool areStatsReady = vpStatisticalTestAbstract::updateStatistics(signal); + if (areStatsReady) { + // Computation of the limits + computeDeltaAndLimits(); + + // Initialize first value + m_wt = m_mean; + } + return areStatsReady; +} + +void vpStatisticalTestEWMA::updateTestSignals(const float &signal) +{ + // Update last value + m_wtprev = m_wt; +// w(t) = alpha * s(t) + (1 - alpha) * w(t- 1); + m_wt = m_wtprev + m_alpha * (signal - m_wtprev); +} + +vpStatisticalTestEWMA::vpStatisticalTestEWMA(const float &alpha) + : vpStatisticalTestAbstract() + , m_alpha(0.f) + , m_wt(0.f) + , m_wtprev(0.f) +{ + init(alpha); +} + +void vpStatisticalTestEWMA::init(const float &alpha) +{ + vpStatisticalTestAbstract::init(); + m_alpha = alpha; + unsigned int nbRequiredSamples = static_cast(std::ceil(3.f / m_alpha)); + setNbSamplesForStat(nbRequiredSamples); + m_wt = 0.f; + m_wtprev = 0.f; +} + +void vpStatisticalTestEWMA::init(const float &alpha, const float &mean, const float &stdev) +{ + vpStatisticalTestAbstract::init(); + m_alpha = alpha; + m_mean = mean; + unsigned int nbRequiredSamples = static_cast(std::ceil(3.f / m_alpha)); + setNbSamplesForStat(nbRequiredSamples); + m_stdev = stdev; + m_wt = mean; + m_wtprev = 0.f; + + // Computation of the limits + computeDeltaAndLimits(); + m_areStatisticsComputed = true; +} + +void vpStatisticalTestEWMA::setAlpha(const float &alpha) +{ + init(alpha); +} diff --git a/modules/core/src/math/misc/vpStatisticalTestHinkley.cpp b/modules/core/src/math/misc/vpStatisticalTestHinkley.cpp new file mode 100644 index 0000000000..6d75931952 --- /dev/null +++ b/modules/core/src/math/misc/vpStatisticalTestHinkley.cpp @@ -0,0 +1,251 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +/** +* +* \file vpStatisticalTestHinkley.cpp +* +* \brief Definition of the vpStatisticalTestHinkley class corresponding to the Hinkley's +* cumulative sum test. +*/ + +#include +#include + +#include // std::fabs +#include +#include // numeric_limits +#include +#include + +vpStatisticalTestHinkley::vpStatisticalTestHinkley() + : vpStatisticalTestAbstract() + , m_dmin2(0.1f) + , m_alpha(0.2f) + , m_Sk(0.f) + , m_Mk(0.f) + , m_Tk(0.f) + , m_Nk(0.f) + , m_computeDeltaAndAlpha(false) + , m_h(4.76f) + , m_k(1.f) +{ + init(); +} + +vpStatisticalTestHinkley::vpStatisticalTestHinkley(const float &alpha, const float &delta_val, const unsigned int &nbSamplesForInit) + : vpStatisticalTestAbstract() + , m_dmin2(delta_val / 2.f) + , m_alpha(alpha) + , m_Sk(0.f) + , m_Mk(0.f) + , m_Tk(0.f) + , m_Nk(0.f) + , m_computeDeltaAndAlpha(false) + , m_h(4.76f) + , m_k(1.f) +{ + init(alpha, delta_val, nbSamplesForInit); +} + +vpStatisticalTestHinkley::vpStatisticalTestHinkley(const float &h, const float &k, const bool &computeAlphaDeltaFromStdev, const unsigned int &nbSamplesForInit) + : vpStatisticalTestAbstract() +{ + init(h, k, computeAlphaDeltaFromStdev, nbSamplesForInit); +} + +vpStatisticalTestHinkley::vpStatisticalTestHinkley(const float &h, const float &k, const float &mean, const float &stdev) + : vpStatisticalTestAbstract() +{ + init(h, k, mean, stdev); +} + +void vpStatisticalTestHinkley::init() +{ + vpStatisticalTestAbstract::init(); + setNbSamplesForStat(30); + setAlpha(m_alpha); + + m_Sk = 0.f; + m_Mk = 0.f; + + m_Tk = 0.f; + m_Nk = 0.f; + + m_computeDeltaAndAlpha = false; +} + +void vpStatisticalTestHinkley::init(const float &alpha, const float &delta_val, const unsigned int &nbSamplesForInit) +{ + init(); + setNbSamplesForStat(nbSamplesForInit); + setAlpha(alpha); + setDelta(delta_val); + m_computeDeltaAndAlpha = false; +} + +void vpStatisticalTestHinkley::init(const float &alpha, const float &delta_val, const float &mean) +{ + init(); + setAlpha(alpha); + setDelta(delta_val); + m_mean = mean; + m_computeDeltaAndAlpha = false; + m_areStatisticsComputed = true; +} + +void vpStatisticalTestHinkley::init(const float &h, const float &k, const bool &computeAlphaDeltaFromStdev, const unsigned int &nbSamples) +{ + if (!computeAlphaDeltaFromStdev) { + throw(vpException(vpException::badValue, "computeAlphaDeltaFromStdev must be true, or use another init function")); + } + init(); + setNbSamplesForStat(nbSamples); + m_h = h; + m_k = k; + m_computeDeltaAndAlpha = true; +} + +void vpStatisticalTestHinkley::init(const float &h, const float &k, const float &mean, const float &stdev) +{ + init(); + m_mean = mean; + m_stdev = stdev; + m_h = h; + m_k = k; + m_computeDeltaAndAlpha = true; + computeAlphaDelta(); + m_areStatisticsComputed = true; +} + +vpStatisticalTestHinkley::~vpStatisticalTestHinkley() { } + +void vpStatisticalTestHinkley::setDelta(const float &delta) { m_dmin2 = delta / 2.f; } + +void vpStatisticalTestHinkley::setAlpha(const float &alpha) +{ + this->m_alpha = alpha; + m_limitDown = m_alpha; + m_limitUp = m_alpha; +} + +void vpStatisticalTestHinkley::computeAlphaDelta() +{ + float delta = m_k * m_stdev; + setDelta(delta); + float alpha = m_h * m_stdev; + setAlpha(alpha); +} + +void vpStatisticalTestHinkley::computeMean(double signal) +{ + // When the mean slowly increases or decreases, especially after an abrupt change of the mean, + // the means tends to drift. To reduce the drift of the mean + // it is updated with the current value of the signal only if + // a beginning of a mean drift is not detected, + // i.e. if ( ((m_Mk-m_Sk) == 0) && ((m_Tk-m_Nk) == 0) ) + if ((std::fabs(m_Mk - m_Sk) <= std::fabs(vpMath::maximum(m_Mk, m_Sk)) * std::numeric_limits::epsilon()) && + (std::fabs(m_Tk - m_Nk) <= std::fabs(vpMath::maximum(m_Tk, m_Nk)) * std::numeric_limits::epsilon())) { + m_mean = (m_mean * (m_count - 1) + signal) / (m_count); + } +} + +void vpStatisticalTestHinkley::computeSk(double signal) +{ + m_Sk += signal - m_mean + m_dmin2; +} + +void vpStatisticalTestHinkley::computeMk() +{ + if (m_Sk > m_Mk) { + m_Mk = m_Sk; + } +} + +void vpStatisticalTestHinkley::computeTk(double signal) +{ + m_Tk += signal - m_mean - m_dmin2; +} + +void vpStatisticalTestHinkley::computeNk() +{ + if (m_Tk < m_Nk) { + m_Nk = m_Tk; + } +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestHinkley::detectDownwardMeanDrift() +{ + vpStatisticalTestAbstract::vpMeanDriftType shift = MEAN_DRIFT_NONE; + if ((m_Mk - m_Sk) > m_alpha) { + shift = MEAN_DRIFT_DOWNWARD; + } + return shift; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestHinkley::detectUpwardMeanDrift() +{ + vpStatisticalTestAbstract::vpMeanDriftType shift = MEAN_DRIFT_NONE; + if ((m_Tk - m_Nk) > m_alpha) { + shift = MEAN_DRIFT_UPWARD; + } + return shift; +} + +bool vpStatisticalTestHinkley::updateStatistics(const float &signal) +{ + bool updateStats = vpStatisticalTestAbstract::updateStatistics(signal); + if (m_areStatisticsComputed) { + // If needed, compute alpha and delta + if (m_computeDeltaAndAlpha) { + computeAlphaDelta(); + } + // Initialize the test signals + m_Mk = 0.f; + m_Nk = 0.f; + m_Sk = 0.f; + m_Tk = 0.f; + } + return updateStats; +} + +void vpStatisticalTestHinkley::updateTestSignals(const float &signal) +{ + computeSk(signal); + computeTk(signal); + + computeMk(); + computeNk(); + + ++m_count; + computeMean(signal); +} diff --git a/modules/core/src/math/misc/vpStatisticalTestMeanAdjustedCUSUM.cpp b/modules/core/src/math/misc/vpStatisticalTestMeanAdjustedCUSUM.cpp new file mode 100644 index 0000000000..0007b0a0eb --- /dev/null +++ b/modules/core/src/math/misc/vpStatisticalTestMeanAdjustedCUSUM.cpp @@ -0,0 +1,156 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +/** +* +* \file vpStatisticalTestMeanAdjustedCUSUM.cpp +* +* \brief Definition of the vpStatisticalTestMeanAdjustedCUSUM class that implements mean adjusted Cumulative Sum +* mean drift test. +*/ + +#include + +void vpStatisticalTestMeanAdjustedCUSUM::computeDeltaAndLimits() +{ + setDelta(m_k * m_stdev); + float limitDown = m_h * m_stdev; + float limitUp = limitDown; + setLimits(limitDown, limitUp); +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestMeanAdjustedCUSUM::detectDownwardMeanDrift() +{ + if (m_sminus >= m_limitDown) { + return MEAN_DRIFT_DOWNWARD; + } + else { + return MEAN_DRIFT_NONE; + } +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestMeanAdjustedCUSUM::detectUpwardMeanDrift() +{ + if (m_splus >= m_limitUp) { + return MEAN_DRIFT_UPWARD; + } + else { + return MEAN_DRIFT_NONE; + } +} + +bool vpStatisticalTestMeanAdjustedCUSUM::updateStatistics(const float &signal) +{ + bool areStatsAvailable = vpStatisticalTestAbstract::updateStatistics(signal); + if (areStatsAvailable) { + // Computation of the limits + if ((m_limitDown < 0.f) && (m_limitUp < 0.f)) { + computeDeltaAndLimits(); + } + + // Initialize first values + m_sminus = 0.f; + m_splus = 0.f; + } + return areStatsAvailable; +} + +void vpStatisticalTestMeanAdjustedCUSUM::updateTestSignals(const float &signal) +{ + m_sminus = std::max(0.f, m_sminus - (signal - m_mean) - m_half_delta); + m_splus = std::max(0.f, m_splus + (signal - m_mean) - m_half_delta); +} + +vpStatisticalTestMeanAdjustedCUSUM::vpStatisticalTestMeanAdjustedCUSUM(const float &h, const float &k, const unsigned int &nbPtsForStats) + : vpStatisticalTestAbstract() + , m_delta(-1.f) + , m_h(h) + , m_half_delta(-1.f) + , m_k(k) + , m_sminus(0.f) + , m_splus(0.f) +{ + init(h, k, nbPtsForStats); +} + +void vpStatisticalTestMeanAdjustedCUSUM::init(const float &h, const float &k, const unsigned int &nbPtsForStats) +{ + vpStatisticalTestAbstract::init(); + setNbSamplesForStat(nbPtsForStats); + m_delta = -1.f; + m_half_delta = -1.f; + setLimits(-1.f, -1.f); // To compute automatically the limits once the signal statistics are available. + m_h = h; + m_k = k; + m_mean = 0.f; + m_sminus = 0.f; + m_splus = 0.f; + m_stdev = 0.f; +} + +void vpStatisticalTestMeanAdjustedCUSUM::init(const float &delta, const float &limitDown, const float &limitUp, const unsigned int &nbPtsForStats) +{ + vpStatisticalTestAbstract::init(); + setDelta(delta); + setLimits(limitDown, limitUp); + setNbSamplesForStat(nbPtsForStats); + m_sminus = 0.f; + m_splus = 0.f; +} + +void vpStatisticalTestMeanAdjustedCUSUM::init(const float &h, const float &k, const float &mean, const float &stdev) +{ + vpStatisticalTestAbstract::init(); + setLimits(-1.f, -1.f); // To compute automatically the limits once the signal statistics are available. + m_h = h; + m_k = k; + m_mean = mean; + m_sminus = 0.f; + m_splus = 0.f; + m_stdev = stdev; + // Compute delta and limits from m_h, m_k and m_stdev + computeDeltaAndLimits(); + m_areStatisticsComputed = true; +} + +void vpStatisticalTestMeanAdjustedCUSUM::init(const float &delta, const float &limitDown, const float &limitUp, const float &mean, const float &stdev) +{ + vpStatisticalTestAbstract::init(); + setDelta(delta); + setLimits(limitDown, limitUp); + m_mean = mean; + m_sminus = 0.f; + m_splus = 0.f; + m_stdev = stdev; + m_sumForMean = 0.f; + m_areStatisticsComputed = true; +} diff --git a/modules/core/src/math/misc/vpStatisticalTestShewhart.cpp b/modules/core/src/math/misc/vpStatisticalTestShewhart.cpp new file mode 100644 index 0000000000..b30bcf0a1f --- /dev/null +++ b/modules/core/src/math/misc/vpStatisticalTestShewhart.cpp @@ -0,0 +1,308 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +/** +* +* \file vpStatisticalTestShewhart.cpp +* +* \brief Definition of the vpStatisticalTestShewhart class that implements Shewhart's +* mean drift test. +*/ + +#include + +#include + +#include + +const bool vpStatisticalTestShewhart::CONST_ALL_WECO_ACTIVATED[vpStatisticalTestShewhart::COUNT_WECO - 1] = { true, true, true, true }; + +std::string vpStatisticalTestShewhart::vpWecoRulesAlarmToString(const vpStatisticalTestShewhart::vpWecoRulesAlarm &alarm) +{ + std::string name; + switch (alarm) { + case THREE_SIGMA_WECO: + name = "3-sigma alarm"; + break; + case TWO_SIGMA_WECO: + name = "2-sigma alarm"; + break; + case ONE_SIGMA_WECO: + name = "1-sigma alarm"; + break; + case SAME_SIDE_WECO: + name = "Same-side alarm"; + break; + case NONE_WECO: + name = "No alarm"; + break; + default: + name = "Unknown WECO alarm"; + } + return name; +} + +void vpStatisticalTestShewhart::computeLimits() +{ + float delta = 3.f * m_stdev; + m_limitDown = m_mean - delta; + m_limitUp = m_mean + delta; + m_oneSigmaNegLim = m_mean - m_stdev; + m_oneSigmaPosLim = m_mean + m_stdev; + m_twoSigmaNegLim = m_mean - 2.f * m_stdev; + m_twoSigmaPosLim = m_mean + 2.f * m_stdev; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestShewhart::detectDownwardMeanDrift() +{ + if (m_nbDataInBuffer < NB_DATA_SIGNAL) { + return vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + } + if ((m_signal[m_idCurrentData] <= m_limitDown) && m_activatedWECOrules[THREE_SIGMA_WECO]) { + m_alarm = THREE_SIGMA_WECO; + return vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD; + } + if (!m_activateWECOrules) { + return vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + } + vpStatisticalTestAbstract::vpMeanDriftType result = vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + int id = vpMath::modulo(m_idCurrentData - (NB_DATA_SIGNAL - 1), NB_DATA_SIGNAL); + int i = 0; + unsigned int nbAboveMean = 0; + unsigned int nbAbove2SigmaLimit = 0; + unsigned int nbAbove1SigmaLimit = 0; + while (i < NB_DATA_SIGNAL) { + // Reinit for next iteration + nbAbove2SigmaLimit = 0; + nbAbove1SigmaLimit = 0; + if (m_signal[id] < m_mean && m_activatedWECOrules[SAME_SIDE_WECO]) { + // Single-side test + ++nbAboveMean; + } + if (i > 3 && m_activatedWECOrules[TWO_SIGMA_WECO]) { + // Two sigma test + for (int idPrev = vpMath::modulo(id - 2, NB_DATA_SIGNAL); idPrev != id; idPrev = vpMath::modulo(idPrev + 1, NB_DATA_SIGNAL)) { + if (m_signal[idPrev] <= m_twoSigmaNegLim) { + ++nbAbove2SigmaLimit; + } + } + if (m_signal[id] <= m_twoSigmaNegLim) { + ++nbAbove2SigmaLimit; + } + if (nbAbove2SigmaLimit >= 2) { + break; + } + } + if (i > 5 && m_activatedWECOrules[ONE_SIGMA_WECO]) { + // One sigma test + for (int idPrev = vpMath::modulo(id - 4, NB_DATA_SIGNAL); idPrev != id; idPrev = vpMath::modulo(idPrev + 1, NB_DATA_SIGNAL)) { + if (m_signal[idPrev] <= m_oneSigmaNegLim) { + ++nbAbove1SigmaLimit; + } + } + if (m_signal[id] <= m_oneSigmaNegLim) { + ++nbAbove1SigmaLimit; + } + if (nbAbove1SigmaLimit >= 4) { + break; + } + } + id = vpMath::modulo(id + 1, NB_DATA_SIGNAL); + ++i; + } + if (nbAboveMean == NB_DATA_SIGNAL) { + m_alarm = SAME_SIDE_WECO; + result = MEAN_DRIFT_DOWNWARD; + } + else if (nbAbove2SigmaLimit >= 2) { + m_alarm = TWO_SIGMA_WECO; + result = MEAN_DRIFT_DOWNWARD; + } + else if (nbAbove1SigmaLimit >= 4) { + m_alarm = ONE_SIGMA_WECO; + result = MEAN_DRIFT_DOWNWARD; + } + return result; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestShewhart::detectUpwardMeanDrift() +{ + if (m_nbDataInBuffer < NB_DATA_SIGNAL) { + return vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + } + if ((m_signal[m_idCurrentData] >= m_limitUp) && m_activatedWECOrules[THREE_SIGMA_WECO]) { + m_alarm = THREE_SIGMA_WECO; + return vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD; + } + if (!m_activateWECOrules) { + return vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + } + vpStatisticalTestAbstract::vpMeanDriftType result = vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + int id = vpMath::modulo(m_idCurrentData - (NB_DATA_SIGNAL - 1), NB_DATA_SIGNAL); + int i = 0; + unsigned int nbAboveMean = 0; + unsigned int nbAbove2SigmaLimit = 0; + unsigned int nbAbove1SigmaLimit = 0; + while (i < NB_DATA_SIGNAL) { + // Reinit for next iteration + nbAbove2SigmaLimit = 0; + nbAbove1SigmaLimit = 0; + if (m_signal[id] > m_mean && m_activatedWECOrules[SAME_SIDE_WECO]) { + // Single-side test + ++nbAboveMean; + } + if (i > 3 && m_activatedWECOrules[TWO_SIGMA_WECO]) { + // Two sigma test + for (int idPrev = vpMath::modulo(id - 2, NB_DATA_SIGNAL); idPrev != id; idPrev = vpMath::modulo(idPrev + 1, NB_DATA_SIGNAL)) { + if (m_signal[idPrev] >= m_twoSigmaPosLim) { + ++nbAbove2SigmaLimit; + } + } + if (m_signal[id] >= m_twoSigmaPosLim) { + ++nbAbove2SigmaLimit; + } + if (nbAbove2SigmaLimit >= 2) { + break; + } + } + if (i > 5 && m_activatedWECOrules[ONE_SIGMA_WECO]) { + // One sigma test + for (int idPrev = vpMath::modulo(id - 4, NB_DATA_SIGNAL); idPrev != id; idPrev = vpMath::modulo(idPrev + 1, NB_DATA_SIGNAL)) { + if (m_signal[idPrev] >= m_oneSigmaPosLim) { + ++nbAbove1SigmaLimit; + } + } + if (m_signal[id] >= m_oneSigmaPosLim) { + ++nbAbove1SigmaLimit; + } + if (nbAbove1SigmaLimit >= 4) { + break; + } + } + id = vpMath::modulo(id + 1, NB_DATA_SIGNAL); + ++i; + } + if (nbAboveMean == NB_DATA_SIGNAL) { + m_alarm = SAME_SIDE_WECO; + result = MEAN_DRIFT_UPWARD; + } + else if (nbAbove2SigmaLimit >= 2) { + m_alarm = TWO_SIGMA_WECO; + result = MEAN_DRIFT_UPWARD; + } + else if (nbAbove1SigmaLimit >= 4) { + m_alarm = ONE_SIGMA_WECO; + result = MEAN_DRIFT_UPWARD; + } + return result; +} + +bool vpStatisticalTestShewhart::updateStatistics(const float &signal) +{ + bool areStatsAvailable = vpStatisticalTestAbstract::updateStatistics(signal); + updateTestSignals(signal); // Store the signal in the circular buffer too. + if (areStatsAvailable) { + computeLimits(); + } + return areStatsAvailable; +} + +void vpStatisticalTestShewhart::updateTestSignals(const float &signal) +{ + m_idCurrentData = (m_idCurrentData + 1) % NB_DATA_SIGNAL; + m_signal[m_idCurrentData] = signal; + if (m_nbDataInBuffer < NB_DATA_SIGNAL) { + ++m_nbDataInBuffer; + } +} + +vpStatisticalTestShewhart::vpStatisticalTestShewhart(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1], const unsigned int &nbSamplesForStats) + : vpStatisticalTestSigma(3, nbSamplesForStats) + , m_nbDataInBuffer(0) + , m_activateWECOrules(activateWECOrules) + , m_idCurrentData(0) + , m_alarm(NONE_WECO) + , m_oneSigmaNegLim(0.f) + , m_oneSigmaPosLim(0.f) + , m_twoSigmaNegLim(0.f) + , m_twoSigmaPosLim(0.f) +{ + init(activateWECOrules, activatedRules, nbSamplesForStats); +} + +vpStatisticalTestShewhart::vpStatisticalTestShewhart(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1], const float &mean, const float &stdev) + : vpStatisticalTestSigma(3) + , m_nbDataInBuffer(0) + , m_activateWECOrules(activateWECOrules) + , m_idCurrentData(0) + , m_alarm(NONE_WECO) + , m_oneSigmaNegLim(0.f) + , m_oneSigmaPosLim(0.f) + , m_twoSigmaNegLim(0.f) + , m_twoSigmaPosLim(0.f) +{ + init(activateWECOrules, activatedRules, mean, stdev); +} + +std::vector vpStatisticalTestShewhart::getSignals() const +{ + std::vector signals; + for (int i = 0; i < NB_DATA_SIGNAL; ++i) { + int id = vpMath::modulo(m_idCurrentData - (NB_DATA_SIGNAL - i - 1), NB_DATA_SIGNAL); + signals.push_back(m_signal[id]); + } + return signals; +} + +void vpStatisticalTestShewhart::init(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1], const unsigned int &nbSamplesForStats) +{ + vpStatisticalTestSigma::init(3.f, nbSamplesForStats); + m_nbDataInBuffer = 0; + memset(m_signal, 0, NB_DATA_SIGNAL * sizeof(float)); + m_activateWECOrules = activateWECOrules; + std::memcpy(m_activatedWECOrules, activatedRules, (COUNT_WECO - 1) * sizeof(bool)); + m_idCurrentData = 0; + m_alarm = NONE_WECO; + m_oneSigmaNegLim = 0.f; + m_oneSigmaPosLim = 0.f; + m_twoSigmaNegLim = 0.f; + m_twoSigmaPosLim = 0.f; +} + +void vpStatisticalTestShewhart::init(const bool &activateWECOrules, const bool activatedRules[COUNT_WECO - 1], const float &mean, const float &stdev) +{ + vpStatisticalTestShewhart::init(activateWECOrules, activatedRules, 30); + m_mean = mean; + m_stdev = stdev; + computeLimits(); + m_areStatisticsComputed = true; +} diff --git a/modules/core/src/math/misc/vpStatisticalTestSigma.cpp b/modules/core/src/math/misc/vpStatisticalTestSigma.cpp new file mode 100644 index 0000000000..5a64cacbbc --- /dev/null +++ b/modules/core/src/math/misc/vpStatisticalTestSigma.cpp @@ -0,0 +1,115 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +/** +* +* \file vpStatisticalTestSigma.cpp +* +* \brief Definition of the vpStatisticalTestSigma class that implements sigma +* mean drift test. +*/ + +#include + +void vpStatisticalTestSigma::computeLimits() +{ + float delta = m_h * m_stdev; + m_limitDown = m_mean - delta; + m_limitUp = m_mean + delta; +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestSigma::detectDownwardMeanDrift() +{ + if (m_s <= m_limitDown) { + return vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD; + } + else { + return vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + } +} + +vpStatisticalTestAbstract::vpMeanDriftType vpStatisticalTestSigma::detectUpwardMeanDrift() +{ + if (m_s >= m_limitUp) { + return vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD; + } + else { + return vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + } +} + +bool vpStatisticalTestSigma::updateStatistics(const float &signal) +{ + bool areStatsAvailable = vpStatisticalTestAbstract::updateStatistics(signal); + if (areStatsAvailable) { + computeLimits(); + } + return areStatsAvailable; +} + +void vpStatisticalTestSigma::updateTestSignals(const float &signal) +{ + m_s = signal; +} + +vpStatisticalTestSigma::vpStatisticalTestSigma(const float &h, const unsigned int &nbSamplesForStats) + : vpStatisticalTestAbstract() + , m_h(h) +{ + init(h, nbSamplesForStats); +} + +vpStatisticalTestSigma::vpStatisticalTestSigma(const float &h, const float &mean, const float &stdev) + : vpStatisticalTestAbstract() + , m_h(h) +{ + init(h, mean, stdev); +} + +void vpStatisticalTestSigma::init(const float &h, const unsigned int &nbSamplesForStats) +{ + vpStatisticalTestAbstract::init(); + m_h = h; + setNbSamplesForStat(nbSamplesForStats); + m_s = 0; +} + +void vpStatisticalTestSigma::init(const float &h, const float &mean, const float &stdev) +{ + vpStatisticalTestAbstract::init(); + m_h = h; + m_mean = mean; + m_s = 0; + m_stdev = stdev; + computeLimits(); + m_areStatisticsComputed = true; +} diff --git a/modules/core/test/math/testSPC.cpp b/modules/core/test/math/testSPC.cpp new file mode 100644 index 0000000000..407330060b --- /dev/null +++ b/modules/core/test/math/testSPC.cpp @@ -0,0 +1,809 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * Description: + * Test various Statistical Process Control methods. + */ + +/*! + \example testSPC.cpp + \brief Test various Statistical Process Control methods. +*/ + +#include +#include + +#include +#include +#include +#include +#include + +bool initializeShewhartTest(const float &mean, const float &stdev, const bool &verbose, const std::string &testName, vpStatisticalTestShewhart &tester) +{ + const bool activateWECOrules = true; + tester.init(activateWECOrules, vpStatisticalTestShewhart::CONST_ALL_WECO_ACTIVATED, mean, stdev); + bool isInitOK = true; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT_INIT = vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + vpStatisticalTestAbstract::vpMeanDriftType drift; + unsigned int i = 0; + while ((i < vpStatisticalTestShewhart::NB_DATA_SIGNAL - 1) && isInitOK) { + drift = tester.testDownUpwardMeanDrift(mean); + isInitOK = (drift == EXPECTED_DRIFT_INIT); + if (isInitOK) { + ++i; + } + } + if (!isInitOK) { + if (verbose) { + std::cerr << "\t" << testName << " test initialization failed: " << std::endl; + std::cerr << "\t\ts(t) = " << tester.getSignal() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_- = " << limitDown << std::endl; + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT_INIT) << std::endl; + } + } + return isInitOK; +} + +void usage(const char *name) +{ + std::cout << "SYNOPSIS " << std::endl; + std::cout << name << "[--mean ] [--stdev ] [-v, --verbose] [-h, --help]" << std::endl; + std::cout << "\nOPTIONS" << std::endl; + std::cout << " --mean" << std::endl; + std::cout << " Permits to set the mean of the input signal." << std::endl; + std::cout << std::endl; + std::cout << " --stdev" << std::endl; + std::cout << " Permits to set the standard deviation of the input signal." << std::endl; + std::cout << std::endl; + std::cout << " -v, --verbose" << std::endl; + std::cout << " Activate verbose mode." << std::endl; + std::cout << std::endl; + std::cout << " -h, --help" << std::endl; + std::cout << " Display the help." << std::endl; + std::cout << std::endl; +} + +bool getOptions(int argc, const char **argv, float &opt_mean, float &opt_stdev, bool &opt_verbose) +{ + int i = 1; + while (i < argc) { + std::string argname(argv[i]); + if ((argname == "--mean") && ((i + 1) < argc)) { + opt_mean = std::atof(argv[i + 1]); + ++i; + } + else if ((argname == "--stdev") && ((i + 1) < argc)) { + opt_stdev = std::atof(argv[i + 1]); + ++i; + } + else if ((argname == "-v") || (argname == "--verbose")) { + opt_verbose = true; + } + else if ((argname == "-h") || (argname == "--help")) { + usage(argv[0]); + return false; + } + else if ((argname == "-c") || (argname == "-d")) { + // Arguments given by CTest by default, do nothing + } + else { + usage(argv[0]); + std::cerr << "Error: unrecognised argument \"" << argv[i] << "\"" << std::endl; + return false; + } + ++i; + } + return true; +} + +int main(int argc, const char **argv) +{ + float opt_mean = 0.f; + float opt_stdev = 1.f; + bool opt_verbose = false; + + bool isParsingOk = getOptions(argc, argv, opt_mean, opt_stdev, opt_verbose); + if (!isParsingOk) { + return EXIT_FAILURE; + } + + bool success = true; + + // vpStatisticalTestEWMA tests + { + if (opt_verbose) { + std::cout << "------ vpStatisticalTestEWMA tests ------" << std::endl; + } + const float alpha = 0.1f; + vpStatisticalTestEWMA tester(alpha); + + // ---- Upward drift test ---- + { + tester.init(alpha, opt_mean, opt_stdev); + + // w(t = 1) >= mu + 3 sigma sqrt(alpha / (2 - alpha)) + // <=> alpha s(t=1) + (1 - alpha) mu >= mu + 3 sigma sqrt(alpha / (2 - alpha)) + // <=> s(t=1) >= (1 / alpha) (alpha mu + 3 sigma sqrt(alpha / (2 - alpha)) + float signal = (1.f / alpha) * (alpha * opt_mean + 3.f * opt_stdev * std::sqrt(alpha / (2.f - alpha))); + signal += 0.5f; // To be sure we are greater than the threshold + vpStatisticalTestAbstract::vpMeanDriftType drift = tester.testDownUpwardMeanDrift(signal); + if (drift != vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD) { + success = false; + if (opt_verbose) { + std::cerr << "Upward drift test failed: " << std::endl; + std::cerr << "\tw(t) = " << tester.getWt() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\tlimit_up = " << limitUp << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Upward drift test succeeded." << std::endl; + } + } + + // ---- Downward drift test ---- + { + tester.init(alpha, opt_mean, opt_stdev); + // w(t = 1) <= mu - 3 sigma sqrt(alpha / (2 - alpha)) + // <=> alpha s(t=1) + (1 - alpha) mu <= mu - 3 sigma sqrt(alpha / (2 - alpha)) + // <=> s(t=1) <= (1 / alpha) (alpha mu - 3 sigma sqrt(alpha / (2 - alpha)) + float signal = (1.f / alpha) * (alpha * opt_mean - 3.f * opt_stdev * std::sqrt(alpha / (2.f - alpha))); + signal -= 0.5f; // To be sure we are greater than the threshold + vpStatisticalTestAbstract::vpMeanDriftType drift = tester.testDownUpwardMeanDrift(signal); + if (drift != vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD) { + success = false; + if (opt_verbose) { + std::cerr << "Downward drift test failed: " << std::endl; + std::cerr << "\tw(t) = " << tester.getWt() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\tlimit_down = " << limitDown << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Downward drift test succeeded." << std::endl; + } + } + } + + // vpStatisticalTestHinkley tests + { + if (opt_verbose) { + std::cout << "------ vpStatisticalTestHinkley tests ------" << std::endl; + } + + const float h = 4.76f, k = 1.f; + vpStatisticalTestHinkley tester(h, k, opt_mean, opt_stdev); + const unsigned int HINKLEY_NB_DATA = 4; + const float HINKLEY_SAMPLE = 2.f * opt_stdev; + + // Hinkley's test upward drift + { + const float HINKLEY_DATA[HINKLEY_NB_DATA] = { HINKLEY_SAMPLE, HINKLEY_SAMPLE, HINKLEY_SAMPLE, HINKLEY_SAMPLE }; + const vpStatisticalTestAbstract::vpMeanDriftType HINKLEY_EXPECTED_RES[HINKLEY_NB_DATA] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD }; + bool isTestOk = true; + unsigned int id = 0; + vpStatisticalTestAbstract::vpMeanDriftType drift; + while ((id < HINKLEY_NB_DATA) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(HINKLEY_DATA[id]); + isTestOk = (drift == HINKLEY_EXPECTED_RES[id]); + if (isTestOk) { + ++id; + } + } + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "Upward drift test failed: " << std::endl; + float Tk = tester.getTk(), Nk = tester.getNk(); + std::cerr << "T(k) = " << Tk << " | N(k) = " << Nk << " => S+(k) = " << Tk - Nk << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "lim_+ = " << limitUp << std::endl; + std::cerr << "drift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "expected drift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(HINKLEY_EXPECTED_RES[id]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Upward drift test succeeded." << std::endl; + } + } + + // Hinkley's test downward drift + { + const float HINKLEY_DATA[HINKLEY_NB_DATA] = { -1.f * HINKLEY_SAMPLE, -1.f * HINKLEY_SAMPLE, -1.f * HINKLEY_SAMPLE, -1.f * HINKLEY_SAMPLE }; + const vpStatisticalTestAbstract::vpMeanDriftType HINKLEY_EXPECTED_RES[HINKLEY_NB_DATA] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD }; + bool isTestOk = true; + unsigned int id = 0; + vpStatisticalTestAbstract::vpMeanDriftType drift; + while ((id < HINKLEY_NB_DATA) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(HINKLEY_DATA[id]); + isTestOk = (drift == HINKLEY_EXPECTED_RES[id]); + if (isTestOk) { + ++id; + } + } + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "Downward drift test failed: " << std::endl; + float Sk = tester.getSk(), Mk = tester.getMk(); + std::cerr << "S(k) = " << Sk << " | M(k) = " << Mk << " => S+(k) = " << Mk - Sk << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "lim_- = " << limitDown << std::endl; + std::cerr << "drift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "expected drift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(HINKLEY_EXPECTED_RES[id]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Downward drift test succeeded." << std::endl; + } + } + } + + // vpStatisticalTestMeanAdjustedCUSUM tests + { + if (opt_verbose) { + std::cout << "------ vpStatisticalTestMeanAdjustedCUSUM tests ------" << std::endl; + } + + const float h = 4.76f, k = 1.f; + vpStatisticalTestMeanAdjustedCUSUM tester(h, k); + tester.init(h, k, opt_mean, opt_stdev); + const unsigned int CUSUM_NB_DATA = 4; + const float CUSUM_SAMPLE = 2.f * opt_stdev; + + // Mean adjusted CUSUM test upward drift + { + const float CUSUM_DATA[CUSUM_NB_DATA] = { CUSUM_SAMPLE, CUSUM_SAMPLE, CUSUM_SAMPLE, CUSUM_SAMPLE }; + const vpStatisticalTestAbstract::vpMeanDriftType CUSUM_EXPECTED_RES[CUSUM_NB_DATA] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD }; + bool isTestOk = true; + unsigned int id = 0; + vpStatisticalTestAbstract::vpMeanDriftType drift; + while ((id < CUSUM_NB_DATA) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(CUSUM_DATA[id]); + isTestOk = (drift == CUSUM_EXPECTED_RES[id]); + if (isTestOk) { + ++id; + } + } + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "Upward drift test failed: " << std::endl; + std::cerr << "\tS+(k) = " << tester.getTestSignalPlus() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\tlim_+ = " << limitUp << std::endl; + std::cerr << "\tdrift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\texpected drift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(CUSUM_EXPECTED_RES[id]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Upward drift test succeeded." << std::endl; + } + } + + // Mean adjusted CUSUM test upward drift + { + const float CUSUM_DATA[CUSUM_NB_DATA] = { -1.f * CUSUM_SAMPLE, -1.f * CUSUM_SAMPLE, -1.f * CUSUM_SAMPLE, -1.f * CUSUM_SAMPLE }; + const vpStatisticalTestAbstract::vpMeanDriftType CUSUM_EXPECTED_RES[CUSUM_NB_DATA] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD }; + bool isTestOk = true; + unsigned int id = 0; + vpStatisticalTestAbstract::vpMeanDriftType drift; + while ((id < CUSUM_NB_DATA) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(CUSUM_DATA[id]); + isTestOk = (drift == CUSUM_EXPECTED_RES[id]); + if (isTestOk) { + ++id; + } + } + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "Downward drift test failed: " << std::endl; + std::cerr << "\tS-(k) = " << tester.getTestSignalMinus() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\tlim_- = " << limitDown << std::endl; + std::cerr << "\tdrift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\texpected drift type : " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(CUSUM_EXPECTED_RES[id]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Downward drift test succeeded." << std::endl; + } + } + } + + // vpStatisticalTestShewhart tests + { + if (opt_verbose) { + std::cout << "------ vpStatisticalTestShewhart tests ------" << std::endl; + } + const bool activateWECOrules = true; + vpStatisticalTestShewhart tester(activateWECOrules, vpStatisticalTestShewhart::CONST_ALL_WECO_ACTIVATED, opt_mean, opt_stdev); + + // Upward drift test + { + if (opt_verbose) { + std::cout << "Upward drift tests" << std::endl; + } + + // 3-sigma test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "3-sigma", tester); + if (!isInitOK) { + success = false; + } + else { + const float signal = 3.5f * opt_stdev; + vpStatisticalTestAbstract::vpMeanDriftType drift = tester.testDownUpwardMeanDrift(signal); + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT = vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM = vpStatisticalTestShewhart::THREE_SIGMA_WECO; + if ((drift != EXPECTED_DRIFT) || (alarm != EXPECTED_ALARM)) { + success = false; + if (opt_verbose) { + std::cerr << "\t3-sigma test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\t3-sigma test succeeded " << std::endl; + } + } + } + + // 2-sigma test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "2-sigma", tester); + if (!isInitOK) { + success = false; + } + else { + const unsigned int NB_SAMPLES = 3; + const float DATA[NB_SAMPLES] = { 2.75f * opt_stdev, 1.5f * opt_stdev, 2.5f * opt_stdev }; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT[NB_SAMPLES] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD }; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM[NB_SAMPLES] = { vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::TWO_SIGMA_WECO }; + vpStatisticalTestAbstract::vpMeanDriftType drift; + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + unsigned int i = 0; + bool isTestOk = true; + while ((i < NB_SAMPLES) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(DATA[i]); + alarm = tester.getAlarm(); + isTestOk = ((drift == EXPECTED_DRIFT[i]) && (alarm == EXPECTED_ALARM[i])); + if (isTestOk) { + ++i; + } + } + + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "\t2-sigma test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT[i]) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM[i]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\t2-sigma test succeeded " << std::endl; + } + } + } + + // 1-sigma test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "1-sigma", tester); + if (!isInitOK) { + success = false; + } + else { + const unsigned int NB_SAMPLES = 5; + const float DATA[NB_SAMPLES] = { 2.75f * opt_stdev, 1.5f * opt_stdev, 0.5f * opt_stdev, 1.5f * opt_stdev, 2.5f * opt_stdev }; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT[NB_SAMPLES] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD }; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM[NB_SAMPLES] = { vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::ONE_SIGMA_WECO }; + vpStatisticalTestAbstract::vpMeanDriftType drift; + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + unsigned int i = 0; + bool isTestOk = true; + while ((i < NB_SAMPLES) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(DATA[i]); + alarm = tester.getAlarm(); + isTestOk = ((drift == EXPECTED_DRIFT[i]) && (alarm == EXPECTED_ALARM[i])); + if (isTestOk) { + ++i; + } + } + + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "\t1-sigma test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT[i]) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM[i]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\t1-sigma test succeeded " << std::endl; + } + } + } + + // Same-side test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "Same-side", tester); + if (!isInitOK) { + success = false; + } + else { + const unsigned int NB_SAMPLES = 8; + const float DATA[NB_SAMPLES] = { 2.75f * opt_stdev, 0.5f * opt_stdev, 1.5f * opt_stdev, 0.5f * opt_stdev, 2.75f * opt_stdev, 0.5f * opt_stdev, 1.5f * opt_stdev, 0.5f * opt_stdev }; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT[NB_SAMPLES] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD }; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM[NB_SAMPLES] = { vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::SAME_SIDE_WECO }; + vpStatisticalTestAbstract::vpMeanDriftType drift; + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + unsigned int i = 0; + bool isTestOk = true; + while ((i < NB_SAMPLES) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(DATA[i]); + alarm = tester.getAlarm(); + isTestOk = ((drift == EXPECTED_DRIFT[i]) && (alarm == EXPECTED_ALARM[i])); + if (isTestOk) { + ++i; + } + } + + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "\tSame-side test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT[i]) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM[i]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\tSame-side test succeeded " << std::endl; + } + } + } + } + + // Downward drift test + { + if (opt_verbose) { + std::cout << "Downward drift tests" << std::endl; + } + + // 3-sigma test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "3-sigma", tester); + if (!isInitOK) { + success = false; + } + else { + const float signal = -3.5f * opt_stdev; + vpStatisticalTestAbstract::vpMeanDriftType drift = tester.testDownUpwardMeanDrift(signal); + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT = vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM = vpStatisticalTestShewhart::THREE_SIGMA_WECO; + if ((drift != EXPECTED_DRIFT) || (alarm != EXPECTED_ALARM)) { + success = false; + if (opt_verbose) { + std::cerr << "\t3-sigma test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_- = " << limitDown << std::endl; + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\t3-sigma test succeeded " << std::endl; + } + } + } + + // 2-sigma test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "2-sigma", tester); + if (!isInitOK) { + success = false; + } + else { + const unsigned int NB_SAMPLES = 3; + const float DATA[NB_SAMPLES] = { -2.75f * opt_stdev, -1.5f * opt_stdev, -2.5f * opt_stdev }; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT[NB_SAMPLES] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD }; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM[NB_SAMPLES] = { vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::TWO_SIGMA_WECO }; + vpStatisticalTestAbstract::vpMeanDriftType drift; + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + unsigned int i = 0; + bool isTestOk = true; + while ((i < NB_SAMPLES) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(DATA[i]); + alarm = tester.getAlarm(); + isTestOk = ((drift == EXPECTED_DRIFT[i]) && (alarm == EXPECTED_ALARM[i])); + if (isTestOk) { + ++i; + } + } + + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "\t2-sigma test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_- = " << limitDown << std::endl; + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT[i]) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM[i]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\t2-sigma test succeeded " << std::endl; + } + } + } + + // 1-sigma test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "1-sigma", tester); + if (!isInitOK) { + success = false; + } + else { + const unsigned int NB_SAMPLES = 5; + const float DATA[NB_SAMPLES] = { -2.75f * opt_stdev, -1.5f * opt_stdev, 1.5f * opt_stdev, -1.5f * opt_stdev, -2.5f * opt_stdev }; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT[NB_SAMPLES] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD }; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM[NB_SAMPLES] = { vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::ONE_SIGMA_WECO }; + vpStatisticalTestAbstract::vpMeanDriftType drift; + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + unsigned int i = 0; + bool isTestOk = true; + while ((i < NB_SAMPLES) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(DATA[i]); + alarm = tester.getAlarm(); + isTestOk = ((drift == EXPECTED_DRIFT[i]) && (alarm == EXPECTED_ALARM[i])); + if (isTestOk) { + ++i; + } + } + + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "\t1-sigma test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_- = " << limitDown << std::endl; + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT[i]) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM[i]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\t1-sigma test succeeded " << std::endl; + } + } + } + + // Same-side test + { + bool isInitOK = initializeShewhartTest(opt_mean, opt_stdev, opt_verbose, "Same-side", tester); + if (!isInitOK) { + success = false; + } + else { + const unsigned int NB_SAMPLES = 8; + const float DATA[NB_SAMPLES] = { -2.75f * opt_stdev, -0.5f * opt_stdev, -1.5f * opt_stdev, -0.5f * opt_stdev, -2.75f * opt_stdev, -0.5f * opt_stdev, -1.5f * opt_stdev, -0.5f * opt_stdev }; + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT[NB_SAMPLES] = { vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_NONE, vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD }; + const vpStatisticalTestShewhart::vpWecoRulesAlarm EXPECTED_ALARM[NB_SAMPLES] = { vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::NONE_WECO, vpStatisticalTestShewhart::SAME_SIDE_WECO }; + vpStatisticalTestAbstract::vpMeanDriftType drift; + vpStatisticalTestShewhart::vpWecoRulesAlarm alarm = tester.getAlarm(); + unsigned int i = 0; + bool isTestOk = true; + while ((i < NB_SAMPLES) && isTestOk) { + drift = tester.testDownUpwardMeanDrift(DATA[i]); + alarm = tester.getAlarm(); + isTestOk = ((drift == EXPECTED_DRIFT[i]) && (alarm == EXPECTED_ALARM[i])); + if (isTestOk) { + ++i; + } + } + + if (!isTestOk) { + success = false; + if (opt_verbose) { + std::cerr << "\tSame-side test failed: " << std::endl; + std::vector s = tester.getSignals(); + std::cerr << "\t\ts(t) = [ "; + for (size_t i = 0; i < s.size(); ++i) { + std::cerr << s[i] << " "; + } + std::cerr << "]" << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\t\tlim_+ = " << limitUp << std::endl; + std::cerr << "\t\tlim_- = " << limitDown << std::endl; + std::cerr << "\t\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\t\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT[i]) << std::endl; + std::cerr << "\t\tdetected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(alarm) << std::endl; + std::cerr << "\t\texpected alarm = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(EXPECTED_ALARM[i]) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "\tSame-side test succeeded " << std::endl; + } + } + } + } + } + + // vpStatisticalTestSigma tests + { + if (opt_verbose) { + std::cout << "------ vpStatisticalTestSigma tests ------" << std::endl; + } + const float h = 3.f; + vpStatisticalTestSigma tester(h, opt_mean, opt_stdev); + + // Upward drift test + { + const float signal = 3.5f * opt_stdev; + vpStatisticalTestAbstract::vpMeanDriftType drift = tester.testDownUpwardMeanDrift(signal); + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT = vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD; + if (drift != EXPECTED_DRIFT) { + success = false; + if (opt_verbose) { + std::cerr << "Upward drift test failed: " << std::endl; + std::cerr << "\ts(t) = " << tester.getSignal() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\tlim_+ = " << limitUp << std::endl; + std::cerr << "\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Upward drift test succeeded " << std::endl; + } + } + + // Downward drift test + { + const float signal = -3.5f * opt_stdev; + vpStatisticalTestAbstract::vpMeanDriftType drift = tester.testDownUpwardMeanDrift(signal); + const vpStatisticalTestAbstract::vpMeanDriftType EXPECTED_DRIFT = vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD; + if (drift != EXPECTED_DRIFT) { + success = false; + if (opt_verbose) { + std::cerr << "Downward drift test failed: " << std::endl; + std::cerr << "\ts(t) = " << tester.getSignal() << std::endl; + float limitDown = 0.f, limitUp = 0.f; + tester.getLimits(limitDown, limitUp); + std::cerr << "\tlim_- = " << limitDown << std::endl; + std::cerr << "\tdetected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift) << std::endl; + std::cerr << "\texpected drift = " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(EXPECTED_DRIFT) << std::endl; + } + } + else if (opt_verbose) { + std::cout << "Downward drift test succeeded " << std::endl; + } + } + } + + if (success) { + std::cout << "Test succeed" << std::endl; + } + else { + std::cout << "Test failed" << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/modules/imgproc/src/vpCircleHoughTransform.cpp b/modules/imgproc/src/vpCircleHoughTransform.cpp index 3d1c4660d3..70aeb1933b 100644 --- a/modules/imgproc/src/vpCircleHoughTransform.cpp +++ b/modules/imgproc/src/vpCircleHoughTransform.cpp @@ -315,12 +315,12 @@ vpCircleHoughTransform::detect(const vpImage &I) // the pixelization of the image const float minCenterPositionDiff = 3.f; if ((m_algoParams.m_centerXlimits.second - m_algoParams.m_centerXlimits.first) < minCenterPositionDiff) { - m_algoParams.m_centerXlimits.second += minCenterPositionDiff / 2.f; - m_algoParams.m_centerXlimits.first -= minCenterPositionDiff / 2.f; + m_algoParams.m_centerXlimits.second += static_cast(minCenterPositionDiff / 2.f); + m_algoParams.m_centerXlimits.first -= static_cast(minCenterPositionDiff / 2.f); } if ((m_algoParams.m_centerYlimits.second - m_algoParams.m_centerYlimits.first) < minCenterPositionDiff) { - m_algoParams.m_centerYlimits.second += minCenterPositionDiff / 2.f; - m_algoParams.m_centerYlimits.first -= minCenterPositionDiff / 2.f; + m_algoParams.m_centerYlimits.second += static_cast(minCenterPositionDiff / 2.f); + m_algoParams.m_centerYlimits.first -= static_cast(minCenterPositionDiff / 2.f); } // First thing, we need to apply a Gaussian filter on the image to remove some spurious noise @@ -397,7 +397,7 @@ void vpCircleHoughTransform::computeVotingMask(const vpImage &I, #endif { bool hasFoundSimilarCircle = false; - unsigned int nbPreviouslyDetected = m_finalCircles.size(); + unsigned int nbPreviouslyDetected = static_cast(m_finalCircles.size()); unsigned int id = 0; // Looking for a circle that was detected and is similar to the one given to the function while ((id < nbPreviouslyDetected) && (!hasFoundSimilarCircle)) { @@ -410,7 +410,7 @@ void vpCircleHoughTransform::computeVotingMask(const vpImage &I, { hasFoundSimilarCircle = true; // We found a circle that is similar to the one given to the function => updating the mask - const unsigned int nbVotingPoints = m_finalCirclesVotingPoints[id].size(); + const unsigned int nbVotingPoints = static_cast(m_finalCirclesVotingPoints[id].size()); for (unsigned int idPoint = 0; idPoint < nbVotingPoints; ++idPoint) { const std::pair &votingPoint = m_finalCirclesVotingPoints[id][idPoint]; #if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_17) @@ -855,11 +855,11 @@ vpCircleHoughTransform::computeCenterCandidates() std::sort(merged_peaks_position_votes.begin(), merged_peaks_position_votes.end(), sortingCenters); nbPeaks = static_cast(merged_peaks_position_votes.size()); - int nbPeaksToKeep = (m_algoParams.m_expectedNbCenters > 0 ? m_algoParams.m_expectedNbCenters : nbPeaks); - nbPeaksToKeep = std::min(nbPeaksToKeep, (int)nbPeaks); + int nbPeaksToKeep = (m_algoParams.m_expectedNbCenters > 0 ? m_algoParams.m_expectedNbCenters : static_cast(nbPeaks)); + nbPeaksToKeep = std::min(nbPeaksToKeep, static_cast(nbPeaks)); for (int i = 0; i < nbPeaksToKeep; i++) { m_centerCandidatesList.push_back(merged_peaks_position_votes[i].first); - m_centerVotes.push_back(merged_peaks_position_votes[i].second); + m_centerVotes.push_back(static_cast(merged_peaks_position_votes[i].second)); } } } @@ -1046,17 +1046,17 @@ vpCircleHoughTransform::computeCircleCandidates() } } - unsigned int nbCandidates = v_r_effective.size(); + unsigned int nbCandidates = static_cast(v_r_effective.size()); for (unsigned int idBin = 0; idBin < nbCandidates; ++idBin) { // If the circle of center CeC_i and radius RCB_k has enough votes, it is added to the list // of Circle Candidates float r_effective = v_r_effective[idBin]; vpImageCircle candidateCircle(vpImagePoint(centerCandidate.first, centerCandidate.second), r_effective); - float proba = computeCircleProbability(candidateCircle, v_votes_effective[idBin]); + float proba = computeCircleProbability(candidateCircle, static_cast(v_votes_effective[idBin])); if (proba > m_algoParams.m_circleProbaThresh) { m_circleCandidates.push_back(candidateCircle); m_circleCandidatesProbabilities.push_back(proba); - m_circleCandidatesVotes.push_back(v_votes_effective[idBin]); + m_circleCandidatesVotes.push_back(static_cast(v_votes_effective[idBin])); if (m_algoParams.m_recordVotingPoints) { if (v_hasMerged_effective[idBin]) { // Remove potential duplicated points @@ -1078,7 +1078,7 @@ vpCircleHoughTransform::computeCircleProbability(const vpImageCircle &circle, co float visibleArc(static_cast(nbVotes)); float theoreticalLenght; if (mp_mask != nullptr) { - theoreticalLenght = circle.computePixelsInMask(*mp_mask); + theoreticalLenght = static_cast(circle.computePixelsInMask(*mp_mask)); } else { theoreticalLenght = circle.computeArcLengthInRoI(vpRect(vpImagePoint(0, 0), m_edgeMap.getWidth(), m_edgeMap.getHeight())); diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt index 7ede770bdd..ae1233344b 100644 --- a/tutorial/CMakeLists.txt +++ b/tutorial/CMakeLists.txt @@ -44,6 +44,7 @@ visp_add_subdirectory(imgproc/contrast-sharpening REQUIRED_DEPS visp_co visp_add_subdirectory(imgproc/count-coins REQUIRED_DEPS visp_core visp_io visp_gui visp_imgproc) visp_add_subdirectory(imgproc/flood-fill REQUIRED_DEPS visp_core visp_io visp_gui visp_imgproc) visp_add_subdirectory(imgproc/hough-transform REQUIRED_DEPS visp_core visp_gui visp_imgproc) +visp_add_subdirectory(mean-drift REQUIRED_DEPS visp_core visp_gui) visp_add_subdirectory(munkres REQUIRED_DEPS visp_core visp_gui) visp_add_subdirectory(robot/flir-ptu REQUIRED_DEPS visp_core visp_robot visp_sensor visp_vision visp_gui visp_vs visp_visual_features visp_detection) visp_add_subdirectory(robot/pioneer REQUIRED_DEPS visp_core visp_robot visp_vs visp_gui) diff --git a/tutorial/mean-drift/CMakeLists.txt b/tutorial/mean-drift/CMakeLists.txt new file mode 100644 index 0000000000..fed2f435fd --- /dev/null +++ b/tutorial/mean-drift/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 3.5) + +project(tutorial-meandrift) + +find_package(VISP REQUIRED visp_core visp_gui) + +set(tutorial_cpp) + +list(APPEND tutorial_cpp tutorial-meandrift.cpp) + +foreach(cpp ${tutorial_cpp}) + visp_add_target(${cpp}) + if(COMMAND visp_add_dependency) + visp_add_dependency(${cpp} "tutorials") + endif() +endforeach() diff --git a/tutorial/mean-drift/tutorial-meandrift.cpp b/tutorial/mean-drift/tutorial-meandrift.cpp new file mode 100644 index 0000000000..9ba89e457c --- /dev/null +++ b/tutorial/mean-drift/tutorial-meandrift.cpp @@ -0,0 +1,708 @@ +/**************************************************************************** + * + * ViSP, open source Visual Servoing Platform software. + * Copyright (C) 2005 - 2023 by Inria. All rights reserved. + * + * This software is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * See the file LICENSE.txt at the root directory of this source + * distribution for additional information about the GNU GPL. + * + * For using ViSP with software that can not be combined with the GNU + * GPL, please contact Inria about acquiring a ViSP Professional + * Edition License. + * + * See https://visp.inria.fr for more information. + * + * This software was developed at: + * Inria Rennes - Bretagne Atlantique + * Campus Universitaire de Beaulieu + * 35042 Rennes Cedex + * France + * + * If you have questions regarding the use of this file, please contact + * Inria at visp@inria.fr + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * +*****************************************************************************/ + +//! \example tutorial-meandrift.cpp + +#include //std::memcpy + +#include +#include +#include +#include +#include +#include +#include +#include + +#if defined(VISP_HAVE_DISPLAY) +namespace TutorialMeanDrift +{ + //! [Enum_For_Test_Choice] + /** + * \brief Enumeration that permits to choose which process test to use. + * + */ +typedef enum TypeTest +{ + HINLKEY_TYPE_TEST = 0, /*!< Use Hinkley's test.*/ + EWMA_TYPE_TEST = 1, /*!< Use Exponentially Weighted Moving Average to perform the tests.*/ + MEAN_ADJUSTED_CUSUM_TYPE_TEST = 2, /*!< Use mean adjusted Cumulative Sum to perform the tests.*/ + SHEWHART_TYPE_TEST = 3, /*!< Shewhart's test.*/ + SIGMA_TYPE_TEST = 4, /*!< Simple test based on the comparisong with the standard deviation.*/ + COUNT_TYPE_TEST = 5, /*!< Number of different aavailable methods.*/ + UNKOWN_TYPE_TEST = COUNT_TYPE_TEST /*!< Unknown method.*/ +}TypeTest; +//! [Enum_For_Test_Choice] + +/** + * \brief Permit to cast a \b TypeTest object into a string, for display purpose. + * + * \param[in] type The \b TypeTest object we want to know the name. + * \return std::string The corresponding name. + */ +std::string typeTestToString(const TypeTest &type) +{ + std::string result; + switch (type) { + case HINLKEY_TYPE_TEST: + result = "hinkley"; + break; + case EWMA_TYPE_TEST: + result = "ewma"; + break; + case MEAN_ADJUSTED_CUSUM_TYPE_TEST: + result = "cusum"; + break; + case SHEWHART_TYPE_TEST: + result = "shewhart"; + break; + case SIGMA_TYPE_TEST: + result = "sigma"; + break; + case UNKOWN_TYPE_TEST: + default: + result = "unknown-type-test"; + break; + } + return result; +} + +/** + * \brief Permit to cast a string into a \b TypeTest, to + * cast a command line argument. + * + * \param[in] name The name of the process test the user wants to use. + * \return TypeTest The corresponding \b TypeTest object. + */ +TypeTest typeTestFromString(const std::string &name) +{ + TypeTest result = UNKOWN_TYPE_TEST; + unsigned int count = static_cast(COUNT_TYPE_TEST); + unsigned int id = 0; + bool hasNotFound = true; + while ((id < count) && hasNotFound) { + TypeTest temp = static_cast(id); + if (typeTestToString(temp) == name) { + result = temp; + hasNotFound = false; + } + ++id; + } + return result; +} + +/** + * \brief Get the list of available \b TypeTest objects that are handled. + * + * \param[in] prefix The prefix that should be placed before the list. + * \param[in] sep The separator between each element of the list. + * \param[in] suffix The suffix that should terminate the list. + * \return std::string The list of handled type of process tests, presented as a string. + */ +std::string getAvailableTypeTest(const std::string &prefix = "<", const std::string &sep = " , ", + const std::string &suffix = ">") +{ + std::string msg(prefix); + unsigned int count = static_cast(COUNT_TYPE_TEST); + unsigned int lastId = count - 1; + for (unsigned int i = 0; i < lastId; i++) { + msg += typeTestToString(static_cast(i)) + sep; + } + msg += typeTestToString(static_cast(lastId)) + suffix; + return msg; +} + +/** + * \brief Cast a number type into a string. + * + * \tparam T Type of number. + * \param[in] number The number to cast. + * \return std::string The corresponding string. + */ +template +std::string numberToString(const T &number) +{ + std::stringstream ss; + ss << number; + return ss.str(); +} + +/** + * \brief Cast a boolean into a string. + * + * \param[in] boolean The boolean to cast. + * \return std::string The corresponding string. + */ +std::string boolToString(const bool &boolean) +{ + if (boolean) { + return "true"; + } + else { + return "false"; + } +} + +/** + * \brief Write the WECO's rules used in the Shewhart's test in human readable format. + * + * \param[in] rules The array indicating which WECO's rules are used. + * \param[in] prefix The first character(s) delimiting the array in the string. + * \param[in] suffix The last character(s) delimiting the array in the string. + * \param[in] sep The separator character(s). + * \return std::string The corresponding string. + */ +std::string wecoRulesToString(const bool rules[vpStatisticalTestShewhart::COUNT_WECO - 1], const std::string &prefix = "[", const std::string &suffix = "]", const std::string &sep = " , ") +{ + std::string rulesAsString = prefix; + for (unsigned int i = 0; i < vpStatisticalTestShewhart::COUNT_WECO - 2; ++i) { + if (rules[i]) { + rulesAsString += "ON"; + } + else { + rulesAsString += "OFF"; + } + rulesAsString += sep; + } + if (rules[vpStatisticalTestShewhart::COUNT_WECO - 2]) { + rulesAsString += "ON"; + } + else { + rulesAsString += "OFF"; + } + rulesAsString += suffix; + return rulesAsString; +} + +/** + * \brief Array that sets all the types of mean drift to deactivated. + */ +const bool CONST_ALL_ALARM_OFF[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT] = { false, false, false, false }; + +/** + * \brief Array that sets all the types of mean drift to activated. + */ +const bool CONST_ALL_ALARM_ON[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT] = { true, true, true, true }; + +/** + * \brief Cast a vector of string into an array of boolean activating / deactivating + * the mean drift alarms. + * + * \param[in] names The names of the alarms to set. + * \param[out] array The corresponding array of boolean. + */ +void vectorOfStringToMeanDriftTypeArray(const std::vector &names, bool array[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT]) +{ + std::memcpy(array, CONST_ALL_ALARM_OFF, vpStatisticalTestAbstract::MEAN_DRIFT_COUNT * sizeof(bool)); + size_t nbNames = names.size(); + for (size_t i = 0; i < nbNames; ++i) { + vpStatisticalTestAbstract::vpMeanDriftType alarmToActivate = vpStatisticalTestAbstract::vpMeanDriftTypeFromString(names[i]); + std::cout << "alarm[" << names[i] << "] (i.e. " << static_cast(alarmToActivate) << ") set to true" << std::endl; + array[static_cast(alarmToActivate)] = true; + if (alarmToActivate == vpStatisticalTestAbstract::MEAN_DRIFT_BOTH) { + array[vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD] = true; + array[vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD] = true; + } + } + if (array[vpStatisticalTestAbstract::MEAN_DRIFT_DOWNWARD] || array[vpStatisticalTestAbstract::MEAN_DRIFT_UPWARD]) { + array[vpStatisticalTestAbstract::MEAN_DRIFT_BOTH] = true; + } +} + +/** + * \brief Cast an array of boolean (de)activating the mean drift alarms into + * the corresponding vector of strings. + * + * \param[in] array The array of boolean indicating which alarm are set. + * \return std::vector The corresponding vector of names of alarms. + */ +std::vector meanDriftArrayToVectorOfString(const bool array[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT]) +{ + std::vector listActivatedAlarms; + unsigned int nbTypeTests = static_cast(vpStatisticalTestAbstract::MEAN_DRIFT_COUNT); + for (unsigned int id = 0; id < nbTypeTests; ++id) { + if (array[id]) { + vpStatisticalTestAbstract::vpMeanDriftType test = static_cast(id); + std::string testName = vpStatisticalTestAbstract::vpMeanDriftTypeToString(test); + listActivatedAlarms.push_back(testName); + } + } + return listActivatedAlarms; +} + +/** + * \brief Cast an array of boolean (de)activating the mean drift alarms into + * a single string listing all the alarms. + * + * \param[in] array The array of boolean indicating which alarm are set. + * \param[in] prefix The returned string prefix. + * \param[in] sep The returned string separator. + * \param[in] suffix The returned string suffix. + * \return std::string The corresponding string listing the names of alarms. + */ +std::string meanDriftArrayToString(const bool array[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT], + const std::string &prefix = "[", const std::string &sep = " , ", + const std::string &suffix = "]") +{ + std::vector listActivatedAlarms = meanDriftArrayToVectorOfString(array); + std::string result = prefix; + size_t nbTestActivated = listActivatedAlarms.size(); + if (nbTestActivated == 0) { + return prefix + " " + suffix; + } + for (size_t i = 0; i < nbTestActivated - 1; ++i) { + result += listActivatedAlarms[i] + sep; + } + result += listActivatedAlarms[nbTestActivated - 1] + suffix; + return result; +} + +/** + * \brief Indicate how many alarms are set. + * + * \param[in] array The array of boolean indicating which alarms are set. + * \return unsigned int The number of alarms that are set. + */ +unsigned int meanDriftArrayToNbActivated(const bool array[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT]) +{ + unsigned int nbActivated = 0; + unsigned int nbTypeAlarms = static_cast(vpStatisticalTestAbstract::MEAN_DRIFT_COUNT); + for (unsigned int id = 0; id < nbTypeAlarms; ++id) { + if (array[id]) { + ++nbActivated; + } + } + return nbActivated; +} + +//! [Structure_Parameters] +/** + * \brief Structure that contains the parameters of the different algorithms. + */ +typedef struct ParametersForAlgo +{ + unsigned int m_test_nbsamples; /*!< Number of samples to compute the mean and stdev, common to all the algorithms.*/ + bool m_test_activatedalarms[vpStatisticalTestAbstract::MEAN_DRIFT_COUNT]; /*!< Flag is true for a type of alarm that must be considered, false otherwise.*/ + unsigned int m_test_nbactivatedalarms; /*!< Number of activated alarms.*/ + float m_cusum_h; /*!< Alarm factor for the mean adjusted CUSUM test.*/ + float m_cusum_k; /*!< Detection factor for the mean adjusted CUSUM test.*/ + float m_ewma_alpha; /*!< Forgetting factor for the EWMA test.*/ + float m_hinkley_alpha; /*!< Alarm threshold for the Hinkley's test. */ + float m_hinkley_delta; /*!< Detection threshold for the Hinkley's test. */ + bool m_hinkley_computealphadelta; /*!< If true, compute alpha and delta of the Hinkley's using the stdev of the signal.*/ + float m_hinkley_h; /*!< Alarm factor permitting to compute alpha from the standard deviation of the signal.*/ + float m_hinkley_k; /*!< Detection factor permitting to compute delta from the standard deviation of the signal.*/ + bool m_shewhart_useWECO; /*!< If true, use the WECO rules for additionnal subtests for Shewhart's test.*/ + bool m_shewhart_rules[vpStatisticalTestShewhart::COUNT_WECO - 1]; /*!< Rules for the Shewart's test. True activate a WECO rule, false deactivate it.*/ + float m_sigma_h; /*!< Alarm factor for the sigma test.*/ + + ParametersForAlgo() + : m_test_nbsamples(30) + , m_cusum_h(4.76f) + , m_cusum_k(1.f) + , m_ewma_alpha(0.1f) + , m_hinkley_alpha(4.76f) + , m_hinkley_delta(1.f) + , m_hinkley_computealphadelta(false) + , m_hinkley_h(4.76f) + , m_hinkley_k(1.f) + , m_shewhart_useWECO(false) + , m_sigma_h(3.f) + { + std::memcpy(m_shewhart_rules, vpStatisticalTestShewhart::CONST_ALL_WECO_ACTIVATED, (vpStatisticalTestShewhart::COUNT_WECO - 1) * sizeof(bool)); + memcpy(m_test_activatedalarms, CONST_ALL_ALARM_ON, vpStatisticalTestAbstract::MEAN_DRIFT_COUNT * sizeof(bool)); + m_test_activatedalarms[vpStatisticalTestAbstract::MEAN_DRIFT_NONE] = false; + m_test_nbactivatedalarms = meanDriftArrayToNbActivated(m_test_activatedalarms); + } +}ParametersForAlgo; +//! [Structure_Parameters] +} + +int testOnSynthetic(const TutorialMeanDrift::TypeTest &type, const TutorialMeanDrift::ParametersForAlgo parameters, + const float &mean, const float &mean_drift, const float &stdev) +{ + const float dt = 10.f; // Emulate a 10ms period + + //! [Plot_Init] + vpPlot plotter(1); + plotter.initGraph(0, 1); + plotter.setTitle(0, "Evolution of the signal"); + plotter.setUnitX(0, "Frame ID"); + plotter.setUnitY(0, "No units"); + //! [Plot_Init] + + //! [Test_Creat] + unsigned int idFrame = 0; + vpStatisticalTestAbstract *p_test = nullptr; + switch (type) { + case TutorialMeanDrift::EWMA_TYPE_TEST: + p_test = new vpStatisticalTestEWMA(parameters.m_ewma_alpha); + break; + case TutorialMeanDrift::HINLKEY_TYPE_TEST: + p_test = new vpStatisticalTestHinkley(parameters.m_hinkley_alpha, parameters.m_hinkley_delta, parameters.m_test_nbsamples); + break; + case TutorialMeanDrift::MEAN_ADJUSTED_CUSUM_TYPE_TEST: + p_test = new vpStatisticalTestMeanAdjustedCUSUM(parameters.m_cusum_h, parameters.m_cusum_k, parameters.m_test_nbsamples); + break; + case TutorialMeanDrift::SHEWHART_TYPE_TEST: + p_test = new vpStatisticalTestShewhart(parameters.m_shewhart_useWECO, parameters.m_shewhart_rules, parameters.m_test_nbsamples); + break; + case TutorialMeanDrift::SIGMA_TYPE_TEST: + p_test = new vpStatisticalTestSigma(parameters.m_sigma_h, parameters.m_test_nbsamples); + break; + default: + throw(vpException(vpException::badValue, TutorialMeanDrift::typeTestToString(type) + " is not handled.")); + break; + } + + if ((type == TutorialMeanDrift::HINLKEY_TYPE_TEST) && parameters.m_hinkley_computealphadelta) { + // Initialization of Hinkley's test in automatic mode + delete p_test; + p_test = new vpStatisticalTestHinkley(parameters.m_hinkley_h, parameters.m_hinkley_k, true, parameters.m_test_nbsamples); + } + //! [Test_Creat] + + float signal; + + //! [Test_Init] + // Initial computation of the mean and stdev of the input signal + for (unsigned int i = 0; i < parameters.m_test_nbsamples; ++i) { + vpGaussRand rndGen(stdev, mean, static_cast(idFrame) * dt); + signal = rndGen(); + p_test->testDownUpwardMeanDrift(signal); + ++idFrame; + } + //! [Test_Init] + + std::cout << "Estimated mean of the input signal: " << p_test->getMean() << std::endl; + std::cout << "Estimated stdev of the input signal: " << p_test->getStdev() << std::endl; + + //! [Loop_Monitor] + float mean_eff = mean; + bool hasToRun = true; + vpStatisticalTestAbstract::vpMeanDriftType drift_type = vpStatisticalTestAbstract::MEAN_DRIFT_NONE; + while (hasToRun) { + vpGaussRand rndGen(stdev, mean_eff, static_cast(idFrame) * dt); + signal = rndGen(); + plotter.plot(0, 0, idFrame - parameters.m_test_nbsamples, signal); + drift_type = p_test->testDownUpwardMeanDrift(signal); + if ((drift_type != vpStatisticalTestAbstract::MEAN_DRIFT_NONE) && (parameters.m_test_activatedalarms[drift_type])) { + hasToRun = false; + } + else { + mean_eff += mean_drift; + ++idFrame; + } + } + //! [Loop_Monitor] + + //! [Failure_Debrief] + std::cout << "Test failed at frame: " << idFrame - parameters.m_test_nbsamples << std::endl; + std::cout << "Type of mean drift: " << vpStatisticalTestAbstract::vpMeanDriftTypeToString(drift_type) << std::endl; + std::cout << "Last signal value: " << signal << std::endl; + if (type == TutorialMeanDrift::EWMA_TYPE_TEST) { + vpStatisticalTestEWMA *p_testEwma = dynamic_cast(p_test); + std::cout << "\tw(t) = " << p_testEwma->getWt() << std::endl; + } + else if (type == TutorialMeanDrift::MEAN_ADJUSTED_CUSUM_TYPE_TEST) { + vpStatisticalTestMeanAdjustedCUSUM *p_testCusum = dynamic_cast(p_test); + std::cout << "\tLower cusum = " << p_testCusum->getTestSignalMinus() << std::endl; + std::cout << "\tUpper cusum = " << p_testCusum->getTestSignalPlus() << std::endl; + } + else if (type==TutorialMeanDrift::SHEWHART_TYPE_TEST) { + vpStatisticalTestShewhart *p_testShewhart = dynamic_cast(p_test); + std::vector signal = p_testShewhart->getSignals(); + size_t nbSignal = signal.size(); + std::cout << "Signal history = [ "; + for (size_t i = 0; i < nbSignal; ++i) { + std::cout << signal[i] << " "; + } + std::cout << "]" << std::endl; + std::cout << "\tWECO alarm type = " << vpStatisticalTestShewhart::vpWecoRulesAlarmToString(p_testShewhart->getAlarm()) << std::endl; + } + else if (type == TutorialMeanDrift::HINLKEY_TYPE_TEST) { + vpStatisticalTestHinkley *p_hinkley = dynamic_cast(p_test); + float Mk = p_hinkley->getMk(); + float Nk = p_hinkley->getNk(); + float Sk = p_hinkley->getSk(); + float Tk = p_hinkley->getTk(); + std::cout << "S+(t) = " << Tk - Nk < alarmNames; + bool hasNotFoundNextParams = true; + for (int j = 1; ((i + j) < argc) && hasNotFoundNextParams; ++j) { + std::string candidate(argv[i+j]); + if (candidate.find("--") != std::string::npos) { + // This is the next command line parameter + hasNotFoundNextParams = false; + } + else { + // This is a name + alarmNames.push_back(candidate); + ++nbArguments; + } + } + TutorialMeanDrift::vectorOfStringToMeanDriftTypeArray(alarmNames, parameters.m_test_activatedalarms); + parameters.m_test_nbactivatedalarms = TutorialMeanDrift::meanDriftArrayToNbActivated(parameters.m_test_activatedalarms); + i += nbArguments; + } + else if ((std::string(argv[i]) == "--cusum-h") && ((i + 1) < argc)) { + parameters.m_cusum_h = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--cusum-k") && ((i + 1) < argc)) { + parameters.m_cusum_k = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--ewma-alpha") && ((i + 1) < argc)) { + parameters.m_ewma_alpha = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--hinkley-alpha") && ((i + 1) < argc)) { + parameters.m_hinkley_alpha = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--hinkley-delta") && ((i + 1) < argc)) { + parameters.m_hinkley_delta = std::atof(argv[i + 1]); + ++i; + } + else if (std::string(argv[i]) == "--hinkley-compute") { + parameters.m_hinkley_computealphadelta = true; + } + else if ((std::string(argv[i]) == "--hinkley-h") && ((i + 1) < argc)) { + parameters.m_hinkley_h = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--hinkley-k") && ((i + 1) < argc)) { + parameters.m_hinkley_k = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--shewhart-rules") && (i + vpStatisticalTestShewhart::COUNT_WECO - 1 < argc)) { + for (int j = 0; j < vpStatisticalTestShewhart::COUNT_WECO - 1; ++j) { + std::string argument = std::string(argv[i + 1 + j]); + if ((argument.find("on") != std::string::npos) || (argument.find("ON") != std::string::npos)) { + parameters.m_shewhart_rules[j] = true; + } + else { + parameters.m_shewhart_rules[j] = false; + } + } + i += vpStatisticalTestShewhart::COUNT_WECO - 1; + } + else if (std::string(argv[i]) == "--shewhart-weco") { + parameters.m_shewhart_useWECO = true; + } + else if ((std::string(argv[i]) == "--sigma-h") && ((i + 1) < argc)) { + parameters.m_sigma_h = std::atof(argv[i + 1]); + ++i; + } + else if ((std::string(argv[i]) == "--help") || (std::string(argv[i]) == "-h")) { + std::cout << "\nSYNOPSIS " << std::endl + << argv[0] + << " [--test ]" + << " [--nb-samples ]" + << " [--alarms ]" + << " [--mean ]" + << " [--mean-drift ]" + << " [--stdev ]" + << " [--cusum-h ]" + << " [--cusum-k ]" + << " [--ewma-alpha ]" + << " [--hinkley-alpha <]0; inf[>]" + << " [--hinkley-delta <]0; inf[>]" + << " [--hinkley-compute]" + << " [--hinkley-h <]0; inf[>]" + << " [--hinkley-k <]0; inf[>]" + << " [--shewhart-rules <3-sigma:{on|off} 2-sigma:{on|off} 1-sigma:{on|off} same-side:{on|off}>" + << " [--shewhart-weco]" + << " [--sigma-h ]" + << " [--help,-h]" << std::endl; + std::cout << "\nOPTIONS " << std::endl + << " --test " << std::endl + << " Type of test to perform on the data." << std::endl + << " Available values: " << TutorialMeanDrift::getAvailableTypeTest() << std::endl + << std::endl + << " --nb-samples " << std::endl + << " Number of samples to compute the mean and standard deviation of the monitored signal." << std::endl + << " Default: " << parameters.m_test_nbsamples << std::endl + << std::endl + << " --alarms " << std::endl + << " Set the mean drift alarms to monitor." << std::endl + << " Default: " << TutorialMeanDrift::meanDriftArrayToString(parameters.m_test_activatedalarms) << std::endl + << " Available: " << vpStatisticalTestAbstract::getAvailableMeanDriftType() << std::endl + << std::endl + << " --mean " << std::endl + << " Mean of the signal." << std::endl + << " Default: " << opt_mean<< std::endl + << std::endl + << " --mean-drift " << std::endl + << " Mean drift for the synthetic data." << std::endl + << " Default: " << opt_meandrift << std::endl + << std::endl + << " --stdev " << std::endl + << " Standard deviation of the signal." << std::endl + << " Default: " << opt_stdev << std::endl + << std::endl + << " --cusum-h " << std::endl + << " The alarm factor that permits to the CUSUM test to determine when the process is out of control" << std::endl + << " from the standard deviation of the signal." << std::endl + << " Default: " << parameters.m_cusum_h << std::endl + << std::endl + << " --cusum-k " << std::endl + << " The factor that permits to determine the slack of the CUSUM test, " << std::endl + << " i.e. the minimum value of the jumps we want to detect, from the standard deviation of the signal." << std::endl + << " Default: " << parameters.m_cusum_k << std::endl + << std::endl + << " --ewma-alpha " << std::endl + << " Forgetting factor for the Exponential Weighted Moving Average (EWMA)." << std::endl + << " Default: " << parameters.m_ewma_alpha << std::endl + << std::endl + << " --hinkley-alpha " << std::endl + << " The alarm threshold indicating that a mean drift occurs for the Hinkley's test." << std::endl + << " Default: " << parameters.m_hinkley_alpha << std::endl + << std::endl + << " --hinkley-delta " << std::endl + << " Detection threshold indicating minimal magnitude we want to detect for the Hinkley's test." << std::endl + << " Default: " << parameters.m_hinkley_delta << std::endl + << std::endl + << " --hinkley-compute" << std::endl + << " If set, the Hinkley's test will compute the alarm and detection thresholds" << std::endl + << " from the standard deviation of the input signal." << std::endl + << " Default: disabled" << std::endl + << std::endl + << " --hinkley-h " << std::endl + << " Alarm factor permitting to compute the alarm threshold for the Hinkley's test." << std::endl + << " Default: " << parameters.m_hinkley_h << std::endl + << std::endl + << " --hinkley-k " << std::endl + << " Detection factor permitting to compute the Detection threshold for the Hinkley's test." << std::endl + << " Default: " << parameters.m_hinkley_k << std::endl + << std::endl + << " --shewhart-rules <3-sigma:{on|off} 2-sigma:{on|off} 1-sigma:{on|off} same-side:{on|off}>" << std::endl + << " Choose the WECO additionnal tests for the Shewhart's test to use. To activate them, --shewart-weco must be used." << std::endl + << " Default: ON ON ON ON" << std::endl + << std::endl + << " --shewhart-weco" << std::endl + << " Activate the WECO additionnal tests for the Shewhart's test." << std::endl + << " Default: deactivated" << std::endl + << std::endl + << " --sigma-h " << std::endl + << " The alarm factor of the sigma test." << std::endl + << " Default: " << parameters.m_sigma_h << std::endl + << std::endl + << " --help, -h" << std::endl + << " Display this helper message." << std::endl + << std::endl; + return EXIT_SUCCESS; + } + else { + std::cout << "\nERROR " << std::endl << " Unknown option " << argv[i] << std::endl; + return EXIT_FAILURE; + } + ++i; + } + + if (parameters.m_test_nbactivatedalarms == 0) { + throw(vpException(vpException::badValue, "Error, at least one type of alarm must be monitored. See " + std::string(argv[0]) + " --help")); + return EXIT_FAILURE; + } + + std::cout << " Activated statistical test : " << TutorialMeanDrift::typeTestToString(opt_typeTest) << std::endl; + std::cout << " Activated alarms : " << TutorialMeanDrift::meanDriftArrayToString(parameters.m_test_activatedalarms) << std::endl; + std::cout << " Nb samples for statistics computation: " << parameters.m_test_nbsamples << std::endl; + std::cout << " Alarm factor CUSUM test : " << (opt_typeTest == TutorialMeanDrift::MEAN_ADJUSTED_CUSUM_TYPE_TEST ? TutorialMeanDrift::numberToString(parameters.m_cusum_h) : "N/A") << std::endl; + std::cout << " Detection factor CUSUM test : " << (opt_typeTest == TutorialMeanDrift::MEAN_ADJUSTED_CUSUM_TYPE_TEST ? TutorialMeanDrift::numberToString(parameters.m_cusum_k) : "N/A") << std::endl; + std::cout << " Forgetting factor EWMA test : " << (opt_typeTest == TutorialMeanDrift::EWMA_TYPE_TEST ? TutorialMeanDrift::numberToString(parameters.m_ewma_alpha) : "N/A") << std::endl; + std::cout << " Alarm threshold Hinkley's test : " << ((opt_typeTest == TutorialMeanDrift::HINLKEY_TYPE_TEST) && (!parameters.m_hinkley_computealphadelta) ? TutorialMeanDrift::numberToString(parameters.m_hinkley_alpha) : "N/A") << std::endl; + std::cout << " Detection threshold Hinkley's test : " << ((opt_typeTest == TutorialMeanDrift::HINLKEY_TYPE_TEST) && (!parameters.m_hinkley_computealphadelta) ? TutorialMeanDrift::numberToString(parameters.m_hinkley_delta) : "N/A") << std::endl; + std::cout << " Alarm factor Hinkley's test : " << ((opt_typeTest == TutorialMeanDrift::HINLKEY_TYPE_TEST) && parameters.m_hinkley_computealphadelta ? TutorialMeanDrift::numberToString(parameters.m_hinkley_h) : "N/A") << std::endl; + std::cout << " Detection factor Hinkley's test : " << ((opt_typeTest == TutorialMeanDrift::HINLKEY_TYPE_TEST) && parameters.m_hinkley_computealphadelta ? TutorialMeanDrift::numberToString(parameters.m_hinkley_k) : "N/A") << std::endl; + std::cout << " Shewhart's test set of WECO rules : " << (parameters.m_shewhart_useWECO && (opt_typeTest == TutorialMeanDrift::SHEWHART_TYPE_TEST) ? TutorialMeanDrift::wecoRulesToString(parameters.m_shewhart_rules) : "N/A") << std::endl; + std::cout << " Shewhart's test use WECO rules : " << (opt_typeTest == TutorialMeanDrift::SHEWHART_TYPE_TEST ? TutorialMeanDrift::boolToString(parameters.m_shewhart_useWECO && (opt_typeTest == TutorialMeanDrift::SHEWHART_TYPE_TEST)) : "N/A") << std::endl; + std::cout << " Alarm factor Sigma test : " << (opt_typeTest == TutorialMeanDrift::SIGMA_TYPE_TEST ? TutorialMeanDrift::numberToString(parameters.m_sigma_h) : "N/A") << std::endl; + std::cout << " Actual mean of the input signal: " << opt_mean << std::endl; + std::cout << " Actual stdev of the input signal: " << opt_stdev << std::endl; + std::cout << " Mean drift of the input signal: " << opt_meandrift << std::endl; + + return testOnSynthetic(opt_typeTest, parameters, opt_mean, opt_meandrift, opt_stdev); +} +#else +int main() +{ + std::cerr << "Recompile ViSP with display capabilities in order to use this tutorial." << std::endl; + return EXIT_FAILURE; +} +#endif