/** @file
 *  @brief @f$\rho@f$ values for 2-layer neural network model.
 *
 *  $Id: nnrho.cpp 2077 2005-05-18 03:25:54Z ling $
 *
 *  I guess the @f$\rho@f$ value would be biased for complex models.
 *  That is, too many ``useless'' hypotheses in the model will make
 *  @f$\rho@f$ really weak in indicating bad examples.
 *  This is why I avoid to have a general @f$\rho@f$ code.
 */

#include <iostream>
#include <fstream>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <boost/random.hpp>
#include <lemga/quickfun.h>
#include <lemga/learnmodel.h>

#define GAUSSIAN_WEIGHT
#define USE_QUICK_TANH
#define INNER_LOOP_SIZE  1000
#define OUTPUT_PACE      10     // # of inner loops between outputs

typedef boost::mt19937 random_type;

int main (unsigned int argc, char* argv[]) {
    if (argc < 6) {
        std::cerr << argv[0] << " dat #i #h wr #net(1K)\n";
        return -1;
    }

    const UINT n1 = std::atoi(argv[2]);
    const UINT n2 = std::atoi(argv[3]);
    const REAL wr = std::atof(argv[4]);
    const UINT ns = std::atoi(argv[5]);
    const UINT wsize = n2*(n1+1) + 1*(n2+1);

    // get the data
    std::ifstream fd(argv[1]);
    if (!fd.is_open()) {
        std::cerr << "test: data file open error\n"; return -2;
    }
    lemga::pDataSet dat = lemga::load_data(fd, (1L<<30)-1, n1, 1);
    const UINT nx = dat->size();
    if (nx > UINT(-1)/2 / INNER_LOOP_SIZE / OUTPUT_PACE)
        std::cerr << "Warning: sum_eipi might overflow\n";

#ifdef USE_QUICK_TANH
    quick_tanh_setup();
    const REAL tanh_error_max = 0.0005;
#   ifndef NDEBUG
    UINT unsafe_count = 0;
#   endif
#endif

    random_type irnd;
    random_type::result_type rnd_seed = std::time(0);
    std::cerr << "seed = " << rnd_seed << std::endl;
    irnd.seed(rnd_seed);

#ifdef GAUSSIAN_WEIGHT
#   define RND_WGT(r)        (randn()*(r))
    boost::normal_distribution<random_type> randn(irnd);
#else
#   define RND_WGT(r)        (uni()*(r+r)-(r))
    boost::uniform_01<random_type> uni(irnd);
#endif

    typedef std::vector<REAL>::iterator VRI;
    std::vector<REAL> w(wsize);
    std::vector<UINT> sum_ei(nx, 0), sum_eipi(nx, 0);
    std::vector<UINT> hist(nx+1, 0);

    for (UINT i = 0; i < ns; ++i) {
    for (UINT il = 0; il < INNER_LOOP_SIZE; ++il) {

        for (VRI pw = w.begin(); pw != w.end(); ++pw)
            *pw = RND_WGT(wr);

#ifdef USE_QUICK_TANH
        REAL ow_sum = 0;  // "absolute sum" of output neuron weights
        for (UINT j = 0; j < n2; ++j)
            ow_sum += std::fabs(w[wsize-1-j]);
        const REAL ow_err = tanh_error_max * 1.1 * ow_sum;
#endif

        UINT err = 0;
        static std::vector<bool> ei(nx);
        for (UINT j = 0; j < nx; ++j) {
            /* get input & out */
            const lemga::Input& inp = dat->x(j);
            const bool out = (dat->y(j)[0] >= 0);

            static std::vector<REAL> sum1(n2);
            /* clear sum1 */
            VRI pw = w.begin();
            for (VRI ps = sum1.begin(); ps != sum1.end(); ++ps, ++pw)
                *ps = *pw;

            /* generate sum1 */
            for (lemga::Input::const_iterator px = inp.begin(); px != inp.end(); ++px) {
                for (VRI ps = sum1.begin(); ps != sum1.end(); ++ps, ++pw)
                    *ps += *pw * *px;
            }

            /* generate the output */
#ifdef USE_QUICK_TANH
            const VRI pw_backup = pw;
#endif
            REAL sum2 = *pw; ++pw;
#ifdef USE_QUICK_TANH
            for (VRI ps = sum1.begin(); ps != sum1.end(); ++ps, ++pw)
                sum2 += *pw * quick_tanh(*ps);
            const bool value_safe = (sum2 > ow_err | sum2 < -ow_err);

            if (!value_safe) // use tanh to recompute the result
#   ifndef NDEBUG
                ++unsafe_count;
            const REAL quick_sum2 = sum2;
#   endif
            {pw = pw_backup; sum2 = *pw; ++pw;
#endif // !USE_QUICK_TANH

            for (VRI ps = sum1.begin(); ps != sum1.end(); ++ps, ++pw)
                sum2 += *pw * std::tanh(*ps);
            assert(pw == w.end());

#ifdef USE_QUICK_TANH
            }
            assert(!value_safe | ((quick_sum2 >= 0) ^ (sum2 < 0)));
#endif

            ei[j] = (sum2 >= 0) ^ out;
            if (ei[j]) ++err;
            // since 0<= uni() < 1, we make the balance by using >= 0
        }

        ++hist[err];
        for (UINT j = 0; j < nx; ++j) if (ei[j]) {
            ++sum_ei[j];
            sum_eipi[j] += err; // INNER_LOOP_SIZE*avg(err) should <= 2e9
        }
        } // Inner loop

        if ((i+1) % OUTPUT_PACE == 0 || i+1 == ns) {
            for (UINT j = 0; j < nx; ++j) std::cout << sum_ei[j] << ' ';
            for (UINT j = 0; j < nx; ++j) std::cout << sum_eipi[j] << ' ';
            for (UINT j = 0; j < nx; ++j) std::cout << hist[j] << ' ';
            std::cout << hist[nx] << std::endl;
            std::fill(sum_ei.begin(), sum_ei.end(), 0);
            std::fill(sum_eipi.begin(), sum_eipi.end(), 0);
            std::fill(hist.begin(), hist.end(), 0);
        }
    }

#if defined(USE_QUICK_TANH) && !defined(NDEBUG)
    std::cerr << "unsafe_count = " << unsafe_count << std::endl;
#endif
    return 0;
}