/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//            Intel License Agreement
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "precomp.hpp"

namespace cv {
namespace ml {

NormalBayesClassifier::Params::Params() {}

class NormalBayesClassifierImpl : public NormalBayesClassifier
{
public:
    NormalBayesClassifierImpl()
    {
        nallvars = 0;
    }

    void setParams(const Params&) {}
    Params getParams() const { return Params(); }

    bool train( const Ptr<TrainData>& trainData, int flags )
    {
        const float min_variation = FLT_EPSILON;
        Mat responses = trainData->getNormCatResponses();
        Mat __cls_labels = trainData->getClassLabels();
        Mat __var_idx = trainData->getVarIdx();
        Mat samples = trainData->getTrainSamples();
        int nclasses = (int)__cls_labels.total();

        int nvars = trainData->getNVars();
        int s, c1, c2, cls;

        int __nallvars = trainData->getNAllVars();
        bool update = (flags & UPDATE_MODEL) != 0;

        if( !update )
        {
            nallvars = __nallvars;
            count.resize(nclasses);
            sum.resize(nclasses);
            productsum.resize(nclasses);
            avg.resize(nclasses);
            inv_eigen_values.resize(nclasses);
            cov_rotate_mats.resize(nclasses);

            for( cls = 0; cls < nclasses; cls++ )
            {
                count[cls]            = Mat::zeros( 1, nvars, CV_32SC1 );
                sum[cls]              = Mat::zeros( 1, nvars, CV_64FC1 );
                productsum[cls]       = Mat::zeros( nvars, nvars, CV_64FC1 );
                avg[cls]              = Mat::zeros( 1, nvars, CV_64FC1 );
                inv_eigen_values[cls] = Mat::zeros( 1, nvars, CV_64FC1 );
                cov_rotate_mats[cls]  = Mat::zeros( nvars, nvars, CV_64FC1 );
            }

            var_idx = __var_idx;
            cls_labels = __cls_labels;

            c.create(1, nclasses, CV_64FC1);
        }
        else
        {
            // check that the new training data has the same dimensionality etc.
            if( nallvars != __nallvars ||
                var_idx.size() != __var_idx.size() ||
                norm(var_idx, __var_idx, NORM_INF) != 0 ||
                cls_labels.size() != __cls_labels.size() ||
                norm(cls_labels, __cls_labels, NORM_INF) != 0 )
                CV_Error( CV_StsBadArg,
                "The new training data is inconsistent with the original training data; varIdx and the class labels should be the same" );
        }

        Mat cov( nvars, nvars, CV_64FC1 );
        int nsamples = samples.rows;

        // process train data (count, sum , productsum)
        for( s = 0; s < nsamples; s++ )
        {
            cls = responses.at<int>(s);
            int* count_data = count[cls].ptr<int>();
            double* sum_data = sum[cls].ptr<double>();
            double* prod_data = productsum[cls].ptr<double>();
            const float* train_vec = samples.ptr<float>(s);

            for( c1 = 0; c1 < nvars; c1++, prod_data += nvars )
            {
                double val1 = train_vec[c1];
                sum_data[c1] += val1;
                count_data[c1]++;
                for( c2 = c1; c2 < nvars; c2++ )
                    prod_data[c2] += train_vec[c2]*val1;
            }
        }

        Mat vt;

        // calculate avg, covariance matrix, c
        for( cls = 0; cls < nclasses; cls++ )
        {
            double det = 1;
            int i, j;
            Mat& w = inv_eigen_values[cls];
            int* count_data = count[cls].ptr<int>();
            double* avg_data = avg[cls].ptr<double>();
            double* sum1 = sum[cls].ptr<double>();

            completeSymm(productsum[cls], 0);

            for( j = 0; j < nvars; j++ )
            {
                int n = count_data[j];
                avg_data[j] = n ? sum1[j] / n : 0.;
            }

            count_data = count[cls].ptr<int>();
            avg_data = avg[cls].ptr<double>();
            sum1 = sum[cls].ptr<double>();

            for( i = 0; i < nvars; i++ )
            {
                double* avg2_data = avg[cls].ptr<double>();
                double* sum2 = sum[cls].ptr<double>();
                double* prod_data = productsum[cls].ptr<double>(i);
                double* cov_data = cov.ptr<double>(i);
                double s1val = sum1[i];
                double avg1 = avg_data[i];
                int _count = count_data[i];

                for( j = 0; j <= i; j++ )
                {
                    double avg2 = avg2_data[j];
                    double cov_val = prod_data[j] - avg1 * sum2[j] - avg2 * s1val + avg1 * avg2 * _count;
                    cov_val = (_count > 1) ? cov_val / (_count - 1) : cov_val;
                    cov_data[j] = cov_val;
                }
            }

            completeSymm( cov, 1 );

            SVD::compute(cov, w, cov_rotate_mats[cls], noArray());
            transpose(cov_rotate_mats[cls], cov_rotate_mats[cls]);
            cv::max(w, min_variation, w);
            for( j = 0; j < nvars; j++ )
                det *= w.at<double>(j);

            divide(1., w, w);
            c.at<double>(cls) = det > 0 ? log(det) : -700;
        }

        return true;
    }

    class NBPredictBody : public ParallelLoopBody
    {
    public:
        NBPredictBody( const Mat& _c, const vector<Mat>& _cov_rotate_mats,
                       const vector<Mat>& _inv_eigen_values,
                       const vector<Mat>& _avg,
                       const Mat& _samples, const Mat& _vidx, const Mat& _cls_labels,
                       Mat& _results, Mat& _results_prob, bool _rawOutput )
        {
            c = &_c;
            cov_rotate_mats = &_cov_rotate_mats;
            inv_eigen_values = &_inv_eigen_values;
            avg = &_avg;
            samples = &_samples;
            vidx = &_vidx;
            cls_labels = &_cls_labels;
            results = &_results;
            results_prob = _results_prob.data ? &_results_prob : 0;
            rawOutput = _rawOutput;
        }

        const Mat* c;
        const vector<Mat>* cov_rotate_mats;
        const vector<Mat>* inv_eigen_values;
        const vector<Mat>* avg;
        const Mat* samples;
        const Mat* vidx;
        const Mat* cls_labels;

        Mat* results_prob;
        Mat* results;
        float* value;
        bool rawOutput;

        void operator()( const Range& range ) const
        {
            int cls = -1;
            int rtype = 0, rptype = 0;
            size_t rstep = 0, rpstep = 0;
            int nclasses = (int)cls_labels->total();
            int nvars = avg->at(0).cols;
            double probability = 0;
            const int* vptr = vidx && !vidx->empty() ? vidx->ptr<int>() : 0;

            if (results)
            {
                rtype = results->type();
                rstep = results->isContinuous() ? 1 : results->step/results->elemSize();
            }
            if (results_prob)
            {
                rptype = results_prob->type();
                rpstep = results_prob->isContinuous() ? 1 : results_prob->step/results_prob->elemSize();
            }
            // allocate memory and initializing headers for calculating
            cv::AutoBuffer<double> _buffer(nvars*2);
            double* _diffin = _buffer;
            double* _diffout = _buffer + nvars;
            Mat diffin( 1, nvars, CV_64FC1, _diffin );
            Mat diffout( 1, nvars, CV_64FC1, _diffout );

            for(int k = range.start; k < range.end; k++ )
            {
                double opt = FLT_MAX;

                for(int i = 0; i < nclasses; i++ )
                {
                    double cur = c->at<double>(i);
                    const Mat& u = cov_rotate_mats->at(i);
                    const Mat& w = inv_eigen_values->at(i);

                    const double* avg_data = avg->at(i).ptr<double>();
                    const float* x = samples->ptr<float>(k);

                    // cov = u w u'  -->  cov^(-1) = u w^(-1) u'
                    for(int j = 0; j < nvars; j++ )
                        _diffin[j] = avg_data[j] - x[vptr ? vptr[j] : j];

                    gemm( diffin, u, 1, noArray(), 0, diffout, GEMM_2_T );
                    for(int j = 0; j < nvars; j++ )
                    {
                        double d = _diffout[j];
                        cur += d*d*w.ptr<double>()[j];
                    }

                    if( cur < opt )
                    {
                        cls = i;
                        opt = cur;
                    }
                    probability = exp( -0.5 * cur );

                    if( results_prob )
                    {
                        if ( rptype == CV_32FC1 )
                            results_prob->ptr<float>()[k*rpstep + i] = (float)probability;
                        else
                            results_prob->ptr<double>()[k*rpstep + i] = probability;
                    }
                }

                int ival = rawOutput ? cls : cls_labels->at<int>(cls);
                if( results )
                {
                    if( rtype == CV_32SC1 )
                        results->ptr<int>()[k*rstep] = ival;
                    else
                        results->ptr<float>()[k*rstep] = (float)ival;
                }
            }
        }
    };

    float predict( InputArray _samples, OutputArray _results, int flags ) const
    {
        return predictProb(_samples, _results, noArray(), flags);
    }

    float predictProb( InputArray _samples, OutputArray _results, OutputArray _resultsProb, int flags ) const
    {
        int value=0;
        Mat samples = _samples.getMat(), results, resultsProb;
        int nsamples = samples.rows, nclasses = (int)cls_labels.total();
        bool rawOutput = (flags & RAW_OUTPUT) != 0;

        if( samples.type() != CV_32F || samples.cols != nallvars )
            CV_Error( CV_StsBadArg,
                     "The input samples must be 32f matrix with the number of columns = nallvars" );

        if( samples.rows > 1 && _results.needed() )
            CV_Error( CV_StsNullPtr,
                     "When the number of input samples is >1, the output vector of results must be passed" );

        if( _results.needed() )
        {
            _results.create(nsamples, 1, CV_32S);
            results = _results.getMat();
        }
        else
            results = Mat(1, 1, CV_32S, &value);

        if( _resultsProb.needed() )
        {
            _resultsProb.create(nsamples, nclasses, CV_32F);
            resultsProb = _resultsProb.getMat();
        }

        cv::parallel_for_(cv::Range(0, nsamples),
                          NBPredictBody(c, cov_rotate_mats, inv_eigen_values, avg, samples,
                                       var_idx, cls_labels, results, resultsProb, rawOutput));
        
        return (float)value;
    }

    void write( FileStorage& fs ) const
    {
        int nclasses = (int)cls_labels.total(), i;

        fs << "var_count" << (var_idx.empty() ? nallvars : (int)var_idx.total());
        fs << "var_all" << nallvars;

        if( !var_idx.empty() )
            fs << "var_idx" << var_idx;
        fs << "cls_labels" << cls_labels;

        fs << "count" << "[";
        for( i = 0; i < nclasses; i++ )
            fs << count[i];

        fs << "]" << "sum" << "[";
        for( i = 0; i < nclasses; i++ )
            fs << sum[i];

        fs << "]" << "productsum" << "[";
        for( i = 0; i < nclasses; i++ )
            fs << productsum[i];

        fs << "]" << "avg" << "[";
        for( i = 0; i < nclasses; i++ )
            fs << avg[i];

        fs << "]" << "inv_eigen_values" << "[";
        for( i = 0; i < nclasses; i++ )
            fs << inv_eigen_values[i];

        fs << "]" << "cov_rotate_mats" << "[";
        for( i = 0; i < nclasses; i++ )
            fs << cov_rotate_mats[i];

        fs << "]";

        fs << "c" << c;
    }

    void read( const FileNode& fn )
    {
        clear();

        fn["var_all"] >> nallvars;

        if( nallvars <= 0 )
            CV_Error( CV_StsParseError,
                     "The field \"var_count\" of NBayes classifier is missing or non-positive" );

        fn["var_idx"] >> var_idx;
        fn["cls_labels"] >> cls_labels;

        int nclasses = (int)cls_labels.total(), i;

        if( cls_labels.empty() || nclasses < 1 )
            CV_Error( CV_StsParseError, "No or invalid \"cls_labels\" in NBayes classifier" );

        FileNodeIterator
            count_it = fn["count"].begin(),
            sum_it = fn["sum"].begin(),
            productsum_it = fn["productsum"].begin(),
            avg_it = fn["avg"].begin(),
            inv_eigen_values_it = fn["inv_eigen_values"].begin(),
            cov_rotate_mats_it = fn["cov_rotate_mats"].begin();

        count.resize(nclasses);
        sum.resize(nclasses);
        productsum.resize(nclasses);
        avg.resize(nclasses);
        inv_eigen_values.resize(nclasses);
        cov_rotate_mats.resize(nclasses);

        for( i = 0; i < nclasses; i++, ++count_it, ++sum_it, ++productsum_it, ++avg_it,
                                    ++inv_eigen_values_it, ++cov_rotate_mats_it )
        {
            *count_it >> count[i];
            *sum_it >> sum[i];
            *productsum_it >> productsum[i];
            *avg_it >> avg[i];
            *inv_eigen_values_it >> inv_eigen_values[i];
            *cov_rotate_mats_it >> cov_rotate_mats[i];
        }

        fn["c"] >> c;
    }

    void clear()
    {
        count.clear();
        sum.clear();
        productsum.clear();
        avg.clear();
        inv_eigen_values.clear();
        cov_rotate_mats.clear();

        var_idx.release();
        cls_labels.release();
        c.release();
        nallvars = 0;
    }

    bool isTrained() const { return !avg.empty(); }
    bool isClassifier() const { return true; }
    int getVarCount() const { return nallvars; }
    String getDefaultModelName() const { return "opencv_ml_nbayes"; }

    int nallvars;
    Mat var_idx, cls_labels, c;
    vector<Mat> count, sum, productsum, avg, inv_eigen_values, cov_rotate_mats;
};


Ptr<NormalBayesClassifier> NormalBayesClassifier::create(const Params&)
{
    Ptr<NormalBayesClassifierImpl> p = makePtr<NormalBayesClassifierImpl>();
    return p;
}

}
}

/* End of file. */