void Pileups::compute_from_alignment(VG& graph, Alignment& alignment) {
    // if we start reversed
    if (alignment.has_path() && alignment.path().mapping(0).position().is_reverse()) {
        alignment = reverse_alignment(alignment,
                                      (function<int64_t(int64_t)>) ([&graph](int64_t id) {
                                          return graph.get_node(id)->sequence().size();
                                          }));
    }
    const Path& path = alignment.path();
    int64_t read_offset = 0;
    for (int i = 0; i < path.mapping_size(); ++i) {
        const Mapping& mapping = path.mapping(i);
        if (graph.has_node(mapping.position().node_id())) {
            const Node* node = graph.get_node(mapping.position().node_id());
            NodePileup* pileup = get_create(node->id());
            int64_t node_offset = mapping.position().offset();
            for (int j = 0; j < mapping.edit_size(); ++j) {
                const Edit& edit = mapping.edit(j);
                // process all pileups in edit.
                // update the offsets as we go
                compute_from_edit(*pileup, node_offset, read_offset, *node,
                                  alignment, mapping, edit);
            }
        }
    }
    assert(alignment.sequence().empty() ||
           alignment.path().mapping_size() == 0 ||
           read_offset == alignment.sequence().length());
}

int main_concat(int argc, char** argv) {

    if (argc == 2) {
        help_concat(argv);
        return 1;
    }

    int c;
    optind = 2; // force optind past command positional argument
    while (true) {
        static struct option long_options[] =
        {
            {"help", no_argument, 0, 'h'},
            {0, 0, 0, 0}
        };

        int option_index = 0;
        c = getopt_long (argc, argv, "h",
                long_options, &option_index);

        // Detect the end of the options.
        if (c == -1)
            break;

        switch (c)
        {
            case 'h':
            case '?':
                help_concat(argv);
                exit(1);
                break;

            default:
                abort ();
        }
    }

    list<VG*> graphs;

    while (optind < argc) {
        VG* graph;
        get_input_file(optind, argc, argv, [&](istream& in) {
            graph = new VG(in);
        });
        graphs.push_back(graph);
    }

    VG merged;
    for (list<VG*>::iterator g = graphs.begin(); g != graphs.end(); ++g) {
        merged.append(**g);
    }

    // output
    merged.serialize_to_ostream(std::cout);

    return 0;
}

string Sampler::alignment_seq(const Alignment& aln) {
    // get the graph corresponding to the alignment path
    Graph sub;
    for (int i = 0; i < aln.path().mapping_size(); ++ i) {
        auto& m = aln.path().mapping(i);
        if (m.has_position() && m.position().node_id()) {
            auto id = aln.path().mapping(i).position().node_id();
            xgidx->neighborhood(id, 2, sub);
        }
    }
    VG g; g.extend(sub);
    return g.path_string(aln.path());
}

void VGset::transform(std::function<void(VG*)> lambda) {
    for (auto& name : filenames) {
        // load
        VG* g = NULL;
        if (name == "-") {
            g = new VG(std::cin, show_progress);
        } else {
            ifstream in(name.c_str());
            g = new VG(in, show_progress);
            in.close();
        }
        g->name = name;
        // apply
        lambda(g);
        // write to the same file
        ofstream out(name.c_str());
        g->serialize_to_ostream(out);
        out.close();
        delete g;
    }
}

	VG handle_to_vg(const HandleGraph* xg) {
		// If xg is a null pointer, throw a runtime error
		if (xg == nullptr) {
			throw runtime_error("There is no xg to convert"); 
		} 
		// Initialize the VG graph
		VG vg;
		// Iterate through each handle in xg and create the same handle in vg
		xg->for_each_handle([&](const handle_t& here) {
			// Get the id of the xg handle
			id_t xg_id = xg->get_id(here);
			// Get the sequence of the xg handle
			string xg_seq = xg->get_sequence(here);
			// Create a handle in vg using the xg id and sequence
			vg.create_handle(xg_seq,xg_id);
		});
		// Iterate through each handle in xg 
		xg->for_each_handle([&](const handle_t& handle) {
			id_t id = xg->get_id(handle);
			bool rev = xg->get_is_reverse(handle);
			// Return a vg handle using the xg handle's id and orientation 
			handle_t current = vg.get_handle(id,rev);
			// Follow the right edges of the xg handle
			xg->follow_edges(handle, false, [&](const handle_t& r) {
				id_t id_r = xg->get_id(r);
				bool rev_r = xg->get_is_reverse(r);
				// Return a vg handle using the xg handle's id and orientation
				handle_t next = vg.get_handle(id_r, rev_r);
				// Create an edge in vg using the handles 
				vg.create_edge(current,next);
			});
			// Follow the left edges of the xg handle
			xg->follow_edges(handle, true, [&](const handle_t& l) {
				id_t id_l = xg->get_id(l);
				bool rev_l = xg->get_is_reverse(l);
				// Return a vg handle using the xg handle's id and orientation
				handle_t prev = vg.get_handle(id_l, rev_l);
				// Use the handles created from following the xg edges to create a vg edge
				vg.create_edge(prev,current); //error here
			});
		});
		return vg;
	}

// add all node traversals that valid walks from this one onto a stack
void stack_up_valid_walks(VG& graph, NodeTraversal walk_head, vector<NodeTraversal>& stack) {
    
    id_t head_id = walk_head.node->id();
    
    if (walk_head.backward) {
        // we are leaving from the start of the node
        
        // get all edges involving this node so we can filter them down to valid walks
        for (Edge* edge : graph.edges_of(walk_head.node)) {
            if (edge->from() == head_id && edge->from_start()) {
                // the edge is part of a valid walk
                Node* next_node = graph.get_node(edge->to());
                bool next_backward = edge->to_end();
                // add the next traversal in the walk to the stack
                stack.push_back(NodeTraversal(next_node, next_backward));
            }
            else if (edge->to() == head_id && !edge->to_end()) {
                // the edge is part of a valid walk in the opposite orientation
                Node* next_node = graph.get_node(edge->from());
                bool next_backward = edge->from_start();
                // add the next traversal in the walk to the stack
                stack.push_back(NodeTraversal(next_node, next_backward));
            }
        }
    }
    else {
        // we are leaving from the end of the node
        
        // get all edges involving this node so we can filter them down to valid walks
        for (Edge* edge : graph.edges_of(walk_head.node)) {
            if (edge->from() == head_id && !edge->from_start()) {
                // the edge is part of a valid walk
                Node* next_node = graph.get_node(edge->to());
                bool next_backward = edge->to_end();
                // add the next traversal in the walk to the stack
                stack.push_back(NodeTraversal(next_node, next_backward));
            }
            else if (edge->to() == head_id && edge->to_end()) {
                // the edge is part of a valid walk in the opposite orientation
                Node* next_node = graph.get_node(edge->from());
                bool next_backward = edge->from_start();
                // add the next traversal in the walk to the stack
                stack.push_back(NodeTraversal(next_node, next_backward));
            }
        }
    }
}

//.........这里部分代码省略.........
            {"nodes", no_argument, 0, 'n'},
            {"edges", no_argument, 0, 'e'},
            {"paths", no_argument, 0, 'o'},
            {"orphans", no_argument, 0, 'p'},
            {"gam", required_argument, 0, 'a'},
            {"xg", required_argument, 0, 'x'},
            {0, 0, 0, 0}
        };

        int option_index = 0;
        c = getopt_long (argc, argv, "hneopa:x:",
                long_options, &option_index);

        // Detect the end of the options.
        if (c == -1)
            break;

        switch (c)
        {

            case 'n':
                check_nodes = true;
                break;

            case 'e':
                check_edges = true;
                break;

            case 'o':
                check_orphans = true;
                break;

            case 'p':
                check_paths = true;
                break;

            case 'a':
                gam_path = optarg;
                break;

            case 'x':
                xg_path= optarg;
                break;

            case 'h':
            case '?':
                help_validate(argv);
                exit(1);
                break;

            default:
                abort ();
        }
    }

    if (!gam_path.empty() || !xg_path.empty()) {
        // GAM validation is its entirely own thing
        if (xg_path.empty()) {
            cerr << "error:[vg validate] xg index (-x) required with (-a)" << endl;
            return 1;
        } else if (gam_path.empty()) {
            cerr << "error:[vg validate] gam alignment (-a) required with (-x)" << endl;
            return 1;
        } else if (check_nodes || check_edges || check_orphans || check_paths) {
            cerr << "error:[vg validate] -n, -e -o, -p cannot be used with -a and -x" << endl;
            return 1;
        }
        ifstream in(xg_path.c_str());
        unique_ptr<xg::XG> xindex = stream::VPKG::load_one<xg::XG>(in);
        in.close();
        get_input_file(gam_path, [&](istream& in) {
                stream::for_each<Alignment>(in, [&](Alignment& aln) {
                        if (!alignment_is_valid(aln, xindex.get())) {
                            exit(1);
                        }
                    });
            });
        return 0;
    } else {

        VG* graph;
        get_input_file(optind, argc, argv, [&](istream& in) {
                graph = new VG(in);
            });

        // if we chose a specific subset, do just them
        if (check_nodes || check_edges || check_orphans || check_paths) {
            if (graph->is_valid(check_nodes, check_edges, check_orphans, check_paths)) {
                return 0;
            } else {
                return 1;
            }
            // otherwise do everything
        } else if (graph->is_valid()) {
            return 0;
        } else {
            return 1;
        }
    }
}

//.........这里部分代码省略.........
            output_buf.push_back(aln);
            stream::write_buffered(cout, output_buf, 100);
        };
        vindex->for_alignment_in_range(start_id, end_id, lambda);
        stream::write_buffered(cout, output_buf, 0);
    }

    if (!aln_on_id_range.empty()) {
        assert(!db_name.empty());
        vector<string> parts = split_delims(aln_on_id_range, ":");
        if (parts.size() == 1) {
            convert(parts.front(), start_id);
            end_id = start_id;
        } else {
            convert(parts.front(), start_id);
            convert(parts.back(), end_id);
        }
        vector<vg::id_t> ids;
        for (auto i = start_id; i <= end_id; ++i) {
            ids.push_back(i);
        }
        vector<Alignment> output_buf;
        auto lambda = [&output_buf](const Alignment& aln) {
            output_buf.push_back(aln);
            stream::write_buffered(cout, output_buf, 100);
        };
        vindex->for_alignment_to_nodes(ids, lambda);
        stream::write_buffered(cout, output_buf, 0);
    }

    if (!to_graph_file.empty()) {
        assert(vindex != nullptr);
        ifstream tgi(to_graph_file);
        VG graph(tgi);
        vector<vg::id_t> ids;
        graph.for_each_node([&](Node* n) { ids.push_back(n->id()); });
        vector<Alignment> output_buf;
        auto lambda = [&output_buf](const Alignment& aln) {
            output_buf.push_back(aln);
            stream::write_buffered(cout, output_buf, 100);
        };
        vindex->for_alignment_to_nodes(ids, lambda);
        stream::write_buffered(cout, output_buf, 0);
    }

    if (!xg_name.empty()) {
        if (!node_ids.empty() && path_name.empty() && !pairwise_distance) {
            // get the context of the node
            vector<Graph> graphs;
            set<vg::id_t> ids;
            for (auto node_id : node_ids) ids.insert(node_id);
            for (auto node_id : node_ids) {
                Graph g;
                xindex.neighborhood(node_id, context_size, g, !use_length);
                if (context_size == 0) {
                    for (auto& edge : xindex.edges_of(node_id)) {
                        // if both ends of the edge are in our targets, keep them
                        if (ids.count(edge.to()) && ids.count(edge.from())) {
                            *g.add_edge() = edge;
                        }
                    }
                }
                graphs.push_back(g);
            }
            VG result_graph;
            for (auto& graph : graphs) {

//.........这里部分代码省略.........
        // Read VG from stdin
        graph = unique_ptr<XG>(new XG());
        graph->from_stream(std::cin, validate_graph, print_graph, store_threads, is_sorted_dag);
    } else if (vg_in.size()) {
        // Read VG from a file
        ifstream in;
        in.open(vg_in.c_str());
        graph = unique_ptr<XG>(new XG());
        graph->from_stream(in, validate_graph, print_graph, store_threads, is_sorted_dag);
    }

    if (in_name.size()) {
        get_input_file(in_name, [&](istream& in) {
            // Load from an XG file or - (stdin)
            graph = stream::VPKG::load_one<XG>(in);
        });
    }

    // Prepare structure tree for serialization
    unique_ptr<sdsl::structure_tree_node> structure;
    
    if (!report_name.empty()) {
        // We need to make a report, so we need the structure. Make a real tree
        // node. The unique_ptr handles deleting.
        structure = unique_ptr<sdsl::structure_tree_node>(new sdsl::structure_tree_node("name", "type"));
    }

    if(!vg_out.empty()) {
        if (graph.get() == nullptr) {
             cerr << "error [vg xg] no xg graph exists to convert; Try: vg xg -i graph.xg -X graph.vg" << endl;
             return 1;
        }
        
        VG converted;
        // Convert the xg graph to vg format
        convert_handle_graph(graph.get(), &converted);
        
        // TODO: The converter doesn't copy circular paths yet.
        // When it does, we can remove all this path copying code.

        // Make a raw Proto Graph to hold Path objects
        Graph path_graph;

        // Since paths are not copied, copy the paths.
        for (size_t rank = 1; rank <= graph->max_path_rank(); rank++) {
            // Extract each path into the path graph
            *path_graph.add_path() = graph->path(graph->path_name(rank));
        }

        // Merge in all the paths
        converted.extend(path_graph);
        
        if (vg_out == "-") {
            converted.serialize_to_ostream(std::cout);
        } else {
            converted.serialize_to_file(vg_out);
        }
    }

    if (!out_name.empty()) {
        // Open a destination file if it is a file we want to write to
        ofstream out_file;
        if (out_name != "-") {
            out_file.open(out_name);
        }
        // Work out where to save to

PathIndex::PathIndex(const list<Mapping>& mappings, VG& vg) {
    // Trace the given path in the given VG graph, collecting sequence
    
    // We're going to build the sequence string
    std::stringstream seq_stream;
    
    // What base are we at in the path?
    size_t path_base = 0;
    
    // What was the last rank? Ranks must always go up.
    int64_t last_rank = -1;
    
    for (auto& mapping : mappings) {
    
        if (!by_id.count(mapping.position().node_id())) {
            // This is the first time we have visited this node in the path.
            
            // Add in a mapping.
            by_id[mapping.position().node_id()] = 
                std::make_pair(path_base, mapping.position().is_reverse());
#ifdef debug
            #pragma omp critical (cerr)
            std::cerr << "Node " << mapping.position().node_id() << " rank " << mapping.rank()
                << " starts at base " << path_base << " with "
                << vg.get_node(mapping.position().node_id())->sequence() << std::endl;
#endif
            
            // Make sure ranks are monotonically increasing along the path, or
            // unset.
            assert(mapping.rank() > last_rank || (mapping.rank() == 0 && last_rank == 0));
            last_rank = mapping.rank();
        }
        
        // Say that this node appears here along the reference in this
        // orientation.
        by_start[path_base] = NodeSide(mapping.position().node_id(), mapping.position().is_reverse());
    
        // Remember that occurrence by node ID.
        node_occurrences[mapping.position().node_id()].push_back(by_start.find(path_base));
    
        // Say this Mapping happens at this base along the path
        mapping_positions[&mapping] = path_base;
    
        // Find the node's sequence
        std::string node_sequence = vg.get_node(mapping.position().node_id())->sequence();
    
        while(path_base == 0 && node_sequence.size() > 0 &&
            (node_sequence[0] != 'A' && node_sequence[0] != 'T' && node_sequence[0] != 'C' &&
            node_sequence[0] != 'G' && node_sequence[0] != 'N')) {
            
            // If the path leads with invalid characters (like "X"), throw them
            // out when computing path positions.
            
            // TODO: this is a hack to deal with the debruijn-brca1-k63 graph,
            // which leads with an X.
            #pragma omp critical (cerr)
            std::cerr << "Warning: dropping invalid leading character "
                << node_sequence[0] << " from node " << mapping.position().node_id()
                << std::endl;
                
            node_sequence.erase(node_sequence.begin());
        }
        
        if (mapping.position().is_reverse()) {
            // Put the reverse sequence in the path
            seq_stream << reverse_complement(node_sequence);
        } else {
            // Put the forward sequence in the path
            seq_stream << node_sequence;
        }
        
        // Whether we found the right place for this node in the reference or
        // not, we still need to advance along the reference path. We assume the
        // whole node (except any leading bogus characters) is included in the
        // path (since it sort of has to be, syntactically, unless it's the
        // first or last node).
        path_base += node_sequence.size();
        
        // TODO: handle leading bogus characters in calls on the first node.
    }
    
    // Record the length of the last mapping's node, since there's no next mapping to work it out from
    last_node_length = mappings.empty() ?
        0 : 
        vg.get_node(mappings.back().position().node_id())->sequence().size();
    
    // Create the actual reference sequence we will use
    sequence = seq_stream.str();
    
#ifdef debug
    // Announce progress.
    #pragma omp critical (cerr)
    std::cerr << "Traced " << path_base << " bp path." << std::endl;
    
    if (sequence.size() < 100) {
        #pragma omp critical (cerr)
        std::cerr << "Sequence: " << sequence << std::endl;
    }
#endif

//.........这里部分代码省略.........

int main_validate(int argc, char** argv) {

    if (argc <= 2) {
        help_validate(argv);
        return 1;
    }

    bool check_nodes = false;
    bool check_edges = false;
    bool check_orphans = false;
    bool check_paths = false;

    int c;
    optind = 2; // force optind past command positional argument
    while (true) {
        static struct option long_options[] =
        {
            {"help", no_argument, 0, 'h'},
            {"nodes", no_argument, 0, 'n'},
            {"edges", no_argument, 0, 'e'},
            {"paths", no_argument, 0, 'o'},
            {"orphans", no_argument, 0, 'p'},
            {0, 0, 0, 0}
        };

        int option_index = 0;
        c = getopt_long (argc, argv, "hneop",
                long_options, &option_index);

        // Detect the end of the options.
        if (c == -1)
            break;

        switch (c)
        {

            case 'n':
                check_nodes = true;
                break;

            case 'e':
                check_edges = true;
                break;

            case 'o':
                check_orphans = true;
                break;

            case 'p':
                check_paths = true;
                break;

            case 'h':
            case '?':
                help_validate(argv);
                exit(1);
                break;

            default:
                abort ();
        }
    }

    VG* graph;
    get_input_file(optind, argc, argv, [&](istream& in) {
        graph = new VG(in);
    });

    // if we chose a specific subset, do just them
    if (check_nodes || check_edges || check_orphans || check_paths) {
        if (graph->is_valid(check_nodes, check_edges, check_orphans, check_paths)) {
            return 0;
        } else {
            return 1;
        }
        // otherwise do everything
    } else if (graph->is_valid()) {
        return 0;
    } else {
        return 1;
    }
}

//.........这里部分代码省略.........
            break;

        case 'y':
            destroy_node_id = parse<int>(optarg);
            break;

        case 'a':
            cactus = true;
            break;

        case 'v':
            vcf_filename = optarg;
            break;
            
        case 'G':
            loci_filename = optarg;
            break;

        case 'M':
            max_degree = parse<int>(optarg);
            break;

        case 'h':
        case '?':
            help_mod(argv);
            exit(1);
            break;

        default:
            abort ();
        }
    }

    VG* graph;
    get_input_file(optind, argc, argv, [&](istream& in) {
        graph = new VG(in);
    });
    
    if (retain_complement) {
        // Compute the actual paths to retain
        set<string> complement;
        graph->paths.for_each_name([&](const string& name) {
            if (!paths_to_retain.count(name)) {
                // Complement the set the user specified by putting in all the
                // paths they didn't mention.
                complement.insert(name);
            }
        });
        
        // Retain the complement of what we were asking for.
        paths_to_retain = complement;
    }

    if (!vcf_filename.empty()) {
        // We need to throw out the parts of the graph that are on alt paths,
        // but not on alt paths for alts used by the first sample in the VCF.

        // This is called with the entire path name string to detect alt
        // paths.
        const function<bool(const string&)>& is_alt = Paths::is_alt;

        // This holds the VCF file we read the variants from. It needs to be the
        // same one used to construct the graph.
        vcflib::VariantCallFile variant_file;
        variant_file.open(vcf_filename);
        if (!variant_file.is_open()) {

CactusSiteFinder::CactusSiteFinder(VG& graph, const string& hint_path_name): graph(graph), hint_path_name(hint_path_name) {
    // Make sure the graph is sorted.
    // cactus needs the nodes to be sorted in order to find a source and sink.
    graph.sort();
}

void renderer::draw_objects(VG& v) {
    VG::iterator it;
    for (it = v.begin(); it != v.end(); ++it) {
        draw(*it);
    }
}

//.........这里部分代码省略.........
    if (gam_in_file_name == "-" && graph_file_name == "-") {
        cerr << "[vg augment] error: graph and gam can't both be from stdin." << endl;
        return 1;
    }
    if (gam_in_file_name == "-" && !gam_out_file_name.empty()) {
        cerr << "[vg augment] error: cannot stream input gam when using -A option (as it requires 2 passes)" << endl;
        return 1;
    }

    if (augmentation_mode != "pileup" && augmentation_mode != "direct") {
        cerr << "[vg augment] error: pileup and direct are currently the only supported augmentation modes (-a)" << endl;
        return 1;
    }

    if (augmentation_mode != "direct" and !gam_out_file_name.empty()) {
        cerr << "[vg augment] error: GAM output only works with \"direct\" augmentation mode" << endl;
        return 1;
    }

    if (augmentation_mode != "pileup" and (!support_file_name.empty() || !pileup_file_name.empty())) {
        cerr << "[vg augment] error: Pileup (-P) and Support (-S) output only work with  \"pileup\" augmentation mode" << endl;
        return 1;
    }

    if (label_paths && (!gam_out_file_name.empty() || !translation_file_name.empty())) {
        cerr << "[vg augment] error: Translation (-Z) and GAM (-A) output do not work with \"label-only\" (-B) mode" << endl;
        return 1;
    }
    
    // read the graph
    if (show_progress) {
        cerr << "Reading input graph" << endl;
    }
    VG* graph;
    get_input_file(graph_file_name, [&](istream& in) {
        graph = new VG(in);
    });
    
    
    Pileups* pileups = nullptr;
    
    if (!pileup_file_name.empty() || augmentation_mode == "pileup") {
        // We will need the computed pileups
        
        // compute the pileups from the graph and gam
        pileups = compute_pileups(graph, gam_in_file_name, thread_count, min_quality, max_mismatches,
                                  window_size, max_depth, use_mapq, show_progress);
    }
        
    if (!pileup_file_name.empty()) {
        // We want to write out pileups.
        if (show_progress) {
            cerr << "Writing pileups" << endl;
        }
        ofstream pileup_file(pileup_file_name);
        if (!pileup_file) {
            cerr << "[vg augment] error: unable to open output pileup file: " << pileup_file_name << endl;
            exit(1);
        }
        pileups->write(pileup_file);
    }

    if (augmentation_mode == "direct" && !gam_in_file_name.empty()) {
        // Augment with the reads
        
        if (!support_file_name.empty()) {

TEST(inequality, Point)
{
    CHECK(Point(1, 2) != Point(3, 4));
}

TEST(equality, Point)
{
    CHECK_EQUAL(Point(1, 2), Point(1, 2));
}