mag.cpp

Go to the documentation of this file.

#include "stdafx.h"
#include "mag.h"

TRnd TMAGNodeSimple::Rnd = TRnd(0);
TRnd TMAGNodeBern::Rnd = TRnd(0);
TRnd TMAGNodeBeta::Rnd = TRnd(0);


// MAG affinity matrix

const double TMAGAffMtx::NInf = -DBL_MAX;

TMAGAffMtx::TMAGAffMtx(const TFltV& SeedMatrix) : SeedMtx(SeedMatrix) {
  MtxDim = (int) sqrt((double)SeedMatrix.Len());
  IAssert(MtxDim*MtxDim == SeedMtx.Len());
}

TMAGAffMtx& TMAGAffMtx::operator = (const TMAGAffMtx& Kronecker) {
  if (this != &Kronecker){
    MtxDim=Kronecker.MtxDim;
    SeedMtx=Kronecker.SeedMtx;
  }
  return *this;
}

bool TMAGAffMtx::IsProbMtx() const {
  for (int i = 0; i < Len(); i++) {
    if (At(i) < 0.0 || At(i) > 1.0) return false;
  }
  return true;
}

void TMAGAffMtx::SetRndMtx(TRnd& Rnd, const int& PrmMtxDim, const double& MinProb) {
  MtxDim = PrmMtxDim;
  SeedMtx.Gen(MtxDim*MtxDim);
  for (int p = 0; p < SeedMtx.Len(); p++) {
    do {
      SeedMtx[p] = Rnd.GetUniDev();
    } while (SeedMtx[p] < MinProb);
  }
}

void TMAGAffMtx::SetEpsMtx(const double& Eps1, const double& Eps0, const int& Eps1Val, const int& Eps0Val) {
  for (int i = 0; i < Len(); i++) {
    double& Val = At(i);
    if (Val == Eps1Val) Val = double(Eps1);
    else if (Val == Eps0Val) Val = double(Eps0);
  }
}

void TMAGAffMtx::AddRndNoise(TRnd& Rnd, const double& SDev) {
  Dump("before");
  double NewVal;
  int c =0;
  for (int i = 0; i < Len(); i++) {
    for(c = 0; ((NewVal = At(i)*Rnd.GetNrmDev(1, SDev, 0.8, 1.2)) < 0.01 || NewVal>0.99) && c <1000; c++) { }
    if (c < 999) { At(i) = NewVal; } else { printf("XXXXX\n"); }
  }
  Dump("after");
}

TStr TMAGAffMtx::GetMtxStr() const {
  TChA ChA("[");
  for (int i = 0; i < Len(); i++) {
    ChA += TStr::Fmt("%g", At(i));
    if ((i+1)%GetDim()==0 && (i+1<Len())) { ChA += "; "; }
    else if (i+1<Len()) { ChA += " "; }
  }
  ChA += "]";
  return TStr(ChA);
}

void TMAGAffMtx::GetLLMtx(TMAGAffMtx& LLMtx) {
  LLMtx.GenMtx(MtxDim);
  for (int i = 0; i < Len(); i++) {
    if (At(i) != 0.0) { LLMtx.At(i) = log(At(i)); }
    else { LLMtx.At(i) = NInf; }
  }
}

void TMAGAffMtx::GetProbMtx(TMAGAffMtx& ProbMtx) {
  ProbMtx.GenMtx(MtxDim);
  for (int i = 0; i < Len(); i++) {
    if (At(i) != NInf) { ProbMtx.At(i) = exp(At(i)); }
    else { ProbMtx.At(i) = 0.0; }
  }
}

void TMAGAffMtx::Swap(TMAGAffMtx& Mtx) {
  ::Swap(MtxDim, Mtx.MtxDim);
  SeedMtx.Swap(Mtx.SeedMtx);
}

double TMAGAffMtx::GetMtxSum() const {
  double Sum = 0;
  for (int i = 0; i < Len(); i++) {
    Sum += At(i); }
  return Sum;
}

double TMAGAffMtx::GetRowSum(const int& RowId) const {
  double Sum = 0;
  for (int c = 0; c < GetDim(); c++) {
    Sum += At(RowId, c); }
  return Sum;
}

double TMAGAffMtx::GetColSum(const int& ColId) const {
  double Sum = 0;
  for (int r = 0; r < GetDim(); r++) {
    Sum += At(r, ColId); }
  return Sum;
}

double TMAGAffMtx::Normalize() {
        double Sum = GetMtxSum();
        if(Sum == 0) {
                return 0;
        }

        for(int i = 0; i < Len(); i++) {
                At(i) = At(i) / Sum;
        }
        return Sum;
}

void TMAGAffMtx::Dump(const TStr& MtxNm, const bool& Sort) const {
  /*printf("%s: %d x %d\n", MtxNm.Empty()?"Mtx":MtxNm.CStr(), GetDim(), GetDim());
  for (int r = 0; r < GetDim(); r++) {
    for (int c = 0; c < GetDim(); c++) { printf("  %8.2g", At(r, c)); }
    printf("\n");
  }*/
  if (! MtxNm.Empty()) printf("%s\n", MtxNm.CStr());
  double Sum=0.0;
  TFltV ValV = SeedMtx;
  if (Sort) { ValV.Sort(false); }
  for (int i = 0; i < ValV.Len(); i++) {
    printf("  %10.4g", ValV[i]());
    Sum += ValV[i];
    if ((i+1) % GetDim() == 0) { printf("\n"); }
  }
  printf(" (sum:%.4f)\n", Sum);
}

// average difference in the parameters
double TMAGAffMtx::GetAvgAbsErr(const TMAGAffMtx& Mtx1, const TMAGAffMtx& Mtx2) {
  TFltV P1 = Mtx1.GetMtx();
  TFltV P2 = Mtx2.GetMtx();
  IAssert(P1.Len() == P2.Len());
  P1.Sort();  P2.Sort();
  double delta = 0.0;
  for (int i = 0; i < P1.Len(); i++) {
    delta += fabs(P1[i] - P2[i]);
  }
  return delta/P1.Len();
}

// average L2 difference in the parameters
double TMAGAffMtx::GetAvgFroErr(const TMAGAffMtx& Mtx1, const TMAGAffMtx& Mtx2) {
  TFltV P1 = Mtx1.GetMtx();
  TFltV P2 = Mtx2.GetMtx();
  IAssert(P1.Len() == P2.Len());
  P1.Sort();  P2.Sort();
  double delta = 0.0;
  for (int i = 0; i < P1.Len(); i++) {
    delta += pow(P1[i] - P2[i], 2);
  }
  return sqrt(delta/P1.Len());
}

// get matrix from matlab matrix notation
TMAGAffMtx TMAGAffMtx::GetMtx(TStr MatlabMtxStr) {
  TStrV RowStrV, ColStrV;
  MatlabMtxStr.ChangeChAll(',', ' ');
  MatlabMtxStr.SplitOnAllCh(';', RowStrV);  IAssert(! RowStrV.Empty());
  RowStrV[0].SplitOnWs(ColStrV);    IAssert(! ColStrV.Empty());
  const int Rows = RowStrV.Len();
  const int Cols = ColStrV.Len();
  IAssert(Rows == Cols);
  TMAGAffMtx Mtx(Rows);
  for (int r = 0; r < Rows; r++) {
    RowStrV[r].SplitOnWs(ColStrV);
    IAssert(ColStrV.Len() == Cols);
    for (int c = 0; c < Cols; c++) {
      Mtx.At(r, c) = (double) ColStrV[c].GetFlt(); }
  }
  return Mtx;
}

TMAGAffMtx TMAGAffMtx::GetRndMtx(TRnd& Rnd, const int& Dim, const double& MinProb) {
  TMAGAffMtx Mtx;
  Mtx.SetRndMtx(Rnd, Dim, MinProb);
  return Mtx;
}


// Simple MAG node attributes (Same Bernoulli for every attribute)

void TMAGNodeSimple::AttrGen(TIntVV& AttrVV, const int& NNodes) {
        IAssert(Dim > 0);
        AttrVV.Gen(NNodes, Dim);
        AttrVV.PutAll(0);
        
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < Dim; l++) {
                        if((TMAGNodeSimple::Rnd).GetUniDev() > Mu) {
                                AttrVV.At(i, l) = 1;
                        }
                }
        }
}

void TMAGNodeSimple::LoadTxt(const TStr& InFNm) {
        FILE *fp = fopen(InFNm.CStr(), "r");
        IAssertR(fp != NULL, "File does not exist: " + InFNm);

        char buf[128];
        char *token;
        TStr TokenStr;
        TFlt Val;

        token = strtok(buf, "&");
        token = strtok(token, " \t");
        TokenStr = TStr(token);
        Mu = TokenStr.GetFlt(Val);

        fclose(fp);
}

void TMAGNodeSimple::SaveTxt(TStrV& OutStrV) const {
        OutStrV.Gen(Dim, 0);

        for(int i = 0; i < Dim; i++) {
                OutStrV.Add(TStr::Fmt("%f", double(Mu)));
        }
}


// MAG node attributes with (different) Bernoulli 

TMAGNodeBern& TMAGNodeBern::operator=(const TMAGNodeBern& Dist) {
        MuV = Dist.MuV;
        Dim = Dist.Dim;
        return (*this);
}

void TMAGNodeBern::AttrGen(TIntVV& AttrVV, const int& NNodes) {
        IAssert(Dim > 0);
        AttrVV.Gen(NNodes, Dim);
        AttrVV.PutAll(0);
        
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < Dim; l++) {
                        if((TMAGNodeBern::Rnd).GetUniDev() > MuV[l]) {
                                AttrVV.At(i, l) = 1;
                        }
                }
        }
}

void TMAGNodeBern::LoadTxt(const TStr& InFNm) {
        FILE *fp = fopen(InFNm.CStr(), "r");
        IAssertR(fp != NULL, "File does not exist: " + InFNm);

        Dim = 0;
        MuV.Gen(10, 0);

        char buf[128];
        char *token;
        TStr TokenStr;
        TFlt Val;

        while(fgets(buf, sizeof(buf), fp) != NULL) {
                token = strtok(buf, "&");
                token = strtok(token, " \t");
                TokenStr = TStr(token);
                MuV.Add(TokenStr.GetFlt(Val));
        }
        Dim = MuV.Len();

        fclose(fp);
}

void TMAGNodeBern::SaveTxt(TStrV& OutStrV) const {
        OutStrV.Gen(Dim, 0);

        for(int i = 0; i < Dim; i++) {
                OutStrV.Add(TStr::Fmt("%f", double(MuV[i])));
        }
}


// MAG node attributes with Beta + Bernoulli family

TMAGNodeBeta& TMAGNodeBeta::operator=(const TMAGNodeBeta& Dist) {
        AlphaV = Dist.AlphaV;
        BetaV = Dist.BetaV;
        Dim = Dist.Dim;
        MuV = Dist.MuV;
        Dirty = Dist.Dirty;
        return (*this);
}

void TMAGNodeBeta::SetBeta(const int& Attr, const double& Alpha, const double& Beta) {
        IAssert(Attr < Dim);
        AlphaV[Attr] = Alpha;
        BetaV[Attr] = Beta;
        Dirty = true;
}
        
void TMAGNodeBeta::SetBetaV(const TFltV& _AlphaV, const TFltV& _BetaV) {
        IAssert(_AlphaV.Len() == _BetaV.Len());
        AlphaV = _AlphaV;
        BetaV = _BetaV;
        Dim = _AlphaV.Len();
        Dirty = true;
}

void TMAGNodeBeta::AttrGen(TIntVV& AttrVV, const int& NNodes) {
        IAssert(Dim > 0);
        AttrVV.Gen(NNodes, Dim);
        AttrVV.PutAll(0);

//      printf("AlphaV = %d, BetaV = %d, MuV = %d\n", AlphaV.Len(), BetaV.Len(), MuV.Len());
        
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < Dim; l++) {
                        double x = TMAGNodeBeta::Rnd.GetGammaDev((int)AlphaV[l]);
                        double y = TMAGNodeBeta::Rnd.GetGammaDev((int)BetaV[l]);
                        MuV[l] = x / (x + y);
                        if((TMAGNodeBeta::Rnd).GetUniDev() > MuV[l]) {
                                AttrVV.At(i, l) = 1;
                        }
                }
        }
        Dirty = false;
}

void TMAGNodeBeta::LoadTxt(const TStr& InFNm) {
        FILE *fp = fopen(InFNm.CStr(), "r");
        IAssertR(fp != NULL, "File does not exist: " + InFNm);

        Dim = 0;
        AlphaV.Gen(10, 0);
        BetaV.Gen(10, 0);

        char buf[128];
        char *token;
        TStr TokenStr;
        TFlt Val;

        while(fgets(buf, sizeof(buf), fp) != NULL) {
                token = strtok(buf, "&");
                
                token = strtok(token, " \t");
                TokenStr = TStr(token);
                AlphaV.Add(TokenStr.GetFlt(Val));

                token = strtok(NULL, " \t");
                TokenStr = TStr(token);
                BetaV.Add(TokenStr.GetFlt(Val));

                Dim++;
        }

        fclose(fp);
}

void TMAGNodeBeta::SaveTxt(TStrV& OutStrV) const {
        OutStrV.Gen(Dim, 0);

        for(int i = 0; i < Dim; i++) {
                OutStrV.Add(TStr::Fmt("%f %f", double(AlphaV[i]), double(BetaV[i])));
        }
}


// MAGFit with Bernoulli node attributes

void TMAGFitBern::SetGraph(const PNGraph& GraphPt) {
        Graph = GraphPt;
        bool NodesOk = true;
        // check that nodes IDs are {0,1,..,Nodes-1}
        for (int nid = 0; nid < Graph->GetNodes(); nid++) {
        if (! Graph->IsNode(nid)) { NodesOk=false; break; } }
        if (! NodesOk) {
        TIntV NIdV;  GraphPt->GetNIdV(NIdV);
        Graph = TSnap::GetSubGraph(GraphPt, NIdV, true);
        for (int nid = 0; nid < Graph->GetNodes(); nid++) {
          IAssert(Graph->IsNode(nid)); }
        }
}

void TMAGFitBern::SetPhiVV(const TIntVV& AttrVV, const int KnownIds) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();

        PhiVV.Gen(NNodes, NAttrs);
        KnownVV.Gen(NNodes, NAttrs);
        
        for(int l = 0; l < NAttrs; l++) {
                for(int i = 0; i < NNodes; i++) {
                        if(int(AttrVV(i, l)) == 0) {
                                PhiVV(i, l) = 0.9999;
                        } else {
                                PhiVV(i, l) = 0.0001;
                        }
                }

                if(l < KnownIds) {
                        KnownVV.PutY(l, true);
                } else {
                        KnownVV.PutY(l, false);
                }
        }
}

void TMAGFitBern::SaveTxt(const TStr& FNm) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        const TFltV MuV = GetMuV();
        TMAGAffMtxV MtxV;
        Param.GetMtxV(MtxV);

        FILE *fp = fopen(FNm.GetCStr(), "w");
        for(int l = 0; l < NAttrs; l++) {
                fprintf(fp, "%.4f\t", double(MuV[l]));
                for(int row = 0; row < 2; row++) {
                        for(int col = 0; col < 2; col++) {
                                fprintf(fp, " %.4f", double(MtxV[l].At(row, col)));
                        }
                        fprintf(fp, (row == 0) ? ";" : "\n");
                }
        }
        fclose(fp);

        fp = fopen((FNm + "f").CStr(), "w");
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        fprintf(fp, "%f ", double(PhiVV(i, l)));
                }
                fprintf(fp, "\n");
        }
        fclose(fp);
}
        
void TMAGFitBern::Init(const TFltV& MuV, const TMAGAffMtxV& AffMtxV) {
        TMAGNodeBern DistParam(MuV);
        Param.SetNodeAttr(DistParam);
        Param.SetMtxV(AffMtxV);

        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();

        PhiVV.Gen(NNodes, NAttrs);
        KnownVV.Gen(NNodes, NAttrs);
        KnownVV.PutAll(false);
}

#if 0
void TMAGFitBern::PerturbInit(const TFltV& MuV, const TMAGAffMtxV& AffMtxV, const double& PerturbRate) {
        IAssert(PerturbRate < 1.0);

        TFltV InitMuV = MuV;    //InitMuV.PutAll(0.5);
        TMAGNodeBern DistParam(InitMuV);
        Param.SetMtxV(AffMtxV);
        TRnd& Rnd = TMAGNodeBern::Rnd;
        TMAGAffMtxV PerturbMtxV = AffMtxV;

        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        
        for(int l = 0; l < NAttrs; l++) {
                double Mu = MuV[l] + PerturbRate * (Rnd.GetUniDev() - 0.5) * 2;
//              double Mu = 0.5;
                if(Mu < 0.01) {  Mu = 0.01;  }
                if(Mu > 0.99) {  Mu = 0.99;  }
                DistParam.SetMu(l, Mu);
//              PhiVV.PutY(l, MuV[l] + Perturb);
                TMAGAffMtx AffMtx(AffMtxV[l]);
                for(int p = 0; p < 4; p++) {
                        AffMtx.At(p) += PerturbRate * (Rnd.GetUniDev() - 0.5) * 2;
                        if(AffMtx.At(p) < 0.05) {  AffMtx.At(p) = 0.05;  }
                        if(AffMtx.At(p) > 0.95) {  AffMtx.At(p) = 0.95;  }
                }
                AffMtx.At(0, 1) = AffMtx.At(1, 0);
                PerturbMtxV[l] = AffMtx;
        }
//      NormalizeAffMtxV(PerturbMtxV);

        printf("\n");
        for(int l = 0; l < NAttrs; l++) {
                printf("Mu = %.3f  ", DistParam.GetMu(l));
                printf("AffMtx = %s\n", PerturbMtxV[l].GetMtxStr().GetCStr());
        }
        Param.SetMtxV(PerturbMtxV);
        Param.SetNodeAttr(DistParam);
        
        PhiVV.Gen(NNodes, NAttrs);
        KnownVV.Gen(NNodes, NAttrs);
        KnownVV.PutAll(false);
}
#endif  //      0

void TMAGFitBern::RandomInit(const TFltV& MuV, const TMAGAffMtxV& AffMtxV, const int& Seed) {
        TRnd& Rnd = TMAGNodeBern::Rnd;
        Rnd.PutSeed(Seed);

        TFltV InitMuV = MuV;    InitMuV.PutAll(0.5);
        TMAGNodeBern DistParam(InitMuV);
        Param.SetMtxV(AffMtxV);
        
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        
        PhiVV.Gen(NNodes, NAttrs);
        KnownVV.Gen(NNodes, NAttrs);
        KnownVV.PutAll(false);

        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        PhiVV.At(i, l) = Rnd.GetUniDev();
//                      PhiVV.At(i, l) = 0.5;
                }
        }
        
        TMAGAffMtxV RndMtxV = AffMtxV;
        for(int l = 0; l < NAttrs; l++) {
                for(int p = 0; p < 4; p++) {
                        RndMtxV[l].At(p) = TMAGNodeBern::Rnd.GetUniDev();
                        if(RndMtxV[l].At(p) < 0.1) {  RndMtxV[l].At(p) = 0.1;  }
                        if(RndMtxV[l].At(p) > 0.9) {  RndMtxV[l].At(p) = 0.9;  }
                }
                RndMtxV[l].At(0, 1) = RndMtxV[l].At(1, 0);
        }
        
        printf("\n");
        for(int l = 0; l < NAttrs; l++) {
                printf("AffMtx = %s\n", RndMtxV[l].GetMtxStr().GetCStr());
        }
        Param.SetMtxV(RndMtxV);
        Param.SetNodeAttr(DistParam);
}

const double TMAGFitBern::GetInCoeff(const int& i, const int& j, const int& l, const int& A, const TMAGAffMtx& Theta) const {
        return (PhiVV.At(i, l) * Theta.At(0, A) + (1.0 - PhiVV.At(i, l)) * Theta.At(1, A));
}

const double TMAGFitBern::GetOutCoeff(const int& i, const int& j, const int& l, const int& A, const TMAGAffMtx& Theta) const {
        return (PhiVV.At(j, l) * Theta.At(A, 0) + (1.0 - PhiVV.At(j, l)) * Theta.At(A, 1));
}

const double TMAGFitBern::GetAvgInCoeff(const int& i, const int& AId, const int& A, const TMAGAffMtx& Theta) const {
        const int NNodes = Param.GetNodes();
        const double Mu_l = AvgPhiV[AId] / double(NNodes);
        return (Mu_l * Theta.At(0, A) + (1.0 - Mu_l) * Theta.At(1, A));
}

const double TMAGFitBern::GetAvgOutCoeff(const int& i, const int& AId, const int& A, const TMAGAffMtx& Theta) const {
        const int NNodes = Param.GetNodes();
        const double Mu_l = AvgPhiV[AId] / double(NNodes);
        return (Mu_l * Theta.At(A, 0) + (1.0 - Mu_l) * Theta.At(A, 1));
}

const double TMAGFitBern::GetProbPhi(const int& NId1, const int& NId2, const int& AId, const int& Attr1, const int& Attr2) const {
        double Prob1 = (Attr1 == 0) ? double(PhiVV.At(NId1, AId)) : (1.0 - PhiVV.At(NId1, AId));
        double Prob2 = (Attr2 == 0) ? double(PhiVV.At(NId2, AId)) : (1.0 - PhiVV.At(NId2, AId));
        return (Prob1 * Prob2);
}

const double TMAGFitBern::GetProbMu(const int& NId1, const int& NId2, const int& AId, const int& Attr1, const int& Attr2, const bool Left, const bool Right) const {
        TMAGNodeBern DistParam = Param.GetNodeAttr();
//      double Mu = DistParam.GetMu(AId);
        double Mu = AvgPhiV[AId] / double(Param.GetNodes());
        double Prob1 = (Left) ? double(PhiVV.At(NId1, AId)) : double(Mu);
        double Prob2 = (Right)? double(PhiVV.At(NId2, AId)) : double(Mu);
        Prob1 = (Attr1 == 0) ? Prob1 : 1.0 - Prob1;
        Prob2 = (Attr2 == 0) ? Prob2 : 1.0 - Prob2;
        return (Prob1 * Prob2);
}

const double TMAGFitBern::GetThetaLL(const int& NId1, const int& NId2, const int& AId) const {
        double LL = 0.0;
        const TMAGAffMtx& Mtx = Param.GetMtx(AId);
        for(int A1 = 0; A1 < 2; A1++) {
                for(int A2 = 0; A2 < 2; A2++) {
                        LL += GetProbPhi(NId1, NId2, AId, A1, A2) * Mtx.At(A1, A2);
                }
        }
        return log(LL);
}

const double TMAGFitBern::GetAvgThetaLL(const int& NId1, const int& NId2, const int& AId, const bool Left, const bool Right) const {
        double LL = 0.0;
        const TMAGAffMtx& Mtx = Param.GetMtx(AId);
        for(int A1 = 0; A1 < 2; A1++) {
                for(int A2 = 0; A2 < 2; A2++) {
                        LL += GetProbMu(NId1, NId2, AId, A1, A2, Left, Right) * Mtx.At(A1, A2);
                }
        }
        return log(LL);
}

const double TMAGFitBern::GetSqThetaLL(const int& NId1, const int& NId2, const int& AId) const {
        double LL = 0.0;
        const TMAGAffMtx& Mtx = Param.GetMtx(AId);
        for(int A1 = 0; A1 < 2; A1++) {
                for(int A2 = 0; A2 < 2; A2++) {
                        LL += GetProbPhi(NId1, NId2, AId, A1, A2) * Mtx.At(A1, A2) * Mtx.At(A1, A2);
                }
        }
        return log(LL);
}

const double TMAGFitBern::GetAvgSqThetaLL(const int& NId1, const int& NId2, const int& AId, const bool Left, const bool Right) const {
        double LL = 0.0;
        const TMAGAffMtx& Mtx = Param.GetMtx(AId);
        for(int A1 = 0; A1 < 2; A1++) {
                for(int A2 = 0; A2 < 2; A2++) {
                        LL += GetProbMu(NId1, NId2, AId, A1, A2, Left, Right) * Mtx.At(A1, A2) * Mtx.At(A1, A2);
                }
        }
        return log(LL);
}

const double TMAGFitBern::GetProdLinWeight(const int& NId1, const int& NId2) const {
        const int NAttrs = Param.GetAttrs();
        double LL = 0.0;

        for(int l = 0; l < NAttrs; l++) {
                LL += GetThetaLL(NId1, NId2, l);
        }
//      return LL;
        return LL + log(NormConst);
}

const double TMAGFitBern::GetAvgProdLinWeight(const int& NId1, const int& NId2, const bool Left, const bool Right) const {
        const int NAttrs = Param.GetAttrs();
        double LL = 0.0;

        for(int l = 0; l < NAttrs; l++) {
                LL += GetAvgThetaLL(NId1, NId2, l, Left, Right);
        }
//      return LL;
        return LL + log(NormConst);
}

const double TMAGFitBern::GetProdSqWeight(const int& NId1, const int& NId2) const {
        const int NAttrs = Param.GetAttrs();
        double LL = 0.0;

        for(int l = 0; l < NAttrs; l++) {
                LL += GetSqThetaLL(NId1, NId2, l);
        }
//      return LL;
        return LL + 2 * log(NormConst);
}

const double TMAGFitBern::GetAvgProdSqWeight(const int& NId1, const int& NId2, const bool Left, const bool Right) const {
        const int NAttrs = Param.GetAttrs();
        double LL = 0.0;

        for(int l = 0; l < NAttrs; l++) {
                LL += GetAvgSqThetaLL(NId1, NId2, l, Left, Right);
        }
//      return LL;
        return LL + 2 * log(NormConst);
}

const double LogSumExp(const double LogVal1, const double LogVal2) {
        double MaxExp = (LogVal1 > LogVal2) ? LogVal1 : LogVal2;
        double Sum = exp(LogVal1 - MaxExp) + exp(LogVal2 - MaxExp);
        return (log(Sum) + MaxExp);
}

const double LogSumExp(const TFltV& LogValV) {
        const int Len = LogValV.Len();
        double MaxExp = -DBL_MAX;
        
        for(int i = 0; i < Len; i++) {
                if(MaxExp < LogValV[i]) {  MaxExp = LogValV[i];  }
        }
        
        double Sum = 0.0;
        for(int i = 0; i < Len; i++) {
                Sum += exp(LogValV[i] - MaxExp);
        }

        return (log(Sum) + MaxExp);
}

const double LogSumExp(const double *LogValArray, const int Len) {
        TFltV TmpV(Len);
        for(int i = 0; i < Len; i++) {  TmpV[i] = LogValArray[i];  }
        return LogSumExp(TmpV);
}

const double TMAGFitBern::GradPhiMI(const double& x, const int& NId, const int& AId, const double& Lambda, const double& DeltaQ, const TFltVV& CntVV) {
        const int NAttrs = CntVV.GetYDim();
        double Grad = DeltaQ - log(x) + log(1.0-x);

        for(int l = 0; l < NAttrs; l++) {
                if(l == AId) {  continue;  }
                const double C0 = PhiVV(NId, l);
                const double C1 = 1.0 - C0;
                Grad -= Lambda * C0 * log(CntVV(0, l) + C0 * x);
                Grad -= Lambda * C1 * log(CntVV(1, l) + C1 * x);
                Grad += Lambda * C0 * log(CntVV(2, l) + C0 * (1-x));
                Grad += Lambda * C1 * log(CntVV(3, l) + C1 * (1-x));
                Grad -= Lambda * log(CntVV(0, l) + CntVV(1, l) + x);
                Grad += Lambda * log(CntVV(2, l) + CntVV(3, l) + (1-x));
        }

        return Grad;
}

const double TMAGFitBern::ObjPhiMI(const double& x, const int& NId, const int& AId, const double& Lambda, const double& Q0, const double& Q1, const TFltVV& CntVV) {
        const int NAttrs = CntVV.GetYDim();
        double Val = x*(Q0 - log(x)) + (1-x)*(Q1 - log(1.0-x));

        for(int l = 0; l < NAttrs; l++) {
                if(l == AId) {  continue;  }
                const double C0 = PhiVV(NId, l);
                const double C1 = 1.0 - C0;
                Val -= Lambda * (CntVV(0, l) + C0 * x) * log(CntVV(0, l) + C0 * x);
                Val -= Lambda * (CntVV(1, l) + C1 * x) * log(CntVV(1, l) + C1 * x);
                Val -= Lambda * (CntVV(2, l) + C0 * (1-x)) * log(CntVV(2, l) + C0 * (1-x));
                Val -= Lambda * (CntVV(3, l) + C1 * (1-x)) * log(CntVV(3, l) + C1 * (1-x));
                Val += Lambda * (CntVV(0, l) + CntVV(1, l) + x) * log(CntVV(0, l) + CntVV(1, l) + x);
                Val += Lambda * (CntVV(2, l) + CntVV(3, l) + 1 - x) * log(CntVV(2, l) + CntVV(3, l) + (1-x));

                if(!(CntVV(0, l) > 0))  printf("CntVV(0, %d) = %.2f\n", l, double(CntVV(0, l)));
                if(!(CntVV(1, l) > 0))  printf("CntVV(1, %d) = %.2f\n", l, double(CntVV(1, l)));
                if(!(CntVV(2, l) > 0))  printf("CntVV(2, %d) = %.2f\n", l, double(CntVV(2, l)));
                if(!(CntVV(3, l) > 0))  printf("CntVV(3, %d) = %.2f\n", l, double(CntVV(3, l)));
        }

        return Val;
}

const double TMAGFitBern::GetEstNoEdgeLL(const int& NId, const int& AId) const {
        // const int NNodes = Param.GetNodes();
        // const int NAttrs = Param.GetAttrs();
        
        TMAGNodeBern DistParam = Param.GetNodeAttr();
        double LL = 0.0;

        return LL;
}

const double TMAGFitBern::UpdatePhi(const int& NId, const int& AId, double& Phi) {
        TMAGAffMtx LLTheta, Theta = Param.GetMtx(AId); 
        TMAGAffMtx SqTheta(Theta);
        const int NNodes = Param.GetNodes();
        // const int NAttrs = Param.GetAttrs();
        Theta.GetLLMtx(LLTheta);
        TMAGNodeBern DistParam = Param.GetNodeAttr();
        const double Mu = DistParam.GetMu(AId);

        for(int i = 0; i < Theta.Len(); i++) {
                SqTheta.At(i) = SqTheta.At(i) * SqTheta.At(i);
        }

        //      Using log-sum-exp trick
        double EdgeQ[2], NonEdgeQ[2];
        double MaxExp[2];
        TFltV NonEdgeLLV[2];
        for(int i = 0; i < 2; i++) {
                EdgeQ[i] = 0.0;
                NonEdgeQ[i] = 0.0;
                MaxExp[i] = -DBL_MAX;
                NonEdgeLLV[i].Gen(4 * NNodes, 0);
        }
        
        for(int j = 0; j < NNodes; j++) {
                if(j == NId) {  continue;       }

                if(Graph->IsEdge(NId, j)) {
                        EdgeQ[0] += GetOutCoeff(NId, j, AId, 0, LLTheta);
                        EdgeQ[1] += GetOutCoeff(NId, j, AId, 1, LLTheta);
                } else {
                        double LinW = GetProdLinWeight(NId, j) - GetThetaLL(NId, j, AId);
                        double SqW = GetProdSqWeight(NId, j) - GetSqThetaLL(NId, j, AId);

                        for(int i = 0; i < 2; i++) {
                                NonEdgeLLV[i].Add(LinW + log(GetOutCoeff(NId, j, AId, i, Theta)));
                                NonEdgeLLV[i].Add(SqW + log(GetOutCoeff(NId, j, AId, i, SqTheta)) + log(0.5));
                        }
                }

                if(Graph->IsEdge(j, NId)) {
                        EdgeQ[0] += GetInCoeff(j, NId, AId, 0, LLTheta);
                        EdgeQ[1] += GetInCoeff(j, NId, AId, 1, LLTheta);
                } else {
                        double LinW = GetProdLinWeight(j, NId) - GetThetaLL(j, NId, AId);
                        double SqW = GetProdSqWeight(j, NId) - GetSqThetaLL(j, NId, AId);

                        for(int i = 0; i < 2; i++) {
                                NonEdgeLLV[i].Add(LinW + log(GetInCoeff(j, NId, AId, i, Theta)));
                                NonEdgeLLV[i].Add(SqW + log(GetInCoeff(j, NId, AId, i, SqTheta)) + log(0.5));
                        }
                }
        }

        NonEdgeQ[0] = LogSumExp(NonEdgeLLV[0]);
        NonEdgeQ[1] = LogSumExp(NonEdgeLLV[1]);
        
        double Q[2];
        Q[0] = log(Mu) + EdgeQ[0] - exp(NonEdgeQ[0]);
        Q[1] = log(1.0 - Mu) + EdgeQ[1] - exp(NonEdgeQ[1]);
//      double Q = Q1 - Q0;
//      printf("  [Phi_{%d}{%d}]  :: Q0 = %f, Q1 = %f\n", NId, AId, Q0, Q1);

//      Phi = 1.0 / (1.0 + exp(Q));
        Phi = Q[0] - LogSumExp(Q, 2);
        Phi = exp(Phi);

        return Phi - PhiVV.At(NId, AId);
}


const double TMAGFitBern::UpdatePhiMI(const double& Lambda, const int& NId, const int& AId, double& Phi) {
        TMAGAffMtx LLTheta, Theta = Param.GetMtx(AId); 
        TMAGAffMtx SqTheta(Theta);
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        Theta.GetLLMtx(LLTheta);
        TMAGNodeBern DistParam = Param.GetNodeAttr();
        const double Mu = DistParam.GetMu(AId);

        for(int i = 0; i < Theta.Len(); i++) {
                SqTheta.At(i) = SqTheta.At(i) * SqTheta.At(i);
        }

        //      Using log-sum-exp trick
        double EdgeQ[2], NonEdgeQ[2];
        double MaxExp[2];
        TFltV NonEdgeLLV[2];
        TFltVV CntVV(4, NAttrs);                CntVV.PutAll(0.0);
        for(int i = 0; i < 2; i++) {
                EdgeQ[i] = 0.0;
                NonEdgeQ[i] = 0.0;
                MaxExp[i] = -DBL_MAX;
                NonEdgeLLV[i].Gen(4 * NNodes, 0);
        }
        
        for(int j = 0; j < NNodes; j++) {
                if(j == NId) {  continue;       }

                for(int l = 0; l < NAttrs; l++) {
                        if(l == AId) {  continue;  }
                        CntVV(0, l) = CntVV(0, l) + PhiVV(j, AId) * PhiVV(j, l);
                        CntVV(1, l) = CntVV(1, l) + PhiVV(j, AId) * (1.0-PhiVV(j, l));
                        CntVV(2, l) = CntVV(2, l) + (1.0-PhiVV(j, AId)) * PhiVV(j, l);
                        CntVV(3, l) = CntVV(3, l) + (1.0-PhiVV(j, AId)) * (1.0-PhiVV(j, l));
                }

                if(Graph->IsEdge(NId, j)) {
                        EdgeQ[0] += GetOutCoeff(NId, j, AId, 0, LLTheta);
                        EdgeQ[1] += GetOutCoeff(NId, j, AId, 1, LLTheta);
                } else {
                        double LinW = GetProdLinWeight(NId, j) - GetThetaLL(NId, j, AId);
                        double SqW = GetProdSqWeight(NId, j) - GetSqThetaLL(NId, j, AId);

                        for(int i = 0; i < 2; i++) {
                                NonEdgeLLV[i].Add(LinW + log(GetOutCoeff(NId, j, AId, i, Theta)));
                                NonEdgeLLV[i].Add(SqW + log(GetOutCoeff(NId, j, AId, i, SqTheta)) + log(0.5));
                        }
                }

                if(Graph->IsEdge(j, NId)) {
                        EdgeQ[0] += GetInCoeff(j, NId, AId, 0, LLTheta);
                        EdgeQ[1] += GetInCoeff(j, NId, AId, 1, LLTheta);
                } else {
                        double LinW = GetProdLinWeight(j, NId) - GetThetaLL(j, NId, AId);
                        double SqW = GetProdSqWeight(j, NId) - GetSqThetaLL(j, NId, AId);

                        for(int i = 0; i < 2; i++) {
                                NonEdgeLLV[i].Add(LinW + log(GetInCoeff(j, NId, AId, i, Theta)));
                                NonEdgeLLV[i].Add(SqW + log(GetInCoeff(j, NId, AId, i, SqTheta)) + log(0.5));
                        }
                }
        }
        
        NonEdgeQ[0] = LogSumExp(NonEdgeLLV[0]);
        NonEdgeQ[1] = LogSumExp(NonEdgeLLV[1]);
        
        double Q[2];
        Q[0] = log(Mu) + EdgeQ[0] - exp(NonEdgeQ[0]);
        Q[1] = log(1.0 - Mu) + EdgeQ[1] - exp(NonEdgeQ[1]);
        double DeltaQ = Q[0] - Q[1];

//      double x[] = {0.1, 0.3, 0.5, 0.7, 0.9};
        double x[] = {PhiVV(NId, AId)};
//      for(int n = 0; n < 5; n++) {
        for(int n = 0; n < 1; n++) {
//              double LrnRate = 0.0002;
                double LrnRate = 0.001;
                for(int step = 0; step < 200; step++) {
                        double Grad = GradPhiMI(x[n], NId, AId, Lambda, DeltaQ, CntVV);
                        if(Grad > 0.0) {  x[n] += LrnRate;  }
                        else {  x[n] -= LrnRate;  }
                        if(x[n] > 0.9999) {  x[n] = 0.9999;  }
                        if(x[n] < 0.0001) {  x[n] = 0.0001;  }
                        LrnRate *= 0.995;
                }
        }

        double MaxVal = -DBL_MAX;
        int MaxX = -1;
//      for(int n = 0; n < 5; n++) {
        for(int n = 0; n < 1; n++) {
                double Val = ObjPhiMI(x[n], NId, AId, Lambda, Q[0], Q[1], CntVV);
                if(Val > MaxVal) {
                        MaxVal = Val;
                        MaxX = n;
                }
        }
        IAssert(MaxX >= 0);

        Phi = x[MaxX];

        return Phi - PhiVV.At(NId, AId);
}


const double TMAGFitBern::UpdateApxPhiMI(const double& Lambda, const int& NId, const int& AId, double& Phi, TFltVV& ProdVV) {
        TMAGAffMtx LLTheta, Theta = Param.GetMtx(AId); 
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        Theta.GetLLMtx(LLTheta);
        TMAGNodeBern DistParam = Param.GetNodeAttr();
        const double Mu = DistParam.GetMu(AId);

        TMAGAffMtx SqTheta(Theta);
        for(int i = 0; i < Theta.Len(); i++) {
                SqTheta.At(i) = SqTheta.At(i) * SqTheta.At(i);
        }

        TFltV ProdV;    ProdVV.GetRow(NId, ProdV);
        ProdV[0] -= GetAvgThetaLL(NId, NId, AId, true, false);
        ProdV[1] -= GetAvgThetaLL(NId, NId, AId, false, true);
        ProdV[2] -= log(2.0) + GetAvgSqThetaLL(NId, NId, AId, true, false);
        ProdV[3] -= log(2.0) + GetAvgSqThetaLL(NId, NId, AId, false, true);

        //      Using log-sum-exp trick
        double EdgeQ[2];
        double MaxExp[2];
        TFltV NonEdgeLLV[2];
        TFltVV CntVV(4, NAttrs);                CntVV.PutAll(0.0);
        for(int i = 0; i < 2; i++) {
                EdgeQ[i] = 0.0;
                MaxExp[i] = -DBL_MAX;
                NonEdgeLLV[i].Gen(4 * NNodes, 0);
        }

        for(int F = 0; F < 2; F++) {
                NonEdgeLLV[F].Add(ProdV[0] + log(GetAvgOutCoeff(NId, AId, F, Theta)));
                NonEdgeLLV[F].Add(ProdV[1] + log(GetAvgInCoeff(NId, AId, F, Theta)));
                NonEdgeLLV[F].Add(ProdV[2] + log(GetAvgOutCoeff(NId, AId, F, SqTheta)));
                NonEdgeLLV[F].Add(ProdV[3] + log(GetAvgInCoeff(NId, AId, F, SqTheta)));
        }
        EdgeQ[0] = -(NNodes - 1) * exp(LogSumExp(NonEdgeLLV[0]));
        EdgeQ[1] = -(NNodes - 1) * exp(LogSumExp(NonEdgeLLV[1]));

        
        for(int l = 0; l < NAttrs; l++) {
                if(l == AId) {  continue;  }
                int BgId = (AId > l) ? AId : l;
                int SmId = (AId + l) - BgId;
                int SmL = (l < AId) ? 1 : 0;
                BgId *= 4;
                CntVV(0, l) = AvgPhiPairVV(SmId, BgId) - PhiVV(NId, AId) * PhiVV(NId, l);
                CntVV(1+SmL, l) = AvgPhiPairVV(SmId, BgId+1+SmL) - PhiVV(NId, AId) * (1.0-PhiVV(NId, l));
                CntVV(2-SmL, l) = AvgPhiPairVV(SmId, BgId+2-SmL) - (1.0-PhiVV(NId, AId)) * PhiVV(NId, l);
                CntVV(3, l) = AvgPhiPairVV(SmId, BgId+3) - (1.0-PhiVV(NId, AId)) * (1.0-PhiVV(NId, l));
        }

        TNGraph::TNodeI NI = Graph->GetNI(NId);
        for(int d = 0; d < NI.GetOutDeg(); d++) {
                int Out = NI.GetOutNId(d);
                if(NId == Out) {  continue;  }
                double LinW = GetProdLinWeight(NId, Out) - GetThetaLL(NId, Out, AId);
                double SqW = GetProdSqWeight(NId, Out) - GetSqThetaLL(NId, Out, AId);

                for(int F = 0; F < 2; F++) {
                        EdgeQ[F] += GetOutCoeff(NId, Out, AId, F, LLTheta);
                        EdgeQ[F] += exp(LinW + log(GetOutCoeff(NId, Out, AId, F, Theta)));
                        EdgeQ[F] += 0.5 * exp(SqW + log(GetOutCoeff(NId, Out, AId, F, SqTheta)));
                }
        }
        for(int d = 0; d < NI.GetInDeg(); d++) {
                int In = NI.GetInNId(d);
                if(NId == In) {  continue;  }
                double LinW = GetProdLinWeight(In, NId) - GetThetaLL(In, NId, AId);
                double SqW = GetProdSqWeight(In, NId) - GetSqThetaLL(In, NId, AId);

                for(int F = 0; F < 2; F++) {
                        EdgeQ[F] += GetInCoeff(In, NId, AId, F, LLTheta);
                        EdgeQ[F] += exp(LinW + log(GetInCoeff(In, NId, AId, F, Theta)));
                        EdgeQ[F] += 0.5 * exp(SqW + log(GetInCoeff(In, NId, AId, F, SqTheta)));
                }
        }

        EdgeQ[0] += log(Mu);
        EdgeQ[1] += log(1.0 - Mu);
        double DeltaQ = EdgeQ[0] - EdgeQ[1];
//      printf("(%d, %d) :: Q[0] = %f, Q[1] = %f\n", NId, AId, EdgeQ[0], EdgeQ[1]);

//      double x[] = {0.1, 0.3, 0.5, 0.7, 0.9};
        double x[] = {PhiVV(NId, AId)};
        TFltV ObjValV;  ObjValV.Gen(60, 0);
//      for(int n = 0; n < 5; n++) {
        for(int n = 0; n < 1; n++) {
//              double LrnRate = 0.0002;
                double LrnRate = 0.001;
                for(int step = 0; step < 50; step++) {
//              for(int step = 0; step < 10; step++) {
                        double Grad = GradPhiMI(x[n], NId, AId, Lambda, DeltaQ, CntVV);
                        if(Grad > 0.0) {  x[n] += LrnRate;  }
                        else {  x[n] -= LrnRate;  }
                        if(x[n] > 0.9999) {  x[n] = 0.9999;  }
                        if(x[n] < 0.0001) {  x[n] = 0.0001;  }
                        if(x[n] == 0.9999 || x[n] == 0.0001) {
                                break;
                        }
                        LrnRate *= 0.995;
                }
                ObjValV.Add(x[n]);
//              ObjValV.Add(PhiVV(NId, AId));
        }

        double MaxVal = -DBL_MAX;
        int MaxX = -1;
//      for(int n = 0; n < 5; n++) {
        for(int n = 0; n < ObjValV.Len(); n++) {
                double Val = ObjPhiMI(ObjValV[n], NId, AId, Lambda, EdgeQ[0], EdgeQ[1], CntVV);
                if(Val > MaxVal) {
                        MaxVal = Val;
                        MaxX = n;
                } else if(MaxX < 0) {
                        printf("(%d, %d) : %f  Q[0] = %f  Q[1] = %f  Val = %f\n", NId, AId, double(x[n]), double(EdgeQ[0]), double(EdgeQ[1]), Val);
                }
        }
        IAssert(MaxX >= 0);

        Phi = ObjValV[MaxX];

        return Phi - PhiVV.At(NId, AId);
}



double TMAGFitBern::DoEStepOneIter(const TFltV& TrueMuV, TFltVV& NewPhiVV, const double& Lambda) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        double MaxDelta = 0, L1 = 0;
        double Val;
        TFltIntIntTrV NewVal;
        int RndCount = 0;
        // double OldMI = 0.0, NewMI = 0.0;
        TFltV MuV(NAttrs);      MuV.PutAll(0.0);
        TIntV NIndV(NNodes), AIndV(NAttrs);

        //      Update Phi
        /*
        for(int i = 0; i < NNodes; i++) {  NIndV[i] = i;  }
        for(int l = 0; l < NAttrs; l++) {  AIndV[l] = l;  }
        if(Randomized) {
                NIndV.Shuffle(TMAGNodeBern::Rnd);
                AIndV.Shuffle(TMAGNodeBern::Rnd);
        }
        */

        NewVal.Gen(NAttrs * 2);
        for(int i = 0; i < NNodes; i++) {
//              const int NId = NIndV[i]%NNodes;
                for(int l = 0; l < NAttrs * 2; l++) {
                        const int NId = TMAGNodeBern::Rnd.GetUniDevInt(NNodes);
                        const int AId = TMAGNodeBern::Rnd.GetUniDevInt(NAttrs);
//                      const int AId = AIndV[l]%NAttrs;
//                      double Delta = UpdatePhi(NId, AId, Val);
                        double Delta = 0.0;
                        if(KnownVV(NId, AId)) {
                                Val = PhiVV.At(NId, AId);
                        } else {
                                Delta = UpdatePhiMI(Lambda, NId, AId, Val);
                        }

//                      PhiVV.At(NId, AId) = Val;
                        NewVal[l] = TFltIntIntTr(Val, NId, AId);
                        
//                      MuV[AId] = MuV[AId] + Val;
                        if(fabs(Delta) > MaxDelta) {
                                MaxDelta = fabs(Delta);
                        }
                        if(Val > 0.3 && Val < 0.7) {    RndCount++;     }
                }

                for(int l = 0; l < NAttrs * 2; l++) {
                        const int NId = NewVal[l].Val2;
                        const int AId = NewVal[l].Val3;
                        PhiVV.At(NId, AId) = NewVal[l].Val1;
                }
        }
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        MuV[l] = MuV[l] + PhiVV.At(i, l);
                }
        }
        for(int l = 0; l < NAttrs; l++) {
                MuV[l] = MuV[l] / double(NNodes);
        }

        TFltV SortMuV = MuV;
        double Avg = 0.0;
        SortMuV.Sort(false);
        for(int l = 0; l < NAttrs; l++) {
                printf("  F[%d] = %.3f", l, double(MuV[l]));
                Avg += SortMuV[l];
                L1 += fabs(TrueMuV[l] - SortMuV[l]);
        }
        printf("\n");
        printf("  Rnd = %d(%.3f)", RndCount, double(RndCount) / double(NNodes * NAttrs));
        printf("  Avg = %.3f\n", Avg / double(NAttrs));
        L1 /= double(NAttrs);

        return L1;
}


double TMAGFitBern::DoEStepApxOneIter(const TFltV& TrueMuV, TFltVV& NewPhiVV, const double& Lambda) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        double MaxDelta = 0, L1 = 0;
        double Val;
        TFltIntIntTrV NewVal;
        int RndCount = 0;
        // double OldMI = 0.0, NewMI = 0.0;
        TFltV MuV(NAttrs);      MuV.PutAll(0.0);
        TFltVV ProdVV(NNodes, 4);       ProdVV.PutAll(0.0);
        TIntV NIndV(NNodes), AIndV(NAttrs);

        //      Update Phi
        /*
        for(int i = 0; i < NNodes; i++) {  NIndV[i] = i;  }
        for(int l = 0; l < NAttrs; l++) {  AIndV[l] = l;  }
        if(Randomized) {
                NIndV.Shuffle(TMAGNodeBern::Rnd);
                AIndV.Shuffle(TMAGNodeBern::Rnd);
        }
        */

        AvgPhiV.Gen(NAttrs);    AvgPhiV.PutAll(0.0);
        AvgPhiPairVV.Gen(NAttrs, 4*NAttrs);             AvgPhiPairVV.PutAll(0.0);
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        for(int p = l+1; p < NAttrs; p++) {
                                int index = 4 * p;
                                AvgPhiPairVV(l, index) += PhiVV(i, l) * PhiVV(i, p);
                                AvgPhiPairVV(l, index+1) += PhiVV(i, l) * (1.0-PhiVV(i, p));
                                AvgPhiPairVV(l, index+2) += (1.0-PhiVV(i, l)) * PhiVV(i, p);
                                AvgPhiPairVV(l, index+3) += (1.0-PhiVV(i, l)) * (1.0-PhiVV(i, p));
                        }
                        AvgPhiV[l] += PhiVV(i, l);
                }
        }
        for(int i = 0; i < NNodes; i++) {
                ProdVV(i, 0) = GetAvgProdLinWeight(i, i, true, false);
                ProdVV(i, 1) = GetAvgProdLinWeight(i, i, false, true);
                ProdVV(i, 2) = GetAvgProdSqWeight(i, i, true, false);
                ProdVV(i, 3) = GetAvgProdSqWeight(i, i, false, true);
        }

        const int Iter = 3;
        int NId;

        NewVal.Gen(NAttrs * Iter);
        for(int i = 0; i < NNodes * Iter; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        const int AId = TMAGNodeBern::Rnd.GetUniDevInt(NAttrs);
                        double Delta = 0.0;
                        if(KnownVV(NId, AId)) {
                                Val = PhiVV.At(NId, AId);
                        } else {
                                Delta = UpdateApxPhiMI(Lambda, NId, AId, Val, ProdVV);
                        }

//                      PhiVV.At(NId, AId) = Val;
                        NewVal[l] = TFltIntIntTr(Val, NId, AId);
                        
//                      MuV[AId] = MuV[AId] + Val;
                        if(fabs(Delta) > MaxDelta) {
                                MaxDelta = fabs(Delta);
                        }
                        if(Val > 0.3 && Val < 0.7) {    RndCount++;     }
                }

                for(int l = 0; l < NAttrs; l++) {
                        const int NId = NewVal[l].Val2;
                        const int AId = NewVal[l].Val3;

                        ProdVV(NId, 0) -= GetAvgThetaLL(NId, NId, AId, true, false);
                        ProdVV(NId, 1) -= GetAvgThetaLL(NId, NId, AId, false, true);
                        ProdVV(NId, 2) -= GetAvgSqThetaLL(NId, NId, AId, true, false);
                        ProdVV(NId, 3) -= GetAvgSqThetaLL(NId, NId, AId, false, true);
                        for(int p = 0; p < NAttrs; p++) {
                                if(p > AId) {
                                        int index = 4 * p;
                                        AvgPhiPairVV(AId, index) -= PhiVV(NId, AId) * PhiVV(NId, p);
                                        AvgPhiPairVV(AId, index+1) -= PhiVV(NId, AId) * (1.0-PhiVV(NId, p));
                                        AvgPhiPairVV(AId, index+2) -= (1.0-PhiVV(NId, AId)) * PhiVV(NId, p);
                                        AvgPhiPairVV(AId, index+3) -= (1.0-PhiVV(NId, AId)) * (1.0-PhiVV(NId, p));
                                } else if (p < AId) {
                                        int index = 4 * AId;
                                        AvgPhiPairVV(p, index) -= PhiVV(NId, p) * PhiVV(NId, AId);
                                        AvgPhiPairVV(p, index+1) -= PhiVV(NId, p) * (1.0-PhiVV(NId, AId));
                                        AvgPhiPairVV(p, index+2) -= (1.0-PhiVV(NId, p)) * PhiVV(NId, AId);
                                        AvgPhiPairVV(p, index+3) -= (1.0-PhiVV(NId, p)) * (1.0-PhiVV(NId, AId));
                                }
                        }
                        AvgPhiV[AId] -= PhiVV(NId, AId);

                        PhiVV.At(NId, AId) = NewVal[l].Val1;
                        
                        ProdVV(NId, 0) += GetAvgThetaLL(NId, NId, AId, true, false);
                        ProdVV(NId, 1) += GetAvgThetaLL(NId, NId, AId, false, true);
                        ProdVV(NId, 2) += GetAvgSqThetaLL(NId, NId, AId, true, false);
                        ProdVV(NId, 3) += GetAvgSqThetaLL(NId, NId, AId, false, true);
                        for(int p = 0; p < NAttrs; p++) {
                                if(p > AId) {
                                        int index = 4 * p;
                                        AvgPhiPairVV(AId, index) += PhiVV(NId, AId) * PhiVV(NId, p);
                                        AvgPhiPairVV(AId, index+1) += PhiVV(NId, AId) * (1.0-PhiVV(NId, p));
                                        AvgPhiPairVV(AId, index+2) += (1.0-PhiVV(NId, AId)) * PhiVV(NId, p);
                                        AvgPhiPairVV(AId, index+3) += (1.0-PhiVV(NId, AId)) * (1.0-PhiVV(NId, p));
                                } else if (p < AId) {
                                        int index = 4 * AId;
                                        AvgPhiPairVV(p, index) += PhiVV(NId, p) * PhiVV(NId, AId);
                                        AvgPhiPairVV(p, index+1) += PhiVV(NId, p) * (1.0-PhiVV(NId, AId));
                                        AvgPhiPairVV(p, index+2) += (1.0-PhiVV(NId, p)) * PhiVV(NId, AId);
                                        AvgPhiPairVV(p, index+3) += (1.0-PhiVV(NId, p)) * (1.0-PhiVV(NId, AId));
                                }
                        }
                        AvgPhiV[AId] += PhiVV(NId, AId);
                }
        }

        for(int l = 0; l < NAttrs; l++) {
                MuV[l] = AvgPhiV[l] / double(NNodes);
        }

        TFltV SortMuV = MuV;
        double Avg = 0.0;
//      SortMuV.Sort(false);
        for(int l = 0; l < NAttrs; l++) {
                printf("  F[%d] = %.3f", l, double(MuV[l]));
                Avg += SortMuV[l];
//              L1 += fabs(TrueMuV[l] - SortMuV[l]);
        }
        printf("\n");
        printf("  Rnd = %d(%.3f)", RndCount, double(RndCount) / double(NNodes * NAttrs));
        printf("  Avg = %.3f\n", Avg / double(NAttrs));
//      printf("  Linf = %f\n", MaxDelta);
//      L1 /= double(NAttrs);

        return L1;
}

double TMAGFitBern::DoEStep(const TFltV& TrueMuV, const int& NIter, double& LL, const double& Lambda) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TFltVV NewPhiVV(NNodes, NAttrs);
        // double MI;

        TFltV Delta(NIter);

        for(int i = 0; i < NIter; i++) {
                TExeTm IterTm;
                printf("EStep iteration : %d\n", (i+1));
                if(ESpeedUp) {
                        Delta[i] = DoEStepApxOneIter(TrueMuV, NewPhiVV, Lambda);
                } else {
                        Delta[i] = DoEStepOneIter(TrueMuV, NewPhiVV, Lambda);
                }
//              PhiVV = NewPhiVV;

                printf("  (Time = %s)\n", IterTm.GetTmStr());
        }
        printf("\n");

        NewPhiVV.Clr();

        return Delta.Last();
}

const double TMAGFitBern::UpdateMu(const int& AId) {
        const int NNodes = Param.GetNodes();
        TMAGNodeBern DistParam = Param.GetNodeAttr();
        const double OldMu = DistParam.GetMu(AId);
        double NewMu = 0.0;

        for(int i = 0; i < NNodes; i++) {
                NewMu += PhiVV.At(i, AId);
        }
        AvgPhiV[AId] = NewMu;
        NewMu /= double(NNodes);

        printf("      [Posterior Mu] = %.4f\n", NewMu);

        double Delta = fabs(NewMu - OldMu);
        DistParam.SetMu(AId, NewMu);
        Param.SetNodeAttr(DistParam);

        return Delta;
}

const void TMAGFitBern::GradAffMtx(const int& AId, const TFltVV& ProdVV, const TFltVV& SqVV, const TMAGAffMtx& CurMtx, TFltV& GradV) {
        const int NNodes = Param.GetNodes();
        // const int NAttrs = Param.GetAttrs();
        GradV.PutAll(0.0);
        
        for(int i = 0; i < NNodes; i++) {
                for(int j = 0; j < NNodes; j++) {
                        double Prod = ProdVV(i, j) - GetThetaLL(i, j, AId);
                        double Sq = SqVV(i, j) - GetSqThetaLL(i, j, AId);

                        for(int p = 0; p < 4; p++) {
                                int Ai = p / 2;
                                int Aj = p % 2;
                                double Prob = GetProbPhi(i, j, AId, Ai, Aj);
                                if(Graph->IsEdge(i, j)) {
                                        GradV[p] += Prob / CurMtx.At(p);
                                } else {
                                        GradV[p] -= Prob * exp(Prod);
                                        GradV[p] -= Prob * exp(Sq) * CurMtx.At(p);
                                }
                        }
                }
        }
}


const void TMAGFitBern::GradApxAffMtx(const int& AId, const TFltVV& ProdVV, const TFltVV& SqVV, const TMAGAffMtx& CurMtx, TFltV& GradV) {
        const int NNodes = Param.GetNodes();
        // const int NAttrs = Param.GetAttrs();
        // const int NSq = NNodes * (NNodes - 1);
        GradV.PutAll(0.0);

        TFltV LogSumV;
        for(int p = 0; p < 4; p++) {
                int Ai = p / 2;
                int Aj = p % 2;
                LogSumV.Gen(NNodes * 4, 0);

                for(int i = 0; i < NNodes; i++) {
                        const double LProd = ProdVV(i, 0) - GetAvgThetaLL(i, i, AId, true, false);
                        const double LSq = SqVV(i, 0) - GetAvgSqThetaLL(i, i, AId, true, false);
                        const double RProd = ProdVV(i, 1) - GetAvgThetaLL(i, i, AId, false, true);
                        const double RSq = SqVV(i, 1) - GetAvgSqThetaLL(i, i, AId, false, true);

                        LogSumV.Add(LProd + log(GetProbMu(i, i, AId, Ai, Aj, true, false)));
                        LogSumV.Add(LSq + log(GetProbMu(i, i, AId, Ai, Aj, true, false)) + log(CurMtx.At(p)));
                        LogSumV.Add(RProd + log(GetProbMu(i, i, AId, Ai, Aj, false, true)));
                        LogSumV.Add(RSq + log(GetProbMu(i, i, AId, Ai, Aj, false, true)) + log(CurMtx.At(p)));
                }
                double LogSum = LogSumExp(LogSumV);
                GradV[p] -= (NNodes - 1) * 0.5 * exp(LogSum);
        }
        
        for(TNGraph::TEdgeI EI = Graph->BegEI(); EI < Graph->EndEI(); EI++) {
                const int NId1 = EI.GetSrcNId();
                const int NId2 = EI.GetDstNId();
                const double ProdOne = GetProdLinWeight(NId1, NId2) - GetThetaLL(NId1, NId2, AId);
                const double SqOne = GetProdSqWeight(NId1, NId2) - GetSqThetaLL(NId1, NId2, AId);

                for(int p = 0; p < 4; p++) {
                        int Ai = p / 2;
                        int Aj = p % 2;
                        double Prob = GetProbPhi(NId1, NId2, AId, Ai, Aj);
                        GradV[p] += Prob / CurMtx.At(p);
                        GradV[p] += Prob * exp(ProdOne);
                        GradV[p] += Prob * exp(SqOne) * CurMtx.At(p);
                }
        }

#if 0
        const double Prod = ProdVV(0, 0) - GetAvgThetaLL(0, 0, AId, false, false);
        const double Sq = SqVV(0, 0) - GetAvgSqThetaLL(0, 0, AId, false, false);
        for(int p = 0; p < 4; p++) {
                int Ai = p / 2;
                int Aj = p % 2;
                GradV[p] -= NSq * exp(Prod) * GetProbMu(0, 0, AId, Ai, Aj, false, false);
                GradV[p] -= NSq * exp(Sq) * GetProbMu(0, 0, AId, Ai, Aj, false, false) * CurMtx.At(p);
        }

        for(TNGraph::TEdgeI EI = Graph->BegEI(); EI < Graph->EndEI(); EI++) {
                const int NId1 = EI.GetSrcNId();
                const int NId2 = EI.GetDstNId();
                const double ProdOne = GetProdLinWeight(NId1, NId2) - GetThetaLL(NId1, NId2, AId);
                const double SqOne = GetProdSqWeight(NId1, NId2) - GetSqThetaLL(NId1, NId2, AId);

                for(int p = 0; p < 4; p++) {
                        int Ai = p / 2;
                        int Aj = p % 2;
                        double Prob = GetProbPhi(NId1, NId2, AId, Ai, Aj);
//                      GradV[p] += Prob / CurMtx.At(p);
//                      GradV[p] += Prob * exp(ProdOne);
//                      GradV[p] += Prob * exp(SqOne) * CurMtx.At(p);
                }
        }
#endif
}


const double TMAGFitBern::UpdateAffMtx(const int& AId, const double& LrnRate, const double& MaxGrad, const double& Lambda, TFltVV& ProdVV, TFltVV& SqVV, TMAGAffMtx& NewMtx) {
        double Delta = 0.0;
        // const int NNodes = Param.GetNodes();
        // const int NAttrs = Param.GetAttrs();
        TMAGAffMtx AffMtx = Param.GetMtx(AId);

        TFltV GradV(4);
        TFltV HessV(4);
        if(MSpeedUp) {
                GradApxAffMtx(AId, ProdVV, SqVV, NewMtx, GradV);
        } else {
                GradAffMtx(AId, ProdVV, SqVV, NewMtx, GradV);
        }

        double Ratio = 1.0;
        for(int p = 0; p < 4; p++) {
                if(fabs(Ratio * LrnRate * GradV[p]) > MaxGrad) {
                        Ratio = MaxGrad / fabs(LrnRate * GradV[p]);
                }
        }

        for(int p = 0; p < 4; p++) {
                GradV[p] *= (Ratio * LrnRate);
                NewMtx.At(p) = AffMtx.At(p) + GradV[p];
//              if(NewMtx.At(p) > 0.9999) {  NewMtx.At(p) = 0.9999;  }
                if(NewMtx.At(p) < 0.0001) {  NewMtx.At(p) = 0.0001;  }
        }

        printf("      [Attr = %d]\n", AId);
    printf("        %s  + [%f, %f; %f %f]  ----->  %s\n", (AffMtx.GetMtxStr()).GetCStr(), double(GradV[0]), double(GradV[1]), double(GradV[2]), double(GradV[3]), (NewMtx.GetMtxStr()).GetCStr());
        
//      Param.SetMtx(AId, NewMtx);
        return Delta;
}


void TMAGFitBern::NormalizeAffMtxV(TMAGAffMtxV& MtxV, const bool UseMu) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = MtxV.Len();
        TFltV MuV = GetMuV();
        double Product = 1.0, ExpEdge = NNodes * (NNodes - 1);
        
        TFltV SumV(NAttrs), EdgeSumV(NAttrs);
        SumV.PutAll(0.0);       EdgeSumV.PutAll(0.0);
        for(int l = 0; l < NAttrs; l++) {
                double Mu = (UseMu) ? double(MuV[l]) : (AvgPhiV[l] / double(NNodes));
                EdgeSumV[l] += Mu * Mu * MtxV[l].At(0, 0);
                EdgeSumV[l] += Mu * (1.0-Mu) * MtxV[l].At(0, 1);
                EdgeSumV[l] += Mu * (1.0-Mu) * MtxV[l].At(1, 0);
                EdgeSumV[l] += (1.0-Mu) * (1.0-Mu) * MtxV[l].At(1, 1);
                SumV[l] = SumV[l] + MtxV[l].At(0, 0);
                SumV[l] = SumV[l] + MtxV[l].At(0, 1);
                SumV[l] = SumV[l] + MtxV[l].At(1, 0);
                SumV[l] = SumV[l] + MtxV[l].At(1, 1);
                Product *= SumV[l];
                ExpEdge *= EdgeSumV[l];
        }
        ExpEdge = Graph->GetEdges() / ExpEdge;
        NormConst *= Product;
//      NormConst = ExpEdge;
        Product = 1.0;
//      Product = pow(Product * ExpEdge, 1.0 / double(NAttrs));
        
        for(int l = 0; l < NAttrs; l++) {
                for(int p = 0; p < 4; p++) {
                        MtxV[l].At(p) = MtxV[l].At(p) * Product / SumV[l];
//                      MtxV[l].At(p) = MtxV[l].At(p) * Product / MtxV[l].At(0, 0);
//                      MtxV[l].At(p) = MtxV[l].At(p) * Product;
//                      if(MtxV[l].At(p) > 0.9999) {  MtxV[l].At(p) = 0.9999;  }
//                      if(MtxV[l].At(p) < 0.0001) {  MtxV[l].At(p) = 0.0001;  }
                }
        }
}

void TMAGFitBern::UnNormalizeAffMtxV(TMAGAffMtxV& MtxV, const bool UseMu) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = MtxV.Len();
        TFltIntPrV MaxEntV(NAttrs);
        TFltV MuV = GetMuV();
        NormalizeAffMtxV(MtxV, UseMu);
        
        double ExpEdge = NNodes * (NNodes - 1);
        for(int l = 0; l < NAttrs; l++) {
                double Mu = MuV[l];
                double EdgeSum = Mu * Mu * MtxV[l].At(0, 0);
                EdgeSum += Mu * (1.0-Mu) * MtxV[l].At(0, 1);
                EdgeSum += Mu * (1.0-Mu) * MtxV[l].At(1, 0);
                EdgeSum += (1.0-Mu) * (1.0-Mu) * MtxV[l].At(1, 1);
                ExpEdge *= EdgeSum;
        }
        NormConst = double(Graph->GetEdges()) / ExpEdge;
//      NormConst *= ExpEdge;
        
        for(int l = 0; l < NAttrs; l++) {
                MaxEntV[l] = TFltIntPr(-1, l);
                for(int p = 0; p < 4; p++) {
                        if(MaxEntV[l].Val1 < MtxV[l].At(p)) {  MaxEntV[l].Val1 = MtxV[l].At(p);  }
                }
        }
        MaxEntV.Sort(false);

        for(int l = 0; l < NAttrs; l++) {
                int CurId = MaxEntV[l].Val2;
                double Factor = pow(NormConst, 1.0 / double(NAttrs - l));
                double MaxFactor = 0.9999 / MaxEntV[l].Val1;
                Factor = (Factor > MaxFactor) ? MaxFactor : Factor;
                NormConst = NormConst / Factor;

                for(int p = 0; p < 4; p++) {
                        MtxV[CurId].At(p) = MtxV[CurId].At(p) * Factor;
                }
        }
}

const void TMAGFitBern::PrepareUpdateAffMtx(TFltVV& ProdVV, TFltVV& SqVV) {
        const int NNodes = Param.GetNodes();
        ProdVV.Gen(NNodes, NNodes);
        SqVV.Gen(NNodes, NNodes);

        for(int i = 0; i < NNodes; i++) {
                for(int j = 0; j < NNodes; j++) {
                        ProdVV(i, j) = GetProdLinWeight(i, j);
                        SqVV(i, j) = GetProdSqWeight(i, j);
                }
        }
}

const void TMAGFitBern::PrepareUpdateApxAffMtx(TFltVV& ProdVV, TFltVV& SqVV) {
        const int NNodes = Param.GetNodes();
        ProdVV.Gen(NNodes, 2);
        SqVV.Gen(NNodes, 2);

        for(int i = 0; i < NNodes; i++) {
                ProdVV(i, 0) = GetAvgProdLinWeight(i, i, true, false);
                ProdVV(i, 1) = GetAvgProdLinWeight(i, i, false, true);
                SqVV(i, 0) = GetAvgProdSqWeight(i, i, true, false);
                SqVV(i, 1) = GetAvgProdSqWeight(i, i, false, true);
        }
}
        
const double TMAGFitBern::UpdateAffMtxV(const int& GradIter, const double& LrnRate, const double& MaxGrad, const double& Lambda, const int& NReal) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        const TMAGNodeBern DistParam = Param.GetNodeAttr();
        const TFltV MuV = DistParam.GetMuV();
        double Delta = 0.0;
        double DecLrnRate = LrnRate, DecMaxGrad = MaxGrad;
        
        TFltVV ProdVV(NNodes, NNodes), SqVV(NNodes, NNodes);
        TMAGAffMtxV NewMtxV, OldMtxV;
        Param.GetMtxV(OldMtxV);
        Param.GetMtxV(NewMtxV);

        for(int g = 0; g < GradIter; g++) {
                if(MSpeedUp) {
                        PrepareUpdateApxAffMtx(ProdVV, SqVV);
                } else {
                        PrepareUpdateAffMtx(ProdVV, SqVV);
                }

                printf("    [Grad step = %d]\n", (g+1));
//              for(int l = 0; l < NAttrs; l++) {
                for(int l = NReal; l < NAttrs; l++) {
                        UpdateAffMtx(l, DecLrnRate, DecMaxGrad, Lambda, ProdVV, SqVV, NewMtxV[l]);
                        Param.SetMtxV(NewMtxV);
                }
                DecLrnRate *= 0.97;
                DecMaxGrad *= 0.97;

                printf("\n");
                NormalizeAffMtxV(NewMtxV, true);
                Param.SetMtxV(NewMtxV);
        }
        NormalizeAffMtxV(NewMtxV, true);
        
        printf( "\nFinal\n");
        for(int l = 0; l < NAttrs; l++) {
                printf("    [");
                for(int p = 0; p < 4; p++) {
//                      NewMtxV[l].At(p) = NewMtxV[l].At(p) * Product / SumV[l];
                        Delta += fabs(OldMtxV[l].At(p) - NewMtxV[l].At(p));
                        printf(" %.4f ", double(NewMtxV[l].At(p)));
                }
                printf("]\n");
        }
        Param.SetMtxV(NewMtxV);
        ProdVV.Clr();           SqVV.Clr();
        return Delta;
}

void TMAGFitBern::DoMStep(const int& GradIter, const double& LrnRate, const double& MaxGrad, const double& Lambda, const int& NReal) {
        // const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        double MuDelta = 0.0, AffMtxDelta = 0.0;

        TExeTm ExeTm;
        
        printf("\n");
        AvgPhiV.Gen(NAttrs);    AvgPhiV.PutAll(0.0);
        for(int l = 0; l < NAttrs; l++) {
//              printf("    [Attr = %d]\n", l);
                MuDelta += UpdateMu(l);
        }
        printf("\n");

        printf("  == Update Theta\n");
        AffMtxDelta += UpdateAffMtxV(GradIter, LrnRate, MaxGrad, Lambda, NReal);
        printf("\n");
        printf("Elpased time = %s\n", ExeTm.GetTmStr());
        printf("\n");
}

void TMAGFitBern::DoEMAlg(const int& NStep, const int& NEstep, const int& NMstep, const double& LrnRate, const double& MaxGrad, const double& Lambda, const double& ReInit, const int& NReal) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TIntV IndexV;
        double LL;
  // double MuDist, MtxDist;

        MuHisV.Gen(NStep + 1, 0);
        MtxHisV.Gen(NStep + 1, 0);
        LLHisV.Gen(NStep + 1, 0);

        printf("--------------------------------------------\n");
        printf("Before EM Iteration\n");
        printf("--------------------------------------------\n");

        TMAGAffMtxV InitMtxV;
        TMAGNodeBern NodeAttr = Param.GetNodeAttr();
        Param.GetMtxV(InitMtxV);
        TFltV InitMuV = NodeAttr.GetMuV();
        
        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        if(! KnownVV(i, l)) {
                                PhiVV.At(i, l) = TMAGNodeBern::Rnd.GetUniDev();
                        }
                }
        }
        
        if(Debug) {
                double LL = ComputeApxLL();
                MuHisV.Add(InitMuV);
                MtxHisV.Add(InitMtxV);
                LLHisV.Add(LL);
        }

        NormalizeAffMtxV(InitMtxV, true);
        Param.SetMtxV(InitMtxV);

        for(int n = 0; n < NStep; n++) {
                printf("--------------------------------------------\n");
                printf("EM Iteration : %d\n", (n+1));
                printf("--------------------------------------------\n");
                
                NodeAttr = Param.GetNodeAttr();
                for(int i = 0; i < NNodes; i++) {
                        for(int l = 0; l < NAttrs; l++) {
                                if(!KnownVV(i, l) && TMAGNodeBern::Rnd.GetUniDev() < ReInit) {
                                        PhiVV.At(i, l) = TMAGNodeBern::Rnd.GetUniDev();
                                }
                        }
                }
                DoEStep(InitMuV, NEstep, LL, Lambda);
                Param.GetMtxV(InitMtxV);
//              NormalizeAffMtxV(InitMtxV);
                Param.SetMtxV(InitMtxV);
                DoMStep(NMstep, LrnRate, MaxGrad, Lambda, NReal);

                printf("\n");

                if(Debug) {
                        double LL = ComputeApxLL();
                        MuHisV.Add(InitMuV);
                        MtxHisV.Add(InitMtxV);
                        LLHisV.Add(LL);
                        printf("    ApxLL = %.2f (Const = %f)\n", LL, double(NormConst));
                }

        }
        Param.GetMtxV(InitMtxV);
        UnNormalizeAffMtxV(InitMtxV, true);
        Param.SetMtxV(InitMtxV);
}

void TMAGFitBern::MakeCCDF(const TFltPrV& RawV, TFltPrV& CcdfV) {
        double Total = 0.0;
        CcdfV.Gen(RawV.Len(), 0);

        for(int i = 0; i < RawV.Len(); i++) {
                if(RawV[i].Val2 <= 0) {  continue;  }
                Total += RawV[i].Val2;
                CcdfV.Add(RawV[i]);
                IAssert(RawV[i].Val2 > 0);
        }
        for(int i = 1; i < CcdfV.Len(); i++) {
                CcdfV[i].Val2 += CcdfV[i-1].Val2;
        }

        for(int i = CcdfV.Len() - 1; i > 0; i--) {
                CcdfV[i].Val2 = (Total - CcdfV[i-1].Val2) ;
                if(CcdfV[i].Val2 <= 0) {  printf("CCDF = %f\n", double(CcdfV[i].Val2));}
                IAssert(CcdfV[i].Val2 > 0);
        }
        CcdfV[0].Val2 = Total;
//      CcdfV[0].Val2 = 1.0;
}


void TMAGFitBern::PlotProperties(const TStr& FNm) {
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TMAGParam<TMAGNodeBern> MAGGen(NNodes, NAttrs);
        TMAGNodeBern MAGNode = Param.GetNodeAttr();
        MAGGen.SetNodeAttr(MAGNode);
        TMAGAffMtxV MtxV;       Param.GetMtxV(MtxV);
        MAGGen.SetMtxV(MtxV);
        
        PNGraph TrG = new TNGraph;
        *TrG = *Graph;

        TIntVV AttrVV(NNodes, NAttrs);
        for(int i = 0; i < NNodes; i++) {
                for(int j = 0; j < NAttrs; j++) {
                        if(PhiVV(i, j) > TMAGNodeBern::Rnd.GetUniDev()) AttrVV(i, j) = 0;
                        else AttrVV(i, j) = 1;
                }
        }
        PNGraph MAG = MAGGen.GenMAG(AttrVV, true, 10000);
//      PNGraph MAG = MAGGen.GenAttrMAG(AttrVV, true, 10000);
        printf("%d edges created for MAG...\n", MAG->GetEdges());
        
        TSnap::DelZeroDegNodes(TrG);
        TSnap::DelZeroDegNodes(MAG);

        TGStatVec GS(tmuNodes, TFSet() | gsdInDeg | gsdOutDeg | gsdWcc | gsdHops | gsdClustCf | gsdSngVec | gsdSngVal | gsdTriadPart);
        
    TGnuPlot InDegP(FNm + "-InDeg"), OutDegP(FNm + "-OutDeg"), SvalP(FNm + "-Sval"), SvecP(FNm + "-Svec"), WccP(FNm + "-Wcc"), HopP(FNm + "-Hop"), TriadP(FNm + "-Triad"), CcfP(FNm + "-Ccf");;

    InDegP.SetXYLabel("Degree", "# of nodes");
    OutDegP.SetXYLabel("Degree", "# of nodes");
    SvalP.SetXYLabel("Rank", "Singular value");
    SvecP.SetXYLabel("Rank", "Primary SngVec component");
    WccP.SetXYLabel("Size of component", "# of components");
    CcfP.SetXYLabel("Degree", "Clustering coefficient");
    HopP.SetXYLabel("Hops", "# of node pairs");
    TriadP.SetXYLabel("# of triads", "# of participating nodes");

    InDegP.SetScale(gpsLog10XY);    InDegP.AddCmd("set key top right");
    OutDegP.SetScale(gpsLog10XY);   OutDegP.AddCmd("set key top right");
    SvalP.SetScale(gpsLog10XY);     SvalP.AddCmd("set key top right");
    SvecP.SetScale(gpsLog10XY);     SvecP.AddCmd("set key top right");
    CcfP.SetScale(gpsLog10XY);      CcfP.AddCmd("set key top right");
    HopP.SetScale(gpsLog10XY);      HopP.AddCmd("set key top right");
    TriadP.SetScale(gpsLog10XY);    TriadP.AddCmd("set key top right");
        InDegP.ShowGrid(false);
        OutDegP.ShowGrid(false);
        SvalP.ShowGrid(false);
        SvecP.ShowGrid(false);
        CcfP.ShowGrid(false);
        HopP.ShowGrid(false);
        TriadP.ShowGrid(false);
        
        const TStr Style[2] = {"lt 1 lw 3 lc rgb 'black'", "lt 2 lw 3 lc rgb 'red'"};
        const TStr Name[2] = {"Real", "MAG"};
        GS.Add(Graph, TSecTm(1), "Real Graph");
        GS.Add(MAG, TSecTm(2), "MAG");

        TFltPrV InDegV, OutDegV, SvalV, SvecV, HopV, WccV, CcfV, TriadV;
        for(int i = 0; i < GS.Len(); i++) {
                MakeCCDF(GS.At(i)->GetDistr(gsdInDeg), InDegV);
                MakeCCDF(GS.At(i)->GetDistr(gsdOutDeg), OutDegV);
                SvalV = GS.At(i)->GetDistr(gsdSngVal);
                SvecV = GS.At(i)->GetDistr(gsdSngVec);
                MakeCCDF(GS.At(i)->GetDistr(gsdClustCf), CcfV);
                HopV = GS.At(i)->GetDistr(gsdHops);
                MakeCCDF(GS.At(i)->GetDistr(gsdTriadPart), TriadV);

                InDegP.AddPlot(InDegV, gpwLines, Name[i], Style[i]);
                OutDegP.AddPlot(OutDegV, gpwLines, Name[i], Style[i]);
                SvalP.AddPlot(SvalV, gpwLines, Name[i], Style[i]);
                SvecP.AddPlot(SvecV, gpwLines, Name[i], Style[i]);
                CcfP.AddPlot(CcfV, gpwLines, Name[i], Style[i]);
                HopP.AddPlot(HopV, gpwLines, Name[i], Style[i]);
                TriadP.AddPlot(TriadV, gpwLines, Name[i], Style[i]);
        }

        InDegP.SaveEps(30);
        OutDegP.SaveEps(30);
        SvalP.SaveEps(30);
        SvecP.SaveEps(30);
        CcfP.SaveEps(30);
        HopP.SaveEps(30);
        TriadP.SaveEps(30);
}

void TMAGFitBern::CountAttr(TFltV& EstMuV) const {
        const int NNodes = PhiVV.GetXDim();
        const int NAttrs = PhiVV.GetYDim();
        EstMuV.Gen(NAttrs);
        EstMuV.PutAll(0.0);

        for(int l = 0; l < NAttrs; l++) {
                for(int i = 0; i < NNodes; i++) {
                        EstMuV[l] = EstMuV[l] + PhiVV(i, l);
                }
                EstMuV[l] = EstMuV[l] / double(NNodes);
        }
}

void TMAGFitBern::SortAttrOrdering(const TFltV& TrueMuV, TIntV& IndexV) const {
        const int NAttrs = TrueMuV.Len();
        // const int NNodes = PhiVV.GetXDim();
        TFltV EstMuV, SortedTrueMuV, SortedEstMuV, TrueIdxV, EstIdxV;
        IndexV.Gen(NAttrs);
        TrueIdxV.Gen(NAttrs);
        EstIdxV.Gen(NAttrs);

        for(int l = 0; l < NAttrs; l++) {
                TrueIdxV[l] = l;
                EstIdxV[l] = l;
        }
        
        CountAttr(EstMuV);
        SortedTrueMuV = TrueMuV;
        SortedEstMuV = EstMuV;
        for(int i = 0; i < NAttrs; i++) {
                if(SortedTrueMuV[i] > 0.5) {  SortedTrueMuV[i] = 1.0 - SortedTrueMuV[i];  }
                if(SortedEstMuV[i] > 0.5) {  SortedEstMuV[i] = 1.0 - SortedEstMuV[i];  }
        }

        for(int i = 0; i < NAttrs; i++) {
                for(int j = i+1; j < NAttrs; j++) {
                        if(SortedTrueMuV[i] < SortedTrueMuV[j]) {
                                SortedTrueMuV.Swap(i, j);
                                TrueIdxV.Swap(i, j);
                        }
                        if(SortedEstMuV[i] < SortedEstMuV[j]) {
                                EstIdxV.Swap((int)SortedEstMuV[i], (int)SortedEstMuV[j]);
                                SortedEstMuV.Swap(i, j);
                        }
                }
        }

        for(int l = 0; l < NAttrs; l++) {
                IndexV[l] = (int)TrueIdxV[(int)EstIdxV[l]];
        }
}

const bool TMAGFitBern::NextPermutation(TIntV& IndexV) const {
        const int NAttrs = IndexV.Len();
        int Pos = NAttrs - 1;
        while(Pos > 0) {
                if(IndexV[Pos-1] < IndexV[Pos]) {
                        break;
                }
                Pos--;
        }
        if(Pos == 0) {
                return false;
        }

        int Val = NAttrs, NewPos = -1;
        for(int i = Pos; i < NAttrs; i++) {
                if(IndexV[i] > IndexV[Pos - 1] && IndexV[i] < Val) {
                        NewPos = i;
                        Val = IndexV[i];
                }
        }
        IndexV[NewPos] = IndexV[Pos - 1];
        IndexV[Pos - 1] = Val;

        TIntV SubIndexV;
    IndexV.GetSubValV(Pos, NAttrs - 1, SubIndexV);
        SubIndexV.Sort(true);
        for(int i = Pos; i < NAttrs; i++) {
                IndexV[i] = SubIndexV[i - Pos];
        }

        return true;
}

const double TMAGFitBern::ComputeJointOneLL(const TIntVV& AttrVV) const {
        double LL = 0.0;
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TMAGAffMtxV MtxV(NAttrs);       Param.GetMtxV(MtxV);
        const TMAGNodeBern NodeAttr = Param.GetNodeAttr();
        const TFltV MuV = NodeAttr.GetMuV();

        for(int l = 0; l < NAttrs; l++) {
                for(int i = 0; i < MtxV[l].Len(); i++) {
                        MtxV[l].At(i) = log(MtxV[l].At(i));
                }
        }

        for(int i = 0; i < NNodes; i++) {
                for(int l = 0; l < NAttrs; l++) {
                        if(AttrVV.At(i, l) == 0) {
                                LL += log(MuV[l]);
                        } else {
                                LL += log(1.0 - MuV[l]);
                        }
                }

                for(int j = 0; j < NNodes; j++) {
                        if(i == j) {  continue;  }

                        double ProbLL = 0.0;
                        for(int l = 0; l < NAttrs; l++) {
                                ProbLL += MtxV[l].At(AttrVV.At(i, l), AttrVV.At(j, l));
                        }

                        if(Graph->IsEdge(i, j)) {
                                LL += ProbLL;
                        } else {
                                LL += log(1-exp(ProbLL));
                        }
                }
        }

        return LL;
}

const double TMAGFitBern::ComputeJointAdjLL(const TIntVV& AttrVV) const {
        double LL = 0.0;
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TMAGAffMtxV MtxV(NAttrs);       Param.GetMtxV(MtxV);
        const TMAGNodeBern NodeAttr = Param.GetNodeAttr();
        const TFltV MuV = NodeAttr.GetMuV();

        for(int l = 0; l < NAttrs; l++) {
                for(int i = 0; i < MtxV[l].Len(); i++) {
                        MtxV[l].At(i) = log(MtxV[l].At(i));
                }
        }

        for(int i = 0; i < NNodes; i++) {
                for(int j = 0; j < NNodes; j++) {
                        if(i == j) {  continue;  }

                        double ProbLL = 0.0;
                        for(int l = 0; l < NAttrs; l++) {
                                ProbLL += MtxV[l].At(AttrVV.At(i, l), AttrVV.At(j, l));
                        }

                        if(Graph->IsEdge(i, j)) {
                                LL += ProbLL;
                        } else {
                                LL += log(1-exp(ProbLL));
                        }
                }
        }

        return LL;
}

const double TMAGFitBern::ComputeJointLL(int NSample) const {
        double LL = 0.0;
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();

        TRnd Rnd(2000);
        TIntVV AttrVV(NNodes, NAttrs);
        int count = 0;
        for(int s = 0; s < NSample; s++) {
                for(int i = 0; i < NNodes; i++) {
                        for(int l = 0; l < NAttrs; l++) {
                        
                                if(Rnd.GetUniDev() <= PhiVV(i, l)) {
                                        AttrVV.At(i, l) = 0;
                                } else {
                                        AttrVV.At(i, l) = 1;
                                }

                                if(PhiVV(i, l) > 0.05 && PhiVV(i, l) < 0.95) count++;
                        }
                }

                LL += ComputeJointOneLL(AttrVV);
        }
        AttrVV.Clr();

        return LL / double(NSample);
}

const double TMAGFitBern::ComputeApxLL() const {
        double LL = 0.0;
  // double LLSelf = 0.0;
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TMAGNodeBern NodeAttr = Param.GetNodeAttr();
        TFltV MuV = NodeAttr.GetMuV();
        TMAGAffMtxV LLMtxV(NAttrs);

        for(int l = 0; l < NAttrs; l++) {
                for(int i = 0; i < NNodes; i++) {
                        LL += PhiVV(i, l) * log(MuV[l]);
                        LL += (1.0 - PhiVV(i, l)) * log(1.0 - MuV[l]);
                        LL -= PhiVV(i, l) * log(PhiVV(i, l));
                        LL -= (1.0 - PhiVV(i, l)) * log(1.0 - PhiVV(i, l));
                }
                TMAGAffMtx Theta = Param.GetMtx(l);
                Theta.GetLLMtx(LLMtxV[l]);
        }

        for(int i = 0; i < NNodes; i++) {
                for(int j = 0; j < NNodes; j++) {
                        if(i == j) {  continue;  }

                        if(Graph->IsEdge(i, j)) {
                                for(int l = 0; l < NAttrs; l++) {
                                        LL += GetProbPhi(i, j, l, 0, 0) * LLMtxV[l].At(0, 0);
                                        LL += GetProbPhi(i, j, l, 0, 1) * LLMtxV[l].At(0, 1);
                                        LL += GetProbPhi(i, j, l, 1, 0) * LLMtxV[l].At(1, 0);
                                        LL += GetProbPhi(i, j, l, 1, 1) * LLMtxV[l].At(1, 1);
                                }
                                LL += log(NormConst);
                        } else {
                                LL += log(1-exp(GetProdLinWeight(i, j)));
                        }
                }
        }

        return LL;
}

const double TMAGFitBern::ComputeApxAdjLL() const {
        double LL = 0.0;
  // double LLSelf = 0.0;
        const int NNodes = Param.GetNodes();
        const int NAttrs = Param.GetAttrs();
        TMAGNodeBern NodeAttr = Param.GetNodeAttr();
        TFltV MuV = NodeAttr.GetMuV();
        MuV.PutAll(0.0);
        TMAGAffMtxV LLMtxV(NAttrs);
        double TotalEdge = 0.0;
        for(int l = 0; l < NAttrs; l++) {
                TMAGAffMtx Theta = Param.GetMtx(l);
                Theta.GetLLMtx(LLMtxV[l]);
        }

        for(int i = 0; i < NNodes; i++) {
                for(int j = 0; j < NNodes; j++) {
                        if(i == j) {  continue;  }

                        if(Graph->IsEdge(i, j)) {
                                for(int l = 0; l < NAttrs; l++) {
                                        LL += GetProbPhi(i, j, l, 0, 0) * LLMtxV[l].At(0, 0);
                                        LL += GetProbPhi(i, j, l, 0, 1) * LLMtxV[l].At(0, 1);
                                        LL += GetProbPhi(i, j, l, 1, 0) * LLMtxV[l].At(1, 0);
                                        LL += GetProbPhi(i, j, l, 1, 1) * LLMtxV[l].At(1, 1);
                                }
                        } else {
                                LL += log(1-exp(GetProdLinWeight(i, j)));
                        }

                        double TempLL = 1.0;
                        for(int l = 0; l < NAttrs; l++) {
                                int Ai = (double(PhiVV(i, l)) > 0.5) ? 0 : 1;
                                int Aj = (double(PhiVV(j, l)) > 0.5) ? 0 : 1;
                                TempLL *= Param.GetMtx(l).At(Ai, Aj);
                        }
                        if(TMAGNodeBern::Rnd.GetUniDev() < TempLL) {
                                TotalEdge += 1.0;
                        }
                }
        }

        return LL;
}

const double TMAGFitBern::ComputeMI(const TIntVV& AttrV, const int AId1, const int AId2) {
        const int NNodes = AttrV.GetXDim();
        double MI = 0.0;
        double Cor = 0.0;

        TFltVV Pxy(2,2);
        TFltV Px(2), Py(2);
        Pxy.PutAll(0.0);
        Px.PutAll(0.0);
        Py.PutAll(0.0);

        for(int i = 0; i < NNodes; i++) {
                int X = AttrV(i, AId1);
                int Y = AttrV(i, AId2);
                Pxy(X, Y) = Pxy(X, Y) + 1;
                Px[X] = Px[X] + 1;
                Py[Y] = Py[Y] + 1;
                Cor += double(X * Y);
        }

        for(int x = 0; x < 2; x++) {
                for(int y = 0; y < 2; y++) {
      MI += Pxy(x, y) / double(NNodes) * (log(Pxy(x, y).Val) - log(Px[x].Val) - log(Py[y].Val) + log((double)NNodes));
                }
        }
        
        return MI;
}

const double TMAGFitBern::ComputeMI(const TFltVV& AttrV, const int AId1, const int AId2) {
        const int NNodes = AttrV.GetXDim();
        double MI = 0.0;
        double Cor = 0.0;

        TFltVV Pxy(2,2);
        TFltV Px(2), Py(2);
        Pxy.PutAll(0.0);
        Px.PutAll(0.0);
        Py.PutAll(0.0);
        
        for(int i = 0; i < NNodes; i++) {
                double X = AttrV(i, AId1);
                double Y = AttrV(i, AId2);
                Pxy(0, 0) = Pxy(0, 0) + X * Y;
                Pxy(0, 1) = Pxy(0, 1) + X * (1 - Y);
                Pxy(1, 0) = Pxy(1, 0) + (1 - X) * Y;
                Pxy(1, 1) = (i+1) - Pxy(0, 0) - Pxy(0, 1) - Pxy(1, 0);
                Px[0] = Px[0] + X;
                Py[0] = Py[0] + Y;
                Cor += double((1-X) * (1-Y));
        }
        Px[1] = NNodes - Px[0];
        Py[1] = NNodes - Py[0];
        
        for(int x = 0; x < 2; x++) {
                for(int y = 0; y < 2; y++) {
                        MI += Pxy(x, y) / double(NNodes) * (log(Pxy(x, y)) - log(Px[x]) - log(Py[y]) + log(double(NNodes)));
                }
        }
        
        return MI;
}

const double TMAGFitBern::ComputeMI(const TIntVV& AttrV) {
        // const int NNodes = AttrV.GetXDim();
        const int NAttrs = AttrV.GetYDim();
        double MI = 0.0;

        for(int l = 0; l < NAttrs; l++) {
                for(int k = l+1; k < NAttrs; k++) {
                        MI += ComputeMI(AttrV, l, k);
                }
        }

        return MI;
}

const double TMAGFitBern::ComputeMI(const TFltVV& AttrV) {
        // const int NNodes = AttrV.GetXDim();
        const int NAttrs = AttrV.GetYDim();
        double MI = 0.0;

        for(int l = 0; l < NAttrs; l++) {
                for(int k = l+1; k < NAttrs; k++) {
                        MI += ComputeMI(AttrV, l, k);
                }
        }

        return MI;
}