statrand.h

Go to the documentation of this file.

#ifndef statrand_h
#define statrand_h
/*+
________________________________________________________________________

 (C) dGB Beheer B.V.; (LICENSE) http://opendtect.org/OpendTect_license.txt
 Author:        Bert Bril
 Date:          Sep 2006
 RCS:           $Id$
________________________________________________________________________

-*/

#include "algomod.h"
#include "gendefs.h"
#include "odset.h"

namespace Stats
{

mExpClass(Algo) RandomGenerator
{
public:

    virtual             ~RandomGenerator()      {}
    virtual void        init(int seed);
                                // seed != 0 gives repeatable numbers

    virtual double      get() const             = 0;    // returns in [0,1]

};


mExpClass(Algo) RandGen : public RandomGenerator
{
public:

                        RandGen();

    virtual double      get() const;
    int                 getInt() const;
    int                 getIndex(int sz) const;
    int                 getIndexFast(int sz,int seed) const;
    od_int64            getIndex(od_int64 sz) const;
    od_int64            getIndexFast(od_int64 sz,od_int64 seed) const;

    template <class T,class SzTp>
    void                subselect(T*,SzTp sz,SzTp targetsz) const;
    template <class T>
    void                subselect(T&,od_int64 targetsz) const;

private:

    int                 seed_;

};


mGlobal(Algo) RandGen randGen();


mExpClass(Algo) NormalRandGen : public RandomGenerator
{
public:

                        NormalRandGen();

    virtual double      get() const;
    float               get(float expect,float stdev) const;
    double              get(double expect,double stdev) const;

protected:

    mutable bool        useothdrawn_;
    mutable double      othdrawn_;

};


template <class T,class SzTp>
inline void Stats::RandGen::subselect( T* arr, SzTp sz, SzTp targetsz ) const
{
    for ( SzTp idx=sz-1; idx>=targetsz; idx-- )
    {
        const SzTp notselidx = getIndex( idx );
        if ( notselidx != idx )
            Swap( arr[notselidx], arr[idx] );
    }
}


template <class ODSET>
inline void Stats::RandGen::subselect( ODSET& ods, od_int64 targetsz ) const
{
    typedef typename ODSET::size_type size_type;
    const size_type sz = ods.size();
    if ( sz <= targetsz ) return;

    for ( size_type idx=sz-1; idx>=targetsz; idx-- )
    {
        const size_type notselidx = getIndex( idx );
        if ( notselidx != idx )
            ods.swap( notselidx, idx );
    }

    removeRange( ods, (size_type)targetsz, sz-1 );
}

}; // namespace Stats


#endif