// -*- C++ -*-

// Copyright 2006-2007 Deutsches Forschungszentrum fuer Kuenstliche Intelligenz
// or its licensors, as applicable.
//
// You may not use this file except under the terms of the accompanying license.
//
// Licensed under the Apache License, Version 2.0 (the "License"); you
// may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Project:  ocr-bpnet - neural network classifier
// File: bpnet.cc
// Purpose: neural network classifier
// Responsible: Hagen Kaprykowsky (kapry@iupr.net)
// Reviewer: Yves Rangoni (rangoni@iupr.net)
// Primary Repository:
// Web Sites: www.iupr.org, www.dfki.de

#include <float.h>
#include "bpnet.h"
#include "ocr-utils.h"
#include "narray-ops.h"
#include "logger.h"

using namespace ocropus;
using namespace colib;
using namespace narray_ops;

namespace iupr_bpnet {

#define MIN_SCORE 1e-6
#define epsilon_stdev 1e-04

    // log the confusion matrix for the training set for each epoch
    Logger logger_confusion_train("bpnet.conf_train");
    // log the confusion matrix for the testing set for each epoch
    Logger logger_confusion_test("bpnet.conf_test");

    // log error rates for each epoch
    Logger logger_errors("bpnet.error_rates");

    param_int   freq_error_rates("freq_error_rates",1,"frequency of error_rates graph plotting");
    param_float graph_scale("graph_scale",1.,"scale factor for graph ploting if error_rates is enabled");
    param_bool  merge_lr("merge_lr",true,"plot learning rate with error rates, otherwise in another plot");

    template<class T>
    bool valid(T &v) {
        for(int i=0;i<v.length1d();i++)
            if(isnan(v.at1d(i))) {
                return false;
            }
        return true;
    }

    void shuffle_feat(objlist<floatarray> &v,intarray &c) {
        floatarray v_tmp;
        v_tmp.resize(v[0].length());
        int c_tmp;

        int n = v.length();
        for(int i=0;i<n-1;i++) {
            int target = rand()%(n-i)+i;
            copy(v_tmp,v[target]);
            copy(v[target],v[i]);
            copy(v[i],v_tmp);
            c_tmp = c[target];
            c[target] = c[i];
            c[i] = c_tmp;
        }
    }

    void write(FILE *stream,floatarray &input) {
        for(int i=0;i<input.length1d();i++) {
            fprintf(stream,"%.10g\n",double(input.at1d(i)));
        }
    }
    void read(floatarray &output,FILE *stream) {
        float value;
        for(int i=0;i<output.length1d();i++) {
            if(fscanf(stream,"%g",&value)!=1) {
                throw "bad file";
                return;
            }
            output.at1d(i) = value;
        }
    }

    // not to be used outside this scope
    //#define NSIGMOID 10000
    //#define SIGMOID_RANGE 15.0
    static const int NSIGMOID=10000;
    static const float SIGMOID_RANGE=15.0;

    // table of the sigmoid function
    static floatarray sigmoid_table(NSIGMOID);

    // slow exponential sigmoid function, not tabulated
    static float slow_sigmoid(float x) {
        return 1.0/(1.0+exp(-x));
    }

    // initialize the simgoid table
    static void init_sigmoid_table() {
        for(int i=0;i<NSIGMOID;i++) {
            sigmoid_table(i) = slow_sigmoid(i*SIGMOID_RANGE/NSIGMOID);
        }
    }

    // tabulated sigmoid function
    static inline float sigmoid(float x,bool &inited) {
        float abs_x = (x<0.0)?-x:x;
        int index;
        float abs_result;

        if(!inited) {
            init_sigmoid_table();
            inited = true;
        }

        if(abs_x>=SIGMOID_RANGE) {
            abs_result = 1.0;
        } else {
            index = int(NSIGMOID*abs_x/SIGMOID_RANGE);
            abs_result = (index>=NSIGMOID) ? 1.0 : sigmoid_table(index);
        }
        return (x<0.0) ? 1.0-abs_result:abs_result;
    }

    // float random number between low and high
    static float random_range(float low,float high) {
        //drand48 is obsolete, replaced by rand():
        float rnd = float(rand())/float(RAND_MAX);
        return (rnd * (high-low) + low);
    }

    // push activations one layer up using 2d weight matrix
    void bp_propagate_activations(floatarray &activations_input,int ninput,
                                  floatarray &activations_output,int noutput,
                                  floatarray &weights,floatarray &offsets,
                                  bool &sigmoid_inited) {
        float total;
        for(int i=0;i<noutput;i++) {
            total = offsets(i);
            for(int j=0;j<ninput;j++) {
                total += weights(i,j)*activations_input(j);
            }
            activations_output(i) = sigmoid(total,sigmoid_inited);
        }
    }

    // determine hidden layer error from error at the output units.
    void bp_propagate_deltas(floatarray &deltas_input,int noutput,
                             floatarray &activations_input,
                             floatarray &delta_output,int ninput,
                             floatarray &weights) {
        for(int j=0;j<ninput;j++) {
            float deriv = activations_input(j)*(1.0-activations_input(j));
            float total = 0.0;
            if(deriv<1e-5) deriv=1e-5;
            for(int i=0;i<noutput;i++)
                total += delta_output(i)*weights(i,j);
            deltas_input(j) = deriv*total;
        }
    }

    // weight update using backpropagation formula
    void bp_update_weights(floatarray &offsets,floatarray &delta_output,
                            int noutput,floatarray &activation_input,int ninput,
                            floatarray &weights,float eta) {
        for(int i=0;i<noutput;i++) {
            for(int j=0;j<ninput;j++) {
                weights(i,j) += eta*delta_output(i)*activation_input(j);
            }
            offsets(i) += eta*delta_output(i);
        }
    }

    // scale every dimension of the input
    // down to mean=0, std_dev=1
    void normalize_input_train(objlist<floatarray> &vectors,doublearray
                               &stdev,doublearray &m_x) {
        CHECK_ARG(stdev.length()==m_x.length());
        int nsamples = vectors.length();
        CHECK_CONDITION(nsamples>0);
        int ninput = m_x.length();
        doublearray m_xx;
        m_xx.resize(ninput);
        fill(m_xx,0.0f);
        fill(m_x,0.0f);
        fill(stdev,0.0f);

        for(int d=0;d<ninput;d++) {
            // calc mean and empirical variance
            for(int n=0;n<nsamples;n++) {
                m_x(d)  += vectors[n](d);
            }
            m_x(d) /= nsamples;
            for(int n=0;n<nsamples;n++) {
                double t = vectors[n](d) - m_x(d);
                m_xx(d) += t * t;
            }
            m_xx(d) /= nsamples;
            double sqr_stdev = m_xx(d);
            if(sqr_stdev < 0.)
                sqr_stdev = 0.;
            stdev(d) = sqrt(sqr_stdev);
            // normalize
            for(int n=0;n<nsamples;n++) {
                if(stdev(d)>epsilon_stdev) {
                    vectors[n](d) = (vectors[n](d)-m_x(d))/stdev(d);
                } else {    // var = 0: all the same;
                    vectors[n](d) = vectors[n](d)-m_x(d);
                }
            }

        } // end dim loop
        ASSERT(valid(m_x));
        ASSERT(valid(stdev));
    } // end normalize_input

    void normalize_input_retrain(objlist<floatarray> &vectors,doublearray
                                 &stdev,doublearray &m_x) {
        CHECK_ARG(stdev.length()==m_x.length());
        int ninput = m_x.length();
        int nsamples = vectors.length();
        for(int d=0;d<ninput;d++) {
            for(int n=0;n<nsamples;n++) {
                if(stdev(d)>epsilon_stdev) {
                    vectors[n](d) = (vectors[n](d)-m_x(d))/stdev(d);
                } else {
                    vectors[n](d) = vectors[n](d)-m_x(d); //var=0: all the same;
                }
            }
        }
    }


#undef NSIGMOID
#undef SIGMOID_RANGE

    void BpnetClassifier::common_presets() {

        init = false;
        sigmoid_inited = false;

        norm = true;
        shuffle = true;
        ninput = -1;
        nhidden = -1;
        noutput = -1;
        epochs = -1;
        learningrate = -1.0f;
        testportion = 0.0f;

        eta_plus = 1.10;
        eta_minus = 0.6;
        eta_tol = 1.02;
        vlr = false;
    }

    BpnetClassifier::BpnetClassifier() {
        common_presets();
        filedump = false;
    }

    BpnetClassifier::BpnetClassifier(const char* path_tmp_mlp) {
        common_presets();
        buf = path_tmp_mlp;
        //strncpy(buf,path_tmp_mlp,sizeof(buf));
        filedump = true;
    }

    void BpnetClassifier::set(const char* name,double value) {
        if(!strcmp(name,"ninput")) ninput = int(value);
        else if(!strcmp(name,"nhidden")) nhidden = int(value);
        else if(!strcmp(name,"noutput")) noutput = int(value);
        else if(!strcmp(name,"epochs")) epochs = int(value);
        else if(!strcmp(name,"learningrate")) learningrate = float(value);
        else if(!strcmp(name,"testportion")) testportion = float(value);
        else if(!strcmp(name,"normalize")) norm = bool(value);
        else if(!strcmp(name,"shuffle")) shuffle = bool(value);
        else if(!strcmp(name,"filedump")) filedump = bool(value);

        else if(!strcmp(name,"vlr")) vlr = bool(value);
        else if(!strcmp(name,"eta_plus")) eta_plus = float(value);
        else if(!strcmp(name,"eta_minus")) eta_minus= float(value);
        else if(!strcmp(name,"eta_tol")) eta_tol = float(value);
        else throw "unknown parameter name";
    }

    void BpnetClassifier::add(floatarray &v,int c) {
        ASSERT(valid(v));
        if(ninput<0)
            ninput = v.dim(0);
        else
            CHECK_CONDITION(v.dim(0)==ninput);

        copy(vectors.push(),v);
        classes.push(c);
        ASSERT(vectors.length()==classes.length());
    }

    void BpnetClassifier::train() {
            if(noutput == -1)
                noutput = max(classes) + 1;
            create();
            if(norm) {
                printf("Normalizing... ");fflush(stdout);logger_errors("Normalizing");

                if(init) {
                    normalize_input_retrain(vectors,stdev,m_x);
                } else {
                    normalize_input_train(vectors,stdev,m_x);
                }
            }
            printf("end\n");
            if(shuffle) {
                printf("Shuffling... ");fflush(stdout);logger_errors("Shuffling");
                shuffle_feat(vectors,classes);
                printf("end\n");
            }
            init_backprop(0.001);
            finishTraining();
            dealloc_train();
            init = true;
    }

    void BpnetClassifier::seal() {
        vectors.dealloc();
        classes.dealloc();
    }

    void BpnetClassifier::discriminant(floatarray &result,floatarray &v) {
        if(v.length() != ninput) {
            throw_fmt("trained with input dimension %d, but got %d",
                       ninput, v.length());
        }
        result.resize(noutput);
        if(norm) {
            normalize_input_classify(v,stdev,m_x);
        }
        copy(input,v);
        forward();
        // Copy output layer to result avoiding scores less than MIN_SCORE
        for(int i=0;i<output.length();i++) {
            result(i) = max(MIN_SCORE, output(i));
        }
    }

    void BpnetClassifier::save(FILE *stream) {
        if(!stream) {
            throw "cannot open output file for bp3 for writing";
        }
        fprintf(stream,"bp3-net %d %d %d %d\n",ninput,nhidden,noutput,norm);
        if(noutput > 0) {
            if(norm) {
                // write mean and variance
                for(int d=0;d<ninput;d++) {
                    fprintf(stream,"%f %f\n",m_x(d),stdev(d));
                }
                ASSERT(valid(m_x));
                ASSERT(valid(stdev));
            }
            write(stream,weights_hidden_input);
            write(stream,hidden_offsets);
            write(stream,weights_output_hidden);
            write(stream,output_offsets);
        }
    }

    void BpnetClassifier::load(FILE *stream) {
        double m_x_tmp,stdev_tmp;
        int norm_tmp;
        if(!stream) {
            throw "bad input format";
        }
        if(fscanf(stream,"bp3-net %d %d %d %d",&ninput,&nhidden,&noutput,&norm_tmp)!=4 ||
            ninput<1||ninput>1000000||nhidden<1||nhidden>1000000||noutput<1||
            noutput>1000000) {
            throw "bad input format";
        }
        norm = bool(norm_tmp);
        create();
        if(norm) {  // read mean and variance
            for(int d=0;d<ninput;d++) {
                CHECK_CONDITION(fscanf(stream,"%lf %lf\n",&m_x_tmp,&stdev_tmp)==2);
                m_x(d) = m_x_tmp;
                stdev(d) = stdev_tmp;
            }
            ASSERT(valid(m_x));
            ASSERT(valid(stdev));
        }
        read(weights_hidden_input,stream);
        read(hidden_offsets,stream);
        read(weights_output_hidden,stream);
        read(output_offsets,stream);
        init = true;
    }

    void BpnetClassifier::create() {
        if(!init) {
            weights_hidden_input.resize(nhidden,ninput);
            weights_output_hidden.resize(noutput,nhidden);
            hidden_offsets.resize(nhidden);
            output_offsets.resize(noutput);
            input.resize(ninput);
            hidden.resize(nhidden);
            output.resize(noutput);
            confusion_train = make_ConfusionMatrix(noutput);
            confusion_test = make_ConfusionMatrix(noutput);
            best_confusion_train = make_ConfusionMatrix(noutput);
            best_confusion_test = make_ConfusionMatrix(noutput);
            m_x.resize(ninput);
            stdev.resize(ninput);
        }
    }

    void BpnetClassifier::forward_one_input(int sample_index, int &target, int &predicted) {

        copy(input,vectors[sample_index]);
        target = classes(sample_index);

        forward();
        predicted = argmax(output);
        ASSERT(predicted>=0 && predicted<noutput);
    }

    void BpnetClassifier::test_on_db(int n1, int n2, float &_error, float &cls_error,
        ConfusionMatrix& cm) {
        ASSERT(n2>n1);
        int cls, predicted;
        _error = 0.0;
        cls_error = 0.0;
        cm.clear();
        for(int sample_index=n1;sample_index<n2;sample_index++) {
            forward_one_input(sample_index, cls, predicted);
            compute_error(cls);
            _error += norm2squared(error);
            cls_error += (predicted!=cls);
            cm.increment(cls,predicted);
        }
        cls_error /= (float(n2-n1));
    }

    void logging_train_test(const char* str,float _error,float cls_error,
        ConfusionMatrix& cm,int ep,float best_cls_error,float best_error) {
        strbuf temp_string;
        if(strcmp(str,"Train")==0) {
            strbuf_format(temp_string, "Train cls error: %10.4f%%\t"
                    "Train error: %10.4f",cls_error*100.,_error);

            logger_confusion_train.format(  "Train confusion matrix"
                                            " for epoch %d",ep+1);
            strbuf conf_str;
            cm.logMatrix(conf_str);
            logger_confusion_train.html((char*)conf_str);

        } else {
            strbuf_format(temp_string, "Test cls error:     %7.4f%%\t"
                                    "Test error:  %10.4f\t"
                                    "Best test cls error: %7.4f%%\t"
                                    "Best test error: %10.4f",
                                    cls_error*100.,_error,
                                    best_cls_error*100., best_error);

            logger_confusion_test.format(   "Test confusion matrix"
                                            " for epoch %d",ep+1);
            strbuf conf_str;
            cm.logMatrix(conf_str);
            logger_confusion_test.html((char*)conf_str);

        }
        puts((char*)temp_string);
        logger_errors(temp_string);
    }

    void write_data_in_R(strbuf &s,const char* v,floatarray &d,int e) {
        s += v;
        strbuf aux;
        for(int i=0;i<e;i++) {
            strbuf_format(aux,"%f,",d[i]);  s += aux;
        }
        strbuf_format(aux,"%f);\n",d[e]);   s += aux;
    }

    void construct_R_file(strbuf& R,floatarray &h_train_error,floatarray &h_train_cls_error,
                        floatarray &h_test_error,floatarray &h_test_cls_error,
                        floatarray &h_lr,int epochs,
                        const char* log_dir,int image_counter) {
        R = "";
        strbuf aux;
        strbuf buf_PDF;

        float graph_dim_x = 14;
        float graph_dim_y = 8;

        graph_dim_x *= graph_scale;
        graph_dim_y *= graph_scale;

        write_data_in_R(R,"tr_err<-c(",h_train_error,epochs);
        write_data_in_R(R,"tr_cls_err<-c(",h_train_cls_error,epochs);
        write_data_in_R(R,"te_err<-c(",h_test_error,epochs);
        write_data_in_R(R,"te_cls_err<-c(",h_test_cls_error,epochs);
        write_data_in_R(R,"lr<-c(",h_lr,epochs);

        strbuf_format(buf_PDF,"%s/ocropus-log-%d-",log_dir,image_counter);

        strbuf_format(aux,"\npdf('%str_te_err.png',width=%g,height=%g);\n",
                    (char*)buf_PDF,graph_dim_x,graph_dim_y);    R += aux;

        if(merge_lr) {
            R += "plot.new();\n";
            strbuf_format(aux,
                        "\tplot.window(xlim=c(%d,%d),ylim=c(min(tr_err[2:%d],"
                        "te_err[2:%d]),max(tr_err[2:%d],te_err[2:%d])));\n",
                        (epochs==0)?0:1,(epochs==0)?1:epochs,epochs+1,
                        epochs+1,epochs+1,epochs+1);    R += aux;
            if(epochs==1) {
                R += "\tpoints(c(1),tr_err[2],col='red',lwd=2,lty=1);\n";
                R += "\tpoints(c(1),te_err[2],col='blue',lwd=2,lty=3);\n";
            } else {
                strbuf_format(aux,"\tlines(c(1:%d),tr_err[2:%d],col='red',"
                            "lwd=2,lty=1);\n",epochs,epochs+1);  R += aux;
                strbuf_format(aux,"\tlines(c(1:%d),te_err[2:%d],col='blue',"
                            "lwd=2,lty=3);\n",epochs,epochs+1); R += aux;
            }
            R += "\taxis(1);\n";
            R += "\taxis(2);\n";
            strbuf_format(aux,"\tlegend('topright',c('Train','Test','LR'),"
                    "lwd=2,lty=c(1,3,1),col=c('red','blue','darkgreen'));\n");
            R += aux;
            strbuf_format(aux,"\ttitle(xlab='Epochs',ylab='Mean error',"
                        "main='Train error, test error, learning rate');\n");
            R += aux;
            strbuf_format(aux,"\tplot.window(xlim=c(1,%d),c(min(lr[2:%d]),"
                        "max(lr[2:%d])));\n",epochs,epochs+1,epochs+1);
            R += aux;
            if(epochs==1) {
                R += "\tpoints(c(1),lr[2],col='darkgreen',lwd=2,lty=1);\n";
            } else {
                strbuf_format(aux,"\tlines(c(1:%d),lr[2:%d],col='darkgreen',"
                            "lwd=2,lty=1);\n",epochs,epochs+1);
                R += aux;
            }
            R += "\taxis(4,col='darkgreen',col.axis='darkgreen');\n";
            R += "\tbox();\n";
        } else {
            if(epochs==1) {
                R += "\tplot(c(1),c(tr_err[2]),col='red',lwd=2,lty=1";
            } else {
                strbuf_format(aux,"\tplot(c(1:%d),tr_err[2:%d],type='l',"
                            "col='red',lwd=2,lty=1",epochs,epochs+1);
                R += aux;
            }
            strbuf_format(aux,",\n\t\tylim=c(min(tr_err[2:%d],te_err[2:%d]),"
                        "max(tr_err[2:%d],te_err[2:%d]))",
                        epochs+1,epochs+1,epochs+1,epochs+1);
            R += aux;
            strbuf_format(aux,",\n\t\txlab='Epochs',ylab='Mean error'"
                        ",\n\t\tmain='Train and test error');\n");
            R += aux;
            if(epochs==1) {
                strbuf_format(aux,"\tpoints(c(1),c(te_err[2]),col='blue',"
                            "lwd=2,lty=3);\n",epochs,epochs+1);
                R += aux;
            } else {
                strbuf_format(aux,"\tlines(c(1:%d),te_err[2:%d],type='l',"
                            "col='blue',lwd=2,lty=3);\n",epochs,epochs+1);
                R += aux;
            }
            R += "\tlegend('topright',c('Train','Test'),lwd=2,lty=c(1,3),"
                    "col=c('red','blue'));\n";
        }
        R += "dev.off();\n\n";
        if(!merge_lr) {
            strbuf_format(buf_PDF,"%s/ocropus-log-%d-",
                        log_dir,image_counter+2);
            strbuf_format(aux,"pdf(\"%slr.png\",width=%g,height=%g);\n",
                        (char*)buf_PDF,graph_dim_x,graph_dim_y);
            R += aux;

            if(epochs==1) {
                strbuf_format(aux,"\tplot(c(1),c(lr[2]),col='darkgreen',"
                            "lwd=2,lty=1",epochs,epochs+1);    R += aux;
            } else {
                strbuf_format(aux,"\tplot(c(1:%d),lr[2:%d],type='l',"
                    "col='darkgreen',lwd=2,lty=1",epochs,epochs+1);  R += aux;
            }
            strbuf_format(aux,",\n\tylim=c(min(lr[2:%d]),max(lr[2:%d])),\n\t",
                        epochs+1,epochs+1);  R += aux;
            R += "xlab='Epochs',ylab='Value',main='Learning rate');\n";
            R += "dev.off();\n\n";
        }
        strbuf_format(buf_PDF,  "%s/ocropus-log-%d-",
                                log_dir,image_counter+1);
        strbuf_format(aux,"pdf('%str_te_cls_err.png',width=%g,height=%g);\n",
                    (char*)buf_PDF,graph_dim_x,graph_dim_y);    R += aux;
        if(epochs==1) {
            R += "\tplot(c(1),c(tr_cls_err[2]),col='red',lwd=2,lty=1";
        } else {
            strbuf_format(aux,"\tplot(c(1:%d),tr_cls_err[2:%d],"
                "type='l',col='red',lwd=2,lty=1",epochs,epochs+1);  R += aux;
        }
        strbuf_format(aux,",\n\t\tylim=c(min(tr_cls_err[2:%d],te_cls_err[2:%d]),"
                    "max(tr_cls_err[2:%d],te_cls_err[2:%d]))",
                    epochs+1,epochs+1,epochs+1,epochs+1);   R += aux;
        R += ",\n\t\txlab=\"Epochs\",ylab='Classification %% error',"
                    "\n\t\tmain='Train and test classification error');\n";
        if(epochs==1) {
            strbuf_format(aux,"\tpoints(c(1),c(te_cls_err[2]),col='blue',"
                        "lwd=2,lty=3);\n",epochs,epochs+1); R += aux;
        } else {
            strbuf_format(aux,"\tlines(c(1:%d),te_cls_err[2:%d],col='blue',"
                        "lwd=2,lty=3);\n",epochs,epochs+1); R += aux;
        }
        R += "\tlegend('topright',c('Train','Test'),lwd=2,"
                    "lty=c(1,3),col=c('red','blue'));\n";
        R += "dev.off();\n\n";
    }

    /// Plot learningrate, training, testing errors and classification errors
    /// up to the given number of epochs.
    /// It requires R, otherwise nothing is added to the logger
    void train_test_curves(
                    floatarray &h_train_error,floatarray &h_train_cls_error,
                    floatarray &h_test_error,floatarray &h_test_cls_error,
                    floatarray &h_lr,int epochs) {
        if(!logger_errors.enabled) return;
        if(epochs == 0) return;
        if(system("R --version > /dev/null") != 0) return;
        if(system("convert > /dev/null") != 0) return;

        const char* log_dir = get_logger_directory();
        int image_counter = get_image_counter();

        strbuf aux;
        strbuf buf_R_filename;
        strbuf R;
        strbuf buf_PDF_filename;

        construct_R_file(R,h_train_error,h_train_cls_error,
                    h_test_error,h_test_cls_error,
                    h_lr,epochs,log_dir,image_counter);

        strbuf_format(buf_R_filename,"%s/curves-%d-%d.R",log_dir,
                        epochs,image_counter);

        stdio R_file = stdio(buf_R_filename,"wt");
        fputs(R,R_file);

        logger_errors.putIndent();
        if(merge_lr) {
            logger_errors.html("Train and test error, learning rate<BR>");
        } else {
            logger_errors.html("Train, test errors<BR>");
        }

        strbuf_format(aux,"<IMG SRC=\"ocropus-log-%d-tr_te_err.png\">\n"
                            "<BR>\n",image_counter);
        logger_errors.html(aux);

        if(!merge_lr) {
            logger_errors.putIndent();
            logger_errors.html("Learning rate<BR>");
            strbuf_format(aux,"<IMG SRC=\"ocropus-log-%d-lr.png\">\n<BR>\n",
                                    image_counter+2);
            logger_errors.html(aux);
        }

        logger_errors.putIndent();
        logger_errors.html("Train, test classification error<BR>");
        strbuf_format(aux,"<IMG SRC=\"ocropus-log-%d-tr_te_cls_err.png\">"
                                "\n<BR>\n",image_counter+1);
        logger_errors.html(aux);


        strbuf cmd;

        strbuf_format(cmd, "(nice -n19 R --quiet --slave CMD BATCH %s ; ",
                            (char*)buf_R_filename);
        strbuf_format(  buf_PDF_filename,"%s/ocropus-log-%d-tr_te_err.png",
                        log_dir,image_counter);
        strbuf_format(  aux,"nice -n19 convert -antialias %s png:%s ; ",
                        (char*)buf_PDF_filename,(char*)buf_PDF_filename);
        cmd += aux;

        if(!merge_lr) {
            strbuf_format(buf_PDF_filename,"%s/ocropus-log-%d-lr.png",
                                    log_dir,image_counter+2);
            strbuf_format(aux,  "nice -n19 convert -antialias %s png:%s ; ",
                                (char*)buf_PDF_filename,(char*)buf_PDF_filename);
            cmd += aux;
        }

        strbuf_format(buf_PDF_filename,"%s/ocropus-log-%d-tr_te_cls_err.png",
                                log_dir,image_counter+1);
        strbuf_format(aux,  "nice -n19 convert -antialias %s png:%s) &",
                            (char*)buf_PDF_filename,(char*)buf_PDF_filename);
        cmd += aux;
        //printf("%s\n",(char*)cmd);
        system(cmd);

        set_image_counter(image_counter+2+(merge_lr?0:1));
    }

    void BpnetClassifier::finishTraining() {
        CHECK_CONDITION(learningrate>0.0f);
        ASSERT(testportion>=0.0f&&testportion<=1.0f);

        floatarray histo_train_error;histo_train_error.resize(epochs+1);
        floatarray histo_train_cls_error;histo_train_cls_error.resize(epochs+1);
        floatarray histo_test_error;histo_test_error.resize(epochs+1);
        floatarray histo_test_cls_error;histo_test_cls_error.resize(epochs+1);
        floatarray histo_lr;histo_lr.resize(epochs+1);

        strbuf temp_string;
        int predicted = -1;
        int cls = -1;
        float train_cls_error;
        //float best_train_cls_error = FLT_MAX;
        float train_error;
        float best_train_error = FLT_MAX;

        float test_cls_error = 0.0f;
        float best_test_cls_error = FLT_MAX;
        float test_error = 0.0f;
        float best_test_error = FLT_MAX;

        int n = vectors.length();
        int ntrain = int(float(n)*(1.0f-testportion));
        int ntest = n-ntrain;

        strbuf_format(temp_string, "ep:%d nh:%d lr:%g tp:%g np:%d i:%d o:%d\n",
                epochs,nhidden,learningrate,testportion,n,ninput,noutput);
        puts((char*)temp_string);
        logger_errors(temp_string);

        strbuf_format(temp_string, "=== Start training on %d samples "
                "(testing on %d) for %d epochs ===\n",ntrain,ntest,epochs);
        puts((char*)temp_string);
        logger_errors(temp_string);

        //local copy of the current mlp
        floatarray weights_hidden_input_temp;
        floatarray hidden_offsets_temp;
        floatarray weights_output_hidden_temp;
        floatarray output_offsets_temp;

        float old_train_error = 0.0f;
        strbuf_format(temp_string, "Epoch: %d\tLR: %f\n",0,learningrate);
        puts((char*)temp_string);
        logger_errors(temp_string);

        test_on_db(0,ntrain,train_error,train_cls_error,*confusion_train);
        old_train_error = best_train_error = train_error;
        logging_train_test( "Train",train_error,train_cls_error,
                            *confusion_train,-1,0,0);

        test_on_db(ntrain,n,test_error,test_cls_error,*confusion_test);
        best_test_error = test_error;
        best_test_cls_error = test_cls_error;
        logging_train_test( "Test ",test_error,test_cls_error,
                            *confusion_test,-1,
                            best_test_cls_error,best_test_error);

        histo_train_error[0] = train_error/ntrain;
        histo_train_cls_error[0] = train_cls_error*100.;
        histo_test_error[0] = test_error/ntest;
        histo_test_cls_error[0] = test_cls_error*100.;
        histo_lr[0] = learningrate;

        for(int epoch=0;epoch<epochs;epoch++) {
            strbuf_format(temp_string, "Epoch: %d\tLR: %f",epoch+1,learningrate);
            puts((char*)temp_string);
            logger_errors(temp_string);
            //confusion_test->printReduced(stdout);
            //keep current network in memory
            copy(weights_hidden_input_temp, weights_hidden_input);
            copy(hidden_offsets_temp,       hidden_offsets);
            copy(weights_output_hidden_temp,weights_output_hidden);
            copy(output_offsets_temp,       output_offsets);

            //train
            for(int sample_index=0;sample_index<ntrain;sample_index++) {
                forward_one_input(sample_index, cls, predicted);
                compute_error(cls);
                backward();
                update();
                if(sample_index % 1000 == 0) {
                    printf("%d/%d\r", sample_index, ntrain);fflush(stdout);
                }
            }
            printf("\n");

            test_on_db(0,ntrain,train_error,train_cls_error,*confusion_train);
            logging_train_test( "Train",train_error,train_cls_error,
                                *confusion_train,epoch,0,0);

            //test
            if(testportion>0.0f) {
                test_on_db(ntrain,n,test_error,test_cls_error,*confusion_test);

                if(test_cls_error < best_test_cls_error) {
                    best_test_cls_error = test_cls_error;
                }
                if(test_error < best_test_error) {
                    if(filedump) {
                        copy(weights_hidden_input_best,weights_hidden_input);
                        copy(hidden_offsets_best,hidden_offsets);
                        copy(weights_output_hidden_best,weights_output_hidden);
                        copy(output_offsets_best,output_offsets);
                        save(stdio((char*)buf,"w"));
                    }
                    best_test_error = test_error;
                    best_confusion_train->set(*confusion_train);
                    best_confusion_test->set(*confusion_test);
                }
                logging_train_test( "Test",test_error,test_cls_error,
                                    *confusion_test,epoch,
                                    best_test_cls_error,best_test_error);
                histo_test_error[epoch+1] = test_error/ntest;
                histo_test_cls_error[epoch+1] = test_cls_error*100;
            } else {
                if(train_error < best_train_error) {
                    if(filedump) {
                        copy(weights_hidden_input_best,weights_hidden_input);
                        copy(hidden_offsets_best,hidden_offsets);
                        copy(weights_output_hidden_best,weights_output_hidden);
                        copy(output_offsets_best,output_offsets);
                        save(stdio((char*)buf,"w"));
                    }
                    best_train_error = train_error;
                    best_confusion_train->set(*confusion_train);
                }
            }
            histo_train_error[epoch+1] = train_error/ntrain;
            histo_train_cls_error[epoch+1] = train_cls_error*100.;
            histo_lr[epoch+1] = learningrate;

            if(vlr) {
                if(train_error / old_train_error >= eta_tol) {
                    learningrate *= eta_minus;
                    printf("eta -\n");
                    logger_errors("eta -");
                    copy(weights_hidden_input,  weights_hidden_input_temp);
                    copy(hidden_offsets,        hidden_offsets_temp);
                    copy(weights_output_hidden, weights_output_hidden_temp);
                    copy(output_offsets,        output_offsets_temp);
                } else {
                    if(train_error < old_train_error) {
                        learningrate *= eta_plus;
                        printf("eta +\n");
                        logger_errors("eta +");
                    } else {
                        printf("eta =\n");
                        logger_errors("eta =");
                    }
                    old_train_error = train_error;
                }
            }

            if((freq_error_rates>0) && ((epoch%freq_error_rates)==0 || (epoch==0)))
                train_test_curves(
                    histo_train_error,histo_train_cls_error,
                    histo_test_error,histo_test_cls_error,
                    histo_lr,epoch+1);
            //confusion_test->printReduced(stdout);
        }
        if(filedump) {  // get the best weights
            copy(weights_hidden_input,weights_hidden_input_best);
            copy(hidden_offsets,hidden_offsets_best);
            copy(weights_output_hidden,weights_output_hidden_best);
            copy(output_offsets,output_offsets_best);
        }
        fflush(stdout);

        train_test_curves(
                histo_train_error,histo_train_cls_error,
                histo_test_error,histo_test_cls_error,
                histo_lr,epochs);
    }

    void BpnetClassifier::init_backprop(float range) {
        error.resize(noutput);
        hidden_deltas.resize(nhidden);
        output_deltas.resize(noutput);
        fill(error,0.0f);
        fill(hidden_deltas,0.0f);
        fill(output_deltas,0.0f);
        if(!init) {
            for(int i=0;i<nhidden;i++) {
                for(int j=0;j<ninput;j++) {
                    weights_hidden_input(i,j) = random_range(-range,range);
                }
            }
            for(int i=0;i<nhidden;i++) {
                hidden_offsets(i) = random_range(-range,range);
            }
            for(int i=0;i<noutput;i++) {
                for(int j=0;j<nhidden;j++) {
                    weights_output_hidden(i,j) = random_range(-range,range);
                }
            }
            for(int i=0;i<noutput;i++) {
                output_offsets(i) = random_range(-range,range);
            }
        }
    }

    void BpnetClassifier::dealloc_train() {
        error.dealloc();
        hidden_deltas.dealloc();
        output_deltas.dealloc();
        classes.clear();
        vectors.clear();
    }

    void BpnetClassifier::forward() {
        bp_propagate_activations(input,ninput,hidden,nhidden,
                                 weights_hidden_input,hidden_offsets,
                                 sigmoid_inited);
        bp_propagate_activations(hidden,nhidden,output,noutput,
                                 weights_output_hidden,output_offsets,
                                 sigmoid_inited);
    }

    void BpnetClassifier::compute_error(int cls) {
        copy(error,output);
        error(cls) -= 1;
    }

    void BpnetClassifier::backward() {
        int i;
        for(i=0;i<noutput;i++) {
            float deriv = output(i)*(1.0-output(i));
            output_deltas(i) = -deriv*error(i);
        }
        bp_propagate_deltas(hidden_deltas,noutput,
                            hidden,output_deltas,nhidden,
                            weights_output_hidden);
    }

    void BpnetClassifier::update() {
        bp_update_weights(output_offsets,output_deltas,noutput,
                          hidden,nhidden,weights_output_hidden,
                          learningrate);
        bp_update_weights(hidden_offsets,hidden_deltas,nhidden,
                          input,ninput,weights_hidden_input,
                          learningrate);
    }

}

namespace ocropus {
    IClassifier *make_BpnetClassifier() {
        using namespace iupr_bpnet;
        return new BpnetClassifier();
    }
    IClassifier *make_BpnetClassifierDumpIntoFile(const char *path) {
        using namespace iupr_bpnet;
        return new BpnetClassifier(path);
    }
};