#include "include/MMVII_all.h"
#include "IndexBinaire.h"
#include "include/MMVII_Tpl_Images.h"

// #include "include/MMVII_2Include_Serial_Tpl.h"
// #include<map>

/** \file cMethodCalIndex.cpp
    \brief Contains class for comuting one bit of info

*/


namespace MMVII
{


/* ==================================== */
/*         cIB_LinearFoncBool           */
/* ==================================== */

cIB_LinearFoncBool::cIB_LinearFoncBool
(
     cAppli_ComputeParamIndexBinaire & anAppli, 
     const cDenseVect<double>&   aVect, 
     double                aThresh 
)  :
   mAppli    (anAppli),
   mVect     (aVect),
   mThresh   (aThresh)
{
}

cIB_LinearFoncBool::cIB_LinearFoncBool
(
     cAppli_ComputeParamIndexBinaire & anAppli, 
     const cDenseVect<double>&   aVect
)  :
    cIB_LinearFoncBool
    (
        anAppli,
        aVect,
        aVect.DotProduct(anAppli.Stat2().Moy())
    )
{
}
   

cIB_LinearFoncBool::cIB_LinearFoncBool(cAppli_ComputeParamIndexBinaire & anAppli,int aK) :
    cIB_LinearFoncBool
    (
        anAppli,
        anAppli.Eigen().EigenVectors().ReadCol(aK)
    )
/*
   mAppli    (anAppli),
   mK        (aK),
   mVect     (mAppli.Eigen().EigenVectors().ReadCol(aK)),
   mThresh   (mVect.DotProduct(mAppli.Stat2().Moy()))
*/
{
}

bool cIB_LinearFoncBool::Calc(const cVecInvRad & aVIR) const
{
    return RCalc(aVIR) > 0;
}

double cIB_LinearFoncBool::RCalc(const cVecInvRad & aVIR) const 
{
    cDenseVect<tREAL8> & aTmp = mAppli.TmpVect();
    CopyIn(aTmp.DIm(),aVIR.mVec.DIm());

    return mVect.DotProduct(aTmp) - mThresh;
}

const cDenseVect<double>&  cIB_LinearFoncBool::Vect() const 
{
   return mVect;
}



/* ==================================== */
/*         cIB_LinearFoncBool           */
/* ==================================== */

template <class Type> void cAppli_ComputeParamIndexBinaire::AddOneEqParamLin(double aPds,const cDataIm1D<Type> & aCdg,double aVal,int aNb)
{
   if (! aNb) return;
   // aCdg  .  aVect - mThr = 0
   cDenseVect<tREAL8> aV(1+mNbValByP);
   for (int aK=0 ; aK<mNbValByP ; aK++)
      aV(aK) = aCdg.GetV(aK);
   aV(mNbValByP) = -1;
   mLSQOpt.AddObservation(aPds,aV,aVal);

}

void cAppli_ComputeParamIndexBinaire::AddVecOneEqParamLin(double aPds,cVecInvRad & aCdg,double aVal)
{
    // const cDataIm1D<tU_INT1> & aDIm = aCdg.mVec.DIm();
    // AddOneEqParamLin(aPds,aDIm,aVal,1);
    AddOneEqParamLin(aPds,aCdg.mVec.DIm(),aVal,1);
}



void cAppli_ComputeParamIndexBinaire::OptimiseScoreAPriori(tPtVBool  aVBOld,const tVPtVIR & aVIR)  
{
   bool Cont = true;
   int aCpt = 0;
   tPtVBool  aV000 = aVBOld;
   while (Cont)
   {
      aCpt++;
      std::vector<double>  aVScore;
      std::vector<double>  aVSAbs;
      std::vector<double>  aVBadMin;
      std::vector<double>  aVBadMax;

      for (const auto & aV : aVIR)
      {
          double aSc = aVBOld->FB().RCalc(*aV);
          aVScore.push_back(aSc);
          aVSAbs.push_back(std::abs(aSc));
      }
      for (int aK=0 ; aK<int(aVIR.size()) ; aK+=2)
      {
         if ((aVScore[aK]>0) != (aVScore[aK+1]>0))
         {
             aVBadMin.push_back(std::min(aVSAbs[aK],aVSAbs[aK+1]));
             aVBadMax.push_back(std::max(aVSAbs[aK],aVSAbs[aK+1]));
         }
      }

      double aMedAbs = ConstMediane(aVSAbs);
      double aMedSigned = ConstMediane(aVScore);
      double aMedBadMin = ConstMediane(aVBadMin);
      double aMedBadMax = ConstMediane(aVBadMax);
      StdOut() << "MedAbs " << aMedAbs  
               << " MedS " << aMedSigned 
               << " MedBad " << aMedBadMin   <<  " " << aMedBadMax
               << " Sc " << aVBOld->Score() << "\n";

      double aSigma = aMedAbs / 2;

   //  Fill least square

      mLSQOpt.Reset();
      double aPdsClose000 = 10;
      double aPdsCloseCur = 1000;
      double aPdsAver0 = 1e6;

      double aSomPds = 0.0;
      std::vector<double> aVPds;

      for (int aK=0 ; aK<int(aVIR.size()) ; aK+=2)
      {
         aVPds.push_back(exp(-Square(aVSAbs[aK]/aSigma)));
         aVPds.push_back(exp(-Square(aVSAbs[aK+1]/aSigma)));
         if ((aVScore[aK]>0) != (aVScore[aK+1]>0))
         {
             int aKToChange =  (aVSAbs[aK] < aVSAbs[aK+1]) ? aK : (aK+1);
             int aKCompl = (2*aK+1) - aKToChange;
             double aPdsCh = aVPds[aKToChange];
             double aPdsCompl = aVPds[aKCompl];
             double aPds = (aPdsCh - aPdsCompl);
             double aVal =  aVScore[aKCompl];
             double aSign = (aVal >0) ? 1 : -1;
             aVal += aSign * aMedBadMax;
             aSomPds += aPds*aVal;
             AddVecOneEqParamLin(aPds,*(aVIR[aKToChange]),aVal);
         }
         else
         {
             for (int aK2=aK; aK2<aK+2 ; aK2++)
             {
                 // Perte de temps, gain pas evident
                 // AddVecOneEqParamLin(aVPds[aK]/10,*(aVIR[aK]),aVScore[aK]);
             }
         }
      }
      AddEqProxCur(aPdsCloseCur*aSomPds,aVBOld);
      AddEqProxCur(aPdsClose000*aSomPds,aV000);
      AddOneEqParamLin(aPdsAver0*aSomPds,mStat2.Moy().DIm(),0.0,1);
      cDenseVect<double>  aSol = mLSQOpt.Solve();

      tPtVBool aVBNew = VecBoolFromSol(aSol,aVBOld->Index());
      double aDelta = aVBOld->Score() - aVBNew->Score() ;

      static std::vector<double> aVDelta ;
      aVDelta.push_back(aDelta);

      StdOut()  << " L200 " << aV000->FB().Vect().L2Dist(aVBNew->FB().Vect())
                << " L2Cur " << aVBOld->FB().Vect().L2Dist(aVBNew->FB().Vect())
                << " SC priori " <<  aDelta 
                << " Med-Delta " << ConstMediane(aVDelta) << "\n";

      if (aDelta < 0)
      {
         aVBOld = aVBNew;
         mVVBool.at(aVBOld->Index()) = aVBOld;
      }
      Cont = (aDelta < 0) && (aCpt<10);
   }
   StdOut() << "CPT=" << aCpt << "\n";
}

void cAppli_ComputeParamIndexBinaire::AddEqProxCur(double aPdsCloseCur,tPtVBool aVB0)
{
   const cDenseVect<double>& aVK0 = aVB0->FB().Vect() ;
   for (int aKP=0 ; aKP<mNbValByP ; aKP++)
   {
       mLSQOpt.AddObsFixVar(aPdsCloseCur,aKP,aVK0(aKP));
   }
}

tPtVBool cAppli_ComputeParamIndexBinaire::VecBoolFromSol(const  cDenseVect<double> & aSol,int aIndex)
{
   cDenseVect<double>  aSF(mNbValByP);
   for (int aKP=0 ; aKP<mNbValByP ; aKP++)
   {
       aSF(aKP) = aSol(aKP);
   }
   
   return tPtVBool (new cVecBool(aIndex,false,new cIB_LinearFoncBool(*this,aSF,aSol(mNbValByP)),mVIR));
}


void cAppli_ComputeParamIndexBinaire::TestNewParamLinear(const std::vector<tPtVBool>& aOldVB,int aK0Vec)
{
   double aPdsCloseCur = 1e-1;
   double aPdsEq       = 1;
   mLSQOpt.Reset();

   AddEqProxCur(aPdsCloseCur,aOldVB[aK0Vec]);
/*
   tPtVBool aVB0 = aOldVB[aK0Vec];
   const cDenseVect<double>& aVK0 = aVB0->FB().Vect() ;
   for (int aKP=0 ; aKP<mNbValByP ; aKP++)
   {
       mLSQOpt.AddObsFixVar(aPdsCloseCur,aKP,aVK0(aKP));
   }
*/
   for (int aKV=0 ; aKV<int(mVVBool.size()) ; aKV++)
   {
      cVecBool & aVB = *(mVVBool[aKV].get());
      double aPds = aPdsEq *  RandUnif_0_1();
      AddOneEqParamLin(aPds,aVB.Cdg0().DIm(),0.0,aVB.Nb0());
      AddOneEqParamLin(aPds,aVB.Cdg1().DIm(),0.0,aVB.Nb1());
   }
   cDenseVect<double>  aSol = mLSQOpt.Solve();

 
   tPtVBool aVB = VecBoolFromSol(aSol,aOldVB[aK0Vec]->Index());

   std::vector<tPtVBool> aNewV = aOldVB;
   aNewV[aK0Vec] = aVB;
   int aK;
   double aSc = ScoreSol(aK,aNewV);

   if (aSc > mBestSc)
   {
      StdOut() << "D2222 " << aVB->FB().Vect().L2Dist(aOldVB[aK0Vec]->FB().Vect()) << mBestSc << " " << aSc << "\n";
      ChangeVB(aSc,aVB,aK0Vec);
   }
   
/*
*/

   // mVVBool = aVSave;
}


const cDenseVect<tREAL4>&  cVecBool::Cdg0() const {return mCdg0;}
int                        cVecBool::Nb0()  const {return mNb0;}
const cDenseVect<tREAL4>&  cVecBool::Cdg1() const {return mCdg1;}
int                        cVecBool::Nb1()  const {return mNb1;}
int                        cVecBool::Index()  const {return mIndex;}



cVecBool::cVecBool(int Index,bool Med,cIB_LinearFoncBool * aFB,const tVPtVIR & aVIR)  :
    mIndex (Index),
    mFB (aFB),
    mCdg0 (aVIR.at(0)->mVec.DIm().Sz()),
    mNb0  (0),
    mCdg1 (aVIR.at(0)->mVec.DIm().Sz()),
    mNb1  (0)
{
   mVB.reserve(aVIR.size());
   std::vector<double>  aVScore;
   aVScore.reserve(aVIR.size());

   for (const auto & aV : aVIR)
   {
       aVScore.push_back(mFB->RCalc(*aV));
   }
   if (1)
   {
       double aMed = ConstMediane(aVScore);

       // StdOut() << "Meeeddd  " << aMed << "\n";
       if (!Med)
       {
          aMed=0;
       }

       for (auto & aV : aVScore)
       {
           aV -= aMed;
       }
   }
   for (const auto & aV : aVScore)
   {
       mVB.push_back(aV>0);
   }

       //mVB.push_back(aC>0);
   CalcCdgScoreAPriori(aVIR);

   MMVII_INTERNAL_ASSERT_medium((mVB.size()%2)==0,"cVecBool odd size");
}



void cVecBool::CalcCdgScoreAPriori(const tVPtVIR & aVIR)  
{
   mCdg0.DIm().InitNull();
   mCdg1.DIm().InitNull();
   
   double aProp1  = 0;
   double aPropEq = 0;

   for (int aK=0 ; aK<int(mVB.size()) ; aK++)
   {
       const auto & aVec =  aVIR.at(aK)->mVec.DIm();
       if (mVB.at(aK))
       {
          AddIn(mCdg1.DIm(),aVec);
          mNb1++;
          aProp1++;
       }
       else
       {
          AddIn(mCdg0.DIm(),aVec);
          mNb0++;
       }
       if ((aK%2==0) && (mVB.at(aK)== mVB.at(aK+1)))
       {
           aPropEq++;
       }

   }
   if (mNb0) mCdg0.DIm() *=  1/double(mNb0);
   if (mNb1) mCdg1.DIm() *=  1/double(mNb1);

   aProp1 /= mVB.size();
   aPropEq /= (mVB.size() / 2.0);

   double aPropThEq =  Square(aProp1) + Square(1-aProp1);
   mScore = aPropEq/aPropThEq;

    // StdOut()  << "Eq=" << aPropEq << " Sc=" << aPropEq/aPropThEq << " P1=" << aProp1 << "\n";
}

cIB_LinearFoncBool&   cVecBool::FB()
{
   return *(mFB.get());
}


/* ==================================== */
/*          cStatDifBits                */
/* ==================================== */

int NbbBitDif(const std::vector<tPtVBool> & aVVB,const cPt2di & aPair)
{
   int aNbBitDif = 0;
   for (const auto & aVB : aVVB)
       if (aVB->KBit(aPair.x()) != aVB->KBit(aPair.y()))
          aNbBitDif++;
   return aNbBitDif;
}

cStatDifBits::cStatDifBits(const std::vector<cPt2di> & aVPair,const std::vector<tPtVBool> & aVVB) :
   mHNbBitDif (aVVB.size()+1,0),
   mStatR     (mHNbBitDif.size(),0.0),
   mStatRCum  (mHNbBitDif.size(),0.0)
{
    for (const auto & aPair : aVPair)
    {
        mHNbBitDif.at(NbbBitDif(aVVB,aPair)) ++;
/*
        int aNbBitDif = 0;
        for (const auto & aVB : aVVB)
            if (aVB->KBit(aP.x()) != aVB->KBit(aP.y()))
                aNbBitDif++;
        mHNbBitDif.at(aNbBitDif) ++;
*/
    }
    for (int aK=0 ; aK<=int(aVVB.size()) ; aK++)
    {
        mStatR.at(aK) = mHNbBitDif[aK] / double(aVPair.size());
    }
    mStatRCum.at(0) = mStatR.at(0);
    for (int aK=1 ; aK<=int(aVVB.size()) ; aK++)
    {
        mStatRCum.at(aK) = mStatR.at(aK) + mStatRCum.at(aK-1);
    }
}

double cStatDifBits::Score(const cStatDifBits & aStatFalse,double aPdsFalse,int &aKMax) const
{
   double aRes = -1e10;
   for (int aK=0 ; aK<int(mStatRCum.size()) ; aK++)
   {
       double aSc = mStatRCum.at(aK) - aPdsFalse * aStatFalse.mStatRCum.at(aK) ;
       if (aSc > aRes)
       {
           aRes = aSc;
           aKMax = aK;
       }
   }
   return aRes;
}

void  cStatDifBits::Show(const cStatDifBits & aStatFalse,int aK1,int aK2,double aVMax) const
{
   aK1 = std::max(aK1,0);
   aK2 = std::min(aK2,int(mStatRCum.size()-1));
   for (int aK=aK1 ; aK<=aK2 ; aK++)
   {
      StdOut() << "Cum[" << aK << "]=" 
                << " " << mStatRCum.at(aK) 
                << " " << aStatFalse.mStatRCum.at(aK) 
                << " " << mStatRCum.at(aK)/ std::max(1e-8,aStatFalse.mStatRCum.at(aK))
                << "\n";
      if (mStatRCum.at(aK)>aVMax)
         return;
   }
}


};

