From 4cec0f0f4ef084aef20d41e86501a5ac00b70d74 Mon Sep 17 00:00:00 2001 From: ashwinsawant17 Date: Mon, 24 Nov 2025 07:38:49 -0500 Subject: [PATCH 1/6] Building segment_steps, and some basic get fns --- flatgfa/src/flatgfa.rs | 92 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) diff --git a/flatgfa/src/flatgfa.rs b/flatgfa/src/flatgfa.rs index b27e0872..beee01f6 100644 --- a/flatgfa/src/flatgfa.rs +++ b/flatgfa/src/flatgfa.rs @@ -550,3 +550,95 @@ pub type FixedGFAStore<'a> = GFAStore<'a, FixedFamily>; /// `FlatGFA`. It exposes an API for building up a GFA data structure, so it is /// useful for creating new ones from scratch. pub type HeapGFAStore = GFAStore<'static, HeapFamily>; + + +pub struct StepRef { + path: Id, + offset: u32 +} + +pub struct StepsBySegIndex { + steps: Vec, + segment_steps: Vec> +} + + +impl StepsBySegIndex { + + pub fn new(fgfa: &FlatGFA) -> Self{ + + + /// helper to extract the segment index from the stepref + fn segment_of_step(fgfa: &FlatGFA, step: &StepRef) -> usize { + let path = &fgfa.paths[step.path]; + let step_span = path.steps; + let step_slice = &fgfa.steps[step_span]; + step_slice[step.offset as usize].segment().index() + } + + // will be our `steps` vector that contains all steprefs + let mut all_steps = Vec::new(); + + // build the vector of all the steprefs, by iterating over each path, so we can construct using a path id and offset + for (path_id, path) in fgfa.paths.items() { + for (offset, _) in fgfa.get_path_steps(path).enumerate() { + all_steps.push(StepRef { + path: path_id, + offset: offset as u32, + }); + } + } + + // organize by the index of the segment in the segment pool + all_steps.sort_by_key(|a| { + segment_of_step(fgfa, a) + }); + + // the working segment's index + let mut seg_ind: usize = segment_of_step(fgfa, &all_steps[0]); + + // the working start of the span of StepRefs + let mut span_start: usize = 0; + + // the vector of spans of step refs that will act as our index + let mut segment_steps: Vec> = Vec::new(); + + + + for (i, step_ref) in all_steps.iter().enumerate() { + let curr_seg_ind: usize = segment_of_step(fgfa, step_ref); + + // if we have moved onto a new segment add to the segment_steps vec + if curr_seg_ind > seg_ind { + // create and push the new span of step refs for this segment + // it will be left inclusive right exclusive + let new_seg: Span = Span::new(Id::new(span_start), Id::new(i)); + segment_steps.push(new_seg); + + // update the working seg_ind and range_start variables to reflect the new segment/range + seg_ind = curr_seg_ind; + span_start = i; + } + + } + + + Self { + steps: all_steps, + segment_steps: segment_steps + } + } + + + /// Returns a slice of StepRefs that cross over the given segment + pub fn get_steps_slice(&self, segment: Id) -> &[StepRef] { + let range: Range = Range::from(&self.segment_steps[segment.index()]); + &(self.steps[range]) + } + + /// Returns the number of steps that cross over this segment + pub fn get_num_steps(&self, segment: Id) -> usize { + self.segment_steps[segment.index()].len() + } + +} \ No newline at end of file From 0eb441595228b90dd4ec375ee3933091af87f183 Mon Sep 17 00:00:00 2001 From: WinashAntsaw Date: Mon, 2 Feb 2026 16:01:58 -0500 Subject: [PATCH 2/6] moved stepsegindex to new index.rs, added to lib.rs, more thoroughly commented --- flatgfa/src/flatgfa.rs | 33 +++++--------- flatgfa/src/index.rs | 100 +++++++++++++++++++++++++++++++++++++++++ flatgfa/src/lib.rs | 1 + 3 files changed, 111 insertions(+), 23 deletions(-) create mode 100644 flatgfa/src/index.rs diff --git a/flatgfa/src/flatgfa.rs b/flatgfa/src/flatgfa.rs index beee01f6..4264073e 100644 --- a/flatgfa/src/flatgfa.rs +++ b/flatgfa/src/flatgfa.rs @@ -551,23 +551,18 @@ pub type FixedGFAStore<'a> = GFAStore<'a, FixedFamily>; /// useful for creating new ones from scratch. pub type HeapGFAStore = GFAStore<'static, HeapFamily>; - pub struct StepRef { path: Id, - offset: u32 + offset: u32, } pub struct StepsBySegIndex { steps: Vec, - segment_steps: Vec> + segment_steps: Vec>, } - impl StepsBySegIndex { - - pub fn new(fgfa: &FlatGFA) -> Self{ - - + pub fn new(fgfa: &FlatGFA) -> Self { /// helper to extract the segment index from the stepref fn segment_of_step(fgfa: &FlatGFA, step: &StepRef) -> usize { let path = &fgfa.paths[step.path]; @@ -575,7 +570,7 @@ impl StepsBySegIndex { let step_slice = &fgfa.steps[step_span]; step_slice[step.offset as usize].segment().index() } - + // will be our `steps` vector that contains all steprefs let mut all_steps = Vec::new(); @@ -588,23 +583,19 @@ impl StepsBySegIndex { }); } } - + // organize by the index of the segment in the segment pool - all_steps.sort_by_key(|a| { - segment_of_step(fgfa, a) - }); + all_steps.sort_by_key(|a| segment_of_step(fgfa, a)); // the working segment's index let mut seg_ind: usize = segment_of_step(fgfa, &all_steps[0]); // the working start of the span of StepRefs let mut span_start: usize = 0; - + // the vector of spans of step refs that will act as our index let mut segment_steps: Vec> = Vec::new(); - - for (i, step_ref) in all_steps.iter().enumerate() { let curr_seg_ind: usize = segment_of_step(fgfa, step_ref); @@ -617,19 +608,16 @@ impl StepsBySegIndex { // update the working seg_ind and range_start variables to reflect the new segment/range seg_ind = curr_seg_ind; - span_start = i; + span_start = i; } - } - Self { steps: all_steps, - segment_steps: segment_steps + segment_steps, } } - /// Returns a slice of StepRefs that cross over the given segment pub fn get_steps_slice(&self, segment: Id) -> &[StepRef] { let range: Range = Range::from(&self.segment_steps[segment.index()]); @@ -640,5 +628,4 @@ impl StepsBySegIndex { pub fn get_num_steps(&self, segment: Id) -> usize { self.segment_steps[segment.index()].len() } - -} \ No newline at end of file +} diff --git a/flatgfa/src/index.rs b/flatgfa/src/index.rs new file mode 100644 index 00000000..82ebc53b --- /dev/null +++ b/flatgfa/src/index.rs @@ -0,0 +1,100 @@ +use std::ops::Range; + +use crate::flatgfa::{FlatGFA, Path, Segment}; +use crate::pool::{Id, Span}; + +/// A reference to a specific Step within a Path +pub struct StepRef { + /// The ID of the Path + path: Id, + + /// The offset of the Step within that Path's Span of Steps + offset: u32, +} + +/// An in memory index that provides all Steps over a given Segment +pub struct StepsBySegIndex { + /// All Steps that are accessible by the index + /// StepRefs that are cross over the same segment are contiguous + steps: Vec, + + /// Spans into the steps vector. All Steps in the Span cross over one Segment + /// A Span's index maps to a Segment's index in a corresponding FlatGFA + segment_steps: Vec>, +} + +// TODO: make this return an Id instead of a usize +/// Extract the segment index from a StepRef in a given FlatGFA +fn segment_of_step(fgfa: &FlatGFA, step: &StepRef) -> usize { + let path = &fgfa.paths[step.path]; + let step_span = path.steps; + let step_slice = &fgfa.steps[step_span]; + step_slice[step.offset as usize].segment().index() +} + +impl StepsBySegIndex { + pub fn new(fgfa: &FlatGFA) -> Self { + // will be our `steps` vector that contains all steprefs + let mut all_steps = Vec::new(); + + // build the vector of all the steprefs, by iterating over each path, so we can construct using a path id and offset + for (path_id, path) in fgfa.paths.items() { + for (offset, _) in fgfa.get_path_steps(path).enumerate() { + all_steps.push(StepRef { + path: path_id, + offset: offset as u32, + }); + } + } + + // sort the steprefs by the index of the segment in the segment pool + all_steps.sort_by_key(|a| segment_of_step(fgfa, a)); + + // once sorted, steps that cross the same segment will be contiguous in all_steps + // this allows us to generate spans over steps that cross over the same segment, + // storing them in segment_steps + + // the working segment's index + let mut seg_ind: usize = segment_of_step(fgfa, &all_steps[0]); + + // the working start of the span of StepRefs + let mut span_start: usize = 0; + + // the vector of spans of step refs that will act as our index + let mut segment_steps: Vec> = Vec::new(); + + // enumerate over each step ref in all_steps to build our spans + for (i, step_ref) in all_steps.iter().enumerate() { + let curr_seg_ind: usize = segment_of_step(fgfa, step_ref); + + // if we have moved onto a new segment add to the segment_steps vec + if curr_seg_ind > seg_ind { + // create and push the new span of step refs for this segment + // it will be left inclusive right exclusive + let new_seg: Span = Span::new(Id::new(span_start), Id::new(i)); + segment_steps.push(new_seg); + + // update the working seg_ind and range_start variables to reflect the new segment/range + seg_ind = curr_seg_ind; + span_start = i; + } + } + + Self { + steps: all_steps, + segment_steps, + } + } + + /// Returns a slice of StepRefs that cross over the given segment + pub fn get_steps_slice(&self, segment: Id) -> &[StepRef] { + let range: Range = Range::from(&self.segment_steps[segment.index()]); + + &(self.steps[range]) + } + + /// Returns the number of steps that cross over this segment + pub fn get_num_steps(&self, segment: Id) -> usize { + self.segment_steps[segment.index()].len() + } +} diff --git a/flatgfa/src/lib.rs b/flatgfa/src/lib.rs index 7b5eeafe..c132ddd9 100644 --- a/flatgfa/src/lib.rs +++ b/flatgfa/src/lib.rs @@ -3,6 +3,7 @@ pub mod file; pub mod flatbed; pub mod flatgfa; pub mod gfaline; +pub mod index; pub mod memfile; pub mod namemap; pub mod ops; From 523f6906fe472409b456db69d905f4a4c94de152 Mon Sep 17 00:00:00 2001 From: WinashAntsaw Date: Mon, 2 Feb 2026 16:11:10 -0500 Subject: [PATCH 3/6] removing code in index.rs from flatgfa.rs --- flatgfa/src/flatgfa.rs | 79 ------------------------------------------ 1 file changed, 79 deletions(-) diff --git a/flatgfa/src/flatgfa.rs b/flatgfa/src/flatgfa.rs index 4264073e..b27e0872 100644 --- a/flatgfa/src/flatgfa.rs +++ b/flatgfa/src/flatgfa.rs @@ -550,82 +550,3 @@ pub type FixedGFAStore<'a> = GFAStore<'a, FixedFamily>; /// `FlatGFA`. It exposes an API for building up a GFA data structure, so it is /// useful for creating new ones from scratch. pub type HeapGFAStore = GFAStore<'static, HeapFamily>; - -pub struct StepRef { - path: Id, - offset: u32, -} - -pub struct StepsBySegIndex { - steps: Vec, - segment_steps: Vec>, -} - -impl StepsBySegIndex { - pub fn new(fgfa: &FlatGFA) -> Self { - /// helper to extract the segment index from the stepref - fn segment_of_step(fgfa: &FlatGFA, step: &StepRef) -> usize { - let path = &fgfa.paths[step.path]; - let step_span = path.steps; - let step_slice = &fgfa.steps[step_span]; - step_slice[step.offset as usize].segment().index() - } - - // will be our `steps` vector that contains all steprefs - let mut all_steps = Vec::new(); - - // build the vector of all the steprefs, by iterating over each path, so we can construct using a path id and offset - for (path_id, path) in fgfa.paths.items() { - for (offset, _) in fgfa.get_path_steps(path).enumerate() { - all_steps.push(StepRef { - path: path_id, - offset: offset as u32, - }); - } - } - - // organize by the index of the segment in the segment pool - all_steps.sort_by_key(|a| segment_of_step(fgfa, a)); - - // the working segment's index - let mut seg_ind: usize = segment_of_step(fgfa, &all_steps[0]); - - // the working start of the span of StepRefs - let mut span_start: usize = 0; - - // the vector of spans of step refs that will act as our index - let mut segment_steps: Vec> = Vec::new(); - - for (i, step_ref) in all_steps.iter().enumerate() { - let curr_seg_ind: usize = segment_of_step(fgfa, step_ref); - - // if we have moved onto a new segment add to the segment_steps vec - if curr_seg_ind > seg_ind { - // create and push the new span of step refs for this segment - // it will be left inclusive right exclusive - let new_seg: Span = Span::new(Id::new(span_start), Id::new(i)); - segment_steps.push(new_seg); - - // update the working seg_ind and range_start variables to reflect the new segment/range - seg_ind = curr_seg_ind; - span_start = i; - } - } - - Self { - steps: all_steps, - segment_steps, - } - } - - /// Returns a slice of StepRefs that cross over the given segment - pub fn get_steps_slice(&self, segment: Id) -> &[StepRef] { - let range: Range = Range::from(&self.segment_steps[segment.index()]); - &(self.steps[range]) - } - - /// Returns the number of steps that cross over this segment - pub fn get_num_steps(&self, segment: Id) -> usize { - self.segment_steps[segment.index()].len() - } -} From 26fb424782570170e5e7640f2189aeab2b0b9d68 Mon Sep 17 00:00:00 2001 From: ashwinsawant17 Date: Tue, 3 Feb 2026 02:15:02 -0500 Subject: [PATCH 4/6] added -n or --index flag to depth command, pushed bug fix that didnt get pushed earlier --- flatgfa/src/cli/cmds.rs | 10 +++++++--- flatgfa/src/cli/main.rs | 4 ++-- flatgfa/src/index.rs | 8 ++++++++ flatgfa/src/ops/depth.rs | 36 ++++++++++++++++++++++++++---------- 4 files changed, 43 insertions(+), 15 deletions(-) diff --git a/flatgfa/src/cli/cmds.rs b/flatgfa/src/cli/cmds.rs index 3bf22cd4..9b55fd9c 100644 --- a/flatgfa/src/cli/cmds.rs +++ b/flatgfa/src/cli/cmds.rs @@ -196,10 +196,14 @@ pub fn extract( /// compute node depth, the number of times paths cross a node #[derive(FromArgs, PartialEq, Debug)] #[argh(subcommand, name = "depth")] -pub struct Depth {} +pub struct Depth { + /// enable usage of StepsBySegIndex + #[argh(switch, short = 'n')] + index: bool, +} -pub fn depth(gfa: &flatgfa::FlatGFA) { - let (depths, uniq_paths) = ops::depth::depth(gfa); +pub fn depth(gfa: &flatgfa::FlatGFA, args: Depth) { + let (depths, uniq_paths) = ops::depth::depth(gfa, args.index); println!("#node.id\tdepth\tdepth.uniq"); for (id, seg) in gfa.segs.items() { diff --git a/flatgfa/src/cli/main.rs b/flatgfa/src/cli/main.rs index 03ee8357..f24226ba 100644 --- a/flatgfa/src/cli/main.rs +++ b/flatgfa/src/cli/main.rs @@ -127,8 +127,8 @@ fn main() -> Result<(), &'static str> { let store = cmds::extract(&gfa, sub_args)?; dump(&store.as_ref(), &args.output); } - Some(Command::Depth(_)) => { - cmds::depth(&gfa); + Some(Command::Depth(sub_args)) => { + cmds::depth(&gfa, sub_args); } Some(Command::Chop(sub_args)) => { let store = cmds::chop(&gfa, sub_args)?; diff --git a/flatgfa/src/index.rs b/flatgfa/src/index.rs index 82ebc53b..64069c4c 100644 --- a/flatgfa/src/index.rs +++ b/flatgfa/src/index.rs @@ -80,6 +80,11 @@ impl StepsBySegIndex { } } + // because we've only pushed a new span whenever we move onto a new segment + // we will need to push the final span after iterating + let new_seg = Span::new(Id::new(span_start), Id::new(all_steps.len())); + segment_steps.push(new_seg); + Self { steps: all_steps, segment_steps, @@ -95,6 +100,9 @@ impl StepsBySegIndex { /// Returns the number of steps that cross over this segment pub fn get_num_steps(&self, segment: Id) -> usize { + let vec_length = self.segment_steps.len(); + + println!("segment_steps length = {}", vec_length); self.segment_steps[segment.index()].len() } } diff --git a/flatgfa/src/ops/depth.rs b/flatgfa/src/ops/depth.rs index b46070d5..1243ee52 100644 --- a/flatgfa/src/ops/depth.rs +++ b/flatgfa/src/ops/depth.rs @@ -1,4 +1,5 @@ use crate::flatgfa; +use crate::index; use bit_vec::BitVec; /// Compute the *depth* of each segment in the variation graph. @@ -8,7 +9,7 @@ use bit_vec::BitVec; /// which only counts each path that tarverses a given segment once. /// /// Both outputs are depth values indexed by segment ID. -pub fn depth(gfa: &flatgfa::FlatGFA) -> (Vec, Vec) { +pub fn depth(gfa: &flatgfa::FlatGFA, use_index: bool) -> (Vec, Vec) { // Our output vectors: the ordinary and unique depths of each segment. let mut depths = vec![0; gfa.segs.len()]; let mut uniq_depths = vec![0; gfa.segs.len()]; @@ -18,15 +19,30 @@ pub fn depth(gfa: &flatgfa::FlatGFA) -> (Vec, Vec) { // subsequent traversals (for the purpose of counting unique depth). let mut seen = BitVec::from_elem(gfa.segs.len(), false); - for path in gfa.paths.all().iter() { - seen.clear(); // All segments are unseen. - for step in &gfa.steps[path.steps] { - let seg_id = step.segment().index(); - depths[seg_id] += 1; - if !seen[seg_id] { - // The first traversal of this path over this segment. - uniq_depths[seg_id] += 1; - seen.set(seg_id, true); + if use_index { + // Build the index + let step_seg_index = index::StepsBySegIndex::new(gfa); + + // iterate over each segment and populate the output vectors + for (id, _) in gfa.segs.items() { + // get the actual offset of the id of thesegment + let ind = id.index(); + + // use that offset to directly modify the depths vector + depths[ind] = step_seg_index.get_num_steps(id); + uniq_depths[ind] = step_seg_index.get_num_steps(id); + } + } else { + for path in gfa.paths.all().iter() { + seen.clear(); // All segments are unseen. + for step in &gfa.steps[path.steps] { + let seg_id = step.segment().index(); + depths[seg_id] += 1; + if !seen[seg_id] { + // The first traversal of this path over this segment. + uniq_depths[seg_id] += 1; + seen.set(seg_id, true); + } } } } From 5ed77de5339420a9212526076abf5f2ff5f5d281 Mon Sep 17 00:00:00 2001 From: ashwinsawant17 Date: Tue, 3 Feb 2026 02:33:36 -0500 Subject: [PATCH 5/6] removed debug print statement from index.rs --- flatgfa/src/index.rs | 3 --- 1 file changed, 3 deletions(-) diff --git a/flatgfa/src/index.rs b/flatgfa/src/index.rs index 64069c4c..6a2063ee 100644 --- a/flatgfa/src/index.rs +++ b/flatgfa/src/index.rs @@ -100,9 +100,6 @@ impl StepsBySegIndex { /// Returns the number of steps that cross over this segment pub fn get_num_steps(&self, segment: Id) -> usize { - let vec_length = self.segment_steps.len(); - - println!("segment_steps length = {}", vec_length); self.segment_steps[segment.index()].len() } } From 887d72a3ece0f9e1c1d33fc7cecd873068d75afb Mon Sep 17 00:00:00 2001 From: ashwinsawant17 Date: Sun, 22 Feb 2026 23:48:23 -0500 Subject: [PATCH 6/6] fixed bug with 0 depth segments, added unique depth functionality --- flatgfa/src/index.rs | 67 ++++++++++++++++++++++++++-------------- flatgfa/src/ops/depth.rs | 31 ++++++++++++++++--- 2 files changed, 70 insertions(+), 28 deletions(-) diff --git a/flatgfa/src/index.rs b/flatgfa/src/index.rs index 6a2063ee..90c42591 100644 --- a/flatgfa/src/index.rs +++ b/flatgfa/src/index.rs @@ -6,30 +6,30 @@ use crate::pool::{Id, Span}; /// A reference to a specific Step within a Path pub struct StepRef { /// The ID of the Path - path: Id, + pub path: Id, /// The offset of the Step within that Path's Span of Steps - offset: u32, + pub offset: u32, } /// An in memory index that provides all Steps over a given Segment pub struct StepsBySegIndex { - /// All Steps that are accessible by the index + /// All Steps that are accessible by the index. /// StepRefs that are cross over the same segment are contiguous - steps: Vec, + pub steps: Vec, - /// Spans into the steps vector. All Steps in the Span cross over one Segment + /// Spans into the steps vector. All Steps in the Span cross over one Segment. /// A Span's index maps to a Segment's index in a corresponding FlatGFA - segment_steps: Vec>, + pub segment_steps: Vec>, } // TODO: make this return an Id instead of a usize /// Extract the segment index from a StepRef in a given FlatGFA -fn segment_of_step(fgfa: &FlatGFA, step: &StepRef) -> usize { +fn segment_of_step(fgfa: &FlatGFA, step: &StepRef) -> Id { let path = &fgfa.paths[step.path]; let step_span = path.steps; let step_slice = &fgfa.steps[step_span]; - step_slice[step.offset as usize].segment().index() + step_slice[step.offset as usize].segment() } impl StepsBySegIndex { @@ -38,24 +38,33 @@ impl StepsBySegIndex { let mut all_steps = Vec::new(); // build the vector of all the steprefs, by iterating over each path, so we can construct using a path id and offset + // package the Id with the StepRef's into tuples so sorting by it is easier for (path_id, path) in fgfa.paths.items() { for (offset, _) in fgfa.get_path_steps(path).enumerate() { - all_steps.push(StepRef { + // the constructed StepRef based on the path and step offset + let step = StepRef { path: path_id, offset: offset as u32, - }); + }; + + // the segment that this step passes over in the form of a Id + let seg = segment_of_step(fgfa, &step); + + // push the tuple of the Id and StepRef into the steps vector + all_steps.push((seg, step)); } } // sort the steprefs by the index of the segment in the segment pool - all_steps.sort_by_key(|a| segment_of_step(fgfa, a)); + // by extracting the actual numeric index from the Id + all_steps.sort_by_key(|a| a.0.index()); // once sorted, steps that cross the same segment will be contiguous in all_steps // this allows us to generate spans over steps that cross over the same segment, // storing them in segment_steps // the working segment's index - let mut seg_ind: usize = segment_of_step(fgfa, &all_steps[0]); + let mut seg_ind: &Id = &all_steps[0].0; // the working start of the span of StepRefs let mut span_start: usize = 0; @@ -63,16 +72,23 @@ impl StepsBySegIndex { // the vector of spans of step refs that will act as our index let mut segment_steps: Vec> = Vec::new(); - // enumerate over each step ref in all_steps to build our spans - for (i, step_ref) in all_steps.iter().enumerate() { - let curr_seg_ind: usize = segment_of_step(fgfa, step_ref); + // fill the segment_steps with empty spans + // TODO: we definitely don't need to do another iteration to fill this with empty spans + // It's likely more efficient to push empty spans as needed + for _ in 0..fgfa.segs.len() { + segment_steps.push(Span::new_empty()); + } + // enumerate over each step ref in all_steps to build our spans + for (i, (curr_seg_ind, _)) in all_steps.iter().enumerate() { // if we have moved onto a new segment add to the segment_steps vec - if curr_seg_ind > seg_ind { - // create and push the new span of step refs for this segment + if curr_seg_ind.index() > seg_ind.index() { + // create and set the new span of step refs for this segment // it will be left inclusive right exclusive - let new_seg: Span = Span::new(Id::new(span_start), Id::new(i)); - segment_steps.push(new_seg); + let new_span: Span = Span::new(Id::new(span_start), Id::new(i)); + + // assign the span to the index in segment_steps that maps to the index of the segment in the FlatGFA segment pool + segment_steps[seg_ind.index()] = new_span; // update the working seg_ind and range_start variables to reflect the new segment/range seg_ind = curr_seg_ind; @@ -80,13 +96,16 @@ impl StepsBySegIndex { } } - // because we've only pushed a new span whenever we move onto a new segment - // we will need to push the final span after iterating - let new_seg = Span::new(Id::new(span_start), Id::new(all_steps.len())); - segment_steps.push(new_seg); + // because we've only set a new span whenever we move onto a new segment + // we will need to set the final span after iterating + let new_span: Span<_> = Span::new(Id::new(span_start), Id::new(all_steps.len())); + segment_steps[seg_ind.index()] = new_span; + + // unpackage all_steps into a vector with just spans + let steps: Vec = all_steps.into_iter().map(|x| x.1).collect(); Self { - steps: all_steps, + steps, segment_steps, } } diff --git a/flatgfa/src/ops/depth.rs b/flatgfa/src/ops/depth.rs index 1243ee52..1467e16e 100644 --- a/flatgfa/src/ops/depth.rs +++ b/flatgfa/src/ops/depth.rs @@ -23,14 +23,37 @@ pub fn depth(gfa: &flatgfa::FlatGFA, use_index: bool) -> (Vec, Vec // Build the index let step_seg_index = index::StepsBySegIndex::new(gfa); + // This bit vector keeps track of whether the current *path* has already been seen + // by the current working *segment*. A reverse approach to the non-indexed option + let mut path_seen: BitVec = BitVec::from_elem(gfa.paths.len(), false); + // iterate over each segment and populate the output vectors - for (id, _) in gfa.segs.items() { + for (seg_id, _) in gfa.segs.items() { + // clear the path_seen vector, since all paths should be unseen + path_seen.clear(); + // get the actual offset of the id of thesegment - let ind = id.index(); + let ind = seg_id.index(); + + // get the span of StepRefs for this segment + let span = step_seg_index.get_steps_slice(seg_id); // use that offset to directly modify the depths vector - depths[ind] = step_seg_index.get_num_steps(id); - uniq_depths[ind] = step_seg_index.get_num_steps(id); + depths[ind] = span.len(); + + // iterate over the span to populate the uniq_depths vector + for stepref in span.iter() { + // extract the path index + let path_id = stepref.path.index(); + + // if the path has not been seen, increment the unique depth + if !path_seen[path_id] { + uniq_depths[ind] += 1; + + // set the path to seen in the path_seen bitvec + path_seen.set(path_id, true); + } + } } } else { for path in gfa.paths.all().iter() {