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
ff.cpp
Go to the documentation of this file.
00001 
00002 // Forest Fire
00003 void TForestFire::InfectAll() {
00004   InfectNIdV.Gen(Graph->GetNodes());
00005   for (TNGraph::TNodeI NI = Graph->BegNI(); NI < Graph->EndNI(); NI++) {
00006     InfectNIdV.Add(NI.GetId()); }
00007 }
00008 
00009 void TForestFire::InfectRnd(const int& NInfect) {
00010   IAssert(NInfect < Graph->GetNodes());
00011   TIntV NIdV(Graph->GetNodes(), 0);
00012   for (TNGraph::TNodeI NI = Graph->BegNI(); NI < Graph->EndNI(); NI++) {
00013     NIdV.Add(NI.GetId()); }
00014   NIdV.Shuffle(Rnd);
00015   InfectNIdV.Gen(NInfect, 0);
00016   for (int i = 0; i < NInfect; i++) {
00017     InfectNIdV.Add(NIdV[i]); }
00018 }
00019 
00020 // burn each link independently (forward with FwdBurnProb, backward with BckBurnProb)
00021 void TForestFire::BurnExpFire() {
00022   const double OldFwdBurnProb = FwdBurnProb;
00023   const double OldBckBurnProb = BckBurnProb;
00024   const int NInfect = InfectNIdV.Len();
00025   const TNGraph& G = *Graph;
00026   TIntH BurnedNIdH;               // burned nodes
00027   TIntV BurningNIdV = InfectNIdV; // currently burning nodes
00028   TIntV NewBurnedNIdV;            // nodes newly burned in current step
00029   bool HasAliveNbrs;              // has unburned neighbors
00030   int NBurned = NInfect, NDiedFire=0;
00031   for (int i = 0; i < InfectNIdV.Len(); i++) {
00032     BurnedNIdH.AddDat(InfectNIdV[i]); }
00033   NBurnedTmV.Clr(false);  NBurningTmV.Clr(false);  NewBurnedTmV.Clr(false);
00034   for (int time = 0; ; time++) {
00035     NewBurnedNIdV.Clr(false);
00036     // for each burning node
00037     for (int node = 0; node < BurningNIdV.Len(); node++) {
00038       const int& BurningNId = BurningNIdV[node];
00039       const TNGraph::TNodeI Node = G.GetNI(BurningNId);
00040       HasAliveNbrs = false;
00041       NDiedFire = 0;
00042       // burn forward links  (out-links)
00043       for (int e = 0; e < Node.GetOutDeg(); e++) {
00044         const int OutNId = Node.GetOutNId(e);
00045         if (! BurnedNIdH.IsKey(OutNId)) { // not yet burned
00046           HasAliveNbrs = true;
00047           if (Rnd.GetUniDev() < FwdBurnProb) {
00048             BurnedNIdH.AddDat(OutNId);  NewBurnedNIdV.Add(OutNId);  NBurned++; }
00049         }
00050       }
00051       // burn backward links (in-links)
00052       if (BckBurnProb > 0.0) {
00053         for (int e = 0; e < Node.GetInDeg(); e++) {
00054           const int InNId = Node.GetInNId(e);
00055           if (! BurnedNIdH.IsKey(InNId)) { // not yet burned
00056             HasAliveNbrs = true;
00057             if (Rnd.GetUniDev() < BckBurnProb) {
00058               BurnedNIdH.AddDat(InNId);  NewBurnedNIdV.Add(InNId);  NBurned++; }
00059           }
00060         }
00061       }
00062       if (! HasAliveNbrs) { NDiedFire++; }
00063     }
00064     NBurnedTmV.Add(NBurned);
00065     NBurningTmV.Add(BurningNIdV.Len() - NDiedFire);
00066     NewBurnedTmV.Add(NewBurnedNIdV.Len());
00067     //BurningNIdV.AddV(NewBurnedNIdV);   // node is burning eternally
00068     BurningNIdV.Swap(NewBurnedNIdV);    // node is burning just 1 time step
00069     if (BurningNIdV.Empty()) break;
00070     FwdBurnProb = FwdBurnProb * ProbDecay;
00071     BckBurnProb = BckBurnProb * ProbDecay;
00072   }
00073   BurnedNIdV.Gen(BurnedNIdH.Len(), 0);
00074   for (int i = 0; i < BurnedNIdH.Len(); i++) {
00075     BurnedNIdV.Add(BurnedNIdH.GetKey(i)); }
00076   FwdBurnProb = OldFwdBurnProb;
00077   BckBurnProb = OldBckBurnProb;
00078 }
00079 
00080 // Node selects N~geometric(1.0-FwdBurnProb)-1 out-links and burns them. Then same for in-links.
00081 // geometirc(p) has mean 1/(p), so for given FwdBurnProb, we burn 1/(1-FwdBurnProb)
00082 void TForestFire::BurnGeoFire() {
00083   const double OldFwdBurnProb=FwdBurnProb;
00084   const double OldBckBurnProb=BckBurnProb;
00085   const int& NInfect = InfectNIdV.Len();
00086   const TNGraph& G = *Graph;
00087   TIntH BurnedNIdH;               // burned nodes
00088   TIntV BurningNIdV = InfectNIdV; // currently burning nodes
00089   TIntV NewBurnedNIdV;            // nodes newly burned in current step
00090   bool HasAliveInNbrs, HasAliveOutNbrs; // has unburned neighbors
00091   TIntV AliveNIdV;                // NIds of alive neighbors
00092   int NBurned = NInfect, time;
00093   for (int i = 0; i < InfectNIdV.Len(); i++) {
00094     BurnedNIdH.AddDat(InfectNIdV[i]); }
00095   NBurnedTmV.Clr(false);  NBurningTmV.Clr(false);  NewBurnedTmV.Clr(false);
00096   for (time = 0; ; time++) {
00097     NewBurnedNIdV.Clr(false);
00098     for (int node = 0; node < BurningNIdV.Len(); node++) {
00099       const int& BurningNId = BurningNIdV[node];
00100       const TNGraph::TNodeI Node = G.GetNI(BurningNId);
00101       // find unburned links
00102       HasAliveOutNbrs = false;
00103       AliveNIdV.Clr(false); // unburned links
00104       for (int e = 0; e < Node.GetOutDeg(); e++) {
00105         const int OutNId = Node.GetOutNId(e);
00106         if (! BurnedNIdH.IsKey(OutNId)) {
00107           HasAliveOutNbrs = true;  AliveNIdV.Add(OutNId); }
00108       }
00109       // number of links to burn (geometric coin). Can also burn 0 links
00110       const int BurnNFwdLinks = Rnd.GetGeoDev(1.0-FwdBurnProb) - 1;
00111       if (HasAliveOutNbrs && BurnNFwdLinks > 0) {
00112         AliveNIdV.Shuffle(Rnd);
00113         for (int i = 0; i < TMath::Mn(BurnNFwdLinks, AliveNIdV.Len()); i++) {
00114           BurnedNIdH.AddDat(AliveNIdV[i]);
00115           NewBurnedNIdV.Add(AliveNIdV[i]);  NBurned++; }
00116       }
00117       // backward links
00118       if (BckBurnProb > 0.0) {
00119         // find unburned links
00120         HasAliveInNbrs = false;
00121         AliveNIdV.Clr(false);
00122         for (int e = 0; e < Node.GetInDeg(); e++) {
00123           const int InNId = Node.GetInNId(e);
00124           if (! BurnedNIdH.IsKey(InNId)) {
00125             HasAliveInNbrs = true;  AliveNIdV.Add(InNId); }
00126         }
00127          // number of links to burn (geometric coin). Can also burn 0 links
00128         const int BurnNBckLinks = Rnd.GetGeoDev(1.0-BckBurnProb) - 1;
00129         if (HasAliveInNbrs && BurnNBckLinks > 0) {
00130           AliveNIdV.Shuffle(Rnd);
00131           for (int i = 0; i < TMath::Mn(BurnNBckLinks, AliveNIdV.Len()); i++) {
00132             BurnedNIdH.AddDat(AliveNIdV[i]);
00133             NewBurnedNIdV.Add(AliveNIdV[i]);  NBurned++; }
00134         }
00135       }
00136     }
00137     NBurnedTmV.Add(NBurned);  NBurningTmV.Add(BurningNIdV.Len());  NewBurnedTmV.Add(NewBurnedNIdV.Len());
00138     // BurningNIdV.AddV(NewBurnedNIdV);   // node is burning eternally
00139     BurningNIdV.Swap(NewBurnedNIdV);   // node is burning just 1 time step
00140     if (BurningNIdV.Empty()) break;
00141     FwdBurnProb = FwdBurnProb * ProbDecay;
00142     BckBurnProb = BckBurnProb * ProbDecay;
00143   }
00144   BurnedNIdV.Gen(BurnedNIdH.Len(), 0);
00145   for (int i = 0; i < BurnedNIdH.Len(); i++) {
00146     BurnedNIdV.Add(BurnedNIdH.GetKey(i)); }
00147   FwdBurnProb = OldFwdBurnProb;
00148   BckBurnProb = OldBckBurnProb;
00149 }
00150 
00151 /*// exact implementation of the algorithm described in KDD '05 paper
00152 void TForestFire::BurnGeoFire() {
00153   const double OldFwdBurnProb=FwdBurnProb;
00154   const double OldBckBurnProb=BckBurnProb;
00155   const double ProbRatio = BckBurnProb/FwdBurnProb;
00156   const int NInfect = InfectNIdV.Len();
00157   const TNGraph& G = *Graph;
00158   TIntH BurnedNIdH;               // burned nodes
00159   TIntV BurningNIdV = InfectNIdV; // currently burning nodes
00160   TIntV NewBurnedNIdV;            // nodes newly burned in current step
00161   bool HasAliveInNbrs, HasAliveOutNbrs; // has unburned neighbors
00162   TIntV AliveNIdV;                // NIds of alive neighbors
00163   int NBurned = NInfect, time;
00164   for (int i = 0; i < InfectNIdV.Len(); i++) {
00165     BurnedNIdH.AddDat(InfectNIdV[i]); }
00166   NBurnedTmV.Clr(false);  NBurningTmV.Clr(false);  NewBurnedTmV.Clr(false);
00167   for (time = 0; ; time++) {
00168     NewBurnedNIdV.Clr(false);
00169     for (int node = 0; node < BurningNIdV.Len(); node++) {
00170       const int& BurningNId = BurningNIdV[node];
00171       const int BurnNLinks = Rnd.GetGeoDev(1.0-(FwdBurnProb)) - 1;
00172       const TNGraph::TNode& Node = G.GetNode(BurningNId);
00173       // find unburned links
00174       if (BurnNLinks > 0) {
00175         AliveNIdV.Clr(false);
00176         for (int e = 0; e < Node.GetOutDeg(); e++) {
00177           const int OutNId = Node.GetOutNId(e);
00178           if (! BurnedNIdH.IsKey(OutNId)) { HasAliveOutNbrs=true;  AliveNIdV.Add(OutNId); }
00179         }
00180         for (int e = 0; e < Node.GetInDeg(); e++) {
00181           const int InNId = Node.GetInNId(e);
00182           if (! BurnedNIdH.IsKey(InNId)) { HasAliveInNbrs=true;  AliveNIdV.Add(-InNId); }
00183         }
00184         AliveNIdV.Shuffle(Rnd);
00185         // add links
00186         for (int e = i; i < AliveNIdV.Len(); i++) {
00187           if (AliveNIdV[i] > 0) BurnedNIdH.AddDat(AliveNIdV[i]);
00188           NewBurnedNIdV.Add(AliveNIdV[i]);  NBurned++;
00189         }
00190       }
00191       HasAliveOutNbrs = false;
00192       AliveNIdV.Clr(false); // unburned links
00193       for (int e = 0; e < Node.GetOutDeg(); e++) {
00194         const int OutNId = Node.GetOutNId(e);
00195         if (! BurnedNIdH.IsKey(OutNId)) {
00196           HasAliveOutNbrs = true;  AliveNIdV.Add(OutNId); }
00197       }
00198       // number of links to burn (geometric coin). Can also burn 0 links
00199       const int BurnNFwdLinks = Rnd.GetGeoDev(1.0-FwdBurnProb) - 1;
00200       if (HasAliveOutNbrs && BurnNFwdLinks > 0) {
00201         AliveNIdV.Shuffle(Rnd);
00202         for (int i = 0; i < TMath::Mn(BurnNFwdLinks, AliveNIdV.Len()); i++) {
00203           BurnedNIdH.AddDat(AliveNIdV[i]);
00204           NewBurnedNIdV.Add(AliveNIdV[i]);  NBurned++; }
00205       }
00206       // backward links
00207       if (BckBurnProb > 0.0) {
00208         // find unburned links
00209         HasAliveInNbrs = false;
00210         AliveNIdV.Clr(false);
00211         for (int e = 0; e < Node.GetInDeg(); e++) {
00212           const int InNId = Node.GetInNId(e);
00213           if (! BurnedNIdH.IsKey(InNId)) {
00214             HasAliveInNbrs = true;  AliveNIdV.Add(InNId); }
00215         }
00216          // number of links to burn (geometric coin). Can also burn 0 links
00217         const int BurnNBckLinks = Rnd.GetGeoDev(1.0-BckBurnProb) - 1;
00218         if (HasAliveInNbrs && BurnNBckLinks > 0) {
00219           AliveNIdV.Shuffle(Rnd);
00220           for (int i = 0; i < TMath::Mn(BurnNBckLinks, AliveNIdV.Len()); i++) {
00221             BurnedNIdH.AddDat(AliveNIdV[i]);
00222             NewBurnedNIdV.Add(AliveNIdV[i]);  NBurned++; }
00223         }
00224       }
00225     }
00226     NBurnedTmV.Add(NBurned);  NBurningTmV.Add(BurningNIdV.Len());  NewBurnedTmV.Add(NewBurnedNIdV.Len());
00227     // BurningNIdV.AddV(NewBurnedNIdV);   // node is burning eternally
00228     BurningNIdV.Swap(NewBurnedNIdV);   // node is burning just 1 time step
00229     if (BurningNIdV.Empty()) break;
00230     FwdBurnProb = FwdBurnProb * ProbDecay;
00231     BckBurnProb = BckBurnProb * ProbDecay;
00232   }
00233   BurnedNIdV.Gen(BurnedNIdH.Len(), 0);
00234   for (int i = 0; i < BurnedNIdH.Len(); i++) {
00235     BurnedNIdV.Add(BurnedNIdH.GetKey(i)); }
00236   FwdBurnProb = OldFwdBurnProb;
00237   BckBurnProb = OldBckBurnProb;
00238 }//*/
00239 
00240 void TForestFire::PlotFire(const TStr& FNmPref, const TStr& Desc, const bool& PlotAllBurned) {
00241   TGnuPlot GnuPlot(FNmPref, TStr::Fmt("%s. ForestFire. G(%d, %d). Fwd:%g  Bck:%g  NInfect:%d",
00242     Desc.CStr(), Graph->GetNodes(), Graph->GetEdges(), FwdBurnProb(), BckBurnProb(), InfectNIdV.Len()));
00243   GnuPlot.SetXYLabel("TIME EPOCH", "Number of NODES");
00244   if (PlotAllBurned) GnuPlot.AddPlot(NBurnedTmV, gpwLinesPoints, "All burned nodes till time");
00245   GnuPlot.AddPlot(NBurningTmV, gpwLinesPoints, "Burning nodes at time");
00246   GnuPlot.AddPlot(NewBurnedTmV, gpwLinesPoints, "Newly burned nodes at time");
00247   GnuPlot.SavePng(TFile::GetUniqueFNm(TStr::Fmt("fireSz.%s_#.png", FNmPref.CStr())));
00248 }
00249 
00250 PNGraph TForestFire::GenGraph(const int& Nodes, const double& FwdProb, const double& BckProb) {
00251   TFfGGen Ff(false, 1, FwdProb, BckProb, 1.0, 0.0, 0.0);
00252   Ff.GenGraph(Nodes);
00253   return Ff.GetGraph();
00254 }
00255 
00257 // Forest Fire Graph Generator
00258 int TFfGGen::TimeLimitSec = 30*60; // 30 minutes
00259 
00260 TFfGGen::TFfGGen(const bool& BurnExpFireP, const int& StartNNodes, const double& ForwBurnProb,
00261                  const double& BackBurnProb, const double& DecayProb, const double& Take2AmbasPrb, const double& OrphanPrb) :
00262  Graph(), BurnExpFire(BurnExpFireP), StartNodes(StartNNodes), FwdBurnProb(ForwBurnProb),
00263  BckBurnProb(BackBurnProb), ProbDecay(DecayProb), Take2AmbProb(Take2AmbasPrb), OrphanProb(OrphanPrb) {
00264   //IAssert(ProbDecay == 1.0); // do not need Decay any more
00265 }
00266 
00267 TStr TFfGGen::GetParamStr() const {
00268   return TStr::Fmt("%s  FWD:%g  BCK:%g, StartNds:%d, Take2:%g, Orphan:%g, ProbDecay:%g",
00269     BurnExpFire?"EXP":"GEO", FwdBurnProb(), BckBurnProb(), StartNodes(), Take2AmbProb(), OrphanProb(), ProbDecay());
00270 }
00271 
00272 TFfGGen::TStopReason TFfGGen::AddNodes(const int& GraphNodes, const bool& FloodStop) {
00273   printf("\n***ForestFire:  %s  Nodes:%d  StartNodes:%d  Take2AmbProb:%g\n", BurnExpFire?"ExpFire":"GeoFire", GraphNodes, StartNodes(), Take2AmbProb());
00274   printf("                FwdBurnP:%g  BckBurnP:%g  ProbDecay:%g  Orphan:%g\n", FwdBurnProb(), BckBurnProb(), ProbDecay(), OrphanProb());
00275   TExeTm ExeTm;
00276   int Burned1=0, Burned2=0, Burned3=0; // last 3 fire sizes
00277   // create initial set of nodes
00278   if (Graph.Empty()) { Graph = PNGraph::New(); }
00279   if (Graph->GetNodes() == 0) {
00280     for (int n = 0; n < StartNodes; n++) { Graph->AddNode(); }
00281   }
00282   int NEdges = Graph->GetEdges();
00283   // forest fire
00284   TRnd Rnd(0);
00285   TForestFire ForestFire(Graph, FwdBurnProb, BckBurnProb, ProbDecay, 0);
00286   // add nodes
00287   for (int NNodes = Graph->GetNodes()+1; NNodes <= GraphNodes; NNodes++) {
00288     const int NewNId = Graph->AddNode(-1);
00289     IAssert(NewNId == Graph->GetNodes()-1); // node ids have to be 0...N
00290     // not an Orphan (burn fire)
00291     if (OrphanProb == 0.0 || Rnd.GetUniDev() > OrphanProb) {
00292       // infect ambassadors
00293       if (Take2AmbProb == 0.0 || Rnd.GetUniDev() > Take2AmbProb || NewNId < 2) {
00294         ForestFire.Infect(Rnd.GetUniDevInt(NewNId)); // take 1 ambassador
00295       } else {
00296         const int AmbassadorNId1 = Rnd.GetUniDevInt(NewNId);
00297         int AmbassadorNId2 = Rnd.GetUniDevInt(NewNId);
00298         while (AmbassadorNId1 == AmbassadorNId2) {
00299           AmbassadorNId2 = Rnd.GetUniDevInt(NewNId); }
00300         ForestFire.Infect(TIntV::GetV(AmbassadorNId1, AmbassadorNId2)); // take 2 ambassadors
00301       }
00302       // burn fire
00303       if (BurnExpFire) { ForestFire.BurnExpFire(); }
00304       else { ForestFire.BurnGeoFire(); }
00305       // add edges to burned nodes
00306       for (int e = 0; e < ForestFire.GetBurned(); e++) {
00307         Graph->AddEdge(NewNId, ForestFire.GetBurnedNId(e));
00308         NEdges++;
00309       }
00310       Burned1=Burned2;  Burned2=Burned3;  Burned3=ForestFire.GetBurned();
00311     } else {
00312       // Orphan (zero out-links)
00313       Burned1=Burned2;  Burned2=Burned3;  Burned3=0;
00314     }
00315     if (NNodes % Kilo(1) == 0) {
00316       printf("(%d, %d)  burned: [%d,%d,%d]  [%s]\n", NNodes, NEdges, Burned1, Burned2, Burned3, ExeTm.GetStr()); }
00317     if (FloodStop && NEdges>GraphNodes && (NEdges/double(NNodes)>1000.0)) { // average node degree is more than 500
00318       printf(". FLOOD. G(%6d, %6d)\n", NNodes, NEdges);  return srFlood; }
00319     if (NNodes % 1000 == 0 && TimeLimitSec > 0 && ExeTm.GetSecs() > TimeLimitSec) {
00320       printf(". TIME LIMIT. G(%d, %d)\n", Graph->GetNodes(), Graph->GetEdges());
00321       return srTimeLimit; }
00322   }
00323   IAssert(Graph->GetEdges() == NEdges);
00324   return srOk;
00325 }
00326 
00327 TFfGGen::TStopReason TFfGGen::GenGraph(const int& GraphNodes, const bool& FloodStop) {
00328   Graph = PNGraph::New();
00329   return AddNodes(GraphNodes, FloodStop);
00330 }
00331 
00332 TFfGGen::TStopReason TFfGGen::GenGraph(const int& GraphNodes, PGStatVec& EvolStat, const bool& FloodStop) {
00333   int GrowthStatNodes = 100;
00334   Graph = PNGraph::New();
00335   AddNodes(StartNodes);
00336   TStopReason SR = srUndef;
00337   while (Graph->GetNodes() < GraphNodes) {
00338     SR = AddNodes(GrowthStatNodes, FloodStop);
00339     if (SR != srOk) { return SR; }
00340     EvolStat->Add(Graph, TSecTm(Graph->GetNodes()));
00341     GrowthStatNodes = int(1.5*GrowthStatNodes);
00342   }
00343   return SR;
00344 }
00345 
00346 void TFfGGen::PlotFireSize(const TStr& FNmPref, const TStr& DescStr) {
00347   TGnuPlot GnuPlot("fs."+FNmPref, TStr::Fmt("%s. Fire size. G(%d, %d)",
00348     DescStr.CStr(), Graph->GetNodes(), Graph->GetEdges()));
00349   GnuPlot.SetXYLabel("Vertex id (iterations)", "Fire size (node out-degree)");
00350   TFltPrV IdToOutDegV;
00351   for (TNGraph::TNodeI NI = Graph->BegNI(); NI < Graph->EndNI(); NI++) {
00352     IdToOutDegV.Add(TFltPr(NI.GetId(), NI.GetOutDeg())); }
00353   IdToOutDegV.Sort();
00354   GnuPlot.AddPlot(IdToOutDegV, gpwImpulses, "Node out-degree");
00355   GnuPlot.SavePng();
00356 }
00357 
00358 void TFfGGen::GenFFGraphs(const double& FProb, const double& BProb, const TStr& FNm) {
00359   const int NRuns = 10;
00360   const int NNodes = 10000;
00361   TGStat::NDiamRuns = 10;
00362   //const double FProb = 0.35, BProb = 0.20;  // ff1
00363   //const double FProb = 0.37, BProb = 0.32;  // ff2
00364   //const double FProb = 0.37, BProb = 0.325; // ff22
00365   //const double FProb = 0.37, BProb = 0.33;  // ff3
00366   //const double FProb = 0.37, BProb = 0.35;  // ff4
00367   //const double FProb = 0.38, BProb = 0.35;  // ff5
00368   TVec<PGStatVec> GAtTmV;
00369   TFfGGen FF(false, 1, FProb, BProb, 1.0, 0, 0);
00370   for (int r = 0; r < NRuns; r++) {
00371     PGStatVec GV = TGStatVec::New(tmuNodes, TGStat::AllStat());
00372     FF.GenGraph(NNodes, GV, true);
00373     for (int i = 0; i < GV->Len(); i++) {
00374       if (i == GAtTmV.Len()) {
00375         GAtTmV.Add(TGStatVec::New(tmuNodes, TGStat::AllStat()));
00376       }
00377       GAtTmV[i]->Add(GV->At(i));
00378     }
00379     IAssert(GAtTmV.Len() == GV->Len());
00380   }
00381   PGStatVec AvgStat = TGStatVec::New(tmuNodes, TGStat::AllStat());
00382   for (int i = 0; i < GAtTmV.Len(); i++) {
00383     AvgStat->Add(GAtTmV[i]->GetAvgGStat(false));
00384   }
00385   AvgStat->PlotAllVsX(gsvNodes, FNm, TStr::Fmt("Forest Fire: F:%g  B:%g (%d runs)", FProb, BProb, NRuns));
00386   AvgStat->Last()->PlotAll(FNm, TStr::Fmt("Forest Fire: F:%g  B:%g (%d runs)", FProb, BProb, NRuns));
00387 }
00388 
00389 /*/////////////////////////////////////////////////
00390 // Forest Fire Phase Transition
00391 TFfPhaseTrans::TFfPhaseTrans(const int& NNds, const int& StartNds,  const double& Take2AmbPr,
00392                              const double& ProbOrphan, const double& ProbIncrement, const int& NRunsPerFB) :
00393   BurExpFire(false), NNodes(NNds), StartNodes(StartNds), NRuns(NRunsPerFB), Take2AmbProb(Take2AmbPr),
00394   OrphanProb(ProbOrphan), ProbInc(ProbIncrement), FBPrGSetH(), FBPrGStatH() {
00395   TakeStatSet = TFSet() | gsEffDiam;
00396 }
00397 
00398 PFfPhaseTrans TFfPhaseTrans::New(const int& NNds, const int& StartNds,  const double& Take2AmbPr,
00399                                  const double& ProbOrphan, const double& ProbIncrement, const int& NRunsPerFB) {
00400   return new TFfPhaseTrans(NNds, StartNds, Take2AmbPr, ProbOrphan, ProbIncrement, NRunsPerFB);
00401 }
00402 
00403 TFfPhaseTrans::TFfPhaseTrans(TSIn& SIn) : BurExpFire(SIn), NNodes(SIn), StartNodes(SIn), NRuns(SIn),
00404   Take2AmbProb(SIn), OrphanProb(SIn), ProbInc(SIn), FBPrGSetH(SIn), FBPrGStatH(SIn), TakeStatSet(SIn) {
00405 }
00406 PFfPhaseTrans TFfPhaseTrans::Load(TSIn& SIn) {
00407   return new TFfPhaseTrans(SIn);
00408 }
00409 
00410 void TFfPhaseTrans::Save(TFOut& SOut) const {
00411   BurExpFire.Save(SOut);  NNodes.Save(SOut);
00412   StartNodes.Save(SOut);  NRuns.Save(SOut);
00413   Take2AmbProb.Save(SOut);  OrphanProb.Save(SOut);
00414   ProbInc.Save(SOut);
00415   FBPrGSetH.Save(SOut);  FBPrGStatH.Save(SOut);
00416   TakeStatSet.Save(SOut);
00417 }
00418 
00419 TStr TFfPhaseTrans::GetTitleStr(const int& ValN) const {
00420   const PGrowthStat AvgStat = GetAvgGStat(ValN);
00421   const double AvgDeg = 2.0*AvgStat->Last()->Edges / double(AvgStat->Last()->Nodes);
00422   const double DiamCf = GetDecDiam(ValN).Val1;
00423   const TFltPr GplCf = GetGplCf(ValN);
00424   return TStr::Fmt("FWD: %.4f  BCK: %.4f  DIAM-INC: %.2f  GPL: %.2f (%.2f) AvgDeg:%.1f",
00425     GetFProb(ValN), GetBProb(ValN), DiamCf, GplCf.Val1(), GplCf.Val2(), AvgDeg);
00426 }
00427 
00428 TStr TFfPhaseTrans::GetFNm(const TStr& FNmPref) const {
00429   return TStr::Fmt("%s.n%dk.s%d", FNmPref.CStr(), int(NNodes/1000.0), StartNodes());
00430 }
00431 
00432 TFltPr TFfPhaseTrans::GetGplCf(const int& ValN) const {
00433   const TGrowthSet& GrowthSet = *GetGSet(ValN);
00434   TMom GplMom;
00435   for (int i = 0; i < GrowthSet.Len(); i++) {
00436     GplMom.Add(GrowthSet[i]->GetGplCf());
00437   }
00438   GplMom.Def();
00439   return TFltPr(GplMom.GetMean(), GplMom.GetSDev());
00440 }
00441 
00442 TFltPr TFfPhaseTrans::GetDecDiam(const int& ValN) const {
00443   const TGrowthSet& GrowthSet = *GetGSet(ValN);
00444   TMom DiamMom;
00445   for (int i = 0; i < GrowthSet.Len(); i++) {
00446     DiamMom.Add(GrowthSet[i]->IsDecEffDiam());
00447   }
00448   DiamMom.Def();
00449   return TFltPr(DiamMom.GetMean(), DiamMom.GetSDev());
00450 }
00451 
00452 TFltPr TFfPhaseTrans::GetEffDiam(const int& ValN, const int& AtTick) const {
00453   const TGrowthSet& GrowthSet = *GetGSet(ValN);
00454   TMom DiamMom;
00455   for (int i = 0; i < GrowthSet.Len(); i++) {
00456     if (AtTick != -1) { DiamMom.Add(GrowthSet[i]->At(AtTick)->EffDiam); }
00457     else { DiamMom.Add(GrowthSet[i]->Last()->EffDiam); }
00458   }
00459   DiamMom.Def();
00460   return TFltPr(DiamMom.GetMean(), DiamMom.GetSDev());
00461 }
00462 
00463 void TFfPhaseTrans::GetFProbV(TFltV& FProbV) const {
00464   THash<TFlt, TInt> FProbH;
00465   for (int i = 0; i < Len(); i++) {
00466     FProbH.AddKey(GetFProb(i)); }
00467   FProbH.GetKeyV(FProbV);
00468   FProbV.Sort();
00469 }
00470 
00471 TFltPr TFfPhaseTrans::RunForestFire(double FwdProb, double BckProb, const bool& Plot) {
00472   FwdProb = TMath::Round(FwdProb, 4);
00473   BckProb = TMath::Round(BckProb, 4);
00474   printf("\n---FwdBurnProb--%g----Back--%g------------------[%s]\n", FwdProb, BckProb, TExeTm::GetCurTm());
00475   if (FBPrGStatH.IsKey(TFltPr(FwdProb, BckProb))) {
00476     printf("*** SKIP -- already have it\n");
00477     const TGrowthStat& Stat = *FBPrGStatH.GetDat(TFltPr(FwdProb, BckProb));
00478     return TFltPr(Stat.GetGplCf(), Stat.IsDecEffDiam());
00479   }
00480   if (BckProb > 1.0) return TFltPr(-1, -1);
00481   PGrowthSet GrowthSet = TGrowthSet::New();
00482   TExeTm ExeTm;
00483   TFfGGen ForestFire(false, StartNodes, FwdProb, BckProb, 1.0, Take2AmbProb, OrphanProb);
00484   ForestFire.TakeStat(TakeStatSet); // gsDiam
00485   GrowthSet->Clr(false);
00486   for (int run = 0; run < NRuns; run++) {
00487     printf("\nRUN %d        [last run %s]\n", run+1, ExeTm.GetTmStr());  ExeTm.Tick();
00488     TFfGGen::TStopReason StopReason =*//* ForestFire.GenGraph(NNodes, true);
00489     if (! ForestFire.GetGrowthStat()->Empty()) {
00490       GrowthSet->Add(ForestFire.GetGrowthStat()); }
00491   }
00492   IAssert(! GrowthSet.Empty());
00493   // add stat
00494   FBPrGSetH.AddDat(TFltPr(FwdProb, BckProb), GrowthSet);
00495   PGrowthStat AvgStat = TGrowthStat::New();
00496   AvgStat->AvgGrowthStat(*GrowthSet);
00497   FBPrGStatH.AddDat(TFltPr(FwdProb, BckProb), AvgStat);
00498   const double DiamCf = LastDecDiam().Val1;
00499   const double GplCf = LastGplCf().Val1;
00500   // plot
00501   if (Plot && ! AvgStat.Empty()) {
00502     const TStr FNm = TStr::Fmt("F%04d.B%04d", int(FwdProb*10000), int(BckProb*10000));
00503     const TStr Title = GetTitleStr(Len()-1);
00504     AvgStat->PlotDiam(FNm, Title);
00505     AvgStat->PlotGpl(FNm, Title, true, false, false, false, false);
00506   }
00507   return TFltPr(GplCf, DiamCf);
00508 }
00509 
00510 void TFfPhaseTrans::FwdProbSteps(const double& MinFwdProb, const double& MaxFwdProb, const double& BckProb) {
00511   const TStr BinFNm = TFile::GetUniqueFNm(TStr::Fmt("phaseFwd.n%dk.s%d_#.bin", int(NNodes/1000.0), StartNodes()));
00512   TFfGGen::TimeLimitSec = 10*60;
00513   for (double FwdProb = MinFwdProb; FwdProb <= MaxFwdProb; FwdProb += ProbInc) {
00514     RunForestFire(FwdProb, BckProb, true);
00515     { TFOut FOut(BinFNm);
00516     Save(FOut); }
00517   }
00518 }
00519 
00520 void TFfPhaseTrans::FwdProbStepsFact(const double& MinFwdProb, const double& MaxFwdProb, const double& BckFact) {
00521   const TStr BinFNm = TFile::GetUniqueFNm(TStr::Fmt("phaseFwd.n%dk.s%d_#.bin", int(NNodes/1000.0), StartNodes()));
00522   TFfGGen::TimeLimitSec = 10*60;
00523   for (double FwdProb = MinFwdProb; FwdProb <= MaxFwdProb; FwdProb += ProbInc) {
00524     RunForestFire(FwdProb, BckFact*FwdProb, true);
00525     { TFOut FOut(BinFNm);
00526     Save(FOut); }
00527   }
00528 }
00529 
00530 void TFfPhaseTrans::FwdBckPhasePlot(const double& MinFwdProb, const double& MaxFwdProb, const double& MinBckFact,
00531                                     const double& MaxBckFact, const int& TimeLimitSec) {
00532   const TStr BinFNm = TFile::GetUniqueFNm(TStr::Fmt("phaseFF.n%dk.s%d_#.bin", int(NNodes/1000.0), StartNodes()));
00533   TFfGGen::TimeLimitSec = TimeLimitSec;
00534   for (double FwdProb = MinFwdProb; FwdProb <= MaxFwdProb; FwdProb += ProbInc) {
00535     for (double BckFact = MinBckFact; BckFact < MaxBckFact+0.001; BckFact += ProbInc) {
00536       const double BckProb = FwdProb * BckFact;
00537       RunForestFire(FwdProb, BckProb, true);
00538       { TFOut FOut(BinFNm);
00539       Save(FOut); }
00540     }
00541   }
00542 }
00543 
00544 void TFfPhaseTrans::FindGplPhaseTr(const double& StartFProb, const double& FollowGplCf, const TStr& FNmPref, const TStr& Desc) {
00545   const TStr FNm = TStr::Fmt("phGPL.%s", GetFNm(FNmPref).CStr());
00546   const double Tolerance = 0.01;
00547   const double MinDelta = 0.001;
00548   const bool Plot = false;
00549   TFfGGen::TimeLimitSec = 10*60;
00550   TGrowthStat::MinNodesEdges = 2*(StartNodes-1)+100;
00551   const int OldNDiamRuns = TGraphStat::NDiamRuns;
00552   TGraphStat::NDiamRuns = 1;  //!!! diameter
00553   TakeStat(TFSet() | gsEffDiam);
00554   FILE *F = fopen((FNm+".txt").CStr(), "wt");
00555   fprintf(F, "FollowGplCf:  %g\n", FollowGplCf);
00556   fprintf(F, "StartNodes:   %d\n", StartNodes());
00557   fprintf(F, "Take2AmbProb: %g\n", Take2AmbProb());
00558   fprintf(F, "OrphanProb:   %g\n", OrphanProb());
00559   fprintf(F, "Tolerance:    %g\n", Tolerance);
00560   double FwdProb = StartFProb, LBckRat=0, RBckRat=1, BckRat, GplCf;
00561   //TFltPr GplPr;
00562   while (FwdProb < 1.0) {
00563     while (true) {
00564       BckRat = (RBckRat+LBckRat) / 2;
00565       fprintf(F, "FWD: %g, (%f -- %f)", FwdProb, LBckRat, RBckRat);
00566       GplCf = RunForestFire(FwdProb, FwdProb*BckRat, Plot).Val1;
00567       IAssert(GplCf != -1);
00568       fprintf(F, "  %f  gpl: %.4f (%.4f)", BckRat, GplCf, LastGplCf().Val2());
00569       if (TMath::IsInEps(GplCf - FollowGplCf, Tolerance)) { fprintf(F, "  ***\n"); break; }
00570       if (RBckRat-LBckRat < MinDelta) { fprintf(F, "  gap\n"); break; }
00571       if (GplCf > FollowGplCf) { RBckRat = BckRat; } else { LBckRat = BckRat; }
00572       fprintf(F, "\n");  fflush(F);
00573     }
00574     FwdProb += ProbInc;
00575     RBckRat = BckRat+0.01;
00576     if (RBckRat > 1.0) RBckRat = 1.0;
00577     LBckRat = RBckRat - 0.2;
00578     if (LBckRat < 0.0) LBckRat = 0.0;
00579     { TFOut FOut(FNm+".bin");
00580      Save(FOut); }
00581     SaveGplPhaseTr(FollowGplCf, FNmPref, Desc);
00582     fprintf(F, "\n");
00583   }
00584   fclose(F);
00585   TGraphStat::NDiamRuns = OldNDiamRuns;
00586 }
00587 
00588 void TFfPhaseTrans::SaveGplPhaseTr(const double& FollowGplCf, const TStr& FNmPref, const TStr& Desc) {
00589   const double Tolerance = 0.02;
00590   THash<TFlt,  TIntFltPr> FProbH;
00591   for (int i = 0; i < Len(); i ++) {
00592     const double FProb = GetFProb(i);
00593     const double GplCf = GetGplCf(i).Val1;
00594     if (TMath::IsInEps(GplCf-FollowGplCf, Tolerance)) {
00595       if (! FProbH.IsKey(FProb)) {
00596         FProbH.AddDat(FProb, TIntFltPr(i, GplCf)); }
00597       else {
00598         const double bestCf = FProbH.GetDat(FProb).Val2;
00599         if (fabs(bestCf - FollowGplCf) > fabs(GplCf - FollowGplCf)) {
00600           FProbH.AddDat(FProb, TIntFltPr(i, GplCf)); }
00601       }
00602     }
00603   }
00604   TVec<TPair<TFlt, TIntFltPr> > FProbV;
00605   FProbH.GetKeyDatPrV(FProbV);  FProbV.Sort();
00606   const bool HasDiam = TakeStatSet.In(gsEffDiam);
00607   FILE *F = fopen(TStr::Fmt("phGPL.%s.tab", GetFNm(FNmPref).CStr()).CStr(), "wt");
00608   if (! Desc.Empty()) fprintf(F, "%s\n", Desc.CStr());
00609   fprintf(F, "FollowGplCf:  %g\n", FollowGplCf);
00610   fprintf(F, "StartNodes:   %d\n", StartNodes());
00611   fprintf(F, "Take2AmbProb: %g\n", Take2AmbProb());
00612   fprintf(F, "OrphanProb:   %g\n", OrphanProb());
00613   fprintf(F, "Tolerance:    %g\n", Tolerance);
00614   fprintf(F, "id\tFProb\tBProb\tBRatio\tGlp\tGlpDev");
00615   if (HasDiam) {  fprintf(F, "\tDiamCf\tDiamDev\tEffDiam"); }
00616   fprintf(F, "\n");
00617   for (int i = 0; i < FProbH.Len(); i++) {
00618     const int Id = FProbV[i].Val2.Val1;
00619     const TFltPr Gpl = GetGplCf(Id);
00620     fprintf(F, "%d\t%f\t%f\t%f\t", Id, GetFProb(Id), GetBProb(Id), GetBProb(Id)/GetFProb(Id));
00621     fprintf(F, "%f\t%f", Gpl.Val1(), Gpl.Val2());
00622     if (HasDiam) {
00623       const TFltPr DiamCf = GetDecDiam(Id);
00624       fprintf(F, "\t%f\t%f\t%f", DiamCf.Val1(), DiamCf.Val2(), GetEffDiam(Id, -1).Val1());
00625     }
00626     fprintf(F, "\n");
00627   }
00628   fclose(F);
00629 }
00630 
00631 void TFfPhaseTrans::FindDiamPhaseTr(const double& StartFProb, const double& FollowDiamCf, const TStr& FNmPref, const TStr& Desc) {
00632   const TStr FNm = TStr::Fmt("phDIAM.%s", GetFNm(FNmPref).CStr());
00633   const double Tolerance = 0.01;
00634   const double MinDelta = 0.001;
00635   const bool Plot = false;
00636   TFfGGen::TimeLimitSec = 10*60;
00637   const int OldNDiamRuns = TGraphStat::NDiamRuns;
00638   TGraphStat::NDiamRuns = 1;
00639   TGrowthStat::MinNodesEdges = 2*(StartNodes-1)+100;
00640   TakeStat(TFSet() | gsEffDiam);
00641   FILE *F = fopen((FNm+".txt").CStr(), "wt");
00642   fprintf(F, "FollowDiamCf: %g\n", FollowDiamCf);
00643   fprintf(F, "StartNodes:   %d\n", StartNodes());
00644   fprintf(F, "Take2AmbProb: %g\n", Take2AmbProb());
00645   fprintf(F, "OrphanProb:   %g\n", OrphanProb());
00646   fprintf(F, "Tolerance:    %g\n", Tolerance);
00647   double FwdProb = StartFProb, LBckRat=0, RBckRat=1, BckRat, DiamCf;
00648   //TFltPr GplPr;
00649   while (FwdProb < 1.0) {
00650     while (true) {
00651       BckRat = (RBckRat+LBckRat) / 2;
00652       fprintf(F, "FWD: %g, (%f -- %f)", FwdProb, LBckRat, RBckRat);
00653       DiamCf = RunForestFire(FwdProb, FwdProb*BckRat, Plot).Val2;
00654       IAssert(DiamCf != -1);
00655       fprintf(F, "  %f  diam: %.4f (%.4f)", BckRat, DiamCf, LastDecDiam().Val2());
00656       if (TMath::IsInEps(DiamCf - FollowDiamCf, Tolerance)) { fprintf(F, "  ***\n"); break; }
00657       if (RBckRat-LBckRat < MinDelta) { fprintf(F, "  gap\n"); break; }
00658       if (DiamCf < FollowDiamCf) { RBckRat = BckRat; } else { LBckRat = BckRat; }
00659       fprintf(F, "\n");  fflush(F);
00660     }
00661     FwdProb += ProbInc;
00662     RBckRat = BckRat+0.05;
00663     if (RBckRat > 1.0) RBckRat = 1.0;
00664     LBckRat = RBckRat - 0.15;
00665     if (LBckRat < 0.0) LBckRat = 0.0;
00666     { TFOut FOut(FNm+".bin");
00667     Save(FOut); }
00668     SaveDiamPhaseTr(FollowDiamCf, FNmPref, Desc);
00669     fprintf(F, "\n");
00670   }
00671   fclose(F);
00672   TGraphStat::NDiamRuns = OldNDiamRuns;
00673 }
00674 
00675 void TFfPhaseTrans::SaveDiamPhaseTr(const double& FollowDiamCf, const TStr& FNmPref, const TStr& Desc) {
00676   const double Tolerance = 0.03;
00677   THash<TFlt,  TIntFltPr> FProbH;
00678   for (int i = 0; i < Len(); i ++) {
00679     const double FProb = GetFProb(i);
00680     const double DiamCf = GetDecDiam(i).Val1;
00681     if (TMath::IsInEps(DiamCf - FollowDiamCf, Tolerance)) {
00682       if (! FProbH.IsKey(FProb)) {
00683         FProbH.AddDat(FProb, TIntFltPr(i, DiamCf)); }
00684       else {
00685         const double bestCf = FProbH.GetDat(FProb).Val2;
00686         if (fabs(bestCf - FollowDiamCf) > fabs(DiamCf - FollowDiamCf)) {
00687           FProbH.AddDat(FProb, TIntFltPr(i, DiamCf)); }
00688       }
00689     }
00690   }
00691   TVec<TPair<TFlt, TIntFltPr> > FProbV;
00692   FProbH.GetKeyDatPrV(FProbV);  FProbV.Sort();
00693   FILE *F = fopen(TStr::Fmt("phDIAM.%s.tab", GetFNm(FNmPref).CStr()).CStr(), "wt");
00694   if (! Desc.Empty()) fprintf(F, "%s\n", Desc.CStr());
00695   fprintf(F, "FollowDiamCf: %g\n", FollowDiamCf);
00696   fprintf(F, "StartNodes:   %d\n", StartNodes());
00697   fprintf(F, "Take2AmbProb: %g\n", Take2AmbProb());
00698   fprintf(F, "OrphanProb:   %g\n", OrphanProb());
00699   fprintf(F, "Tolerance:    %g\n", Tolerance);
00700   fprintf(F, "id\tFProb\tBProb\tBRatio\tDiamCf\tDiamDev\tGplCf\tGplDev\tEffDiam\n");
00701   for (int i = 0; i < FProbV.Len(); i++) {
00702     const int Id = FProbV[i].Val2.Val1;
00703     const TFltPr DiamCf = GetDecDiam(Id);
00704     const TFltPr GplCf = GetGplCf(Id);
00705     const TFltPr EffDiam = GetEffDiam(Id, -1);
00706     fprintf(F, "%d\t%f\t%f\t%f\t", Id, GetFProb(Id), GetBProb(Id), GetBProb(Id)/GetFProb(Id));
00707     fprintf(F, "%f\t%f\t%f\t%f\t%f\n", DiamCf.Val1(), DiamCf.Val2(), GplCf.Val1(), GplCf.Val2(), EffDiam.Val1());
00708   }
00709   fclose(F);
00710 }
00711 
00712 void TFfPhaseTrans::Merge(const PFfPhaseTrans& FfPhaseTrans) {
00713   Merge(*FfPhaseTrans);
00714 }
00715 
00716 void TFfPhaseTrans::Merge(const TFfPhaseTrans& FfPhaseTrans) {
00717   printf("Merging:\n");
00718   printf("  source      %6d  (Fwd,Bck) pairs\n", FfPhaseTrans.Len());
00719   printf("  destination %6d  (Fwd,Bck) pairs\n", Len());
00720   IAssert(BurExpFire == FfPhaseTrans.BurExpFire);
00721   IAssert(NNodes == FfPhaseTrans.NNodes);
00722   IAssert(StartNodes == FfPhaseTrans.StartNodes);
00723   IAssert(Take2AmbProb == FfPhaseTrans.Take2AmbProb);
00724   IAssert(FBPrGSetH.Len() == FBPrGStatH.Len());
00725   IAssert(FfPhaseTrans.FBPrGSetH.Len() == FfPhaseTrans.FBPrGStatH.Len());
00726   for (int i1 = 0; i1 < FfPhaseTrans.FBPrGSetH.Len(); i1++) {
00727     IAssert(FfPhaseTrans.FBPrGSetH.GetKey(i1) == FfPhaseTrans.FBPrGStatH.GetKey(i1));
00728     const TFltPr& Key = FfPhaseTrans.FBPrGSetH.GetKey(i1);
00729     if (! FBPrGStatH.IsKey(Key)) {
00730       const PGrowthStat Stat = FfPhaseTrans.FBPrGStatH[i1];
00731       const PGrowthSet Set = FfPhaseTrans.FBPrGSetH[i1];
00732       FBPrGStatH.AddDat(Key, Stat);
00733       FBPrGSetH.AddDat(Key, Set);
00734     }
00735   }
00736   printf("  ** merged   %6d  (Fwd,Bck) pairs\n", Len());
00737 }
00738 //*/
00739 
00741 // Undirected Forest Fire (does not densify!)
00742 
00743 // Node selects N~geometric(1.0-BurnProb)-1 links and burns them. 
00744 // geometirc(p) has mean 1/(p), so for given BurnProb, we burn 1/(1-BurnProb) links in average
00745 int TUndirFFire::BurnGeoFire(const int& StartNId) {
00746   BurnedSet.Clr(false);
00747   BurningNIdV.Clr(false);  
00748   NewBurnedNIdV.Clr(false);
00749   AliveNIdV.Clr(false);
00750   const TUNGraph& G = *Graph;
00751   int NBurned = 1;
00752   BurnedSet.AddKey(StartNId);
00753   BurningNIdV.Add(StartNId);
00754   while (! BurningNIdV.Empty()) {
00755     for (int node = 0; node < BurningNIdV.Len(); node++) {
00756       const int& BurningNId = BurningNIdV[node];
00757       const TUNGraph::TNodeI& Node = G.GetNI(BurningNId);
00758       // find unburned links
00759       AliveNIdV.Clr(false); // unburned links
00760       for (int e = 0; e < Node.GetOutDeg(); e++) {
00761         const int OutNId = Node.GetOutNId(e);
00762         if (! BurnedSet.IsKey(OutNId)) {
00763           AliveNIdV.Add(OutNId); }
00764       }
00765       // number of links to burn (geometric coin). Can also burn 0 links
00766       const int BurnNLinks = Rnd.GetGeoDev(1.0-BurnProb) - 1;
00767       if (! AliveNIdV.Empty() && BurnNLinks > 0) {
00768         AliveNIdV.Shuffle(Rnd);
00769         for (int i = 0; i < TMath::Mn(BurnNLinks, AliveNIdV.Len()); i++) {
00770           BurnedSet.AddKey(AliveNIdV[i]);
00771           NewBurnedNIdV.Add(AliveNIdV[i]);
00772           NBurned++;
00773         }
00774       }
00775     }
00776     BurningNIdV.Swap(NewBurnedNIdV);   // each node is burning just 1 time step
00777     // BurningNIdV.AddV(NewBurnedNIdV);   // all nodes are burning eternally
00778     NewBurnedNIdV.Clr(false);
00779   }
00780   IAssert(BurnedSet.Len() == NBurned);
00781   return NBurned;
00782 }
00783 
00784 TFfGGen::TStopReason TUndirFFire::AddNodes(const int& GraphNodes, const bool& FloodStop) {
00785   printf("\n***Undirected GEO ForestFire: graph(%d,%d) add %d nodes, burn prob %.3f\n", 
00786     Graph->GetNodes(), Graph->GetEdges(), GraphNodes, BurnProb);
00787   TExeTm ExeTm;
00788   int Burned1=0, Burned2=0, Burned3=0; // last 3 fire sizes
00789   TIntPrV NodesEdgesV;
00790   // create initial set of nodes
00791   if (Graph.Empty()) { Graph = PUNGraph::New(); }
00792   if (Graph->GetNodes() == 0) { Graph->AddNode(); }
00793   int NEdges = Graph->GetEdges();
00794   // forest fire
00795   for (int NNodes = Graph->GetNodes()+1; NNodes <= GraphNodes; NNodes++) {
00796     const int NewNId = Graph->AddNode(-1);
00797     IAssert(NewNId == Graph->GetNodes()-1); // node ids have to be 0...N
00798     const int StartNId = Rnd.GetUniDevInt(NewNId);
00799     const int NBurned = BurnGeoFire(StartNId);
00800     // add edges to burned nodes
00801     for (int e = 0; e < NBurned; e++) {
00802       Graph->AddEdge(NewNId, GetBurnedNId(e)); }
00803     NEdges += NBurned;
00804     Burned1=Burned2;  Burned2=Burned3;  Burned3=NBurned;
00805     if (NNodes % Kilo(1) == 0) {
00806       printf("(%d, %d)    burned: [%d,%d,%d]  [%s]\n", NNodes, NEdges, Burned1, Burned2, Burned3, ExeTm.GetStr()); 
00807       NodesEdgesV.Add(TIntPr(NNodes, NEdges));
00808     }
00809     if (FloodStop && NEdges>1000 && NEdges/double(NNodes)>100.0) { // average node degree is more than 50
00810       printf("!!! FLOOD. G(%6d, %6d)\n", NNodes, NEdges);  return TFfGGen::srFlood; }
00811   }
00812   printf("\n");
00813   IAssert(Graph->GetEdges() == NEdges);
00814   return TFfGGen::srOk;
00815 }