swh:1:snp:327e80ff7d96a6df4667f20ec612357647e9f63b
Raw File
Tip revision: 71cea6b7054d14cabe643a38cd2b88a50ed09da1 authored by Software Heritage on 01 February 2019, 00:00:00 UTC
hal: Deposit 358 in collection hal
Tip revision: 71cea6b
tst_chronomodeltest.cpp
#ifndef UNIT_TEST
#define UNIT_TEST

#include <QString>
#include <QProgressDialog>
#include <QtTest>

#include "ChronoApp.h"
#include <QtWidgets>
#include "MainController.h"
#include "StdUtilities.h"
#include "Functions.h"
#include "MetropolisVariable.h"
#include "Model.h"
#include "Phase.h"
#include "Event.h"

#include "fftw3.h"

#include <iostream>
#include <cmath>
#include <errno.h>
//#include <fenv.h>
#include "Generator.h"
#include "Functions.h"



QString formatFunction(const double valueToFormat, const bool forCSV = false);

/**
 * @brief myFuzzyCompare to replace qFuzzyCompare, because I can change precision
 * @param p1
 * @param p2
 * @return
 */
#define PRECISION 0.0000000001 //default 0.000000000001
static inline bool myFuzzyCompare(double p1, double p2)
{
    return (qAbs(p1 - p2) <= PRECISION * qMin(qAbs(p1), qAbs(p2)));
}

class ChronomodelTest : public QObject
{
    Q_OBJECT

public:
    ChronomodelTest();
   // QString formatFunction(const double valueToFormat, const bool forCSV = false);
    Quartiles quartilesTypeTest(const QVector<double> &trace, const int quartileType, const double p);
    QPair<int, double> gammaQuartile(const QVector<double> &trace, const int quartileType, const double p);

private Q_SLOTS:
    //void testCase1();
    void mersenneTwister();
    void timeRange();
    void trace();
    void fft();
    void hpd();
    void tempo();
    void quantile();
    void quartileForTrace();
};

ChronomodelTest::ChronomodelTest()
{
}

/* void ChronomodelTest::testCase1()
{
    QString str = "Hello";
        QCOMPARE(str.toUpper(), QString("HELLO"));

    QBENCHMARK {
            str.trimmed();
        }
}
*/

/**
  * @brief Debuggage mersenne twister
 *
 */

void ChronomodelTest::mersenneTwister()
{
   Generator::initGenerator(int(111));

   QString text;
   const QString reponse111 ("; -758.484; -396.421; 707.813; -404.568; -244.275; -581.953; 299.128; 203.405; 424.546; 840.891");

   for (int i=0; i<10; ++i)
       text= text+"; "+QString::number( Generator::randomUniform(-1000, 1000));

   QVERIFY(reponse111 == text);

   // 2d test
   text="";
   Generator::initGenerator(int(3649));

   const QString reponse3694 ("; -757.681; -1.8087; 363.35; -982.55; -505.08; 683.574; -629.103; -477.18; 691.793; -829.783");

   for (int i=0; i<10; ++i)
       text= text+"; "+QString::number( Generator::randomUniform(-1000, 1000));

   QCOMPARE(reponse3694, text);

}

/**
  * @brief Debuggage timerange
  *
  */
void ChronomodelTest::timeRange()
{
    QVector<double> trace1;
    trace1<<1.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4<<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8<<3.12;

    QVector<double> trace2;
    trace2<<5.18<<7.5<<6.18<<13.4<<9.3<<11.4<<9.1<<14.15<<13.19<<9.6<<11.3<<9.9<<15.13<<13.20<<14.14<<6.2<<9.2<<8.2<<7.13<<9.11;

     QPair<double, double> range = timeRangeFromTraces(trace1, trace2, 95., "test timeRangeFromTraces 1");
     QPair<double, double> ansRang = qMakePair<double, double> (1.1, 14.199);

    QVERIFY(range.first == ansRang.first);
    QVERIFY(qFuzzyCompare(range.second, ansRang.second));

    QVector<double> traceBeta;
    traceBeta<<1.1<<3.8<<5.6<<1.99<<7.4<<4.2<<5.7<<3.6<<3.4<<2.7<<3.1<<4.3<<3.87<<4.91<<6.11<<5.2<<1.18<<2.4<<6.8<<3.12;
    traceBeta<<1.11<<3.81<<5.61<<1.991<<7.41<<4.21<<5.71<<3.61<<3.41<<2.71<<3.11<<4.31<<3.871<<4.911<<6.111<<5.21<<1.181<<2.41<<6.81<<3.121;
    traceBeta<<1.112<<3.812<<5.612<<1.9912<<7.412<<4.212<<5.712<<3.612<<3.412<<2.712<<3.112<<4.312<<3.8712<<4.9112<<6.1112<<5.212<<1.1812<<2.412<<6.812<<3.1212;

    QVector<double> traceAlpha;
    traceAlpha<<11.9<<7.5<<12.18<<13.4<<9.3<<11.4<<9.1<<14.15<<13.019<<10.6<<11.3<<9.9<<15.13<<13.19<<14.14<<6.2<<9.2<<8.2<<7.13<<9.01;
    traceAlpha<<11.91<<7.51<<12.181<<13.41<<9.31<<11.41<<9.11<<14.151<<13.0191<<10.61<<11.31<<9.91<<15.131<<13.1921<<14.141<<6.21<<9.21<<8.21<<7.131<<9.011;
    traceAlpha<<11.912<<7.512<<12.1812<<13.412<<9.312<<11.412<<9.112<<14.1512<<13.01912<<10.612<<11.312<<9.912<<15.1312<<13.19212<<14.1412<<6.212<<9.212<<8.212<<7.1312<<9.0112;

    QPair<double, double> gapRange = gapRangeFromTraces(traceBeta, traceAlpha, 95., "test gapRangeFromTraces");
    QPair<double, double> ansGap = qMakePair<double, double> ( -INFINITY , INFINITY);

    QCOMPARE(gapRange, ansGap);

    QVector<double>tBeta;
    tBeta<<0.0<<5.0<<5.0<<10.0;

    QVector<double>tAlpha;
    tAlpha<<11.0<<16.0<<16.0<<21.0;

    QPair<double,double> gapRang2 = gapRangeFromTraces(tBeta, tAlpha, 50., "test gapRangeFromTraces 2");

    QPair<double, double> ansGap2 = qMakePair<double, double> (5.005,13.5037481259);
    QVERIFY(myFuzzyCompare(gapRang2.first, ansGap2.first));
    QVERIFY(myFuzzyCompare(gapRang2.second,ansGap2.second));
}

/**
  * @brief Debuggage fullrunTrace
  */

void ChronomodelTest::trace()
{

    ChainSpecs chain1;
    chain1.mNumBurnIter = 3; //used
    chain1.mBurnIterIndex=1;

    chain1.mBatchIterIndex = 2;
    chain1.mBatchIndex = 2; //used
    chain1.mNumBatchIter = 3; //used

    chain1.mNumRunIter = 10; //used
    chain1.mRunIterIndex = 5;

    chain1.mThinningInterval = 1; //used

    MetropolisVariable fullTrace;
    fullTrace.mRawTrace = new QVector<double>();
    // Adding chain1 in fullTrace
    for (auto i=1; i<21; ++i)
        fullTrace.mRawTrace->append(i);
    //----------chain2
    ChainSpecs chain2;
    chain2.mNumBurnIter = 3; //used
    chain2.mBurnIterIndex=1;

    chain2.mBatchIterIndex = 2;
    chain2.mBatchIndex = 2; //used
    chain2.mNumBatchIter = 3; //used

    chain2.mNumRunIter = 10; //used 2 iter more unused
    chain2.mRunIterIndex = 5;

    chain2.mThinningInterval = 1; //used
 // Adding chain2 in fullTrace
    for (auto i=1; i<21; ++i)
        fullTrace.mRawTrace->append(i);

    //----------chain3
    ChainSpecs chain3;
 /*   chain3.mNumBurnIter=2;
    chain3.mBurnIterIndex=1;

    chain3.mBatchIterIndex = 2;
    chain3.mBatchIndex = 1;
    chain3.mNumBatchIter = 1;

    chain3.mBurnIterIndex = 1;
    chain3.mNumBurnIter = 2;

    chain3.mNumRunIter = 5;
    chain3.mRunIterIndex = 5;

    chain3.mThinningInterval = 1;
    */
    chain3 = chain1;
    // Adding chain3 in fullTrace
    for (auto i=1; i<21; ++i)
        fullTrace.mRawTrace->append(i);


    // Create List of chain
    QList<ChainSpecs> chains;
    chains.append(chain1);
    chains.append(chain2);
    chains.append(chain3);

    // Using runRawTraceForChain function
    QVector<double> runTrace1 =  fullTrace.runRawTraceForChain(chains, 0);
  //  qDebug()<<"trace1"<<runTrace1;

    QVector<double> runTrace2 =  fullTrace.runRawTraceForChain(chains, 1);
 //   qDebug()<<"trace2"<<runTrace2;

    // Controle
    QVERIFY(runTrace1 == runTrace2);


    QVector<double> runTrace3 =  fullTrace.runRawTraceForChain(chains, 2);
//    qDebug()<<"trace3"<<runTrace3;

    // Controle
    QVERIFY(runTrace1 == runTrace3);

    fullTrace.mFormatedTrace->resize(fullTrace.mRawTrace->size());
    std::copy(fullTrace.mRawTrace->cbegin(), fullTrace.mRawTrace->cend(), fullTrace.mFormatedTrace->begin());

    // Using fullRunTrace function
    QVector<double> fullRunTrace = fullTrace.fullRunTrace(chains);
//    qDebug()<<"fullRunTrace"<<fullRunTrace;

    // Create good Answer
    QVector<double> fullRunTraceVerif = runTrace1;
    fullRunTraceVerif.append(runTrace2);
    fullRunTraceVerif.append(runTrace3);

    // Controle
    QVERIFY(fullRunTrace == fullRunTraceVerif);
}

QString formatFunction(const double valueToFormat, const bool forCSV)
{
    (void) forCSV;
    QLocale loc = QLocale();
    return loc.toString(valueToFormat);
}

void ChronomodelTest::fft()
{
    int fftLen (32);
    const double bandwidth (1.06);
    const double tmin (-10);
    const double tmax (20);

    ChainSpecs chain1;
    chain1.mNumBurnIter = 1; //used
    chain1.mBurnIterIndex=1;

    chain1.mBatchIterIndex = 2;
    chain1.mBatchIndex = 0; //used
    chain1.mNumBatchIter = 0; //used

    chain1.mNumRunIter = 18; //used
    chain1.mRunIterIndex = 5;

    chain1.mThinningInterval = 1; //used


    QVector<double> trace1;
    trace1<<1.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4<<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8<<3.12;

     MetropolisVariable density1;
     density1.mSupport = MetropolisVariable::eR;
     density1.setName("test unit");
     //density1.mFormatedTrace = &trace1;

     density1.mFormatedTrace->resize(trace1.size());
     std::copy(trace1.cbegin(), trace1.cend(), density1.mFormatedTrace->begin());

     QList<ChainSpecs> chain1Spe;
     chain1Spe.append(chain1);

    density1.generateHistos(chain1Spe, fftLen, bandwidth, tmin, tmax) ;
    //qDebug()<<"histo1"<<density1.mHisto;
    QVERIFY(qFuzzyCompare(density1.mHisto.firstKey(), -5.35157149234));

    ChainSpecs chain2 (chain1);
    ChainSpecs chain3 (chain1);

    QList<ChainSpecs> chains;
    chains.append(chain1);
    chains.append(chain2 );
    chains.append(chain3);

    QVector<double> trace2 (trace1);
    trace2<<1.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4<<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8<<3.12;
    trace2<<1.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4<<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8<<3.12;

    MetropolisVariable density2;
    density2.mSupport = MetropolisVariable::eR;
    density2.setName("test unit2");

    density2.mFormatedTrace->resize(trace2.size());
    std::copy(trace2.cbegin(), trace2.cend(), density2.mFormatedTrace->begin());

    // Using generateHistos() function which use fft
    density2.generateHistos(chains, fftLen, bandwidth, tmin, tmax);
    //qDebug()<<"histo2"<<density2.mHisto;

    QVERIFY(qFuzzyCompare(density2.mHisto.firstKey(), -4.06316390051));
}

void ChronomodelTest::hpd()
{
    //--- test HPD --- Only one chain
    ChainSpecs chainHPD1;
    chainHPD1.mNumBurnIter = 10; //used
    chainHPD1.mBurnIterIndex=0;

    chainHPD1.mBatchIterIndex = 0;
    chainHPD1.mBatchIndex = 10; //used
    chainHPD1.mNumBatchIter = 100; //used

    chainHPD1.mNumRunIter = 200; //used
    chainHPD1.mRunIterIndex = 0;

    chainHPD1.mThinningInterval = 5; //used


    QList<ChainSpecs> chainHPDs;
    chainHPDs.append(chainHPD1);


    QVector<double> traceHPD;
    traceHPD<<414.408053235927<<414.408053235927<<412.436722904401<<412.436722904401<<412.436722904401<<412.436722904401<<411.711834420135<<404.966479053834<<413.69476291588<<413.14348411438;
    traceHPD<<406.619756727929<<413.176346507546<<411.529083778499<<411.227728655378<<407.207172079998<<410.679958716126<<408.971364930604<<414.203568623728<<414.203568623728<<414.203568623728<<411.348234239404<<410.15306043888<<411.14005458578;
    traceHPD<<416.452988344086<<418.839268060008<<418.090686816061<<418.090686816061<<415.531647611625<<409.107137284467<<413.405114566042<<413.405114566042<<411.298588403653<<409.762980960855<<410.324567473625<<407.990387626478<<404.225507494716;
    traceHPD<<402.527803036477<<402.527803036477<<405.633185027145<<400.149771473239<<400.149771473239<<400.149771473239<<400.149771473239<<401.67722548249<<399.114505580032<<402.456807169368<<405.831190357112<<405.08442739979<<410.244747379223<<411.222409464041;
    traceHPD<<414.236725854778<<419.165939997713<<419.352313750215<<421.092859866884<<417.993652667721<<417.979263828908<<420.348942857362<<421.680930187683<<417.091573180055<<413.744198236104<<412.537340193616<<412.346493302029<<412.345338514746<<408.887403784499<<406.604157034121<<403.133063186189<<408.351952584557<<417.573893292046;
    traceHPD<<415.651152836086<<415.651152836086<<415.651152836086<<415.651152836086<<412.223046532598<<412.13383341022<<412.13383341022<<414.250789991845<<411.072839073443<<412.169501994485<<412.588456504881<<400.369216807319<<400.369216807319<<400.369216807319<<400.760555444567<<400.760555444567<<406.820239884521<<402.777128776677<<405.077621444105;
    traceHPD<<401.678868627691<<400.252193795241<<404.500421682083<<405.361935980527<<406.207509954835<<401.406274689322<<401.406274689322<<403.003477765683<<403.489160828872<<405.569274976211<<402.442651467919<<402.77186106237<<402.147304038975<<404.139133966681<<399.012891817584<<398.98217893389<<390.609422003471<<389.548845745389<<388.570589336105<<392.297942975421<<390.594971787731;
    traceHPD<<398.791891393196<<406.333425974158<<399.570839507341<<392.164849428161<<387.897983705244<<390.755840558811<<393.339775262995<<392.851299426197<<399.812418419347<<393.781780514019<<393.781780514019<<390.061932600189<<390.868563595625<<405.399252207496<<406.829388991955<<407.860088954454<<407.860088954454<<415.788371003507<<412.933480296261<<405.633692752962<<401.921164670616<<406.617695439799<<408.630159488386;
    traceHPD<<410.735847971243<<410.735847971243<<414.125022860809<<404.492204474488<<408.090928283779<<402.066684458103<<409.167819877058<<412.171988210103<<406.306972756689<<405.322595798767<<410.381966079171<<399.407933886885<<399.407933886885<<390.774934961717<<394.850923432576<<390.97273447792<<390.97273447792<<396.137341652167<<396.196160011687<<402.85628436897<<402.85628436897<<404.307024942622<<404.807124937587<<403.109841224752<<400.614510957774;
    traceHPD<<398.015655979332<<400.617156911614<<402.461411349569<<402.461411349569<<402.461411349569<<400.701443321604<<399.60984376723<<398.832540084484<<402.414134539381<<401.471651602716<<397.461251206002<<394.808077013267<<394.264255240029<<397.575775243442<<392.61685603766<<393.11004924234<<392.24866461044<<396.614211124405<<399.210660943994<<400.443632093341<<400.443632093341<<404.042670239977<<404.042670239977<<404.042670239977<<404.042670239977<<404.042670239977;
    traceHPD<<402.470042972322<<402.470042972322<<402.470042972322<<402.037722549295<<396.613578626651<<393.206917659911<<397.692432427219<<401.375126519345<<401.375126519345<<401.375126519345<<401.375126519345<<401.375126519345<<403.111521396346<<403.739982501994<<399.63432512217<<399.228241452197<<399.074448898776<<386.819171470435<<386.819171470435<<386.819171470435<<386.819171470435<<389.736981817529<<391.420741877606<<393.764158078232<<391.690193455264<<389.633203710203<<389.294824984334<<389.294824984334;
    traceHPD<<385.718179846552<<390.365889867664<<386.646349612216<<388.351353326969<<388.351353326969<<390.082650456181<<391.986593256303<<391.986593256303<<384.74202130658<<384.74202130658<<392.329318414147<<388.238202126929<<391.547116314306<<391.547116314306<<386.862973576222<<395.746756167884<<395.53547248711<<386.101483286121<<386.801297904029<<386.801297904029<<392.471742364309<<396.130375870542<<388.874908256844<<389.731227717135<<393.293565692474<<386.602496937093<<389.986215162834<<389.080577372964<<389.444920433468<<389.794248527088;
    traceHPD<<395.863656210345<<391.0003506294<<392.561688365552<<395.015454208541<<394.540993255349<<391.794148611399<<391.794148611399<<385.82417443828<<386.227004535155<<386.227004535155<<385.754583469814<<385.754583469814<<385.754583469814<<385.754583469814<<393.104787893047<<394.18588675156<<394.18588675156<<387.992928275005<<387.992928275005<<389.868006656267<<387.154513109008<<387.845527311302<<387.845527311302<<388.737355852481<<396.072266202164<<400.848865961229<<395.177872257826<<395.177872257826<<398.97718975905<<397.984081600766<<397.984081600766<<399.47073982129<<399.172089537439<<398.414893663652<<398.414893663652;
    traceHPD<<398.414893663652<<401.406359701689<<394.774616224662<<394.594378504006<<392.672430887586<<386.395815746491<<386.395815746491<<386.928233396189<<388.513783425642<<393.068067295845<<396.237395395151<<392.711410221192<<398.564026407974<<399.126860670782<<399.126860670782<<396.950365632326<<387.210931275424<<390.597366775789<<398.292153859114<<395.220575836686<<399.210204804223<<399.845745133883<<399.845745133883<<399.845745133883<<394.654422090703<<390.239437635159<<386.323304825494<<386.323304825494<<382.630745285456<<382.630745285456<<382.630745285456<<382.630745285456<<383.842590333824<<379.306123462394<<371.818360309935<<371.741739660982<<370.652223891523<<370.652223891523<<369.065367823052<<369.548111318734;
    traceHPD<<369.548111318734<<372.414667148777<<368.914452092381<<368.914452092381<<368.914452092381<<368.914452092381<<368.425478187388<<374.496551136715<<366.200218760374<<378.161482986371<<375.023933930611<<375.023933930611<<370.162325923627<<369.596205385303<<365.579067194548<<365.579067194548<<367.582081778992<<361.159521610803<<367.529563019744<<371.033972903679<<374.679123943658<<374.679123943658<<374.454223043217<<374.454223043217<<370.892007260271<<374.717524298941<<369.965893441175<<368.084302915927<<365.865728667836<<365.578672900005<<368.621859335405<<364.029729571481<<365.987970001355<<366.611351089861<<366.611351089861<<366.611351089861<<367.798223279835<<369.786098998821<<366.99670601046<<366.99670601046<<368.695546698343;
    traceHPD<<365.694427183407<<365.694427183407<<366.180319319617<<366.180319319617<<368.484685538836<<378.62293632334<<381.273744871393<<382.784272783819<<382.784272783819<<382.784272783819<<383.148608776784<<383.208589618073<<383.208589618073<<383.208589618073<<379.665937804328<<379.657307873144<<379.657307873144<<377.849766180868<<375.817298194941<<376.971593942885<<365.006411893557<<365.006411893557<<364.235837839385<<363.710429475451<<363.710429475451<<365.741424466402<<365.741424466402<<365.741424466402<<369.663337659036<<369.663337659036<<365.136857892888<<368.882794151337<<369.768026757394<<369.079778684447<<369.079778684447<<365.120125351566<<370.202171906735<<368.52098704578<<368.52098704578<<368.657971782146<<368.657971782146;
    traceHPD<<367.211478528171<<363.765851274479<<365.381155201819<<365.381155201819<<365.381155201819<<359.680700330139<<361.610497434523<<368.133859985643<<368.133859985643<<362.020797828876<<362.020797828876<<363.339471826316<<363.92560173929<<355.694306487978<<360.015894833516<<359.39885470578<<366.127147241041<<354.374052370797<<350.325958250723<<351.792460451679<<351.792460451679<<351.792460451679<<351.792460451679<<351.371121760587<<350.874505834257<<342.423738783193<<342.423738783193<<342.423738783193<<342.423738783193<<342.423738783193<<342.423738783193<<342.385533905907<<342.385533905907<<342.385533905907<<343.165678359892<<345.581078053658<<342.633865203896<<342.633865203896<<341.874822505034<<346.948050263357<<342.688652147293<<348.816817592353<<342.451941830833;
    traceHPD<<346.200696217542<<348.137036533064<<354.399163293087<<349.998121376179<<344.781042048741<<344.781042048741<<352.069342317908<<359.585457902732<<355.52929776752<<352.766826864977<<351.791028506215<<355.179683651852<<350.983147169651<<343.719353523664<<343.531389895632<<341.981314392451<<341.981314392451<<345.956603916383<<353.352393954279<<350.660315131508<<345.384258960392<<345.384258960392<<345.384258960392<<341.583372000694<<341.583372000694<<349.854665611846<<352.020426424075<<352.020426424075<<354.941503983649<<354.941503983649<<354.887859504137<<352.28981025569<<346.726822367464<<341.495928283528<<341.495928283528<<342.716332441939<<333.880970823601<<336.544604670112<<329.488733622616<<330.070827292924<<330.746254267742<<328.272940487712<<331.85698763246<<332.31350844828<<332.941130065379<<339.926526547775<<331.763225184<<331.763225184;
    traceHPD<<329.309142371605<<332.732848863886<<323.946995477496<<323.946995477496<<323.946995477496<<325.969508750475<<325.969508750475<<324.713409776352<<324.713409776352<<322.655101356143<<327.779594416258<<322.469301180415<<322.56803650841<<322.56803650841<<322.56803650841<<326.456468175661<<328.536560109085<<329.084298452632<<324.655443653518<<323.878528117165<<323.878528117165<<327.604851108206<<326.216888001571<<331.628866337165<<329.546407206846<<330.830599842408<<335.213072153393<<330.702064627615<<329.962102935424<<329.108182537563<<331.879981321525<<333.5471192983<<333.5471192983<<332.343018744164<<331.522043866973<<327.413682306204<<326.588760348691<<332.643016623344<<332.643016623344<<330.594987293227<<328.795514488995<<322.595149935502<<322.595149935502<<322.595149935502<<322.595149935502<<322.595149935502<<322.595149935502<<322.595149935502<<322.595149935502<<322.595149935502<<327.782988886235;
    traceHPD<<332.666139908848<<329.583387422592<<329.583387422592<<327.513140728314<<326.177304005971<<323.369282592049<<323.286061362282<<327.954275454953<<331.714127287941<<339.046300947955<<339.046300947955<<334.446184793816<<327.69373305492<<337.451146675251<<341.377138392625<<337.300220875688<<337.629585879084<<337.629585879084<<329.266691027197<<324.630876267082<<324.630876267082<<324.630876267082<<324.630876267082<<331.494209230382<<331.494209230382<<336.480710974833<<335.095252718263<<325.922990221272<<325.922990221272<<325.922990221272<<324.657441237863<<324.883436797924<<324.883436797924<<318.951833847438<<311.766822885996<<311.641175925045<<311.969010257988<<312.148579018609<<313.574566640483<<311.99447237871<<311.99447237871<<311.99447237871<<312.066513690812<<310.064192532501<<313.312116868981<<313.312116868981<<315.292643206556<<315.34379346171<<316.429515520286<<315.290254313554<<315.014803496384<<310.824878647417<<310.824878647417<<309.085192064164;
    traceHPD<<312.516163024402<<312.463283351601<<316.777548239368<<309.529815540798<<309.845707542111<<306.714564003614<<306.645957198475<<304.089912162255<<301.003502810254<<302.284954893373<<300.475840629344<<298.082261788312<<303.119585146014<<304.488944430092<<301.900389434372<<294.386603092018<<299.245104009498<<295.886590647187<<299.393131268905<<299.328841274111<<292.885413671976<<292.885413671976<<289.857084771695<<302.764763012081<<302.764763012081<<300.547540574117<<300.547540574117<<300.547540574117<<299.940279661663<<298.909051532571<<301.114204733416<<300.752880664847<<300.752880664847<<296.800838611634<<300.408044754069<<295.447745394605<<299.338382337403<<301.524291490993<<302.099659906604<<298.930779832704<<299.070231656719<<294.031047673528<<296.166913448561<<284.926029512043<<284.926029512043<<287.206336060751<<289.909182814292<<297.159956210352<<300.193885914498<<291.743037017262<<293.768060702068<<293.863096264322<<293.365033895613<<286.932368198093<<282.248060802937;
    traceHPD<<277.155985269899<<277.155985269899<<277.155985269899<<273.646850866735<<270.363458602489<<265.415949946688<<259.373945624711<<259.373945624711<<264.654592517873<<264.654592517873<<267.074093457746<<267.074093457746<<259.615269338226<<264.534744796782<<265.371606843418<<262.895363518842<<262.895363518842<<265.686454869927<<253.377102677762<<253.989962719847<<253.989962719847<<257.928816656231<<250.888778041833<<250.888778041833<<255.02905380704<<263.93905320867<<254.695141600829<<254.695141600829<<260.105801050169<<260.679043677023<<252.848494190081<<253.801098498648<<252.032260941763<<260.519074559369<<263.363217453972<<266.092929513119<<269.168553589354<<267.395505057576<<264.231699109617<<270.977853749993<<270.977853749993<<268.177077879901<<268.177077879901<<271.175813639644<<266.485175035962<<267.33995223714<<268.588564220595<<268.588564220595<<267.723762251019<<267.723762251019<<267.723762251019<<266.887545052931<<262.099959933089<<259.959642518007<<254.462376293934<<254.462376293934<<255.926882474279<<251.868481345864<<262.570882402628<<253.662230681332<<251.920869634849<<255.114690503203<<258.012456432514<<262.409484556021<<262.409484556021<<257.486755880142<<257.878462741281<<264.980464046927<<264.980464046927<<267.009475932311<<267.009475932311<<266.034914343054<<262.695040567259<<256.515003207711<<256.515003207711<<261.962847722377<<263.61965091811<<259.549364969228<<259.146308340014<<263.617423963802<<261.803391228123<<261.964237401697<<264.675365015166<<268.03166954277<<268.03166954277<<268.122344235418<<268.122344235418<<265.707704953959<<266.968564591064<<263.855463901891<<263.837749053731<<263.140503015017<<264.508217158556<<267.715878833149<<267.715878833149<<263.977032885882<<263.977032885882<<259.872180756733<<249.952923028581<<253.827152526701<<253.827152526701<<248.301993774296<<248.301993774296<<248.301993774296<<249.877887737806<<249.877887737806<<252.524984320747<<252.524984320747<<254.083934577489<<252.657283119605<<252.657283119605<<248.914023629009<<248.914023629009<<248.914023629009<<253.92039552459<<251.620025936157<<255.838515826664<<255.106337971193<<256.318477351316<<256.318477351316<<255.297798755043<<255.297798755043<<252.717615354557<<256.118826326062<<257.422499279639<<257.520521080325<<259.835017675119<<256.385015617245<<260.347045393209<<265.408350740382<<266.95671819623<<266.43651827959<<264.357681143242<<260.250354930664<<262.3907505218<<260.407172993096<<255.377098231414<<256.732482629582<<258.080559020474<<253.968566323629<<254.650389670196<<246.723911621437<<248.071683774246<<251.163285203089<<249.493112643474<<252.769702713178<<252.769702713178<<249.767075279164<<255.687939038902<<262.077774259267<<256.863054167274<<258.587469753123<<258.587469753123<<258.587469753123<<251.817277719697<<251.817277719697<<257.834811135464<<255.977877070495<<258.413541904515<<258.413541904515<<252.199602482123<<250.038471286128<<251.914650591038<<261.06021249049<<259.453413148471<<259.453413148471<<247.587287412022<<247.587287412022<<244.976451068199<<244.976451068199<<249.953315932868<<245.412831269243<<243.447269026172<<246.031838577846<<246.031838577846<<246.031838577846<<245.743449473778<<245.743449473778<<245.743449473778<<245.743449473778<<245.743449473778<<246.596619899407<<246.170121694801<<245.269365842877<<250.682419341689<<252.128921155813<<253.172415460557<<248.211788295978<<257.996684023578<<259.294037888242<<263.809482006924<<256.717607147205<<259.766949009648<<263.955527296511<<257.710987939398<<255.376683900369<<258.262285612214<<258.262285612214<<258.262285612214<<256.689275774626<<252.590158605626<<252.590158605626<<256.796586567164<<257.521591713257<<261.47438383992<<265.932975944623<<265.932975944623<<265.932975944623<<260.941864162086<<262.571289717947<<264.608638180013<<261.686449962682<<261.686449962682<<265.220380472048<<266.143840908637<<266.685283146313<<262.841249625358<<262.76481993808<<255.819347341415<<255.958532666183<<257.327663374613<<248.058925522198<<253.671232994211<<257.270486840991<<260.433578868097<<258.959529782016<<255.60001852577<<260.10196710842<<255.519438265507<<254.023320167971<<251.102155672831<<253.679919283937<<245.932286035243<<246.74664010999<<250.006616188891<<257.941292640374<<256.360681874049<<255.853137201907<<254.73726014991<<250.020245118803<<250.020245118803<<248.04441065118<<251.158936999423<<253.456319621353<<250.178089252683<<250.104646811619<<247.482465921464<<248.188561027282<<246.41886597623<<247.516665021434<<251.289396763565<<251.553331136105<<247.434569550706<<247.434569550706<<247.434569550706<<247.434569550706<<247.434569550706<<247.434569550706<<247.054016667867<<247.054016667867<<244.711512847412<<233.161910038371<<235.903824091089<<242.34142140883<<242.34142140883<<242.34142140883<<237.182398297211<<241.892890449042<<237.621543231228<<240.546816448617<<235.489214766872<<236.132205732214<<233.931060642479<<233.931060642479<<233.931060642479<<231.893698241506<<231.893698241506<<234.523631936837<<234.523631936837<<230.29040466299<<230.29040466299<<231.014978717445<<234.080196419359<<232.181453734963<<232.181453734963<<232.181453734963<<229.684217695569<<229.684217695569<<230.987389046955<<230.987389046955<<230.987389046955<<236.367640642823<<239.083272591013<<247.119452748052<<246.143380541954<<241.394786871413<<245.918193911942<<247.065039767476<<247.065039767476<<248.511387542991<<247.793387907078<<239.102453619672<<232.351814761583<<232.351814761583<<232.351814761583<<241.608245720703<<240.809074152846<<240.006789765464<<238.484645620216<<235.581600899624<<231.005111475849<<220.337577185127<<219.237420979786<<221.187141035589<<221.187141035589<<220.207270042869<<225.416023622633<<226.208361917769<<226.75191437163<<226.497085581229<<226.803327965426<<230.177950409008<<225.974994062618<<222.100702382225<<221.125440130974<<219.335023515645<<214.36380076509<<219.773429707547<<220.518183581604<<223.896775187469<<223.896775187469<<223.634101086478<<230.852365498617<<230.852365498617<<230.852365498617<<230.852365498617<<226.198909935915<<220.634004210567<<219.142929448475<<218.228120702827<<224.855617336605<<223.503517145126<<227.775818169041<<228.301222977624<<222.312310821823<<223.929177035601<<225.517364610624<<225.517364610624<<226.605152027705<<226.605152027705<<226.605152027705<<226.605152027705<<229.053772392931<<226.736350877758<<221.086005861626<<217.981420962924<<213.731924864975<<210.838287420016<<214.137105208305<<222.757296232435<<221.413835616176<<223.03951977555<<221.385004201965<<225.139087816097<<201.696259404116<<208.048891173069<<202.444749175845<<209.027322565724<<211.286934338954<<205.753268671814<<197.853452746673<<212.695307925484<<208.350442327874<<199.519568208289<<209.23749463309<<207.380478439051<<208.135023612111<<200.131161298219<<189.457750858686<<181.694979365202<<194.892731389807<<194.530891763125<<207.227846421787<<197.582493401<<204.342416680262<<199.431916785557<<203.185978166691<<207.694888993961<<202.95637604834<<201.8685454529<<200.777648625028<<202.246854363589<<194.805200336976<<195.216568246443<<188.283282543232<<191.690849058651<<189.386139113436<<188.222703936485<<189.553868796024<<192.413238655967<<202.945967647128<<193.491131022231<<199.706461443975;

    MetropolisVariable densityHPD;
    densityHPD.mSupport = MetropolisVariable::eR;
    densityHPD.setName("test HPD");

    densityHPD.mRawTrace = &traceHPD;
    densityHPD.mFormatedTrace = &traceHPD;
    densityHPD.generateHistos(chainHPDs, 1024, 1.06, -1000., 2000.);
    densityHPD.generateHPD(95.);

    QList<QPair<double, QPair<double, double> > > intervals = intervalsForHpd(densityHPD.mHPD, 95.);
  //  qDebug()<<"histo2 HPD"<<intervals[0].first<<intervals[0].second.first<<intervals[0].second.second;

    QPair<double, double> answer = qMakePair<double, double> (181.253726753 , 217.061790521);

   QVERIFY(qFuzzyCompare(intervals[0].first, 94.9814736407));
   QVERIFY(myFuzzyCompare(intervals[0].second.first, answer.first));
    /* obtain with ChronoModel v1.6
     * Posterior event date. :
        MAP : 202,437023 Mean : 200,472803 Std deviation : 9,324248
        Q1 : 194,8052 Q2 (Median) : 201,782402 Q3 : 207,380478
        HPD Region (95%) : [181,253727 : 217,061791] (94,981474%) BC/AD
        Credibility Interval (95%) : [188,222704 : 212,695308] BC/AD
        Acceptation rate (all acquire iterations) : 73,5% (MH : proposal = adapt. Gaussian random walk)
      */


}

void ChronomodelTest::tempo()
{

    Event ev1;
    QVector<double> trace1;
    // i had <<0.0 for the init plot, not used in the RawTrace
    trace1<<0.0  <<0.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4<<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8<<3.12;
    //trace1<<0.0 <<6.6<<7.7<<8.8<<9.9<<9.9 <<1.1<<2.2<<3.3<<4.4<<5.5;//mNumRunIter=10

    ev1.mTheta.mRawTrace = new QVector<double>(trace1.size());
    std::copy(trace1.begin(), trace1.end(), ev1.mTheta.mRawTrace->begin());


    Event ev2;
    QVector<double> trace2;
    trace2<<0.0  <<0.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4<<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8 <<3.12;
    //trace2<<0.0<<1.1<<2.2<<3.3<<4.4<<5.5<<6.6<<7.7<<8.8<<9.9<<9.9; //mNumRunIter=10

    ev2.mTheta.mRawTrace = new QVector<double>(trace2.size());
    std::copy(trace2.begin(), trace2.end(), ev2.mTheta.mRawTrace->begin() );

    Event ev3;
    QVector<double> trace3;
    trace3 <<0.0 <<4.2<<6.1<<4.3<<8.7<<9.1<<11.6<<5.2<<1.18<<2.4<<6.8<<3.12  <<0.1<<3.8<<5.6<<9.9<<7.4<<8.2<<5.7<<8.6<<9.4; // here mNumRunIter=21

    ev3.mTheta.mRawTrace = new QVector<double>(trace3.size());
    std::copy(trace3.begin(), trace3.end(), ev3.mTheta.mRawTrace->begin() );


    Phase phase;
    phase.mEvents.append(&ev1);
    //phase.mEvents.append(&ev2);
    //phase.mEvents.append(&ev3);


    ChainSpecs chain1;
    chain1.mNumBurnIter = 0; //not used
    chain1.mBurnIterIndex = 0;

    chain1.mBatchIterIndex = 0;
    chain1.mBatchIndex = 0; //not used
    chain1.mNumBatchIter = 0; //not used

    chain1.mNumRunIter = 20; //used
    chain1.mRunIterIndex = 5;

    chain1.mThinningInterval = 1; //used

    Model model;
    model.mSettings.mTmin = 0.;
    model.mSettings.mTmax = 20.;
    model.mPhases.append(&phase);
    model.mChains.append(chain1);

 /*   model.generateTempo();

    for (auto val : model.mPhases[0]->mTempoCredibilityInf.toStdMap()) {
        qDebug()<<val.first<<" "<<val.second*chain1.mNumRunIter;
    }

    // controle Tempo
    int expected = chain1.mNumRunIter * phase.mEvents.size();
    int total = model.mPhases[0]->mTempo.last()*chain1.mNumRunIter;
qDebug()<<" Total = "<<total<<" Total expected"<<expected;
    QVERIFY(total == expected);
*/
}

Quartiles ChronomodelTest::quartilesTypeTest(const QVector<double>& trace, const int quartileType, const double p)

{
    Q_ASSERT(&trace);
    Quartiles Q;
    QVector<double> traceSorted (trace);

    QPair<int, double> parQ1 = gammaQuartile(trace, quartileType, p); // first is j and second is gamma
    QPair<int, double> parQ2 = gammaQuartile(trace, quartileType, 0.5);
    QPair<int, double> parQ3 = gammaQuartile(trace, quartileType, 1-p);

    std::sort(traceSorted.begin(), traceSorted.end());

    // Q1 determination
    if (parQ1.first<=0)
       Q.Q1 = (double)traceSorted.first();

    else if (parQ1.first < traceSorted.size())
            Q.Q1 = (1.- parQ1.second)*(double)traceSorted.at(parQ1.first-1) + parQ1.second*(double)traceSorted.at(parQ1.first);
    else
        Q.Q1 = (double)traceSorted.last();

    // Q2 determination
    if (parQ2.first<=0)
       Q.Q2 = (double)traceSorted.first();

    else if (parQ2.first < traceSorted.size())
            Q.Q2 = (1.- parQ2.second)*(double)traceSorted.at(parQ2.first-1) + parQ2.second*(double)traceSorted.at(parQ2.first);
    else
        Q.Q2 = (double)traceSorted.last();

    // Q3 determination
    if (parQ3.first<=0)
       Q.Q3 = (double)traceSorted.first();

    else if (parQ3.first < traceSorted.size())
            Q.Q3 = (1.- parQ3.second)*(double)traceSorted.at(parQ3.first-1) + parQ3.second*(double)traceSorted.at(parQ3.first);
    else
        Q.Q3 = (double)traceSorted.last();

    return Q;
}

/**
 * @brief ChronomodelTest::gammaQuartile used with quantile, find the gamma coef corresponding to the
 * type of R Calcul, to use with the general formula
 * @param trace
 * @param quartileType
 * @param p
 * @return
 */
QPair<int, double> ChronomodelTest::gammaQuartile(const QVector<double> &trace, const int quartileType, const double p)
{
    const int n (trace.size());
    int j (0);
    // We use jFloor which is the floor value of j but in the original double type
    // because when we cacul g in the 3 first cases we need the double format
    double jFloor(0.);

    double m (0.);
    double g (0.);
    double gamma (0.);
    double k (0.);


    switch (quartileType) {
    /* Case 1 to 3 are discontinuous case */

    /* It is different to R but it is identique to QuantileType1 in the article "Quantile calculations in R"
     * http://tolstoy.newcastle.edu.au/R/e17/help/att-1067/Quartiles_in_R.pdf
     */
    case 1:
        m = 0.;
        jFloor = floor((n * p) + m);
        j = (int)jFloor;
        g = n*p + m - jFloor;

        gamma = (g<1e-10 ? 0 : 1.) ;

        break;

    case 2: // same probleme as type 1
        m = 0.;
        jFloor = floor((n * p) + m);
        j = (int)jFloor;
        g = n*p + m - jFloor;
        gamma = (g==0. ? 0.5 : 1.) ;
        break;

    case 3: // OK with R
        m = -0.5;
        jFloor = floor((n * p) + m);
        j = (int)jFloor;
        g = n*p + m - jFloor;
        gamma = (g==0. && isEven(j) ? 0. : 1.);
        break;

    /* Case 4 to 9 are continuous case */
    case 4: // OK with R
        m = 0.;
        k = p * n;
        g = k - floor(k);
        j = (int) floor(k);
        gamma = g ;
        break;
    case 5: // OK with R
        m = 0.;
        k = (p * n) + 0.5;
        g = k - floor(k);
        j = (int) floor(k);
        gamma = g ;
        break;
    case 6:
        k = p * (n+1);
        g = k - floor(k);
        j = (int) floor(k);
        gamma = g;
        break;
    case 7: // OK with R, this is the default type in R software
        k = (p*(n-1) + 1);
        g = k - floor(k);
        j = (int) floor(k);
        gamma = g ;
        break;

    /* OK with R, it is the formula with Bos-Levenbach (1953) parameter
     *  http://www.barringer1.com/wa_files/The-plotting-of-observations-on-probability-paper.pdf
     */
    case 8:
        k = p * (n + 0.4) + 0.3;
        g = k - floor(k);
        j = (int) floor(k);
        gamma = g ;
        break;
    case 9:
        k = p * (n + 2./8.) + 3./8.;
        g = k - floor(k);
        j = (int) floor(k);
        gamma = g;
        break;

    default:
        gamma = 0.;
        break;
    }

    return qMakePair(j, gamma);
}

/**
 * @brief ChronomodelTest::quantile The comparaison is done with the R function quantile()
 * It's only a inner test to develop other function
 * serie <- c(0.1, 3.8, 5.6 ,9.9 ,7.4 ,8.2 ,5.7 ,8.6 ,9.4, 4.2 ,6.1, 4.3, 8.7 ,9.1 ,11.6 ,5.2 ,1.18, 2.4, 6.8, 3.12)
 * quantile(serie, type = 9)
 */
void ChronomodelTest::quantile()
{
    QVector<double> serie1;
    serie1 <<0.1<< 3.8<< 5.6 <<9.9 <<7.4 <<8.2 <<5.7 <<8.6 <<9.4<< 4.2 <<6.1<< 4.3<< 8.7 <<9.1 <<11.6 <<5.2 <<1.18<< 2.4<< 6.8<< 3.12;
  for (double i(0.); i<=100.; i=i+25.) {
        Quartiles q = quartilesTypeTest(serie1, 9, i/100. );
qDebug()<< "serie1 "<<i << q.Q1;
    }

/*     QVector<double> suite;
    suite<<0<<1<<2<<3<<4<<5<<6<<7<<8<<9<<10;

    for (double i(0.); i<=100.; i=i+10.) {
        Quartiles q = quartilesTypeTest(suite, 1, i/100. );
qDebug()<< "suite"<<i << q.Q1;
    }
*/
    QVector<double> serie2;
    serie2<< 0<<1<<2<<3<<4<<5<<6<<7<<8<<9<<10<<11<<12<<13<<14<<15<<16<<17<<18<<19<<20;
      for (double i(0); i<=100.; i=i+25.)
    {
          Quartiles q = quartilesTypeTest(serie2, 9, i/100. );
    qDebug()<< "serie2 "<<i << q.Q1;
     }

qDebug()<<"---------------------------------------";
/*     QVector<double> data12;
      data12 <<2<<3<<5<<7<<11<<13<<17<<19<<23<<29<<31<<37;
      int quantTy=2;
      Quartiles q = quartilesTypeTest(data12, quantTy, 0.25 );
      qDebug()<< "data12.Q1"<< q.Q1;
      q = quartilesTypeTest(data12, quantTy, 0.5 );
      qDebug()<< "data12.Q2"<< q.Q1;
      q = quartilesTypeTest(data12, quantTy, 0.75 );
      qDebug()<< "data12.Q3"<< q.Q1;
*/
    QVERIFY(true == true);
}


void ChronomodelTest::quartileForTrace()
{
    QVector<double> serie1;
    serie1 <<0.1<< 3.8<< 5.6 <<9.9 <<7.4 <<8.2 <<5.7 <<8.6 <<9.4<< 4.2 <<6.1<< 4.3<< 8.7 <<9.1 <<11.6 <<5.2 <<1.18<< 2.4<< 6.8<< 3.12;

    // used in MetropolisVariable::generateNumericalResults
    Quartiles Q1 = quartilesForTrace(serie1);
        qDebug()<< "serie1 0.25"<< Q1.Q1;
        qDebug()<< "serie1 0.50"<< Q1.Q2;
        qDebug()<< "serie1 0.75"<< Q1.Q3;

    // used in Model::generateTempo Type=8
    Quartiles cred = quartilesType(serie1, 8, 0.025);
    /* Result with R
     * 2.5%   50% 97.5%
     * 0.1   5.9  11.6
    */
    qDebug()<< "cred 0.025"<< cred.Q1;
    qDebug()<< "cred 0.50"<< cred.Q2;
    qDebug()<< "cred 0.975"<< cred.Q3;

    QVERIFY(cred.Q1 == 0.1);
    QVERIFY(cred.Q2 == 5.9);
    QVERIFY(cred.Q3 == 11.6);

    cred = quartilesType(serie1, 8, 0.25);
    /* Result with R
     * 25%      50%      75%
     * 3.966667 5.900000 8.658333
     */
    qDebug()<< "cred 0.25"<< cred.Q1;
    qDebug()<< "cred 0.50"<< cred.Q2;
    qDebug()<< "cred 0.75"<< cred.Q3;

    QVERIFY(cred.Q1 == 3.96);
    QVERIFY(cred.Q2 == 5.9);
    QVERIFY(cred.Q3 == 8.66);
}




QTEST_APPLESS_MAIN(ChronomodelTest)

#include "tst_chronomodeltest.moc"

#endif
back to top