SNAP Library 2.0, User Reference  2013-05-13 16:33:57
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Public Member Functions | Static Public Member Functions | Private Attributes | Friends
TKroneckerLL Class Reference

!!!!! MYUNGHWAN, CHECK! More...

#include <kronecker.h>

List of all members.

Public Member Functions

 TKroneckerLL ()
 TKroneckerLL (const PNGraph &GraphPt, const TFltV &ParamV, const double &PermPSwapNd=0.2)
 TKroneckerLL (const PNGraph &GraphPt, const TKronMtx &ParamMtx, const double &PermPSwapNd=0.2)
 TKroneckerLL (const PNGraph &GraphPt, const TKronMtx &ParamMtx, const TIntV &NodeIdPermV, const double &PermPSwapNd=0.2)
int GetNodes () const
int GetKronIters () const
PNGraph GetGraph () const
void SetGraph (const PNGraph &GraphPt)
const TKronMtx & GetProbMtx () const
const TKronMtx & GetLLMtx () const
int GetParams () const
int GetDim () const
void SetDebug (const bool Debug)
const TFltV & GetLLHist () const
const TVec< TKronMtx > & GetParamHist () const
bool IsObsNode (const int &NId) const
bool IsObsEdge (const int &NId1, const int &NId2) const
bool IsLatentNode (const int &NId) const
bool IsLatentEdge (const int &NId1, const int &NId2) const
void SetPerm (const char &PermId)
void SetOrderPerm ()
void SetRndPerm ()
void SetDegPerm ()
void SetBestDegPerm ()
 !!!!! MYUNGHWAN, CHECK!
void SetPerm (const TIntV &NodePermV)
void SetIPerm (const TIntV &Perm)
 !!!!! MYUNGHWAN, CHECK!
const TIntV & GetPermV () const
void AppendIsoNodes ()
void RestoreGraph (const bool RestoreNodes=true)
 !!!!! MYUNGHWAN, CHECK!
double GetFullGraphLL () const
double GetFullRowLL (int RowId) const
double GetFullColLL (int ColId) const
double GetEmptyGraphLL () const
double GetApxEmptyGraphLL () const
void InitLL (const TFltV &ParamV)
void InitLL (const TKronMtx &ParamMtx)
void InitLL (const PNGraph &GraphPt, const TKronMtx &ParamMtx)
double CalcGraphLL ()
double CalcApxGraphLL ()
double GetLL () const
double GetAbsErr () const
double NodeLLDelta (const int &NId) const
double SwapNodesLL (const int &NId1, const int &NId2)
bool SampleNextPerm (int &NId1, int &NId2)
double GetEmptyGraphDLL (const int &ParamId) const
double GetApxEmptyGraphDLL (const int &ParamId) const
const TFltV & CalcGraphDLL ()
const TFltV & CalcFullApxGraphDLL ()
const TFltV & CalcApxGraphDLL ()
double NodeDLLDelta (const int ParamId, const int &NId) const
void UpdateGraphDLL (const int &SwapNId1, const int &SwapNId2)
const TFltV & GetDLL () const
double GetDLL (const int &ParamId) const
void SampleGradient (const int &WarmUp, const int &NSamples, double &AvgLL, TFltV &GradV)
double GradDescent (const int &NIter, const double &LrnRate, double MnStep, double MxStep, const int &WarmUp, const int &NSamples)
double GradDescent2 (const int &NIter, const double &LrnRate, double MnStep, double MxStep, const int &WarmUp, const int &NSamples)
void SetRandomEdges (const int &NEdges, const bool isDir=true)
 !!!!! MYUNGHWAN, CHECK!
void MetroGibbsSampleSetup (const int &WarmUp)
void MetroGibbsSampleNext (const int &WarmUp, const bool DLLUpdate=false)
void RunEStep (const int &GibbsWarmUp, const int &WarmUp, const int &NSamples, TFltV &LLV, TVec< TFltV > &DLLV)
double RunMStep (const TFltV &LLV, const TVec< TFltV > &DLLV, const int &GradIter, const double &LrnRate, double MnStep, double MxStep)
void RunKronEM (const int &EMIter, const int &GradIter, double LrnRate, double MnStep, double MxStep, const int &GibbsWarmUp, const int &WarmUp, const int &NSamples, const TKronEMType &Type=kronNodeMiss, const int &NMissing=-1)
TFltV TestSamplePerm (const TStr &OutFNm, const int &WarmUp, const int &NSamples, const TKronMtx &TrueMtx, const bool &DoPlot=true)
TFltQu TestKronDescent (const bool &DoExact, const bool &TruePerm, double LearnRate, const int &WarmUp, const int &NSamples, const TKronMtx &TrueParam)
void GradDescentConvergence (const TStr &OutFNm, const TStr &Desc1, const bool &SamplePerm, const int &NIters, double LearnRate, const int &WarmUp, const int &NSamples, const int &AvgKGraphs, const TKronMtx &TrueParam)

Static Public Member Functions

static PKroneckerLL New ()
static PKroneckerLL New (const PNGraph &GraphPt, const TKronMtx &ParamMtx, const double &PermPSwapNd=0.1)
static PKroneckerLL New (const PNGraph &GraphPt, const TKronMtx &ParamMtx, const TIntV &NodeIdPermV, const double &PermPSwapNd=0.2)
static double CalcChainR2 (const TVec< TFltV > &ChainLLV)
static void ChainGelmapRubinPlot (const TVec< TFltV > &ChainLLV, const TStr &OutFNm, const TStr &Desc)
static void TestBicCriterion (const TStr &OutFNm, const TStr &Desc1, const PNGraph &G, const int &GradIters, double LearnRate, const int &WarmUp, const int &NSamples, const int &TrueN0)
static void TestGradDescent (const int &KronIters, const int &KiloSamples, const TStr &Permutation)

Private Attributes

TCRef CRef
PNGraph Graph
TInt Nodes
TInt KronIters
TFlt PermSwapNodeProb
TIntTrV GEdgeV
TIntTrV LEdgeV
TInt LSelfEdge
TIntV NodePerm
TIntV InvertPerm
TInt RealNodes
TInt RealEdges
TKronMtx ProbMtx
TKronMtx LLMtx
TFlt LogLike
TFltV GradV
TKronEMType EMType
TInt MissEdges
TBool DebugMode
TFltV LLV
TVec< TKronMtx > MtxV

Friends

class TPt< TKroneckerLL >

Detailed Description

!!!!! MYUNGHWAN, CHECK!

Definition at line 116 of file kronecker.h.

Constructor & Destructor Documentation

TKroneckerLL::TKroneckerLL ( ) [inline]

Definition at line 150 of file kronecker.h.

: Nodes(-1), KronIters(-1), PermSwapNodeProb(0.2), RealNodes(-1), RealEdges(-1), LogLike(TKronMtx::NInf), EMType(kronNodeMiss), MissEdges(-1), DebugMode(false) { }
TKroneckerLL::TKroneckerLL ( const PNGraph &  GraphPt,
const TFltV &  ParamV,
const double &  PermPSwapNd = 0.2 
)

Definition at line 783 of file kronecker.cpp.

                                                                                                : PermSwapNodeProb(PermPSwapNd) {
  InitLL(GraphPt, TKronMtx(ParamV));
}
TKroneckerLL::TKroneckerLL ( const PNGraph &  GraphPt,
const TKronMtx &  ParamMtx,
const double &  PermPSwapNd = 0.2 
)

Definition at line 787 of file kronecker.cpp.

                                                                                                      : PermSwapNodeProb(PermPSwapNd) {
  InitLL(GraphPt, ParamMtx);
}
TKroneckerLL::TKroneckerLL ( const PNGraph &  GraphPt,
const TKronMtx &  ParamMtx,
const TIntV &  NodeIdPermV,
const double &  PermPSwapNd = 0.2 
)

Definition at line 791 of file kronecker.cpp.

                                                                                                                                : PermSwapNodeProb(PermPSwapNd) {
  InitLL(GraphPt, ParamMtx);
  NodePerm = NodeIdPermV;
  SetIPerm(NodePerm);
}

Member Function Documentation

void TKroneckerLL::AppendIsoNodes ( )

Definition at line 914 of file kronecker.cpp.

                                  {
  Nodes = (int) pow((double)ProbMtx.GetDim(), KronIters);
  // add nodes until filling the Kronecker graph model
  for (int nid = Graph->GetNodes(); nid < Nodes; nid++) {
          Graph->AddNode(nid);
  }
}
const TFltV & TKroneckerLL::CalcApxGraphDLL ( )

Definition at line 1194 of file kronecker.cpp.

                                           {
  for (int ParamId = 0; ParamId < LLMtx.Len(); ParamId++) {
    double DLL = GetApxEmptyGraphDLL(ParamId);
    for (int nid = 0; nid < Nodes; nid++) {
      const TNGraph::TNodeI Node = Graph->GetNI(nid);
      const int SrcNId = NodePerm[nid];
      for (int e = 0; e < Node.GetOutDeg(); e++) {
        const int DstNId = NodePerm[Node.GetOutNId(e)];
        DLL = DLL - LLMtx.GetApxNoEdgeDLL(ParamId, SrcNId, DstNId, KronIters)
          + LLMtx.GetEdgeDLL(ParamId, SrcNId, DstNId, KronIters);
      }
    }
    GradV[ParamId] = DLL;
  }
  return GradV;
}
double TKroneckerLL::CalcApxGraphLL ( )

Definition at line 1037 of file kronecker.cpp.

                                    {
  LogLike = GetApxEmptyGraphLL(); // O(N_0)
  for (int nid = 0; nid < Nodes; nid++) {
    const TNGraph::TNodeI Node = Graph->GetNI(nid);
    const int SrcNId = NodePerm[nid];
    for (int e = 0; e < Node.GetOutDeg(); e++) {
      const int DstNId = NodePerm[Node.GetOutNId(e)];
      LogLike = LogLike - LLMtx.GetApxNoEdgeLL(SrcNId, DstNId, KronIters)
        + LLMtx.GetEdgeLL(SrcNId, DstNId, KronIters);
    }
  }
  return LogLike;
}
double TKroneckerLL::CalcChainR2 ( const TVec< TFltV > &  ChainLLV) [static]

Definition at line 1895 of file kronecker.cpp.

                                                            {
  const double J = ChainLLV.Len();
  const double K = ChainLLV[0].Len();
  TFltV AvgJV;    McMcGetAvgJ(ChainLLV, AvgJV);
  double AvgAvg;  McMcGetAvgAvg(AvgJV, AvgAvg);
  IAssert(AvgJV.Len() == ChainLLV.Len());
  double InChainVar=0, OutChainVar=0;
  // between chain var
  for (int j = 0; j < AvgJV.Len(); j++) {
    OutChainVar += TMath::Sqr(AvgJV[j] - AvgAvg); }
  OutChainVar = OutChainVar * (K/double(J-1));
  printf("*** %g chains of len %g\n", J, K);
  printf("  ** between chain var: %f\n", OutChainVar);
  //within chain variance
  for (int j = 0; j < AvgJV.Len(); j++) {
    const TFltV& ChainV = ChainLLV[j];
    for (int k = 0; k < ChainV.Len(); k++) {
      InChainVar += TMath::Sqr(ChainV[k] - AvgJV[j]); }
  }
  InChainVar = InChainVar * 1.0/double(J*(K-1));
  printf("  ** within chain var: %f\n", InChainVar);
  const double PostVar = (K-1)/K * InChainVar + 1.0/K * OutChainVar;
  printf("  ** posterior var: %f\n", PostVar);
  const double ScaleRed = sqrt(PostVar/InChainVar);
  printf("  ** scale reduction (< 1.2): %f\n\n", ScaleRed);
  return ScaleRed;
}
const TFltV & TKroneckerLL::CalcFullApxGraphDLL ( )

Definition at line 1176 of file kronecker.cpp.

                                               {
  for (int ParamId = 0; ParamId < LLMtx.Len(); ParamId++) {
    double DLL = 0.0;
    for (int NId1 = 0; NId1 < Nodes; NId1++) {
      for (int NId2 = 0; NId2 < Nodes; NId2++) {
        if (Graph->IsEdge(NId1, NId2)) {
          DLL += LLMtx.GetEdgeDLL(ParamId, NodePerm[NId1], NodePerm[NId2], KronIters);
        } else {
          DLL += LLMtx.GetApxNoEdgeDLL(ParamId, NodePerm[NId1], NodePerm[NId2], KronIters);
        }
      }
    }
    GradV[ParamId] = DLL;
  }
  return GradV;
}
const TFltV & TKroneckerLL::CalcGraphDLL ( )

Definition at line 1158 of file kronecker.cpp.

                                        {
  for (int ParamId = 0; ParamId < LLMtx.Len(); ParamId++) {
    double DLL = 0.0;
    for (int NId1 = 0; NId1 < Nodes; NId1++) {
      for (int NId2 = 0; NId2 < Nodes; NId2++) {
        if (Graph->IsEdge(NId1, NId2)) {
          DLL += LLMtx.GetEdgeDLL(ParamId, NodePerm[NId1], NodePerm[NId2], KronIters);
        } else {
          DLL += LLMtx.GetNoEdgeDLL(ParamId, NodePerm[NId1], NodePerm[NId2], KronIters);
        }
      }
    }
    GradV[ParamId] = DLL;
  }
  return GradV;
}
double TKroneckerLL::CalcGraphLL ( )

Definition at line 1022 of file kronecker.cpp.

                                 {
  LogLike = GetEmptyGraphLL(); // takes O(N^2)
  for (int nid = 0; nid < Nodes; nid++) {
    const TNGraph::TNodeI Node = Graph->GetNI(nid);
    const int SrcNId = NodePerm[nid];
    for (int e = 0; e < Node.GetOutDeg(); e++) {
      const int DstNId = NodePerm[Node.GetOutNId(e)];
      LogLike = LogLike - LLMtx.GetNoEdgeLL(SrcNId, DstNId, KronIters)
        + LLMtx.GetEdgeLL(SrcNId, DstNId, KronIters);
    }
  }
  return LogLike;
}
void TKroneckerLL::ChainGelmapRubinPlot ( const TVec< TFltV > &  ChainLLV,
const TStr &  OutFNm,
const TStr &  Desc 
) [static]

Definition at line 1924 of file kronecker.cpp.

                                                                                                         {
  TFltPrV LenR2V; // how does potential scale reduction chainge with chain length
  TVec<TFltV> SmallLLV(ChainLLV.Len());
  const int K = ChainLLV[0].Len();
  const int Buckets=1000;
  const int BucketSz = K/Buckets;
  for (int b = 1; b < Buckets; b++) {
    const int End = TMath::Mn(BucketSz*b, K-1);
    for (int c = 0; c < ChainLLV.Len(); c++) {
      ChainLLV[c].GetSubValV(0, End, SmallLLV[c]); }
    LenR2V.Add(TFltPr(End, TKroneckerLL::CalcChainR2(SmallLLV)));
  }
  LenR2V.Add(TFltPr(K, TKroneckerLL::CalcChainR2(ChainLLV)));
  TGnuPlot::PlotValV(LenR2V, TStr::Fmt("gelman-%s", OutFNm.CStr()), TStr::Fmt("%s. %d chains of len %d. BucketSz: %d.",
    Desc.CStr(), ChainLLV.Len(), ChainLLV[0].Len(), BucketSz), "Chain length", "Potential scale reduction");
}
double TKroneckerLL::GetAbsErr ( ) const [inline]

Definition at line 205 of file kronecker.h.

{ return fabs(pow((double)Graph->GetEdges(), 1.0/double(KronIters)) - ProbMtx.GetMtxSum()); }
double TKroneckerLL::GetApxEmptyGraphDLL ( const int &  ParamId) const

Definition at line 1147 of file kronecker.cpp.

                                                                 {
  double Sum=0.0, SumSq=0.0;
  for (int i = 0; i < ProbMtx.Len(); i++) {
    Sum += ProbMtx.At(i);
    SumSq += TMath::Sqr(ProbMtx.At(i));
  }
  // d/dx -sum(x_i) - 0.5sum(x_i^2) = d/dx sum(theta)^k - 0.5 sum(theta^2)^k
  return -KronIters*pow(Sum, KronIters-1) - KronIters*pow(SumSq, KronIters-1)*ProbMtx.At(ParamId);
}
double TKroneckerLL::GetApxEmptyGraphLL ( ) const

Definition at line 987 of file kronecker.cpp.

                                              {
  double Sum=0.0, SumSq=0.0;
  for (int i = 0; i < ProbMtx.Len(); i++) {
    Sum += ProbMtx.At(i);
    SumSq += TMath::Sqr(ProbMtx.At(i));
  }
  return -pow(Sum, KronIters) - 0.5*pow(SumSq, KronIters);
}
int TKroneckerLL::GetDim ( ) const [inline]

Definition at line 165 of file kronecker.h.

{ return ProbMtx.GetDim(); }
const TFltV& TKroneckerLL::GetDLL ( ) const [inline]

Definition at line 218 of file kronecker.h.

{ return GradV; }
double TKroneckerLL::GetDLL ( const int &  ParamId) const [inline]

Definition at line 219 of file kronecker.h.

{ return GradV[ParamId]; }
double TKroneckerLL::GetEmptyGraphDLL ( const int &  ParamId) const

Definition at line 1136 of file kronecker.cpp.

                                                              {
  double DLL = 0.0;
  for (int NId1 = 0; NId1 < Nodes; NId1++) {
    for (int NId2 = 0; NId2 < Nodes; NId2++) {
      DLL += LLMtx.GetNoEdgeDLL(ParamId, NodePerm[NId1], NodePerm[NId2], KronIters);
    }
  }
  return DLL;
}
double TKroneckerLL::GetEmptyGraphLL ( ) const

Definition at line 976 of file kronecker.cpp.

                                           {
  double LL = 0;
  for (int NId1 = 0; NId1 < LLMtx.GetNodes(KronIters); NId1++) {
    for (int NId2 = 0; NId2 < LLMtx.GetNodes(KronIters); NId2++) {
      LL = LL + LLMtx.GetNoEdgeLL(NId1, NId2, KronIters);
    }
  }
  return LL;
}
double TKroneckerLL::GetFullColLL ( int  ColId) const

Definition at line 966 of file kronecker.cpp.

                                                 {
  double ColLL = 0.0;
  const int MtxDim = LLMtx.GetDim();
  for (int level = 0; level < KronIters; level++) {
    ColLL += LLMtx.GetColSum(ColId % MtxDim);
    ColId /= MtxDim;
  }
  return ColLL;
}
double TKroneckerLL::GetFullGraphLL ( ) const

Definition at line 943 of file kronecker.cpp.

                                          {
  // the number of times a seed matrix element appears in
  // the full kronecker adjacency matrix after KronIter
  // kronecker multiplications
  double ElemCnt = 1;
  const double dim = LLMtx.GetDim();
  // count number of times x appears in the full kronecker matrix
  for (int i = 1; i < KronIters; i++) {
    ElemCnt = dim*dim*ElemCnt + TMath::Power(dim, 2*i);
  }
  return ElemCnt * LLMtx.GetMtxSum();
}
double TKroneckerLL::GetFullRowLL ( int  RowId) const

Definition at line 956 of file kronecker.cpp.

                                                 {
  double RowLL = 0.0;
  const int MtxDim = LLMtx.GetDim();
  for (int level = 0; level < KronIters; level++) {
    RowLL += LLMtx.GetRowSum(RowId % MtxDim);
    RowId /= MtxDim;
  }
  return RowLL;
}
PNGraph TKroneckerLL::GetGraph ( ) const [inline]

Definition at line 160 of file kronecker.h.

{ return Graph; }
int TKroneckerLL::GetKronIters ( ) const [inline]

Definition at line 159 of file kronecker.h.

{ return KronIters; }
double TKroneckerLL::GetLL ( ) const [inline]

Definition at line 204 of file kronecker.h.

{ return LogLike; }
const TFltV& TKroneckerLL::GetLLHist ( ) const [inline]

Definition at line 168 of file kronecker.h.

{ return LLV; }
const TKronMtx& TKroneckerLL::GetLLMtx ( ) const [inline]

Definition at line 163 of file kronecker.h.

{ return LLMtx; }
int TKroneckerLL::GetNodes ( ) const [inline]

Definition at line 158 of file kronecker.h.

{ return Nodes; }
const TVec<TKronMtx>& TKroneckerLL::GetParamHist ( ) const [inline]

Definition at line 169 of file kronecker.h.

{ return MtxV; }
int TKroneckerLL::GetParams ( ) const [inline]

Definition at line 164 of file kronecker.h.

{ return ProbMtx.Len(); }
const TIntV& TKroneckerLL::GetPermV ( ) const [inline]

Definition at line 185 of file kronecker.h.

{ return NodePerm; }
const TKronMtx& TKroneckerLL::GetProbMtx ( ) const [inline]

Definition at line 162 of file kronecker.h.

{ return ProbMtx; }
double TKroneckerLL::GradDescent ( const int &  NIter,
const double &  LrnRate,
double  MnStep,
double  MxStep,
const int &  WarmUp,
const int &  NSamples 
)

!!!!! MYUNGHWAN, CHECK!

!!!!! MYUNGHWAN, CHECK!

Definition at line 1299 of file kronecker.cpp.

                                                                                                                                              {
  printf("\n----------------------------------------------------------------------\n");
  printf("Fitting graph on %d nodes, %d edges\n", Graph->GetNodes(), Graph->GetEdges());
  printf("Kron iters:  %d (== %d nodes)\n\n", KronIters(), ProbMtx.GetNodes(KronIters()));
  TExeTm IterTm, TotalTm;
  double OldLL=-1e10, CurLL=0;
  const double EZero = pow((double) Graph->GetEdges(), 1.0/double(KronIters));
  TFltV CurGradV, LearnRateV(GetParams()), LastStep(GetParams());
  LearnRateV.PutAll(LrnRate);
  TKronMtx NewProbMtx = ProbMtx;

  if(DebugMode) {  
          LLV.Gen(NIter, 0);
          MtxV.Gen(NIter, 0);
  }

  for (int Iter = 0; Iter < NIter; Iter++) {
    printf("%03d] ", Iter);
    SampleGradient(WarmUp, NSamples, CurLL, CurGradV);
    for (int p = 0; p < GetParams(); p++) {
      LearnRateV[p] *= 0.95;
      if (Iter < 1) {
        while (fabs(LearnRateV[p]*CurGradV[p]) > MxStep) { LearnRateV[p] *= 0.95; }
        while (fabs(LearnRateV[p]*CurGradV[p]) < 0.02) { LearnRateV[p] *= (1.0/0.95); } // move more
      } else {
        // set learn rate so that move for each parameter is inside the [MnStep, MxStep]
        while (fabs(LearnRateV[p]*CurGradV[p]) > MxStep) { LearnRateV[p] *= 0.95; printf(".");}
        while (fabs(LearnRateV[p]*CurGradV[p]) < MnStep) { LearnRateV[p] *= (1.0/0.95); printf("*");}
        if (MxStep > 3*MnStep) { MxStep *= 0.95; }
      }
      NewProbMtx.At(p) = ProbMtx.At(p) + LearnRateV[p]*CurGradV[p];
      if (NewProbMtx.At(p) > 0.9999) { NewProbMtx.At(p)=0.9999; }
      if (NewProbMtx.At(p) < 0.0001) { NewProbMtx.At(p)=0.0001; }
    }
    printf("  trueE0: %.2f (%d),  estE0: %.2f (%d),  ERR: %f\n", EZero, Graph->GetEdges(),
      ProbMtx.GetMtxSum(), ProbMtx.GetEdges(KronIters), fabs(EZero-ProbMtx.GetMtxSum()));
    printf("  currLL: %.4f, deltaLL: %.4f\n", CurLL, CurLL-OldLL); // positive is good
    for (int p = 0; p < GetParams(); p++) {
      printf("    %d]  %f  <--  %f + %9f   Grad: %9.1f   Rate: %g\n", p, NewProbMtx.At(p),
        ProbMtx.At(p), (double)(LearnRateV[p]*CurGradV[p]), CurGradV[p](), LearnRateV[p]());
    }
    if (Iter+1 < NIter) { // skip last update
      ProbMtx = NewProbMtx;  ProbMtx.GetLLMtx(LLMtx); }
    OldLL=CurLL;
    printf("\n");  fflush(stdout);

        if(DebugMode) {  
                LLV.Add(CurLL);
                MtxV.Add(NewProbMtx);
        }
  }
  printf("TotalExeTm: %s %g\n", TotalTm.GetStr(), TotalTm.GetSecs());
  ProbMtx.Dump("FITTED PARAMS", false);
  return CurLL;
}
double TKroneckerLL::GradDescent2 ( const int &  NIter,
const double &  LrnRate,
double  MnStep,
double  MxStep,
const int &  WarmUp,
const int &  NSamples 
)

Definition at line 1355 of file kronecker.cpp.

                                                                                                                                               {
  printf("\n----------------------------------------------------------------------\n");
  printf("GradDescent2\n");
  printf("Fitting graph on %d nodes, %d edges\n", Graph->GetNodes(), Graph->GetEdges());
  printf("Skip moves that make likelihood smaller\n");
  printf("Kron iters:  %d (== %d nodes)\n\n", KronIters(), ProbMtx.GetNodes(KronIters()));
  TExeTm IterTm, TotalTm;
  double CurLL=0, NewLL=0;
  const double EZero = pow((double) Graph->GetEdges(), 1.0/double(KronIters));
  TFltV CurGradV, NewGradV, LearnRateV(GetParams()), LastStep(GetParams());
  LearnRateV.PutAll(LrnRate);
  TKronMtx NewProbMtx=ProbMtx, CurProbMtx=ProbMtx;
  bool GoodMove = false;
  // Start
  for (int Iter = 0; Iter < NIter; Iter++) {
    printf("%03d] ", Iter);
    if (! GoodMove) { SampleGradient(WarmUp, NSamples, CurLL, CurGradV); }
    CurProbMtx = ProbMtx;
    // update parameters
    for (int p = 0; p < GetParams(); p++) {
      while (fabs(LearnRateV[p]*CurGradV[p]) > MxStep) { LearnRateV[p] *= 0.95; printf(".");}
      while (fabs(LearnRateV[p]*CurGradV[p]) < MnStep) { LearnRateV[p] *= (1.0/0.95); printf("*");}
      NewProbMtx.At(p) = CurProbMtx.At(p) + LearnRateV[p]*CurGradV[p];
      if (NewProbMtx.At(p) > 0.9999) { NewProbMtx.At(p)=0.9999; }
      if (NewProbMtx.At(p) < 0.0001) { NewProbMtx.At(p)=0.0001; }
      LearnRateV[p] *= 0.95;
    }
    printf("  ");
    ProbMtx=NewProbMtx;  ProbMtx.GetLLMtx(LLMtx);
    SampleGradient(WarmUp, NSamples, NewLL, NewGradV);
    if (NewLL > CurLL) { // accept the move
      printf("== Good move:\n");
      printf("  trueE0: %.2f (%d),  estE0: %.2f (%d),  ERR: %f\n", EZero, Graph->GetEdges(),
        ProbMtx.GetMtxSum(), ProbMtx.GetEdges(KronIters), fabs(EZero-ProbMtx.GetMtxSum()));
      printf("  currLL: %.4f  deltaLL: %.4f\n", CurLL, NewLL-CurLL); // positive is good
      for (int p = 0; p < GetParams(); p++) {
        printf("    %d]  %f  <--  %f + %9f   Grad: %9.1f   Rate: %g\n", p, NewProbMtx.At(p),
          CurProbMtx.At(p), (double)(LearnRateV[p]*CurGradV[p]), CurGradV[p](), LearnRateV[p]()); }
      CurLL = NewLL;
      CurGradV = NewGradV;
      GoodMove = true;
    } else {
      printf("** BAD move:\n");
      printf("  *trueE0: %.2f (%d),  estE0: %.2f (%d),  ERR: %f\n", EZero, Graph->GetEdges(),
        ProbMtx.GetMtxSum(), ProbMtx.GetEdges(KronIters), fabs(EZero-ProbMtx.GetMtxSum()));
      printf("  *curLL:  %.4f  deltaLL: %.4f\n", CurLL, NewLL-CurLL); // positive is good
      for (int p = 0; p < GetParams(); p++) {
        printf("   b%d]  %f  <--  %f + %9f   Grad: %9.1f   Rate: %g\n", p, NewProbMtx.At(p),
          CurProbMtx.At(p), (double)(LearnRateV[p]*CurGradV[p]), CurGradV[p](), LearnRateV[p]()); }
      // move to old position
      ProbMtx = CurProbMtx;  ProbMtx.GetLLMtx(LLMtx);
      GoodMove = false;
    }
    printf("\n");  fflush(stdout);
  }
  printf("TotalExeTm: %s %g\n", TotalTm.GetStr(), TotalTm.GetSecs());
  ProbMtx.Dump("FITTED PARAMS\n", false);
  return CurLL;
}
void TKroneckerLL::GradDescentConvergence ( const TStr &  OutFNm,
const TStr &  Desc1,
const bool &  SamplePerm,
const int &  NIters,
double  LearnRate,
const int &  WarmUp,
const int &  NSamples,
const int &  AvgKGraphs,
const TKronMtx &  TrueParam 
)

Definition at line 2017 of file kronecker.cpp.

                                                                                                             {
  TExeTm IterTm;
  int Iter;
  double OldLL=0, MyLL=0, AvgAbsErr=0, AbsSumErr=0;
  TFltV MyGradV, SDevV;
  TFltV LearnRateV(GetParams());  LearnRateV.PutAll(LearnRate);
  TFltPrV EZeroV, DiamV, Lambda1V, Lambda2V, AvgAbsErrV, AvgLLV;
  TFltPrV TrueEZeroV, TrueDiamV, TrueLambda1V, TrueLambda2V, TrueLLV;
  TFltV SngValV;  TSnap::GetSngVals(Graph, 2, SngValV);  SngValV.Sort(false);
  const double TrueEZero = pow((double) Graph->GetEdges(), 1.0/double(KronIters));
  const double TrueEffDiam = TSnap::GetAnfEffDiam(Graph, false, 10);
  const double TrueLambda1 = SngValV[0];
  const double TrueLambda2 = SngValV[1];
  if (! TrueParam.Empty()) {
    const TKronMtx CurParam = ProbMtx;  ProbMtx.Dump();
    InitLL(TrueParam);  SetOrderPerm();  CalcApxGraphLL(); printf("TrueLL: %f\n", LogLike());
    OldLL = LogLike;  InitLL(CurParam);
  }
  const double TrueLL = OldLL;
  if (! SamplePerm) { SetOrderPerm(); } else { SetDegPerm(); }
  for (Iter = 0; Iter < NIters; Iter++) {
    if (! SamplePerm) {
      // don't sample over permutations
      CalcApxGraphDLL();  CalcApxGraphLL();   // fast
      MyLL = LogLike;  MyGradV = GradV;
    } else {
      // sample over permutations (approximate calculations)
      SampleGradient(WarmUp, NSamples, MyLL, MyGradV);
    }
    double SumDiam=0, SumSngVal1=0, SumSngVal2=0;
    for (int trial = 0; trial < AvgKGraphs; trial++) {
      // generate kronecker graph
      PNGraph KronGraph = TKronMtx::GenFastKronecker(ProbMtx, KronIters, true, 0); // approx
      //PNGraph KronGraph = TKronMtx::GenKronecker(ProbMtx, KronIters, true, 0); // true
      SngValV.Clr(true);  TSnap::GetSngVals(KronGraph, 2, SngValV);  SngValV.Sort(false);
      SumDiam += TSnap::GetAnfEffDiam(KronGraph, false, 10);
      SumSngVal1 += SngValV[0];  SumSngVal2 += SngValV[1];
    }
    // how good is the current fit
    AvgLLV.Add(TFltPr(Iter, MyLL));
    EZeroV.Add(TFltPr(Iter, ProbMtx.GetMtxSum()));
    DiamV.Add(TFltPr(Iter, SumDiam/double(AvgKGraphs)));
    Lambda1V.Add(TFltPr(Iter, SumSngVal1/double(AvgKGraphs)));
    Lambda2V.Add(TFltPr(Iter, SumSngVal2/double(AvgKGraphs)));
    TrueLLV.Add(TFltPr(Iter, TrueLL));
    TrueEZeroV.Add(TFltPr(Iter, TrueEZero));
    TrueDiamV.Add(TFltPr(Iter, TrueEffDiam));
    TrueLambda1V.Add(TFltPr(Iter, TrueLambda1));
    TrueLambda2V.Add(TFltPr(Iter, TrueLambda2));
    if (Iter % 10 == 0) {
      const TStr Desc = TStr::Fmt("%s. Iter: %d, G(%d, %d)  K(%d, %d)", Desc1.Empty()?OutFNm.CStr():Desc1.CStr(),
        Iter, Graph->GetNodes(), Graph->GetEdges(), ProbMtx.GetNodes(KronIters), ProbMtx.GetEdges(KronIters));
      PlotTrueAndEst("LL."+OutFNm, Desc, "Average LL", AvgLLV, TrueLLV);
      PlotTrueAndEst("E0."+OutFNm, Desc, "E0 (expected number of edges)", EZeroV, TrueEZeroV);
      PlotTrueAndEst("Diam."+OutFNm+"-Diam", Desc, "Effective diameter", DiamV, TrueDiamV);
      PlotTrueAndEst("Lambda1."+OutFNm, Desc, "Lambda 1", Lambda1V, TrueLambda1V);
      PlotTrueAndEst("Lambda2."+OutFNm, Desc, "Lambda 2", Lambda2V, TrueLambda2V);
      if (! TrueParam.Empty()) {
        PlotTrueAndEst("AbsErr."+OutFNm, Desc, "Average Absolute Error", AvgAbsErrV, TFltPrV()); }
    }
    if (! TrueParam.Empty()) {
      AvgAbsErr = TKronMtx::GetAvgAbsErr(ProbMtx, TrueParam);
      AvgAbsErrV.Add(TFltPr(Iter, AvgAbsErr));
    } else { AvgAbsErr = 1.0; }
    // update parameters
    AbsSumErr = fabs(ProbMtx.GetMtxSum() - TrueEZero);
    // update parameters
    for (int p = 0; p < GetParams(); p++) {
      LearnRateV[p] *= 0.99;
      while (fabs(LearnRateV[p]*MyGradV[p]) > 0.1) { LearnRateV[p] *= 0.99; printf(".");}
      while (fabs(LearnRateV[p]*MyGradV[p]) < 0.002) { LearnRateV[p] *= (1.0/0.95); printf("*");}
      printf("    %d]  %f  <--  %f + %9f   Grad: %9.1f,  Rate:%g\n", p, ProbMtx.At(p) + LearnRateV[p]*MyGradV[p],
        ProbMtx.At(p), (double)(LearnRateV[p]*MyGradV[p]), MyGradV[p](), LearnRateV[p]());
      ProbMtx.At(p) = ProbMtx.At(p) + LearnRateV[p]*MyGradV[p];
      // box constraints
      if (ProbMtx.At(p) > 1.0) { ProbMtx.At(p)=1.0; }
      if (ProbMtx.At(p) < 0.001) { ProbMtx.At(p)=0.001; }
    }
    printf("%d] LL: %g, ", Iter, MyLL);
    printf("  avgAbsErr: %.4f, absSumErr: %.4f, newLL: %.2f, deltaLL: %.2f\n", AvgAbsErr, AbsSumErr, MyLL, OldLL-MyLL);
    if (AvgAbsErr < 0.001) { printf("***CONVERGED!\n");  break; }
    printf("\n");  fflush(stdout);
    ProbMtx.GetLLMtx(LLMtx);  OldLL = MyLL;
  }
  TrueParam.Dump("True  Thetas", true);
  ProbMtx.Dump("Final Thetas", true);
  printf("  AvgAbsErr: %f\n  AbsSumErr: %f\n  Iterations: %d\n", AvgAbsErr, AbsSumErr, Iter);
  printf("Iteration run time: %s, sec: %g\n\n", IterTm.GetTmStr(), IterTm.GetSecs());
}
void TKroneckerLL::InitLL ( const TFltV &  ParamV)

Definition at line 996 of file kronecker.cpp.

                                             {
  InitLL(TKronMtx(ParamV));
}
void TKroneckerLL::InitLL ( const TKronMtx &  ParamMtx)

Definition at line 1000 of file kronecker.cpp.

                                                  {
  IAssert(ParamMtx.IsProbMtx());
  ProbMtx = ParamMtx;
  ProbMtx.GetLLMtx(LLMtx);
  LogLike = TKronMtx::NInf;
  if (GradV.Len() != ProbMtx.Len()) {
    GradV.Gen(ProbMtx.Len()); }
  GradV.PutAll(0.0);
}
void TKroneckerLL::InitLL ( const PNGraph &  GraphPt,
const TKronMtx &  ParamMtx 
)

Definition at line 1010 of file kronecker.cpp.

                                                                          {
  IAssert(ParamMtx.IsProbMtx());
  ProbMtx = ParamMtx;
  ProbMtx.GetLLMtx(LLMtx);
  SetGraph(GraphPt);
  LogLike = TKronMtx::NInf;
  if (GradV.Len() != ProbMtx.Len()) {
    GradV.Gen(ProbMtx.Len()); }
  GradV.PutAll(0.0);
}
bool TKroneckerLL::IsLatentEdge ( const int &  NId1,
const int &  NId2 
) const [inline]

Definition at line 175 of file kronecker.h.

{ return !IsObsEdge(NId1, NId2);        }
bool TKroneckerLL::IsLatentNode ( const int &  NId) const [inline]

Definition at line 174 of file kronecker.h.

{ return !IsObsNode(NId);       }
bool TKroneckerLL::IsObsEdge ( const int &  NId1,
const int &  NId2 
) const [inline]

Definition at line 173 of file kronecker.h.

{ IAssert(RealNodes > 0);       return ((NId1 < RealNodes) && (NId2 < RealNodes));      }
bool TKroneckerLL::IsObsNode ( const int &  NId) const [inline]

Definition at line 172 of file kronecker.h.

{ IAssert(RealNodes > 0);       return (NId < RealNodes);       }
void TKroneckerLL::MetroGibbsSampleNext ( const int &  WarmUp,
const bool  DLLUpdate = false 
)

Definition at line 1503 of file kronecker.cpp.

                                                                               {
        int NId1 = 0, NId2 = 0, hit = 0, GId = 0;
        TIntTr EdgeToRemove, NewEdge;
        double RndAccept;

        if(LEdgeV.Len()) {
                for(int i = 0; i < WarmUp; i++) {
                        hit = TKronMtx::Rnd.GetUniDevInt(LEdgeV.Len());

                        NId1 = LEdgeV[hit].Val1;        NId2 = LEdgeV[hit].Val2;
                        GId = LEdgeV[hit].Val3;
                        SetRandomEdges(1, true);
                        NewEdge = LEdgeV.Last();

                        RndAccept = (1.0 - exp(LLMtx.GetEdgeLL(NewEdge.Val1, NewEdge.Val2, KronIters))) / (1.0 - exp(LLMtx.GetEdgeLL(NId1, NId2, KronIters)));
                        RndAccept = (RndAccept > 1.0) ? 1.0 : RndAccept;

                        if(TKronMtx::Rnd.GetUniDev() > RndAccept) { //  reject
                                Graph->DelEdge(NewEdge.Val1, NewEdge.Val2);
                                if(NewEdge.Val1 != NewEdge.Val2) {      GEdgeV.DelLast();       }
                                else {  LSelfEdge--;    }
                                LEdgeV.DelLast();
                        } else {        //      accept
                                Graph->DelEdge(NId1, NId2);
                                LEdgeV[hit] = LEdgeV.Last();
                                LEdgeV.DelLast();
                                if(NId1 == NId2) {
                                        LSelfEdge--;
                                        if(NewEdge.Val1 != NewEdge.Val2) {
                                                GEdgeV[GEdgeV.Len()-1].Val3 = hit;
                                        }
                                } else {
                                        IAssertR(GEdgeV.Last().Val3 >= 0, "Invalid indexing");

                                        GEdgeV[GId] = GEdgeV.Last();
                                        if(NewEdge.Val1 != NewEdge.Val2) {
                                                GEdgeV[GId].Val3 = hit;
                                        }
                                        LEdgeV[GEdgeV[GId].Val3].Val3 = GId;
                                        GEdgeV.DelLast();
                                }

                        LogLike += LLMtx.GetApxNoEdgeLL(EdgeToRemove.Val1, EdgeToRemove.Val2, KronIters) - LLMtx.GetEdgeLL(EdgeToRemove.Val1, EdgeToRemove.Val2, KronIters);
                        LogLike += -LLMtx.GetApxNoEdgeLL(NewEdge.Val1, NewEdge.Val2, KronIters) + LLMtx.GetEdgeLL(NewEdge.Val1, NewEdge.Val2, KronIters);

                                if(DLLUpdate) {
                                        for (int p = 0; p < LLMtx.Len(); p++) {
                                                GradV[p] += LLMtx.GetApxNoEdgeDLL(p, EdgeToRemove.Val1, EdgeToRemove.Val2, KronIters) - LLMtx.GetEdgeDLL(p, EdgeToRemove.Val1, EdgeToRemove.Val2, KronIters);
                                                GradV[p] += -LLMtx.GetApxNoEdgeDLL(p, NewEdge.Val1, NewEdge.Val2, KronIters) + LLMtx.GetEdgeDLL(p, NewEdge.Val1, NewEdge.Val2, KronIters);
                                        }
                                }
                        }
                }
        }

//      CalcApxGraphLL();
        for (int s = 0; s < WarmUp; s++) {
                if(SampleNextPerm(NId1, NId2)) {
                        if(DLLUpdate)   UpdateGraphDLL(NId1, NId2);
                }
        }
}
void TKroneckerLL::MetroGibbsSampleSetup ( const int &  WarmUp)

Definition at line 1476 of file kronecker.cpp.

                                                          {
        double alpha = log(ProbMtx.GetMtxSum()) / log(double(ProbMtx.GetDim()));
        int NId1 = 0, NId2 = 0;
        int NMissing;

        RestoreGraph(false);
        if(EMType == kronEdgeMiss) {
                CalcApxGraphLL();
                for (int s = 0; s < WarmUp; s++)        SampleNextPerm(NId1, NId2);
        }

        if(EMType == kronFutureLink) {
                NMissing = (int) (pow(ProbMtx.GetMtxSum(), KronIters) - pow(double(RealNodes), alpha));
        } else if(EMType == kronEdgeMiss) {
                NMissing = MissEdges;
        } else {
                NMissing = (int) (pow(ProbMtx.GetMtxSum(), KronIters) * (1.0 - pow(double(RealNodes) / double(Nodes), 2)));
        }
        NMissing = (NMissing < 1) ? 1 : NMissing;

        SetRandomEdges(NMissing, true);

        CalcApxGraphLL();
        for (int s = 0; s < WarmUp; s++)        SampleNextPerm(NId1, NId2);
}
static PKroneckerLL TKroneckerLL::New ( ) [inline, static]

Definition at line 154 of file kronecker.h.

{ return new TKroneckerLL(); }
PKroneckerLL TKroneckerLL::New ( const PNGraph &  GraphPt,
const TKronMtx &  ParamMtx,
const double &  PermPSwapNd = 0.1 
) [static]

Definition at line 797 of file kronecker.cpp.

                                                                                                          {
  return new TKroneckerLL(GraphPt, ParamMtx, PermPSwapNd);
}
PKroneckerLL TKroneckerLL::New ( const PNGraph &  GraphPt,
const TKronMtx &  ParamMtx,
const TIntV &  NodeIdPermV,
const double &  PermPSwapNd = 0.2 
) [static]

Definition at line 801 of file kronecker.cpp.

                                                                                                                                    {
  return new TKroneckerLL(GraphPt, ParamMtx, NodeIdPermV, PermPSwapNd);
}
double TKroneckerLL::NodeDLLDelta ( const int  ParamId,
const int &  NId 
) const

Definition at line 1214 of file kronecker.cpp.

                                                                         {
  if (! Graph->IsNode(NId)) { return 0.0; } // zero degree node
  double Delta = 0.0;
  const TNGraph::TNodeI Node = Graph->GetNI(NId);
  const int SrcRow = NodePerm[NId];
  for (int e = 0; e < Node.GetOutDeg(); e++) {
    const int DstCol = NodePerm[Node.GetOutNId(e)];
    Delta += - LLMtx.GetApxNoEdgeDLL(ParamId, SrcRow, DstCol, KronIters)
      + LLMtx.GetEdgeDLL(ParamId, SrcRow, DstCol, KronIters);
  }
  const int SrcCol = NodePerm[NId];
  for (int e = 0; e < Node.GetInDeg(); e++) {
    const int DstRow = NodePerm[Node.GetInNId(e)];
    Delta += - LLMtx.GetApxNoEdgeDLL(ParamId, DstRow, SrcCol, KronIters)
      + LLMtx.GetEdgeDLL(ParamId, DstRow, SrcCol, KronIters);
  }
  // double counter self-edge
  if (Graph->IsEdge(NId, NId)) {
    Delta += + LLMtx.GetApxNoEdgeDLL(ParamId, SrcRow, SrcCol, KronIters)
      - LLMtx.GetEdgeDLL(ParamId, SrcRow, SrcCol, KronIters);
    IAssert(SrcRow == SrcCol);
  }
  return Delta;
}
double TKroneckerLL::NodeLLDelta ( const int &  NId) const

Definition at line 1055 of file kronecker.cpp.

                                                     {
  if (! Graph->IsNode(NId)) { return 0.0; } // zero degree node
  double Delta = 0.0;
  const TNGraph::TNodeI Node = Graph->GetNI(NId);
  // out-edges
  const int SrcRow = NodePerm[NId];
  for (int e = 0; e < Node.GetOutDeg(); e++) {
    const int DstCol = NodePerm[Node.GetOutNId(e)];
    Delta += - LLMtx.GetApxNoEdgeLL(SrcRow, DstCol, KronIters)
      + LLMtx.GetEdgeLL(SrcRow, DstCol, KronIters);
  }
  //in-edges
  const int SrcCol = NodePerm[NId];
  for (int e = 0; e < Node.GetInDeg(); e++) {
    const int DstRow = NodePerm[Node.GetInNId(e)];
    Delta += - LLMtx.GetApxNoEdgeLL(DstRow, SrcCol, KronIters)
      + LLMtx.GetEdgeLL(DstRow, SrcCol, KronIters);
  }
  // double counted self-edge
  if (Graph->IsEdge(NId, NId)) {
    Delta += + LLMtx.GetApxNoEdgeLL(SrcRow, SrcCol, KronIters)
      - LLMtx.GetEdgeLL(SrcRow, SrcCol, KronIters);
    IAssert(SrcRow == SrcCol);
  }
  return Delta;
}
void TKroneckerLL::RestoreGraph ( const bool  RestoreNodes = true)

!!!!! MYUNGHWAN, CHECK!

Definition at line 923 of file kronecker.cpp.

                                                       {
        //      remove from Graph
        int NId1, NId2;
        for (int e = 0; e < LEdgeV.Len(); e++) {
        NId1 = LEdgeV[e].Val1;  NId2 = LEdgeV[e].Val2;
                Graph->DelEdge(NId1, NId2);
//              GEdgeV.DelIfIn(LEdgeV[e]);
        }
        if(LEdgeV.Len() - LSelfEdge)
                GEdgeV.Del(GEdgeV.Len() - LEdgeV.Len() + LSelfEdge, GEdgeV.Len() - 1);
        LEdgeV.Clr();
        LSelfEdge = 0;

        if(RestoreNodes) {
                for(int i = Graph->GetNodes()-1; i >= RealNodes; i--) {
                        Graph->DelNode(i);
                }
        }
}
void TKroneckerLL::RunEStep ( const int &  GibbsWarmUp,
const int &  WarmUp,
const int &  NSamples,
TFltV &  LLV,
TVec< TFltV > &  DLLV 
)

Definition at line 1567 of file kronecker.cpp.

                                                                                                                         {
        TExeTm ExeTm, TotalTm;
        LLV.Gen(NSamples, 0);
        DLLV.Gen(NSamples, 0);

        ExeTm.Tick();
        for(int i = 0; i < 2; i++)      MetroGibbsSampleSetup(WarmUp);
        printf("  Warm-Up [%u] : %s\n", WarmUp, ExeTm.GetTmStr());
        CalcApxGraphLL();
        for(int i = 0; i < GibbsWarmUp; i++)    MetroGibbsSampleNext(10, false);
        printf("  Gibbs Warm-Up [%u] : %s\n", GibbsWarmUp, ExeTm.GetTmStr());

        ExeTm.Tick();
        CalcApxGraphLL();
        CalcApxGraphDLL();
        for(int i = 0; i < NSamples; i++) {
                MetroGibbsSampleNext(50, false);

                LLV.Add(LogLike);
                DLLV.Add(GradV);

                int OnePercent = (i+1) % (NSamples / 10);
                if(OnePercent == 0) {
                        int TenPercent = ((i+1) / (NSamples / 10)) * 10;
                        printf("  %3u%% done : %s\n", TenPercent, ExeTm.GetTmStr());
                }
        }
}
void TKroneckerLL::RunKronEM ( const int &  EMIter,
const int &  GradIter,
double  LrnRate,
double  MnStep,
double  MxStep,
const int &  GibbsWarmUp,
const int &  WarmUp,
const int &  NSamples,
const TKronEMType &  Type = kronNodeMiss,
const int &  NMissing = -1 
)

Definition at line 1692 of file kronecker.cpp.

                                                                                                                                                                                                                               {
        printf("\n----------------------------------------------------------------------\n");
        printf("Fitting graph on %d nodes, %d edges\n", int(RealNodes), int(RealEdges));
        printf("Kron iters:  %d (== %d nodes)\n\n", KronIters(), ProbMtx.GetNodes(KronIters()));

        TFltV LLV(NSamples);
        TVec<TFltV> DLLV(NSamples);
        //int count = 0;

        EMType = Type;
        MissEdges = NMissing;
        AppendIsoNodes();
        SetRndPerm();

        if(DebugMode) {
                LLV.Gen(EMIter, 0);
                MtxV.Gen(EMIter, 0);
        }

        for(int i = 0; i < EMIter; i++) {
                printf("\n----------------------------------------------------------------------\n");
                printf("%03d EM-iter] E-Step\n", i+1);
                RunEStep(GibbsWarmUp, WarmUp, NSamples, LLV, DLLV);
                printf("\n\n");

                printf("%03d EM-iter] M-Step\n", i+1);
                double CurLL = RunMStep(LLV, DLLV, GradIter, LrnRate, MnStep, MxStep);
                printf("\n\n");

                if(DebugMode) {
                        LLV.Add(CurLL);
                        MtxV.Add(ProbMtx);
                }
        }

        RestoreGraph();
}
double TKroneckerLL::RunMStep ( const TFltV &  LLV,
const TVec< TFltV > &  DLLV,
const int &  GradIter,
const double &  LrnRate,
double  MnStep,
double  MxStep 
)

Definition at line 1597 of file kronecker.cpp.

                                                                                                                                                 {
        TExeTm IterTm, TotalTm;
        double OldLL=LogLike, CurLL=0;
        const double alpha = log(double(RealEdges)) / log(double(RealNodes));
        const double EZero = pow(double(ProbMtx.GetDim()), alpha);

        TFltV CurGradV(GetParams()), LearnRateV(GetParams()), LastStep(GetParams());
        LearnRateV.PutAll(LrnRate);

        TKronMtx NewProbMtx = ProbMtx;
        const int NSamples = LLV.Len();
        const int ReCalcLen = NSamples / 10;

        for (int s = 0; s < LLV.Len(); s++) {
                CurLL += LLV[s];
                for(int p = 0; p < GetParams(); p++) { CurGradV[p] += DLLV[s][p]; }
        }
        CurLL /= NSamples;
        for(int p = 0; p < GetParams(); p++) { CurGradV[p] /= NSamples; }

        double MaxLL = CurLL;
        TKronMtx MaxProbMtx = ProbMtx;
        TKronMtx OldProbMtx = ProbMtx;

        for (int Iter = 0; Iter < GradIter; Iter++) {
                printf("    %03d] ", Iter+1);
                IterTm.Tick();

                for (int p = 0; p < GetParams(); p++) {
                        if (Iter < 1) {
                                while (fabs(LearnRateV[p]*CurGradV[p]) > MxStep) { LearnRateV[p] *= 0.95; }
                                while (fabs(LearnRateV[p]*CurGradV[p]) < 5 * MnStep) { LearnRateV[p] *= (1.0/0.95); } // move more
                        } else {
                        // set learn rate so that move for each parameter is inside the [MnStep, MxStep]
                                while (fabs(LearnRateV[p]*CurGradV[p]) > MxStep) { LearnRateV[p] *= 0.95; printf(".");}
                                while (fabs(LearnRateV[p]*CurGradV[p]) < MnStep) { LearnRateV[p] *= (1.0/0.95); printf("*");}
                                if (MxStep > 3*MnStep) { MxStep *= 0.95; }
                        }
                        NewProbMtx.At(p) = ProbMtx.At(p) + LearnRateV[p]*CurGradV[p];
                        if (NewProbMtx.At(p) > 0.9999) { NewProbMtx.At(p)=0.9999; }
                        if (NewProbMtx.At(p) < 0.0001) { NewProbMtx.At(p)=0.0001; }
                        LearnRateV[p] *= 0.95;
                }
                printf("  trueE0: %.2f (%u from %u),  estE0: %.2f (%u from %u),  ERR: %f\n", EZero, RealEdges(), RealNodes(), ProbMtx.GetMtxSum(), Graph->GetEdges(), Graph->GetNodes(), fabs(EZero-ProbMtx.GetMtxSum()));
                printf("      currLL: %.4f, deltaLL: %.4f\n", CurLL, CurLL-OldLL); // positive is good
                for (int p = 0; p < GetParams(); p++) {
                        printf("      %d]  %f  <--  %f + %9f   Grad: %9.1f   Rate: %g\n", p, NewProbMtx.At(p),
                        ProbMtx.At(p), (double)(LearnRateV[p]*CurGradV[p]), CurGradV[p](), LearnRateV[p]());
                }

                OldLL=CurLL;
                if(Iter == GradIter - 1) {
                        break;
                }

                CurLL = 0;
                CurGradV.PutAll(0.0);
                TFltV OneDLL;

                CalcApxGraphLL();
                CalcApxGraphDLL();

                for(int s = 0; s < NSamples; s++) {
                        ProbMtx = OldProbMtx;  ProbMtx.GetLLMtx(LLMtx);
                        MetroGibbsSampleNext(10, true);
                        ProbMtx = NewProbMtx;  ProbMtx.GetLLMtx(LLMtx);
                        if(s % ReCalcLen == ReCalcLen/2) {
                                CurLL += CalcApxGraphLL();
                                OneDLL = CalcApxGraphDLL();
                        } else {
                                CurLL += LogLike;
                                OneDLL = GradV;
                        }
                        for(int p = 0; p < GetParams(); p++) {
                                CurGradV[p] += OneDLL[p];
                        }
                }
                CurLL /= NSamples;

                if(MaxLL < CurLL) {
                        MaxLL = CurLL;  MaxProbMtx = ProbMtx;
                }

                printf("    Time: %s\n", IterTm.GetTmStr());
                printf("\n");  fflush(stdout);
        }
        ProbMtx = MaxProbMtx;   ProbMtx.GetLLMtx(LLMtx);

        printf("    FinalLL : %f,   TotalExeTm: %s\n", MaxLL, TotalTm.GetTmStr());
        ProbMtx.Dump("    FITTED PARAMS", false);

        return MaxLL;
}
void TKroneckerLL::SampleGradient ( const int &  WarmUp,
const int &  NSamples,
double &  AvgLL,
TFltV &  GradV 
)

Definition at line 1271 of file kronecker.cpp.

                                                                                                        {
  printf("SampleGradient: %s (%s warm-up):", TInt::GetMegaStr(NSamples).CStr(), TInt::GetMegaStr(WarmUp).CStr());
  int NId1=0, NId2=0, NAccept=0;
  TExeTm ExeTm1;
  if (WarmUp > 0) {
    CalcApxGraphLL();
    for (int s = 0; s < WarmUp; s++) { SampleNextPerm(NId1, NId2); }
    printf("  warm-up:%s,", ExeTm1.GetTmStr());  ExeTm1.Tick();
  }
  CalcApxGraphLL(); // re-calculate LL (due to numerical errors)
  CalcApxGraphDLL();
  AvgLL = 0;
  AvgGradV.Gen(LLMtx.Len());  AvgGradV.PutAll(0.0);
  printf("  sampl");
  for (int s = 0; s < NSamples; s++) {
    if (SampleNextPerm(NId1, NId2)) { // new permutation
      UpdateGraphDLL(NId1, NId2);  NAccept++; }
    for (int m = 0; m < LLMtx.Len(); m++) { AvgGradV[m] += GradV[m]; }
    AvgLL += GetLL();
  }
  printf("ing");
  AvgLL = AvgLL / double(NSamples);
  for (int m = 0; m < LLMtx.Len(); m++) {
    AvgGradV[m] = AvgGradV[m] / double(NSamples); }
  printf(":%s (%.0f/s), accept %.1f%%\n", ExeTm1.GetTmStr(), double(NSamples)/ExeTm1.GetSecs(),
    double(100*NAccept)/double(NSamples));
}
bool TKroneckerLL::SampleNextPerm ( int &  NId1,
int &  NId2 
)

Definition at line 1111 of file kronecker.cpp.

                                                      {
  // pick 2 uniform nodes and swap
  if (TKronMtx::Rnd.GetUniDev() < PermSwapNodeProb) {
    NId1 = TKronMtx::Rnd.GetUniDevInt(Nodes);
    NId2 = TKronMtx::Rnd.GetUniDevInt(Nodes);
    while (NId2 == NId1) { NId2 = TKronMtx::Rnd.GetUniDevInt(Nodes); }
  } else {
    // pick uniform edge and swap endpoints (slow as it moves around high degree nodes)
    const int e = TKronMtx::Rnd.GetUniDevInt(GEdgeV.Len());
    NId1 = GEdgeV[e].Val1;  NId2 = GEdgeV[e].Val2;
  }
  const double U = TKronMtx::Rnd.GetUniDev();
  const double OldLL = LogLike;
  const double NewLL = SwapNodesLL(NId1, NId2);
  const double LogU = log(U);
  if (LogU > NewLL - OldLL) { // reject
    LogLike = OldLL;
    NodePerm.Swap(NId2, NId1); //swap back
        InvertPerm.Swap(NodePerm[NId2], NodePerm[NId1]); // swap back
    return false;
  }
  return true; // accept new sample
}
void TKroneckerLL::SetBestDegPerm ( )

!!!!! MYUNGHWAN, CHECK!

Definition at line 842 of file kronecker.cpp.

                                  {
  NodePerm.Gen(Nodes);
  const int NZero = ProbMtx.GetDim();
  TFltIntPrV DegV(Nodes), CDegV(Nodes);
  TFltV Row(NZero);
  TFltV Col(NZero);
  for(int i = 0; i < NZero; i++) {
          for(int j = 0; j < NZero; j++) {
                  Row[i] += ProbMtx.At(i, j);
                  Col[i] += ProbMtx.At(j, i);
          }
  }

  for(int i = 0; i < Nodes; i++) {
          TNGraph::TNodeI NodeI = Graph->GetNI(i);
          int NId = i;
          double RowP = 1.0, ColP = 1.0;
          for(int j = 0; j < KronIters; j++) {
                  int Bit = NId % NZero;
                  RowP *= Row[Bit];             ColP *= Col[Bit];
                  NId /= NZero;
          }
          CDegV[i] = TFltIntPr(RowP + ColP, i);
          DegV[i] = TFltIntPr(NodeI.GetDeg(), i);
  }
  DegV.Sort(false);             CDegV.Sort(false);
  for(int i = 0; i < Nodes; i++) {
          NodePerm[DegV[i].Val2] = CDegV[i].Val2;
  }
  SetIPerm(NodePerm);
}
void TKroneckerLL::SetDebug ( const bool  Debug) [inline]

Definition at line 167 of file kronecker.h.

{ DebugMode = Debug; }
void TKroneckerLL::SetDegPerm ( )

Definition at line 828 of file kronecker.cpp.

                              {
  TIntPrV DegNIdV;
  for (TNGraph::TNodeI NI = Graph->BegNI(); NI < Graph->EndNI(); NI++) {
    DegNIdV.Add(TIntPr(NI.GetDeg(), NI.GetId()));
  }
  DegNIdV.Sort(false);
  NodePerm.Gen(DegNIdV.Len(), 0);
  for (int i = 0; i < DegNIdV.Len(); i++) {
    NodePerm.Add(DegNIdV[i].Val2);
  }
  SetIPerm(NodePerm);
}
void TKroneckerLL::SetGraph ( const PNGraph &  GraphPt)

Definition at line 882 of file kronecker.cpp.

                                                  {
  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)); }
  }
  Nodes = Graph->GetNodes();
  IAssert(LLMtx.GetDim() > 1 && LLMtx.Len() == ProbMtx.Len());
  KronIters = (int) ceil(log(double(Nodes)) / log(double(ProbMtx.GetDim())));
  // edge vector (for swap-edge permutation proposal)
//  if (PermSwapNodeProb < 1.0) { /// !!!!! MYUNGHWAN, CHECK! WHY IS THIS COMMENTED OUT
    GEdgeV.Gen(Graph->GetEdges(), 0);
    for (TNGraph::TEdgeI EI = Graph->BegEI(); EI < Graph->EndEI(); EI++) {
      if (EI.GetSrcNId() != EI.GetDstNId()) {
        GEdgeV.Add(TIntTr(EI.GetSrcNId(), EI.GetDstNId(), -1));
      }
    }
//  }

  RealNodes = Nodes;
  RealEdges = Graph->GetEdges();
  LEdgeV = TIntTrV();
  LSelfEdge = 0;
}
void TKroneckerLL::SetIPerm ( const TIntV &  Perm)

!!!!! MYUNGHWAN, CHECK!

Definition at line 875 of file kronecker.cpp.

                                             {
        InvertPerm.Gen(Perm.Len());
        for (int i = 0; i < Perm.Len(); i++) {
                InvertPerm[Perm[i]] = i;
        }
}
void TKroneckerLL::SetOrderPerm ( )

Definition at line 813 of file kronecker.cpp.

                                {
  NodePerm.Gen(Nodes, 0);
  for (int i = 0; i < Graph->GetNodes(); i++) {
    NodePerm.Add(i); }
  SetIPerm(NodePerm);
}
void TKroneckerLL::SetPerm ( const char &  PermId)

Definition at line 805 of file kronecker.cpp.

                                             {
  if (PermId == 'o') { SetOrderPerm(); }
  else if (PermId == 'd') { SetDegPerm(); }
  else if (PermId == 'r') { SetRndPerm(); }
  else if (PermId == 'b') { SetBestDegPerm(); }
  else FailR("Unknown permutation type (o,d,r)");
}
void TKroneckerLL::SetPerm ( const TIntV &  NodePermV) [inline]

Definition at line 183 of file kronecker.h.

{ NodePerm = NodePermV; SetIPerm(NodePerm); }
void TKroneckerLL::SetRandomEdges ( const int &  NEdges,
const bool  isDir = true 
)

!!!!! MYUNGHWAN, CHECK!

Definition at line 1417 of file kronecker.cpp.

                                                                     {
        int count = 0, added = 0, collision = 0;
        const int MtxDim = ProbMtx.GetDim();
        const double MtxSum = ProbMtx.GetMtxSum();
        TVec<TFltIntIntTr> ProbToRCPosV; // row, col position
        double CumProb = 0.0;

        for(int r = 0; r < MtxDim; r++) {
                for(int c = 0; c < MtxDim; c++) {
                        const double Prob = ProbMtx.At(r, c);
                        if (Prob > 0.0) {
                                CumProb += Prob;
                                ProbToRCPosV.Add(TFltIntIntTr(CumProb/MtxSum, r, c));
                        }
                }
        }

        int Rng, Row, Col, n, NId1, NId2;
        while(added < NEdges) {
                Rng = Nodes;    Row = 0;        Col = 0;
                for (int iter = 0; iter < KronIters; iter++) {
                        const double& Prob = TKronMtx::Rnd.GetUniDev();
                        n = 0; while(Prob > ProbToRCPosV[n].Val1) { n++; }
                        const int MtxRow = ProbToRCPosV[n].Val2;
                        const int MtxCol = ProbToRCPosV[n].Val3;
                        Rng /= MtxDim;
                        Row += MtxRow * Rng;
                        Col += MtxCol * Rng;
                }

                count++;

                NId1 = InvertPerm[Row]; NId2 = InvertPerm[Col];

                //      Check conflicts
                if(EMType != kronEdgeMiss && IsObsEdge(NId1, NId2)) {
                        continue;
                }

                if (! Graph->IsEdge(NId1, NId2)) {
                        Graph->AddEdge(NId1, NId2);
                        if(NId1 != NId2)        { GEdgeV.Add(TIntTr(NId1, NId2, LEdgeV.Len())); }
                        else { LSelfEdge++; }
                        LEdgeV.Add(TIntTr(NId1, NId2, GEdgeV.Len()-1));
                        added++;
                        if (! isDir) {
                                if (NId1 != NId2) {
                                   Graph->AddEdge(NId2, NId1);
                                   GEdgeV.Add(TIntTr(NId2, NId1, LEdgeV.Len()));
                                   LEdgeV.Add(TIntTr(NId2, NId1, GEdgeV.Len()-1));
                                   added++;
                                }
                        }
                } else { collision ++; }
        }
//      printf("total = %d / added = %d / collision = %d\n", count, added, collision);
}
void TKroneckerLL::SetRndPerm ( )

Definition at line 820 of file kronecker.cpp.

                              {
  NodePerm.Gen(Nodes, 0);
  for (int i = 0; i < Graph->GetNodes(); i++) {
    NodePerm.Add(i); }
  NodePerm.Shuffle(TKronMtx::Rnd);
  SetIPerm(NodePerm);
}
double TKroneckerLL::SwapNodesLL ( const int &  NId1,
const int &  NId2 
)

Definition at line 1083 of file kronecker.cpp.

                                                                 {
  // subtract old LL (remove nodes)
  LogLike = LogLike - NodeLLDelta(NId1) - NodeLLDelta(NId2);
  const int PrevId1 = NodePerm[NId1], PrevId2 = NodePerm[NId2];
  // double-counted edges
  if (Graph->IsEdge(NId1, NId2)) {
    LogLike += - LLMtx.GetApxNoEdgeLL(PrevId1, PrevId2, KronIters)
      + LLMtx.GetEdgeLL(PrevId1, PrevId2, KronIters); }
  if (Graph->IsEdge(NId2, NId1)) {
    LogLike += - LLMtx.GetApxNoEdgeLL(PrevId2, PrevId1, KronIters)
      + LLMtx.GetEdgeLL(PrevId2, PrevId1, KronIters); }
  // swap
  NodePerm.Swap(NId1, NId2);
  InvertPerm.Swap(NodePerm[NId1], NodePerm[NId2]);
  // add new LL (add nodes)
  LogLike = LogLike + NodeLLDelta(NId1) + NodeLLDelta(NId2);
  const int NewId1 = NodePerm[NId1], NewId2 = NodePerm[NId2];
  // correct for double-counted edges
  if (Graph->IsEdge(NId1, NId2)) {
    LogLike += + LLMtx.GetApxNoEdgeLL(NewId1, NewId2, KronIters)
      - LLMtx.GetEdgeLL(NewId1, NewId2, KronIters); }
  if (Graph->IsEdge(NId2, NId1)) {
    LogLike += + LLMtx.GetApxNoEdgeLL(NewId2, NewId1, KronIters)
      - LLMtx.GetEdgeLL(NewId2, NewId1, KronIters); }
  return LogLike;
}
void TKroneckerLL::TestBicCriterion ( const TStr &  OutFNm,
const TStr &  Desc1,
const PNGraph &  G,
const int &  GradIters,
double  LearnRate,
const int &  WarmUp,
const int &  NSamples,
const int &  TrueN0 
) [static]

Definition at line 2109 of file kronecker.cpp.

                                                                              {
  TFltPrV BicV, MdlV, LLV;
  const double rndGP = G->GetEdges()/TMath::Sqr(double(G->GetNodes()));
  const double RndGLL = G->GetEdges()*log(rndGP )+ (TMath::Sqr(double(G->GetNodes()))-G->GetEdges())*log(1-rndGP);
  LLV.Add(TFltPr(1, RndGLL));
  BicV.Add(TFltPr(1, -RndGLL + 0.5*TMath::Sqr(1)*log(TMath::Sqr(G->GetNodes()))));
  MdlV.Add(TFltPr(1, -RndGLL + 32*TMath::Sqr(1)+2*(log((double)1)+log((double)G->GetNodes()))));
  for (int NZero = 2; NZero < 10; NZero++) {
    const TKronMtx InitKronMtx = TKronMtx::GetInitMtx(NZero, G->GetNodes(), G->GetEdges());
    InitKronMtx.Dump("INIT PARAM", true);
    TKroneckerLL KronLL(G, InitKronMtx);
    KronLL.SetPerm('d'); // degree perm
    const double LastLL = KronLL.GradDescent(GradIters, LearnRate, 0.001, 0.01, WarmUp, NSamples);
    LLV.Add(TFltPr(NZero, LastLL));
    BicV.Add(TFltPr(NZero, -LastLL + 0.5*TMath::Sqr(NZero)*log(TMath::Sqr(G->GetNodes()))));
    MdlV.Add(TFltPr(NZero, -LastLL + 32*TMath::Sqr(NZero)+2*(log((double)NZero)+log((double)KronLL.GetKronIters()))));
    { TGnuPlot GP("LL-"+OutFNm, Desc1);
    GP.AddPlot(LLV, gpwLinesPoints, "Log-likelihood", "linewidth 1 pointtype 6 pointsize 2");
    GP.SetXYLabel("NZero", "Log-Likelihood");  GP.SavePng(); }
    { TGnuPlot GP("BIC-"+OutFNm, Desc1);
    GP.AddPlot(BicV, gpwLinesPoints, "BIC", "linewidth 1 pointtype 6 pointsize 2");
    GP.SetXYLabel("NZero", "BIC");  GP.SavePng(); }
    { TGnuPlot GP("MDL-"+OutFNm, Desc1);
    GP.AddPlot(MdlV, gpwLinesPoints, "MDL", "linewidth 1 pointtype 6 pointsize 2");
    GP.SetXYLabel("NZero", "MDL");  GP.SavePng(); }
  }
}
void TKroneckerLL::TestGradDescent ( const int &  KronIters,
const int &  KiloSamples,
const TStr &  Permutation 
) [static]

Definition at line 2138 of file kronecker.cpp.

                                                                                                        {
  const TStr OutFNm = TStr::Fmt("grad-%s%d-%dk", Permutation.CStr(), KronIters, KiloSamples);
  TKronMtx KronParam = TKronMtx::GetMtx("0.8 0.6; 0.6 0.4");
  PNGraph Graph  = TKronMtx::GenFastKronecker(KronParam, KronIters, true, 0);
  TKroneckerLL KronLL(Graph, KronParam);
  TVec<TFltV> GradVV(4), SDevVV(4);  TFltV XValV;
  int NId1 = 0, NId2 = 0, NAccept = 0;
  TVec<TMom> GradMomV(4);
  TExeTm ExeTm;
  if (Permutation == "r") KronLL.SetRndPerm();
  else if (Permutation == "d") KronLL.SetDegPerm();
  else if (Permutation == "o") KronLL.SetOrderPerm();
  else FailR("Unknown permutation (r,d,o)");
  KronLL.CalcApxGraphLL();
  KronLL.CalcApxGraphDLL();
  for (int s = 0; s < 1000*KiloSamples; s++) {
    if (KronLL.SampleNextPerm(NId1, NId2)) { // new permutation
      KronLL.UpdateGraphDLL(NId1, NId2);  NAccept++; }
    if (s > 50000) { //warm up period
      for (int m = 0; m < 4; m++) { GradVV[m].Add(KronLL.GradV[m]); }
      if ((s+1) % 1000 == 0) {
        printf(".");
        for (int m = 0; m < 4; m++) { GradVV[m].Add(KronLL.GradV[m]); }
        XValV.Add((s+1));
        if ((s+1) % 100000 == 0) {
          TGnuPlot GP(OutFNm, TStr::Fmt("Gradient vs. samples. %d nodes, %d edges", Graph->GetNodes(), Graph->GetEdges()), true);
          for (int g = 0; g < GradVV.Len(); g++) {
            GP.AddPlot(XValV, GradVV[g], gpwLines, TStr::Fmt("grad %d", g)); }
          GP.SetXYLabel("sample index","log Gradient");
          GP.SavePng();
        }
      }
    }
  }
  printf("\n");
  for (int m = 0; m < 4; m++) {
    GradMomV[m].Def();
    printf("grad %d: mean: %12f  sDev: %12f  median: %12f\n", m,
      GradMomV[m].GetMean(), GradMomV[m].GetSDev(), GradMomV[m].GetMedian());
  }
}
TFltQu TKroneckerLL::TestKronDescent ( const bool &  DoExact,
const bool &  TruePerm,
double  LearnRate,
const int &  WarmUp,
const int &  NSamples,
const TKronMtx &  TrueParam 
)

Definition at line 1942 of file kronecker.cpp.

                                                                                                                                                                   {
  printf("Test gradient descent on a synthetic kronecker graphs:\n");
  if (DoExact) { printf("  -- Exact gradient calculations\n"); }
  else { printf("  -- Approximate gradient calculations\n"); }
  if (TruePerm) { printf("  -- No permutation sampling (use true permutation)\n"); }
  else { printf("  -- Sample permutations (start with degree permutation)\n"); }
  TExeTm IterTm;
  int Iter;
  double OldLL=0, MyLL=0, AvgAbsErr, AbsSumErr;
  TFltV MyGradV, SDevV;
  TFltV LearnRateV(GetParams());  LearnRateV.PutAll(LearnRate);
  if (TruePerm) {
    SetOrderPerm();
  }
  else {
    /*printf("SET EIGEN VECTOR PERMUTATIONS\n");
    TFltV LeftSV, RightSV;
    TGSvd::GetSngVec(Graph, LeftSV, RightSV);
    TFltIntPrV V;
    for (int v=0; v<LeftSV.Len();v++) { V.Add(TFltIntPr(LeftSV[v], v)); }
    V.Sort(false);
    NodePerm.Gen(Nodes, 0);
    for (int v=0; v < V.Len();v++) { NodePerm.Add(V[v].Val2); } //*/
    //printf("RANDOM PERMUTATION\n");   SetRndPerm();
    printf("DEGREE  PERMUTATION\n");  SetDegPerm();
  }
  for (Iter = 0; Iter < 100; Iter++) {
    if (TruePerm) {
      // don't sample over permutations
      if (DoExact) { CalcGraphDLL();  CalcGraphLL(); } // slow, O(N^2)
      else { CalcApxGraphDLL();  CalcApxGraphLL(); }   // fast
      MyLL = LogLike;  MyGradV = GradV;
    } else {
      printf(".");
      // sample over permutations (approximate calculations)
      SampleGradient(WarmUp, NSamples, MyLL, MyGradV);
    }
    printf("%d] LL: %g, ", Iter, MyLL);
    AvgAbsErr = TKronMtx::GetAvgAbsErr(ProbMtx, TrueParam);
    AbsSumErr = fabs(ProbMtx.GetMtxSum() - TrueParam.GetMtxSum());
    printf("  avgAbsErr: %.4f, absSumErr: %.4f, newLL: %.2f, deltaLL: %.2f\n", AvgAbsErr, AbsSumErr, MyLL, OldLL-MyLL);
    for (int p = 0; p < GetParams(); p++) {
      // set learn rate so that move for each parameter is inside the [0.01, 0.1]
      LearnRateV[p] *= 0.9;
      //printf("%d: rate: %f   delta:%f\n", p, LearnRateV[p], fabs(LearnRateV[p]*MyGradV[p]));
      while (fabs(LearnRateV[p]*MyGradV[p]) > 0.1) { LearnRateV[p] *= 0.9; }
      //printf("   rate: %f   delta:%f\n", LearnRateV[p], fabs(LearnRateV[p]*MyGradV[p]));
      while (fabs(LearnRateV[p]*MyGradV[p]) < 0.001) { LearnRateV[p] *= (1.0/0.9); }
      //printf("   rate: %f   delta:%f\n", LearnRateV[p], fabs(LearnRateV[p]*MyGradV[p]));
      printf("    %d]  %f  <--  %f + %f    lrnRate:%g\n", p, ProbMtx.At(p) + LearnRateV[p]*MyGradV[p],
        ProbMtx.At(p), (double)(LearnRateV[p]*MyGradV[p]), LearnRateV[p]());
      ProbMtx.At(p) = ProbMtx.At(p) + LearnRateV[p]*MyGradV[p];
      // box constraints
      if (ProbMtx.At(p) > 0.99) { ProbMtx.At(p)=0.99; }
      if (ProbMtx.At(p) < 0.01) { ProbMtx.At(p)=0.01; }
    }
    ProbMtx.GetLLMtx(LLMtx);  OldLL = MyLL;
    if (AvgAbsErr < 0.01) { printf("***CONVERGED!\n");  break; }
    printf("\n");  fflush(stdout);
  }
  TrueParam.Dump("True  Thetas", true);
  ProbMtx.Dump("Final Thetas", true);
  printf("  AvgAbsErr: %f\n  AbsSumErr: %f\n  Iterations: %d\n", AvgAbsErr, AbsSumErr, Iter);
  printf("Iteration run time: %s, sec: %g\n\n", IterTm.GetTmStr(), IterTm.GetSecs());
  return TFltQu(AvgAbsErr, AbsSumErr, Iter, IterTm.GetSecs());
}
TFltV TKroneckerLL::TestSamplePerm ( const TStr &  OutFNm,
const int &  WarmUp,
const int &  NSamples,
const TKronMtx &  TrueMtx,
const bool &  DoPlot = true 
)

Definition at line 1795 of file kronecker.cpp.

                                                                                                                                          {
  printf("Sample permutations: %s (warm-up: %s)\n", TInt::GetMegaStr(NSamples).CStr(), TInt::GetMegaStr(WarmUp).CStr());
  int NId1=0, NId2=0, NAccept=0;
  TExeTm ExeTm;
  const int PlotLen = NSamples/1000+1;
  double TrueLL=-1, AvgLL=0.0;
  TFltV AvgGradV(GetParams());
  TFltPrV TrueLLV(PlotLen, 0); // true log-likelihood (under the correct permutation)
  TFltPrV EstLLV(PlotLen, 0);  // estiamted log-likelihood (averaged over last 1k permutation)
  TFltPrV AcceptV;             // sample acceptance ratio
  TFltV SampleLLV(NSamples, 0);
  TVec<TFltPrV> GradVV(GetParams());
  for (int g = 0; g < GetParams(); g++) {
    GradVV[g].Gen(PlotLen, 0); }
  if (! TrueMtx.Empty()) {
    TIntV PermV=NodePerm;  TKronMtx CurMtx=ProbMtx;  ProbMtx.Dump();
    InitLL(TrueMtx);  SetOrderPerm();  CalcApxGraphLL();  printf("TrueLL: %f\n", LogLike());
    TrueLL=LogLike;  InitLL(CurMtx); NodePerm=PermV;
  }
  CalcApxGraphLL();
  printf("LogLike at start:       %f\n", LogLike());
  if (WarmUp > 0) {
    EstLLV.Add(TFltPr(0, LogLike));
    if (TrueLL != -1) { TrueLLV.Add(TFltPr(0, TrueLL)); }
    for (int s = 0; s < WarmUp; s++) { SampleNextPerm(NId1, NId2); }
    printf("  warm-up:%s,", ExeTm.GetTmStr());  ExeTm.Tick();
  }
  printf("LogLike afterm warm-up: %f\n", LogLike());
  CalcApxGraphLL(); // re-calculate LL (due to numerical errors)
  CalcApxGraphDLL();
  EstLLV.Add(TFltPr(WarmUp, LogLike));
  if (TrueLL != -1) { TrueLLV.Add(TFltPr(WarmUp, TrueLL)); }
  printf("  recalculated:         %f\n", LogLike());
  // start sampling
  printf("  sampling (average per 1000 samples)\n");
  TVec<TFltV> SamplVV(5);
  for (int s = 0; s < NSamples; s++) {
    if (SampleNextPerm(NId1, NId2)) { // new permutation
      UpdateGraphDLL(NId1, NId2);  NAccept++; }
    for (int m = 0; m < AvgGradV.Len(); m++) { AvgGradV[m] += GradV[m]; }
    AvgLL += GetLL();
    SampleLLV.Add(GetLL());
    /*SamplVV[0].Add(GetLL()); // gives worse autocoreelation than the avg below
    SamplVV[1].Add(GradV[0]);
    SamplVV[2].Add(GradV[1]);
    SamplVV[3].Add(GradV[2]);
    SamplVV[4].Add(GradV[3]);*/
    if (s > 0 && s % 1000 == 0) {
      printf(".");
      for (int g = 0; g < AvgGradV.Len(); g++) {
        GradVV[g].Add(TFltPr(WarmUp+s, AvgGradV[g] / 1000.0)); }
      EstLLV.Add(TFltPr(WarmUp+s, AvgLL / 1000.0));
      if (TrueLL != -1) { TrueLLV.Add(TFltPr(WarmUp+s, TrueLL)); }
      AcceptV.Add(TFltPr(WarmUp+s, NAccept/1000.0));
      // better (faster decaying) autocorrelation when one takes avg. of 1000 consecutive samples
      /*SamplVV[0].Add(AvgLL);
      SamplVV[1].Add(AvgGradV[0]);
      SamplVV[2].Add(AvgGradV[1]);
      SamplVV[3].Add(AvgGradV[2]);
      SamplVV[4].Add(AvgGradV[3]); //*/
      if (s % 100000 == 0 && DoPlot) {
        const TStr Desc = TStr::Fmt("P(NodeSwap)=%g. Nodes: %d, Edges: %d, Params: %d, WarmUp: %s, Samples: %s", PermSwapNodeProb(),
          Graph->GetNodes(), Graph->GetEdges(), GetParams(), TInt::GetMegaStr(WarmUp).CStr(), TInt::GetMegaStr(NSamples).CStr());
        PlotGrad(EstLLV, TrueLLV, GradVV, AcceptV, OutFNm, Desc);
        for (int n = 0; n < SamplVV.Len(); n++) {
          PlotAutoCorrelation(SamplVV[n], 1000, TStr::Fmt("%s-n%d", OutFNm.CStr(), n), Desc); }
        printf("  samples: %d, time: %s, samples/s: %.1f\n", s, ExeTm.GetTmStr(), double(s+1)/ExeTm.GetSecs());
      }
      AvgLL = 0;  AvgGradV.PutAll(0);  NAccept=0;
    }
  }
  if (DoPlot) {
    const TStr Desc = TStr::Fmt("P(NodeSwap)=%g. Nodes: %d, Edges: %d, Params: %d, WarmUp: %s, Samples: %s", PermSwapNodeProb(),
      Graph->GetNodes(), Graph->GetEdges(), GetParams(), TInt::GetMegaStr(WarmUp).CStr(), TInt::GetMegaStr(NSamples).CStr());
    PlotGrad(EstLLV, TrueLLV, GradVV, AcceptV, OutFNm, Desc);
    for (int n = 0; n < SamplVV.Len(); n++) {
      PlotAutoCorrelation(SamplVV[n], 1000, TStr::Fmt("%s-n%d", OutFNm.CStr(), n), Desc); }
  }
  return SampleLLV; // seems to work better for potential scale reduction plot
}
void TKroneckerLL::UpdateGraphDLL ( const int &  SwapNId1,
const int &  SwapNId2 
)

Definition at line 1241 of file kronecker.cpp.

                                                                          {
  for (int ParamId = 0; ParamId < LLMtx.Len(); ParamId++) {
    // permutation before the swap (swap back to previous position)
    NodePerm.Swap(SwapNId1, SwapNId2);
    // subtract old DLL
    TFlt& DLL = GradV[ParamId];
    DLL = DLL - NodeDLLDelta(ParamId, SwapNId1) - NodeDLLDelta(ParamId, SwapNId2);
    // double-counted edges
    const int PrevId1 = NodePerm[SwapNId1], PrevId2 = NodePerm[SwapNId2];
    if (Graph->IsEdge(SwapNId1, SwapNId2)) {
      DLL += - LLMtx.GetApxNoEdgeDLL(ParamId, PrevId1, PrevId2, KronIters)
        + LLMtx.GetEdgeDLL(ParamId, PrevId1, PrevId2, KronIters); }
    if (Graph->IsEdge(SwapNId2, SwapNId1)) {
      DLL += - LLMtx.GetApxNoEdgeDLL(ParamId, PrevId2, PrevId1, KronIters)
        + LLMtx.GetEdgeDLL(ParamId, PrevId2, PrevId1, KronIters); }
    // permutation after the swap (restore the swap)
    NodePerm.Swap(SwapNId1, SwapNId2);
    // add new DLL
    DLL = DLL + NodeDLLDelta(ParamId, SwapNId1) + NodeDLLDelta(ParamId, SwapNId2);
    const int NewId1 = NodePerm[SwapNId1], NewId2 = NodePerm[SwapNId2];
    // double-counted edges
    if (Graph->IsEdge(SwapNId1, SwapNId2)) {
      DLL += + LLMtx.GetApxNoEdgeDLL(ParamId, NewId1, NewId2, KronIters)
        - LLMtx.GetEdgeDLL(ParamId, NewId1, NewId2, KronIters); }
    if (Graph->IsEdge(SwapNId2, SwapNId1)) {
      DLL += + LLMtx.GetApxNoEdgeDLL(ParamId, NewId2, NewId1, KronIters)
        - LLMtx.GetEdgeDLL(ParamId, NewId2, NewId1, KronIters); }
  }
}

Friends And Related Function Documentation

friend class TPt< TKroneckerLL > [friend]

Definition at line 245 of file kronecker.h.

Member Data Documentation

TCRef TKroneckerLL::CRef [private]

Definition at line 119 of file kronecker.h.

TBool TKroneckerLL::DebugMode [private]

Definition at line 141 of file kronecker.h.

TKronEMType TKroneckerLL::EMType [private]

Definition at line 138 of file kronecker.h.

TIntTrV TKroneckerLL::GEdgeV [private]

Definition at line 125 of file kronecker.h.

TFltV TKroneckerLL::GradV [private]

Definition at line 136 of file kronecker.h.

PNGraph TKroneckerLL::Graph [private]

Definition at line 120 of file kronecker.h.

TIntV TKroneckerLL::InvertPerm [private]

Definition at line 129 of file kronecker.h.

TInt TKroneckerLL::KronIters [private]

Definition at line 121 of file kronecker.h.

TIntTrV TKroneckerLL::LEdgeV [private]

Definition at line 126 of file kronecker.h.

TKronMtx TKroneckerLL::LLMtx [private]

Definition at line 134 of file kronecker.h.

TFltV TKroneckerLL::LLV [private]

Definition at line 142 of file kronecker.h.

TFlt TKroneckerLL::LogLike [private]

Definition at line 135 of file kronecker.h.

TInt TKroneckerLL::LSelfEdge [private]

Definition at line 127 of file kronecker.h.

TInt TKroneckerLL::MissEdges [private]

Definition at line 139 of file kronecker.h.

TVec<TKronMtx> TKroneckerLL::MtxV [private]

Definition at line 143 of file kronecker.h.

TIntV TKroneckerLL::NodePerm [private]

Definition at line 128 of file kronecker.h.

TInt TKroneckerLL::Nodes [private]

Definition at line 121 of file kronecker.h.

TFlt TKroneckerLL::PermSwapNodeProb [private]

Definition at line 123 of file kronecker.h.

TKronMtx TKroneckerLL::ProbMtx [private]

Definition at line 134 of file kronecker.h.

TInt TKroneckerLL::RealEdges [private]

Definition at line 132 of file kronecker.h.

TInt TKroneckerLL::RealNodes [private]

Definition at line 131 of file kronecker.h.

The documentation for this class was generated from the following files: