discrepancy.h

/*****************************************************************************
 *   Copyright (C) 2004-2015 The PaGMO development team,                     *
 *   Advanced Concepts Team (ACT), European Space Agency (ESA)               *
 *                                                                           *
 *   https://github.com/esa/pagmo                                            *
 *                                                                           *
 *   act@esa.int                                                             *
 *                                                                           *
 *   This program 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.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program; if not, write to the                           *
 *   Free Software Foundation, Inc.,                                         *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.               *
 *****************************************************************************/

#ifndef PAGMO_UTIL_DISCREPANCY_H
#define PAGMO_UTIL_DISCREPANCY_H

#include <iostream>
#include <vector>
#include <math.h>
#include <algorithm>
#include <boost/shared_ptr.hpp>

#include "../config.h"
#include "../exceptions.h"
#include "../rng.h"

namespace pagmo{ namespace util {


namespace discrepancy {

double van_der_corput(unsigned int n, unsigned int base);
unsigned int prime ( int n );
unsigned int prime_ge ( unsigned int n );
class __PAGMO_VISIBLE project_2_simplex
{
        public:
                project_2_simplex(unsigned int dim) : m_dim(dim) {}
                std::vector<double> operator()(std::vector<double> retval) const;
        private:
                unsigned int m_dim;
};

class base;
typedef boost::shared_ptr<base> base_ptr;

//---------------------------------------------------------
//--------------------BASE CLASS---------------------------
//---------------------------------------------------------


class __PAGMO_VISIBLE base
{
public:

        base(unsigned int dim, unsigned int count = 1) : m_dim(dim), m_count(count) {}

        virtual std::vector<double> operator()() = 0;

        virtual std::vector<double> operator()(unsigned int n) = 0;
        virtual base_ptr clone() const = 0;
        virtual~base();
protected:
        unsigned int m_dim;
        unsigned int m_count;
};

//---------------------------------------------------------
//--------------------DERIVED CLASSES----------------------
//---------------------------------------------------------


class __PAGMO_VISIBLE halton : public base
{
        public:
                halton(unsigned int dim, unsigned int count = 1);
                base_ptr clone() const;
                std::vector<double> operator()();
                std::vector<double> operator()(unsigned int n);
        private:
                std::vector<unsigned int> m_primes;
};


class __PAGMO_VISIBLE faure : public base
{
        public:
        faure(unsigned int dim, unsigned int count = 1);
        base_ptr clone() const;
        std::vector<double> operator()();
        std::vector<double> operator()(unsigned int n);
        private:
                int *binomial_table ( int qs, int m, int n );
                void faure_orig ( unsigned int dim_num, unsigned int *seed, double quasi[] );
                double *faure_generate ( int dim_num, int n, int skip );
                int i4_log_i4 ( int i4, int j4 );
                int i4_min ( int i1, int i2 );
                int i4_power ( int i, int j );
        private:
                int *m_coef;
                int m_hisum_save;
                int m_qs;
                int *m_ytemp;

};


class __PAGMO_VISIBLE simplex : public base
{
public:
        simplex(unsigned int dim, unsigned int count);
        base_ptr clone() const;
        std::vector<double> operator()();
        std::vector<double> operator()(unsigned int n);
private:
        halton m_generator;
        project_2_simplex m_projector;
};



class __PAGMO_VISIBLE sobol : public base
{
        public:
                sobol(unsigned int dim, unsigned int count);
                base_ptr clone() const;
                std::vector<double> operator()();
                std::vector<double> operator()(unsigned int n);
        private:
                int i8_bit_lo0 ( long long int n );
                void i8_sobol ( unsigned int dim_num, long long int *seed, double quasi[ ] );
        private:
                unsigned int m_dim_num_save;
                bool m_initialized;
                long long int m_maxcol;
                long long int m_seed_save;
                double recipd;
                long long int lastq[1111]; //1111 is maximum dimension.
                long long int poly[1111];
                long long int v[1111][62]; //2^62 is approx. limit of points requested.
};


class __PAGMO_VISIBLE lhs : public base
{
        public:
                lhs(unsigned int dim, unsigned int count);
                base_ptr clone() const;
                std::vector<double> operator()();
                std::vector<double> operator()(unsigned int n);
        private:
                std::vector<double> latin_random ( unsigned int dim_num, unsigned int point_num);
                unsigned int *perm_uniform ( unsigned int n);
        private:
                bool m_initialised;
                std::vector<double> m_set;
                unsigned int m_next;
};

}}} //namespace discrepancy

#endif