C++ Digraph::target方法代码示例

static void ILT_Helper( Digraph& dg, ILT_EdgeDataMap2& datas, ArcFilter& ef, NodeFilter& nf, EdgeArray& sortedAirEdges, IntMap& cmap )
{
    // 记录用风地点的始末搜索得到的count
    AcDbIntArray sAirs, tAirs;

    double st = Timer::rdtscSeconds();

    // 显示排序后，用风地点的顺序编号
    //DisplayAirEdgeNum(sortedAirEdges, datas);

    int len = sortedAirEdges.length();
    acutPrintf( _T( "\n用风地点个数:%d" ), len );
    for( int i = 0; i < len; i++ )
    {
        Digraph::Arc e = sortedAirEdges[i];
        Digraph::Node u = dg.source( e );
        Digraph::Node v = dg.target( e );

        sAirs.append( CountInFlowNodes( dg, ef, nf, u, cmap ) );
        tAirs.append( CountOutFlowNodes( dg, ef, nf, v, cmap ) );
        //acutPrintf(_T("\n第%d个用风地点【始节点】=%d  【末节点】=%d"), i, sAirs[i], tAirs[i]);
    }

    /*for(Digraph::NodeIt n(dg); n!=INVALID; ++n)
    {
    	if(nf[n] && imap[n]==0) imap[n] = 1;
    }*/

    double et = Timer::rdtscSeconds();
    acutPrintf( _T( "\n第1次正反dfs查找耗费时间:%.5f" ), et - st );

    // 显示分支的始末节点的cmap值(测试用)
    //DisplayEdgeCMapValue(dg, sortedAirEdges, datas, cmap, sAirs, tAirs);

    //assert(sAirs.length() == len);
    //assert(tAirs.length() == len);

    st = Timer::rdtscSeconds();

    for( int i = len - 1; i >= 0; i-- )
    {
        Digraph::Arc e = sortedAirEdges[i];
        Digraph::Node u = dg.source( e );
        Digraph::Node v = dg.target( e );

        MinusInFlowNodes( dg, ef, nf, u, cmap, sAirs[i] );
        PlusOutFlowNodes( dg, ef, nf, v, cmap, tAirs[i] );
    }

    et = Timer::rdtscSeconds();
    acutPrintf( _T( "\n第2次正反dfs查找耗费时间:%.5f" ), et - st );

    /*
     * 处理完之后，可能出现一些分支的始节点(或末节点)的cmap=0
     * 而末节点(或始节点)的cmap不等于0
     * 这些分支中可能有部分分支也可视为主要进回风之间的联络巷
     */
    // 显示分支的始末节点的cmap值(测试用)
    //DisplayEdgeCMapValue(dg, sortedAirEdges, datas, cmap, sAirs, tAirs);
}

static void PrintVN(Digraph& dg, EdgeIdMap& eid, NodeIdMap& nid, VNO_EdgeMap& ed, VNO_FanMap&fd, VNO_NodeMap& nd, double Q)
{
    cout<<"初始总风量："<<Q<<endl;
    cout<<setiosflags(ios::left)<<setw(12);
    cout<<"ID\t分支(始节点,末节点)\t风阻\t调节风阻\t风量\t始节点压力\t末节点压力"<<endl;
    for(Digraph::ArcIt e(dg); e!=INVALID; ++e)
    {
        int id = eid[e];
        Digraph::Node u = dg.source(e);
        Digraph::Node v = dg.target(e);
        double r = ed[e]->r;
        double q = ed[e]->q;
        double delta_r = ed[e]->delta_r;
        cout<<dg.id(e)<<"\t"<<"e("<<id<<")=("<<nid[u]<<","<<nid[v]<<")"<<"\t"
            <<"r(e"<<id<<")="<<r<<"\t"
            <<"d_r(e"<<id<<")="<<delta_r<<"\t"
            <<"q(e"<<id<<")="<<q<<"\t"
            <<"p(v"<<nid[u]<<")="<<nd[u]->p<<"\t"
            <<"p(v"<<nid[v]<<")="<<nd[v]->p<<endl;
        if(fd[e]!=0)
        {
            cout<<"\t-->分支e"<<id<<"上有风机\n"
                <<fd[e]->a0<<" "
                <<fd[e]->a1<<" "
                <<fd[e]->a2<<" "
                <<endl;
        }
        if(ed[e]->fq > 0)
        {
            cout<<"\t-->分支e"<<id<<"固定风量:"<<ed[e]->fq<<endl;
        }
    }
}

// 使用宽度优先搜索BFS计算最短距离
static void SortAirEdges( Digraph& dg, ArcFilter& ef, Digraph::Node s, Digraph::Node t, const EdgeArray& airEdges, EdgeArray& es )
{
    AFGraph afg( dg, ef );

    Bfs<AFGraph> aBfs( afg );
    aBfs.run( s );
    //if(!aDfs.reached(t)) return; // s->t不可达

    typedef ReverseDigraph<AFGraph> RDigraph;
    RDigraph rdg( afg );
    Bfs<RDigraph> bBfs( rdg );
    bBfs.run( t );
    //if(!bDfs.reached(s)) return;

    typedef std::vector<AirEdgeDist> AirEdgeArray;
    AirEdgeArray aes;

    int len = airEdges.length();
    for( int i = 0; i < len; i++ )
    {
        Digraph::Arc e = airEdges[i];
        Digraph::Node u = dg.source( e );
        Digraph::Node v = dg.target( e );

        AirEdgeDist aed;
        aed.e = e;
        if( aBfs.reached( u ) )
        {
            aed.sd = aBfs.dist( u );
        }
        else
        {
            acutPrintf( _T( "\n第%d个用风地点【始点】不可达" ), i );
        }

        if( bBfs.reached( v ) )
        {
            aed.td = bBfs.dist( v );
        }
        else
        {
            acutPrintf( _T( "\n第%d个用风地点【末点】不可达" ), i );
        }

        aes.push_back( aed );
    }

    // 降序排序
    std::sort( aes.begin(), aes.end() );

    for( AirEdgeArray::iterator itr = aes.begin(); itr != aes.end(); itr++ )
    {
        es.append( itr->e );
    }
}

static void FindStation_Helper( Digraph& dg, Digraph::Arc e, EdgeArray& es, bool searchStart = true, bool searchEnd = true )
{
    if( searchStart ) // 是否需要从分支的始节点搜索
    {
        FindStationFromNode( dg, dg.source( e ), es, true );
    }
    if( searchEnd ) // 是否需要从分支的末节点搜索
    {
        FindStationFromNode( dg, dg.target( e ), es, false );
    }
}

// 标记用风地点的始末节点和网络源汇节点
// 在搜索的时候，不搜索这些节点
static void BuildNodeFilter( Digraph& dg, EdgeArray& airEdges, Digraph::Node s, Digraph::Node t, NodeFilter& nf )
{
    int len = airEdges.length();
    acutPrintf( _T( "\n用风地点个数:%d" ), len );
    for( int i = 0; i < len; i++ )
    {
        Digraph::Arc e = airEdges[i];
        nf[dg.source( e )] = false;
        nf[dg.target( e )] = false;
    }

    nf[s] = false;
    nf[t] = false;

    // 原始网络的源点和汇点也要标记
    for( Digraph::OutArcIt e( dg, s ); e != INVALID; ++e )
    {
        nf[dg.target( e )] = false;
    }
    for( Digraph::InArcIt e( dg, t ); e != INVALID; ++e )
    {
        nf[dg.source( e )] = false;
    }
}

static void FindILTEdges( Digraph& dg, ILT_EdgeDataMap2& datas, IntMap& cmap, EdgeArray& edges )
{
    for( Digraph::ArcIt e( dg ); e != INVALID; ++e )
    {
        if( datas[e]->et == ET_VIRTUAL ) continue;
        Digraph::Node u = dg.source( e );
        Digraph::Node v = dg.target( e );

        // 排除直接相连的分支
        // 联络巷分支的始节点出度必须大于1，末节点的入度必须大于1
        if( countOutArcs( dg, u ) <= 1 || countInArcs( dg, v ) <= 1 ) continue;

        if( cmap[u]*cmap[v] < 0 ) edges.append( e );
    }
}

static void FilterILTEdge( Digraph& dg, IntMap& imap, ArcFilter& ef )
{
    // 查找联络巷
    for( Digraph::ArcIt e( dg ); e != INVALID; ++e )
    {
        Digraph::Node u = dg.source( e );
        Digraph::Node v = dg.target( e );

        // 排除直接相连的分支
        // 联络巷分支的始节点出度必须大于1，末节点的入度必须大于1
        if( countOutArcs( dg, u ) <= 1 || countInArcs( dg, v ) <= 1 ) continue;

        if( imap[u]*imap[v] < 0 ) ef[e] = false;
    }
}

static void DisplayEdgeCMapValue( Digraph& dg, EdgeArray& sortedAirEdges,
                                  ILT_EdgeDataMap2& datas, IntMap& cmap,
                                  AcDbIntArray& sAirs, AcDbIntArray& tAirs )
{
    AcDbObjectIdArray objIds;
    AcDbIntArray ivs, ivs2;
    for( Digraph::ArcIt e( dg ); e != INVALID; ++e )
    {
        objIds.append( datas[e]->objId );
        int pos = sortedAirEdges.find( e );
        if( pos != -1 )
        {
            ivs.append( sAirs[pos] );
            ivs2.append( tAirs[pos] );
        }
        else
        {
            ivs.append( cmap[dg.source( e )] );
            ivs2.append( cmap[dg.target( e )] );
        }
    }
    DisplayIntValue2( objIds, ivs, ivs2 );
}

static void Init_IMap( Digraph& dg, ArcFilter& ef, NodeFilter& nf, EdgeArray& airEdges, IntMap& imap )
{
    AFGraph afg( dg, ef );

    for( Digraph::NodeIt n( dg ); n != INVALID; ++n )
    {
        imap[n] = ( nf[n] ? 1 : 0 ); // 默认所有节点位于进风区
    }

    int len = airEdges.length();
    for( int i = 0; i < len; i++ )
    {
        Digraph::Arc e = airEdges[i];
        Digraph::Node u = dg.source( e );
        Digraph::Node v = dg.target( e );

        Dfs<AFGraph> bDfs( afg );
        bDfs.run( v );
        for( Digraph::NodeIt n( dg ); n != INVALID; ++n )
        {
            if( nf[n] && bDfs.reached( n ) ) imap[n] = -1; // 节点位于回风区
        }
    }
}

//指定一些特殊分支,要求搜索一条路径时尽可能多的通过这些特殊分支
//例如:搜索一条串联通风路径(通过多个用风地点的一条路径)
//采用二分匹配实现失败!!!
static bool MaxKeyEdgePass_Match(Digraph& dg, EdgeArray& airEdges, Digraph::Node s, Digraph::Node t, EdgeArray& mpp)
{
    EdgeArray p;
    if(!GraphUtils::DFS_Path(dg, s, t, p)) return false;

    typedef Digraph::NodeMap<Digraph::Node> DDMap;
    DDMap ddm(dg, INVALID);

    NodeArray left, right;
    left.push_back(s); right.push_back(t);
    for(size_t i=0;i<airEdges.size();i++)
    {
        Digraph::Arc e = airEdges[i];
        Digraph::Node u = dg.source(e);
        Digraph::Node v = dg.target(e);

        p.clear();
        bool ret = GraphUtils::DFS_Path(dg, v, t, p);
        if(!ret) continue;

        p.clear();
        ret = GraphUtils::DFS_Path(dg, s, u, p);
        if(!ret) continue;

        left.push_back(v); right.push_back(u);
    }
    //cout<<"left="<<left.size()<<" right="<<right.size()<<endl;

    ListGraph g;
    ListGraph::NodeMap<Digraph::Node> udm(g);
    typedef std::vector<ListGraph::Node> UNodeArray;
    UNodeArray left_nodes, right_nodes;
    typedef std::map<Digraph::Node, ListGraph::Node> DUMap;
    DUMap dum;

    //添加节点
    for(size_t i=0;i<left.size();i++)
    {
        Digraph::Node u = left[i];
        ListGraph::Node uu = g.addNode();
        udm[uu] = u; left_nodes.push_back(uu);
        //cout<<dg.id(u)<< " ";
    }
    //cout<<endl;
    for(size_t i=0;i<right.size();i++)
    {
        Digraph::Node u = right[i];
        ListGraph::Node uu = g.addNode();
        udm[uu] = u; right_nodes.push_back(uu);
        //cout<<dg.id(u)<<" ";
    }
    //cout<<endl;
    for(size_t i=0;i<right.size();i++)
    {
        Digraph::Node u = right[i];
        ListGraph::Node uu = g.addNode();
        udm[uu] = u; left_nodes.push_back(uu);
        //cout<<dg.id(u)<< " ";
    }
    //cout<<endl;
    for(size_t i=0;i<left.size();i++)
    {
        Digraph::Node u = left[i];
        ListGraph::Node uu = g.addNode();
        udm[uu] = u; right_nodes.push_back(uu);
        //cout<<dg.id(u)<< " ";
    }
    //cout<<endl;
    //添加分支
    for(size_t i=0;i<left_nodes.size();i++)
    {
        ListGraph::Node uv = left_nodes[i];

        for(size_t j=0;j<right_nodes.size();j++)
        {
            ListGraph::Node uu = right_nodes[j];

            Digraph::Node v = udm[uv];
            Digraph::Node u = udm[uu];
            if(u == v) continue;

            p.clear();
            if(GraphUtils::DFS_Path(dg, v, u, p))
            {
                //ListGraph::Node uv = left_nodes[i];
                //ListGraph::Node uu = right_nodes[j];
                ListGraph::Edge e = g.addEdge(uu,uv);
                //cout<<"二分图:"<<dg.id(v)<<"<-->"<<dg.id(u)<<endl;
            }
        }
    }
    //二分图最大匹配
    MaxMatching<ListGraph> mm(g);
    mm.run();

    //cout<<"分支数:"<<countEdges(g)<<" 匹配数:"<<mm.matchingSize()<<endl;
    for(ListGraph::EdgeIt e(g);e!=INVALID;++e)
//.........这里部分代码省略.........

//指定一些特殊分支,要求搜索一条路径时尽可能多的通过这些特殊分支
//例如:搜索一条串联通风路径(通过多个用风地点的一条路径)
//使用关键路径算法实现
static bool MaxKeyEdgePass_CPM(Digraph& dg, EdgeArray& airEdges, Digraph::Node s, Digraph::Node t, EdgeArray& mpp)
{
    EdgeArray p;
    if(!GraphUtils::DFS_Path(dg, s, t, p)) return false;

    NodeArray left, right;
    left.push_back(s); right.push_back(t);
    for(size_t i=0;i<airEdges.size();i++)
    {
        Digraph::Arc e = airEdges[i];
        Digraph::Node u = dg.source(e);
        Digraph::Node v = dg.target(e);

        p.clear();
        bool ret = GraphUtils::DFS_Path(dg, v, t, p);
        if(!ret) continue;

        p.clear();
        ret = GraphUtils::DFS_Path(dg, s, u, p);
        if(!ret) continue;

        left.push_back(v); right.push_back(u);
    }
    //cout<<"left="<<left.size()<<" right="<<right.size()<<endl;

    //构造新的图
    Digraph new_dg;
    //记录新图节点和原始图节点的关系
    Digraph::NodeMap<Digraph::Node> udm(new_dg);
    //记录原始图节点与新构造图节点的关系
    typedef std::map<Digraph::Node, Digraph::Node> DDMap;
    DDMap ddm;
    //左右节点数组
    NodeArray left_nodes, right_nodes;

    //添加节点
    for(size_t i=0;i<left.size();i++)
    {
        Digraph::Node u = left[i], uu = INVALID;
        DDMap::iterator itr = ddm.find(u);
        if(itr == ddm.end())
        {
            uu = new_dg.addNode();
        }
        else
        {
            uu = itr->second;
        }
        udm[uu] = u; ddm[u] = uu; 
        left_nodes.push_back(uu);
        //cout<<dg.id(u)<< " ";
    }
    //cout<<endl;
    for(size_t i=0;i<right.size();i++)
    {
        Digraph::Node u = right[i], uu = INVALID;
        DDMap::iterator itr = ddm.find(u);
        if(itr == ddm.end())
        {
            uu = new_dg.addNode();
        }
        else
        {
            uu = itr->second;
        }
        udm[uu] = u; ddm[u] = uu; 
        right_nodes.push_back(uu);
        //cout<<dg.id(u)<<" ";
    }
    //cout<<endl;
    //添加分支
    for(size_t i=0;i<left_nodes.size();i++)
    {
        Digraph::Node uv = left_nodes[i];

        for(size_t j=0;j<right_nodes.size();j++)
        {
            Digraph::Node uu = right_nodes[j];
            Digraph::Node v = udm[uv];
            Digraph::Node u = udm[uu];
            if(u == v || (v == s && u == t)) continue;

            p.clear();
            if(GraphUtils::DFS_Path(dg, v, u, p))
            {
                Digraph::Arc e = new_dg.addArc(uv,uu);
                //cout<<"新图:"<<dg.id(v)<<"<-->"<<dg.id(u)<<endl;
            }
        }
    }
    for(size_t i=0;i<left_nodes.size();i++)
    {
        Digraph::Node uv = left_nodes[i];
        Digraph::Node uu = right_nodes[i];
        Digraph::Node v = udm[uv];
        Digraph::Node u = udm[uu];
        if(u == v || u == t) continue;
//.........这里部分代码省略.........

int main(int argc, char *argv[]) 
{
  int nt;
  Digraph g;  // graph declaration
  string digraph_steiner_filename;
  DiNodeName vname(g);  // name of graph nodes
  Digraph::NodeMap<double> px(g),py(g);  // xy-coodinates for each node
  Digraph::NodeMap<int> vcolor(g);// color of nodes
  Digraph::ArcMap<int> ecolor(g); // color of edges
  ArcWeight lpvar(g);    // used to obtain the contents of the LP variables
  ArcWeight weight(g);   // edge weights
  Digraph::ArcMap<GRBVar> x(g); // binary variables for each arc
  vector <DiNode> V;
  int seed=0;
  srand48(1);

  // uncomment one of these lines to change default pdf reader, or insert new one
  //set_pdfreader("open");    // pdf reader for Mac OS X
  //set_pdfreader("xpdf");    // pdf reader for Linux
  //set_pdfreader("evince");  // pdf reader for Linux

    
  // double cutoff;   // used to prune non promissing branches (of the B&B tree)
  if (argc!=2) {cout<< endl << "Usage: "<< argv[0]<<"  <digraph_steiner_filename>"<< endl << endl;
    cout << "Examples:      " << argv[0] << " gr_berlin52.steiner" << endl;
    cout << "               " << argv[0] << " gr_usa48.steiner" << endl << endl;
    exit(0);}

  digraph_steiner_filename = argv[1];

  //int time_limit = 3600; // solution must be obtained within time_limit seconds
  GRBEnv env = GRBEnv();
  GRBModel model = GRBModel(env);
  model.getEnv().set(GRB_IntParam_LazyConstraints, 1);
  model.getEnv().set(GRB_IntParam_Seed, seed);
  model.set(GRB_StringAttr_ModelName, "Oriented Steiner Tree with GUROBI"); // prob. name
  model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE); // is a minimization problem

  ReadListDigraphSteiner(digraph_steiner_filename,g,vname,weight,px,py,1,nt,V); 
  Steiner_Instance T(g,vname,px,py,weight,nt,V);
  //for (DiNodeIt v(g);v!=INVALID;++v){ if(v==T.V[0])vcolor[v]=RED; else vcolor[v]=BLUE;}
  //for (int i=1;i<T.nt;i++) vcolor[T.V[i]] = MAGENTA;
  //for (ArcIt e(g); e != INVALID; ++e) ecolor[e] = BLUE;
  //ViewListDigraph(g,vname,px,py,vcolor,ecolor,"Triangulated graph");
  
  // Generate the binary variables and the objective function
  // Add one binary variable for each edge and set its cost in the objective function
  for (Digraph::ArcIt e(g); e != INVALID; ++e) {
    char name[100];
    sprintf(name,"X_%s_%s",vname[g.source(e)].c_str(),vname[g.target(e)].c_str());
    x[e] = model.addVar(0.0, 1.0, weight[e],GRB_BINARY,name); }
  model.update(); // run update to use model inserted variables
  try {
    //if (time_limit >= 0) model.getEnv().set(GRB_DoubleParam_TimeLimit,time_limit);
    //model.getEnv().set(GRB_DoubleParam_ImproveStartTime,10); //try better sol. aft. 10s
    // if (cutoff > 0) model.getEnv().set(GRB_DoubleParam_Cutoff, cutoff );

    ConnectivityCuts cb = ConnectivityCuts(T , x);
    model.setCallback(&cb);
    model.update();
    //model.write("model.lp"); system("cat model.lp");
    model.optimize();

    double soma=0.0;
    for (DiNodeIt v(g);v!=INVALID;++v) vcolor[v]=BLUE; // all nodes BLUE
    for (int i=0;i<T.nt;i++) vcolor[T.V[i]]=MAGENTA; // change terminals to MAGENTA
    vcolor[T.V[0]]=RED; // change root to RED
    for (Digraph::ArcIt e(g); e!=INVALID; ++e) {
      lpvar[e] = x[e].get(GRB_DoubleAttr_X);
      if (lpvar[e] > 1.0 - EPS) { soma += weight[e]; ecolor[e] = RED; }
      else ecolor[e] = NOCOLOR; }
    cout << "Steiner Tree Value = " << soma << endl;
    ViewListDigraph(g,vname,px,py,vcolor,ecolor,
	"Steiner Tree cost in graph with "+IntToString(T.nnodes)+
	" nodes and "+IntToString(T.nt)+" terminals: "+DoubleToString(soma));
  } catch (...) {cout << "Error during callback..." << endl; }
  return 0;
}

Node TarjanHD::hd(const Digraph& g,
                  const DoubleArcMap& w,
                  SubDigraph& subG,
                  NodeNodeMap& mapToOrgG,
                  NodeNodeMap& G2T,
                  const ArcList& sortedArcs,
                  int i)
{
  assert(i >= 0);
  
  int m = static_cast<int>(sortedArcs.size());
  assert(m == lemon::countArcs(subG));
  
  int r = m - i;
  if (r == 0 || r == 1)
  {
    // add to _T a subtree rooted at node r,
    // labeled with w(e_m) and having n children labeled with
    // the vertices of subG
    Node root = _T.addNode();
    _label[root] = w[sortedArcs.back()];
    _T2OrgG[root] = lemon::INVALID;
    
    for (SubNodeIt v(subG); v != lemon::INVALID; ++v)
    {
      Node vv = G2T[v];
      if (vv == lemon::INVALID)
      {
        vv = _T.addNode();
        _label[vv] = -1;
        if (mapToOrgG[v] != lemon::INVALID)
        {
          _orgG2T[mapToOrgG[v]] = vv;
          _T2OrgG[vv] = mapToOrgG[v];
        }
      }
      _T.addArc(vv, root);
    }
    
    return root;
  }
  else
  {
    int j = (i + m) % 2 == 0 ? (i + m) / 2 : (i + m) / 2 + 1;
    // remove arcs j+1 .. m
    ArcListIt arcEndIt, arcIt = sortedArcs.begin();
    for (int k = 1; k <= m; ++k)
    {
      if (k == j + 1)
      {
        arcEndIt = arcIt;
      }
      if (k > j)
      {
        subG.disable(*arcIt);
      }
      ++arcIt;
    }
    
    // compute SCCs
    IntNodeMap comp(g, -1);
    int numSCC = lemon::stronglyConnectedComponents(subG, comp);
    if (numSCC == 1)
    {
      ArcList newSortedArcs(sortedArcs.begin(), arcEndIt);
      
      // _subG is strongly connected
      return hd(g, w, subG, mapToOrgG, G2T, newSortedArcs, i);
    }
    else
    {
      // determine strongly connected components
      NodeHasher<Digraph> hasher(g);
      NodeSetVector components(numSCC, NodeSet(42, hasher));
      for (SubNodeIt v(subG); v != lemon::INVALID; ++v)
      {
        components[comp[v]].insert(v);
      }
      
      NodeVector roots(numSCC, lemon::INVALID);
      double w_i = i > 0 ? w[getArcByRank(sortedArcs, i)] : -std::numeric_limits<double>::max();
      for (int k = 0; k < numSCC; ++k)
      {
        const NodeSet& component = components[k];
        if (component.size() > 1)
        {
          // construct new sorted arc list for component: O(m) time
          ArcList newSortedArcs;
          for (ArcListIt arcIt = sortedArcs.begin(); arcIt != arcEndIt; ++arcIt)
          {
            Node u = g.source(*arcIt);
            Node v = g.target(*arcIt);

            bool u_in_comp = component.find(u) != component.end();
            bool v_in_comp = component.find(v) != component.end();
            
            if (u_in_comp && v_in_comp)
            {
              newSortedArcs.push_back(*arcIt);
            }
//.........这里部分代码省略.........

int TarjanHD::constructCondensedGraph(const Digraph& g,
                                      const DoubleArcMap& w,
                                      const NodeNodeMap& mapToOrgG,
                                      const NodeNodeMap& G2T,
                                      const ArcList& sortedArcs,
                                      const IntNodeMap& comp,
                                      const NodeSetVector& components,
                                      const NodeVector& roots,
                                      const int j,
                                      Digraph& c,
                                      DoubleArcMap& ww,
                                      NodeNodeMap& mapCToOrgG,
                                      NodeNodeMap& C2T,
                                      ArcList& newSortedArcs)
{
  // construct the condensed graph:
  // each strongly connected component is collapsed into a single node, and
  // from the resulting sets of multiple arcs retain only those with minimum weight
  lemon::DynArcLookUp<Digraph> lookUpArc(c);
  
  // first construct the nodes
  const size_t numSCC = components.size();
  NodeVector nodesOfC(numSCC, lemon::INVALID);
  for (size_t k = 0; k < numSCC; ++k)
  {
    Node v = nodesOfC[k] = c.addNode();
    
    if (components[k].size() == 1)
    {
      mapCToOrgG[v] = mapToOrgG[*components[k].begin()];
      C2T[v] = G2T[*components[k].begin()];
    }
    else
    {
      mapCToOrgG[v] = lemon::INVALID;
      C2T[v] = roots[k];
    }
  }
  
  // next construct the arcs: O(m) time
  for (ArcListIt arcIt = sortedArcs.begin(); arcIt != sortedArcs.end(); ++arcIt)
  {
    Arc a = *arcIt;
    Node u = g.source(a);
    Node v = g.target(a);
    
    Node uu = nodesOfC[comp[u]];
    Node vv = nodesOfC[comp[v]];
    if (uu != vv)
    {
      Arc aa = lookUpArc(uu, vv);
      if (aa == lemon::INVALID)
      {
        aa = c.addArc(uu, vv);
        newSortedArcs.push_back(aa);
        ww[aa] = w[a];
      }
#ifdef DEBUG
      else
      {
        assert(w[a] > ww[aa]);
      }
#endif
    }
  }
  
  return get_i(newSortedArcs, ww, w[getArcByRank(sortedArcs, j)]);
}

int main(int argc, char *argv[]) 
{
  int k,found;
  Digraph g;  // graph declaration
  string digraph_kpaths_filename, source_node_name, target_node_name;
  DiNodeName vname(g);  // name of graph nodes
  Digraph::NodeMap<double> px(g),py(g);  // xy-coodinates for each node
  Digraph::NodeMap<int> vcolor(g);// color of nodes
  Digraph::ArcMap<int> ecolor(g); // color of edges
  ArcWeight lpvar(g);    // used to obtain the contents of the LP variables
  ArcWeight weight(g);   // edge weights
  Digraph::ArcMap<GRBVar> x(g); // binary variables for each arc
  vector <DiNode> V;
  DiNode source,target;
  int seed=0;
  srand48(1);

  // uncomment one of these lines to change default pdf reader, or insert new one
  //set_pdfreader("open");    // pdf reader for Mac OS X
  //set_pdfreader("xpdf");    // pdf reader for Linux
  set_pdfreader("evince");  // pdf reader for Linux
  //set_pdfreader("open -a Skim.app");
  // double cutoff;   // used to prune non promissing branches (of the B&B tree)
  if (argc!=5) {cout<<endl<<"Usage: "<< argv[0]<<"  <digraph_kpaths_filename>  <source_node_name>  <target_node_name>  <k>"<< endl << endl;
    cout << "Example:      " << argv[0] << " digr_triang_sparse_100 12 50 5" << endl << endl;
    exit(0);}

  digraph_kpaths_filename = argv[1];
  source_node_name = argv[2];
  target_node_name = argv[3];
  k = atoi(argv[4]);

  GRBEnv env = GRBEnv();
  GRBModel model = GRBModel(env);
  model.getEnv().set(GRB_IntParam_Seed, seed);
  model.set(GRB_StringAttr_ModelName, "Oriented k-Paths with GUROBI"); // prob. name
  model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE); // is a minimization problem

  ReadListDigraph(digraph_kpaths_filename,g,vname,weight,px,py,0);
  found=0;
  for (DiNodeIt v(g);v!=INVALID;++v)
    if(vname[v]==source_node_name){source=v;found=1;break;}
  if (!found) {cout<<"Could not find source node "<<source_node_name<<endl;exit(0);}
  found=0;
  for (DiNodeIt v(g);v!=INVALID;++v)
    if(vname[v]==target_node_name){target=v;found=1;break;}
  if (!found) {cout<<"Could not find target node "<<target_node_name<<endl;exit(0);}
    
  kPaths_Instance T(g,vname,px,py,weight,source,target,k);
  
  //for (DiNodeIt v(g);v!=INVALID;++v){ if(v==T.V[0])vcolor[v]=RED; else vcolor[v]=BLUE;}
  //for (int i=1;i<T.nt;i++) vcolor[T.V[i]] = MAGENTA;
  //for (ArcIt e(g); e != INVALID; ++e) ecolor[e] = BLUE;
  //ViewListDigraph(g,vname,px,py,vcolor,ecolor,"Triangulated graph");
  
  // Generate the binary variables and the objective function
  // Add one binary variable for each edge and set its cost in the objective function
  for (Digraph::ArcIt e(g); e != INVALID; ++e) {
    char name[100];
    sprintf(name,"X_%s_%s",vname[g.source(e)].c_str(),vname[g.target(e)].c_str());
    x[e] = model.addVar(0.0, 1.0, weight[e],GRB_CONTINUOUS,name); }
  model.update(); // run update to use model inserted variables
  try {
    //if (time_limit >= 0) model.getEnv().set(GRB_DoubleParam_TimeLimit,time_limit);
    //model.getEnv().set(GRB_DoubleParam_ImproveStartTime,10); //try better sol. aft. 10s
    // if (cutoff > 0) model.getEnv().set(GRB_DoubleParam_Cutoff, cutoff );
    //model.write("model.lp"); system("cat model.lp");

    // Add degree constraint for each node (sum of solution edges incident to a node is 2)
    for (Digraph::NodeIt v(g); v!=INVALID; ++v) {
      GRBLinExpr exprin, exprout;
      for (Digraph::InArcIt e(g,v); e != INVALID; ++e) exprin += x[e];
      for (Digraph::OutArcIt e(g,v); e != INVALID; ++e) exprout += x[e];

      if (v==source)      model.addConstr(exprout - exprin == k );
      else if (v==target) model.addConstr(exprin - exprout == k );
      else                model.addConstr(exprin - exprout == 0 );
    }

    model.optimize();

    double soma=0.0;
    for (DiNodeIt v(g);v!=INVALID;++v) vcolor[v]=BLUE; // all nodes BLUE
    vcolor[source]=RED; // change the colors of the source node
    vcolor[target]=RED; // and the target node to RED
    for (Digraph::ArcIt e(g); e!=INVALID; ++e) {
      lpvar[e] = x[e].get(GRB_DoubleAttr_X);
      if (lpvar[e] > 1.0 - EPS) soma += weight[e];
    }
    cout << "kPaths Tree Value = " << soma << endl;

    //-----------------------------------------------------------------
    // By Lucas Prado Melo: coloring paths by bfs
    bool ok = true;
    for(int i=0; i < k && ok; i++) {
      queue<DiNode> q;
      Digraph::NodeMap<Arc> pre(g, INVALID);
      q.push(source);
      while (!q.empty()) {
        DiNode cur = q.front();
//.........这里部分代码省略.........