//  Copyright (2013) Cédric Coussinet (cedric.coussinet@nomoseed.net)
//
//  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 3 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, see <http://www.gnu.org/licenses/>

#include <QDir>
#include <QFile>

#include <libxml/parser.h>
#include <libxml/xpath.h>
#include <libxml/xpathInternals.h>

#include <gsl/gsl_sf_lambert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort.h>

#include "discretization.h"

QStringList extractGaussiankernel (const QString scheme, QString intervals, const int n){
    const QStringList temp = intervals.simplified().split(' ');
    QStringList intervalLabel;
    const int borderSize = temp.length()/2;
    double *border = new double[borderSize];
    for (int i=0;i<temp.length();i++){
        if (i%2)
            border[i/2]= temp[i].toDouble();
        else
            intervalLabel.append(temp[i]);
    }
    for (int i=1;i<borderSize;i++){
        if (border[i-1]>= border[i])
            return QStringList("#The order of intervals must be strictly increasing.");
    }
    QStringList labelTest = intervalLabel;
    for (int i=labelTest.length()-1;i>=0;i--)
        if (labelTest.count(labelTest[i])==1)
            labelTest.removeAt(i);
    for (int i=1;i<labelTest.length();i++)
        if (labelTest[i-1]!=labelTest[i])
            if (labelTest.indexOf(labelTest[i-1],i)!=-1)
                return QStringList("#Two items are alternated twice.");
    QStringList label = intervalLabel;
    label.removeDuplicates();
    const int labelSize = label.length();
    double *mean = new double[labelSize];
    double *sd = new double[labelSize];
    double *dx = new double[labelSize];
    double *yl = new double[labelSize];
    double *yr = new double[labelSize];
    double *xl = new double[labelSize];
    double *xr = new double[labelSize];
    for (int i=0;i<labelSize;i++){
        mean[i] = 0.0;
        sd[i] = 0.0;
        dx[i] = 0.0;
        yl[i] = 0.0;
        yr[i] = 0.0;
        xl[i] = 0.0;
        xr[i] = 0.0;
    }
    for (int i=0;i<labelSize;i++){
        if (intervalLabel.count(label[i]) == 1){
            if (i == 0){
                xl[i] = border[0]*2.0-border[1];
                xr[i] = border[0];
            }
            else if (intervalLabel.last() == label[i]){
                xl[i] = border[borderSize-1];
                xr[i] = border[borderSize-1]*2.0 - border[borderSize-2];
            }
            else{
                xl[i] = border[intervalLabel.indexOf(label[i])-1];
                xr[i] = border[intervalLabel.indexOf(label[i])];
            }
            dx [i] = xr[i] - xl[i];
            mean[i] = (xr[i] + xl[i])/2.0;
        }
        else{
            if (intervalLabel.first() == label[i])
                xl[i] = border[0]+(border[0]-border[1]);
            else
                xl[i] = border[intervalLabel.indexOf(label[i])-1];
            if (intervalLabel.last() == label[i])
                xr[i] = border[borderSize-1] + (border[borderSize-1] - border[borderSize-2]);
            else
                xr[i] = border[intervalLabel.lastIndexOf(label[i])];
                dx [i] = xr[i] - xl[i];
                mean[i] = (xr[i] + xl[i])/2.0;
        }
    }
    double dxmax = 0;
    for (int i=0;i<labelSize;i++){
        if (dxmax < dx[i])
            dxmax = dx[i];
    }
    sd[0] = dxmax/2.0;
    int secondIndex = 0;
    int beforeLastIndex = 0;
    const int lastIndex = label.indexOf(intervalLabel.last());
    for (int i=1;i<labelSize;i++){
        if (i != lastIndex)
            if (dxmax == dx[i]){
                sd[i] = dxmax/2.0;
                if (secondIndex == 0)
                    secondIndex = i;
                beforeLastIndex = i;
            }
    }
    sd[lastIndex] = sd[0];
    const double SDLIMIT_2 = -1/(2*sd[0]*sd[0]);
    const double SQRT2PI = sqrt(2*M_PI);
    const double SDLIMIT_SQRT2PI = 1 /(SQRT2PI*sd[0]);
    if (secondIndex != 0){
        int i = secondIndex;
        yl[i] = exp((mean[i]-xl[i])*(mean[i]-xl[i])*SDLIMIT_2)*SDLIMIT_SQRT2PI;
        mean[0] =  xr[0] - M_SQRT2*sd[0]*sqrt(log(SDLIMIT_SQRT2PI/yl[i]));
        i = beforeLastIndex;
        yr[i] = exp((mean[i]-xr[i])*(mean[i]-xr[i])*SDLIMIT_2)*SDLIMIT_SQRT2PI;
        mean[lastIndex] = xl[lastIndex] + M_SQRT2*sd[0]*sqrt(log(SDLIMIT_SQRT2PI/yr[i]));
    }
    double y;

    for (int k=0;k<labelSize;k++){
        for (int i=0;i<labelSize;i++){
            if (sd[i] != 0){
                y = exp((mean[i]-xr[k])*(mean[i]-xr[k])*SDLIMIT_2)*SDLIMIT_SQRT2PI;
                if (yr[k] < y)
                    yr[k] = y;
                y = exp((mean[i]-xl[k])*(mean[i]-xl[k])*SDLIMIT_2)*SDLIMIT_SQRT2PI;
                if (yl[k] < y)
                    yl[k] = y;
            }
        }
    }
    size_t *list = new size_t[labelSize];
    gsl_sort_index(list, dx, 1, labelSize);
    for (int j=labelSize-1;j>=0;j--){
        const int i = list[j];
        if (sd[i] == 0){
            sd[i] = dx[i]/2.0;
            mean[i] = xl[i] + M_SQRT2*sd[i]*sqrt(log(1/(SQRT2PI*sd[i]*yl[i])));
            double p = exp(-(mean[i]-xr[i])*(mean[i]-xr[i])/(2*sd[i]*sd[i]))/(SQRT2PI*sd [i]);
            int k = 1 ;
            double err = 1.0;
            double dsd = sd[i];
            double delta;
            while (k < 1000 && err > 0.000000000001){
                err = abs(yr[i]-p);
                k++;
                dsd = dsd/2;
                if (p>yr[i])
                    sd[i] = sd[i] - dsd;
                else
                    sd[i] = sd[i] + dsd;
                delta = log(1/(SQRT2PI*sd[i]*yl[i]));
                if (delta>0)
                    mean[i] = xl[i] + M_SQRT2*sd[i]*sqrt(delta);
                else
                    mean[i] = xl[i] + M_SQRT2*sd[i]*sqrt(-delta);
                p = exp(-(mean[i]-xr[i])*(mean[i]-xr[i])/(2.0*sd[i]*sd[i]))/(SQRT2PI*sd[i]);
            }
            const double SDI = -1.0/(2.0*sd[i]*sd[i]);
            const double SDI_SQRT2PI = 1.0 /(SQRT2PI*sd[i]);
            for (int k=0;k<labelSize;k++){
                y = exp((mean[i]-xr[k])*(mean[i]-xr[k])*SDI)*SDI_SQRT2PI;
                if (yr[k] < y)
                    yr[k] = y;
                y = exp((mean[i]-xl[k])*(mean[i]-xl[k])*SDI)*SDI_SQRT2PI;
                if (yl[k] < y)
                    yl[k] = y;
            }
        }
    }

    QStringList result;
    QStringList line;
    if (n == 1)
        line.append(scheme+ "_item");
    line.append(scheme+ "_value_"+QString::number(n));
    line.append(scheme+ "_tolerance_"+QString::number(n));
    result.append(line.join(","));
    line.clear();
    for (int i=0;i<labelSize;i++){
        if (n == 1)
            line.append(label[i]);
        line.append(QString::number(mean[i],'E',5));
        line.append(QString::number(sd[i],'E',5));
        result.append(line.join(","));
        line.clear();
    }
    delete [] border;
    delete [] mean;
    delete [] sd;
    delete [] dx;
    delete [] yl;
    delete [] yr;
    delete [] xl;
    delete [] xr;

    return result;
}


QStringList extractGaussiankernelCentred (const QString scheme, QString intervals, const int n){

    const QStringList temp = intervals.simplified().split(' ');
    QStringList label;
    const int borderSize = temp.length()/2;
    double *border = new double[borderSize];
    double *mean = new double[borderSize + 1];
    double *sd = new double[borderSize + 1];
    double *alfInterval = new double[borderSize - 1];
    for (int i=0;i<temp.length();i++){
        if (i%2)
            border[i/2]= temp[i].toDouble();
        else
            label.append(temp[i]);
        if (i < borderSize + 1){
            mean[i]=0;
            sd[i]=0;
        }
    }
    for (int i=1;i<borderSize;i++){
        if (border[i-1]>= border[i])
            return QStringList("#The order of intervals must be strictly increasing.");
    }
    if (label.removeDuplicates()!=0)
        return QStringList("#All items must be unique.");
    double sdinit = 0;
    for (int i=0;i<borderSize-1;i++){
        alfInterval[i]=(border[i+1]-border[i])/2;
        mean[i+1]=border[i] + alfInterval[i];
        if (sdinit < alfInterval[i])
            sdinit = alfInterval[i];
    }
    const double p = exp(-0.5)/(sqrt(2*M_PI)*sdinit);
    double sdmax = 0;
    for (int i=0;i<borderSize-1;i++){
        if (alfInterval[i]<sdinit)
            sd[i+1]=alfInterval[i]/(sqrt(-gsl_sf_lambert_W0(-2*M_PI*p*p*alfInterval[i]*alfInterval[i])));
        else
            sd[i+1]=sdinit;
        if (sdmax < sd[i+1])
            sdmax = sd[i+1];
    }
    mean[0] =  border[0] - M_SQRT2*sdmax *sqrt(log(1/(sqrt(2*M_PI)*sdmax *p)));
    mean[borderSize] = border[borderSize-1] + M_SQRT2*sdmax*sqrt(log(1/(sqrt(2*M_PI)*sdmax *p)));
    sd[0]=sdmax;
    sd[borderSize]=sdmax;

    QStringList result;
    QStringList line;
    if (n == 1)
        line.append(scheme+ "_item");
    line.append(scheme+ "_value_"+QString::number(n));
    line.append(scheme+ "_tolerance_"+QString::number(n));
    result.append(line.join(","));
    line.clear();
    for (int i=0;i<label.length();i++){
        if (n == 1)
            line.append(label[i]);
        line.append(QString::number(mean[i],'E',5));
        line.append(QString::number(sd[i],'E',5));
        result.append(line.join(","));
        line.clear();
    }
    delete [] border;
    delete [] mean;
    delete [] sd;
    delete [] alfInterval;

    return result;
}

void macroActualize(const char* file)
{
    QFile error("macro.log");
    error.open(QFile::WriteOnly);
    xmlSubstituteEntitiesDefault (1);
    xmlLoadExtDtdDefaultValue = 1;
    xmlDocPtr doc;
    doc = xmlReadFile (file, "UTF-8", 0 );
    xmlNsPtr nsdis = NULL;
    xmlNsPtr nssdk = NULL;
    xmlNodePtr node = NULL;
    xmlXPathInit();
    xmlXPathContextPtr ctxt = xmlXPathNewContext(doc);
    xmlXPathRegisterNs(ctxt, (const xmlChar*)"sdk",(const xmlChar*)"http://www.nomoseed.org/sdk");
    xmlXPathRegisterNs(ctxt, (const xmlChar*)"discretization",(const xmlChar*)"http://www.nomoseed.org/discretization");
    xmlXPathObjectPtr xpathResMacros = xmlXPathEvalExpression((const xmlChar*)"/*/*/sdk:macro[@active='true' and discretization:discretization]", ctxt);
    QString data;
    for (int i=0; i < xpathResMacros->nodesetval->nodeNr; i++){
        nssdk = xmlSearchNsByHref(doc, xpathResMacros->nodesetval->nodeTab[i], (const xmlChar*)"http://www.nomoseed.org/sdk");
        ctxt->node = xpathResMacros->nodesetval->nodeTab[i];
        xmlXPathObjectPtr xpathResCSV = xmlXPathEvalExpression((const xmlChar*)"sdk:csv/text()", ctxt);
        if (xmlXPathCastToBoolean(xpathResCSV))
            data = QString((char*) xpathResCSV->nodesetval->nodeTab[0]->content);
        else
            data = "";
        QStringList csv = data.split("\n");
        for (int k=csv.length()-1;k>=0;k--)
            if (csv[k].simplified()=="")
                csv.removeAt(k);
        xmlXPathObjectPtr xpathResCentrer = xmlXPathEvalExpression((const xmlChar*)"discretization:discretization/@center = 'true'", ctxt);
        xmlXPathObjectPtr xpathResInterval = xmlXPathEvalExpression((const xmlChar*)".//discretization:*[@interval]", ctxt);
        const QString name((char*) xmlGetProp(xpathResMacros->nodesetval->nodeTab[i], (const xmlChar*)"name"));
        for (int j=0;j < xpathResInterval->nodesetval->nodeNr; j++){
          node = xpathResInterval->nodesetval->nodeTab[j];
          QStringList result;
          if (xmlXPathCastToBoolean(xpathResCentrer))
              result = extractGaussiankernelCentred (name,QString((char*)xmlGetProp(node,(const xmlChar*)"interval")),j+1);
           else
              result = extractGaussiankernel (name, QString((char*)xmlGetProp(node,(const xmlChar*)"interval")),j+1);
          if (result[0][0] == '#'){
              error.write("Macro discretization error in "+ name.toUtf8() +": " + result[0].toUtf8() + "\n");
              xmlFreeDoc (doc);
              xmlXPathFreeObject(xpathResMacros);
              if (xpathResCSV)
                xmlXPathFreeObject(xpathResCSV);
              if (xpathResCentrer)
                xmlXPathFreeObject(xpathResCentrer);
              xmlXPathFreeObject(xpathResInterval);
              xmlXPathFreeContext(ctxt);
              error.close();
              xmlCleanupParser ();
              return;
          }
          if (csv.isEmpty())
              for (int k = 0;k < result.length();k++){
                csv.append(result[k]);}
          else
            for (int k = 0;k < result.length();k++)
                csv[k] = csv[k]+","+result[k];
          xmlSetProp(node, (const xmlChar*)"tolerance", (const xmlChar*)QString("#"+name+"_tolerance_"+QString::number(j+1)).toUtf8().data());
          xmlSetProp(node, (const xmlChar*)"value", (const xmlChar*)QString("#"+name+"_value_"+QString::number(j+1)).toUtf8().data());
          xmlRemoveProp(xmlSetProp(node, (const xmlChar*)"interval", NULL));
        }
        xmlXPathFreeObject(xpathResInterval);
        xmlXPathObjectPtr xpathResInformation = xmlXPathEvalExpression((const xmlChar*)"discretization:discretization/discretization:rule/discretization:conclusion/discretization:information", ctxt);
        if (xmlXPathCastToBoolean(xpathResInformation))
           xmlSetProp(xpathResInformation->nodesetval->nodeTab[0], (const xmlChar*) "value", (const xmlChar*) QString("#"+name+"_item").toUtf8().data());
        else{
          xmlXPathObjectPtr xpathResTemp = xmlXPathEvalExpression((const xmlChar*) ".//discretization:rule/discretization:conclusion", ctxt);
          nsdis = xmlSearchNsByHref(doc, xpathResTemp->nodesetval->nodeTab[0], (const xmlChar*)"http://www.nomoseed.org/discretization");
          node = xmlNewNode(nsdis, (const xmlChar*) "information");
          xmlSetProp (node, (const xmlChar*) "value", (const xmlChar*) QString("#"+name+"_item").toUtf8().data());
          xmlAddChild(xpathResTemp->nodesetval->nodeTab[0], node);
          xmlXPathFreeObject(xpathResTemp);
        }
        xmlXPathFreeObject(xpathResInformation);
        if (xmlXPathCastToBoolean(xpathResCSV)){
            xmlNodePtr node = xpathResCSV->nodesetval->nodeTab[0];
            xmlReplaceNode(node, xmlNewText((const xmlChar*) csv.join("\n").toUtf8().constData()));
            xmlFreeNode(node);
        }
        else{
          node = xmlNewNode(nssdk, (const xmlChar*) "csv");
          xmlAddChild(node, xmlNewText((const xmlChar*) csv.join("\n").toUtf8().constData()));
          xmlXPathObjectPtr xpathResTemp = xmlXPathEvalExpression((const xmlChar*) "discretization:discretization", ctxt);
          xmlAddPrevSibling(xpathResTemp->nodesetval->nodeTab[0], node);
          xmlXPathFreeObject(xpathResTemp);
        }
        if (xmlXPathCastToBoolean(xpathResCSV))
            xmlXPathFreeObject(xpathResCSV);
        if (xmlXPathCastToBoolean(xpathResCentrer))
            xmlXPathFreeObject(xpathResCentrer);
    }
    xmlSaveFileEnc (file, doc, "UTF-8");
    xmlXPathFreeObject(xpathResMacros);
    xmlXPathFreeContext(ctxt);
    xmlFreeDoc (doc);
    error.close();
    xmlCleanupParser ();
}
