///|
/// A complete sequence analysis result suitable for a report or dashboard.
pub struct SequenceAnalysis {
  sequence : String
  length : Int
  windows : Array[SequenceWindow]
  composition : CompositionReport
  gc_skew : Int
  longest_run : Int
  motifs : Array[String]
} derive(Debug, Eq)

///|
pub struct AlignmentSummary {
  left_length : Int
  right_length : Int
  matches : Int
  mismatches : Int
  gaps : Int
  identity : Double
  edit_distance : Int
} derive(Debug, Eq)

///|
pub struct SequenceRun {
  base : Char
  start : Int
  end : Int
  length : Int
} derive(Debug, Eq)

///|
pub fn analyze_sequence(
  sequence : String,
  window_size : Int,
  step : Int,
) -> SequenceAnalysis {
  let normalized = normalize_sequence(sequence)
  let motifs = common_dna_motifs(normalized)
  {
    sequence: normalized,
    length: normalized.length(),
    windows: sequence_windows(normalized, window_size, step),
    composition: composition_report(normalized),
    gc_skew: sequence_gc_skew(normalized),
    longest_run: longest_homopolymer(normalized),
    motifs,
  }
}

///|
pub fn SequenceAnalysis::to_json_string(self : SequenceAnalysis) -> String {
  let windows = self.windows.map(fn(window) {
    "{\"start\":" +
    window.start.to_string() +
    ",\"end\":" +
    window.end.to_string() +
    ",\"sequence\":" +
    json_string(window.sequence) +
    ",\"gc_percent\":" +
    window.gc_percent.to_string() +
    "}"
  })
  "{\"sequence\":" +
  json_string(self.sequence) +
  ",\"length\":" +
  self.length.to_string() +
  ",\"windows\":" +
  array_json(windows) +
  ",\"composition\":" +
  self.composition.to_json_string() +
  ",\"gc_skew\":" +
  self.gc_skew.to_string() +
  ",\"longest_run\":" +
  self.longest_run.to_string() +
  ",\"motifs\":" +
  array_json(self.motifs.map(fn(item) { json_string(item) })) +
  "}"
}

///|
pub fn sequence_edit_distance(left : String, right : String) -> Int {
  let a = left.to_upper()
  let b = right.to_upper()
  let previous = Array::make(b.length() + 1, 0)
  for j = 0; j <= b.length(); j = j + 1 {
    previous[j] = j
  }
  for i = 1; i <= a.length(); i = i + 1 {
    let current = Array::make(b.length() + 1, 0)
    current[0] = i
    for j = 1; j <= b.length(); j = j + 1 {
      let substitution = previous[j - 1] +
        (if a.get_char(i - 1).unwrap() == b.get_char(j - 1).unwrap() {
          0
        } else {
          1
        })
      let insertion = current[j - 1] + 1
      let deletion = previous[j] + 1
      current[j] = Int::min(substitution, Int::min(insertion, deletion))
    }
    for j = 0; j <= b.length(); j = j + 1 {
      previous[j] = current[j]
    }
  }
  previous[b.length()]
}

///|
pub fn sequence_pairwise_identity(left : String, right : String) -> Double {
  let a = left.to_upper()
  let b = right.to_upper()
  let length = Int::max(a.length(), b.length())
  if length == 0 {
    return 100.0
  }
  let mut matches = 0
  let overlap = Int::min(a.length(), b.length())
  for i in 0.. AlignmentSummary {
  let a = left.to_upper()
  let b = right.to_upper()
  let overlap = Int::min(a.length(), b.length())
  let mut matches = 0
  for i in 0.. String {
  "{\"left_length\":" +
  self.left_length.to_string() +
  ",\"right_length\":" +
  self.right_length.to_string() +
  ",\"matches\":" +
  self.matches.to_string() +
  ",\"mismatches\":" +
  self.mismatches.to_string() +
  ",\"gaps\":" +
  self.gaps.to_string() +
  ",\"identity\":" +
  self.identity.to_string() +
  ",\"edit_distance\":" +
  self.edit_distance.to_string() +
  "}"
}

///|
pub fn sequence_runs(sequence : String) -> Array[SequenceRun] {
  let runs : Array[SequenceRun] = []
  let normalized = sequence.to_upper()
  if normalized == "" {
    return runs
  }
  let mut start = 0
  let mut current = normalized.get_char(0).unwrap()
  for i in 1.. SequenceRun? {
  let runs = sequence_runs(sequence)
  if runs.length() == 0 {
    None
  } else {
    Some(
      runs.fold(init=runs[0], fn(best, current) {
        if current.length > best.length {
          current
        } else {
          best
        }
      }),
    )
  }
}

///|
pub fn gc_skew_profile(sequence : String, window_size : Int) -> Array[Int] {
  let result : Array[Int] = []
  if window_size <= 0 {
    return result
  }
  for window in sequence_windows(sequence, window_size, window_size) {
    result.push(sequence_gc_skew(window.sequence))
  }
  result
}

///|
pub fn gc_profile(sequence : String, window_size : Int) -> Array[Double] {
  sequence_windows(sequence, window_size, window_size).map(fn(window) {
    window.gc_percent
  })
}

///|
pub fn base_profile(
  sequence : String,
  base : Char,
  window_size : Int,
) -> Array[Int] {
  sequence_windows(sequence, window_size, window_size).map(fn(window) {
    count_base(window.sequence, base)
  })
}

///|
pub fn sequence_entropy_profile(
  sequence : String,
  window_size : Int,
) -> Array[Double] {
  sequence_windows(sequence, window_size, window_size).map(fn(window) {
    composition_report(window.sequence).entropy
  })
}

///|
pub fn common_dna_motifs(sequence : String) -> Array[String] {
  let motifs = ["ATG", "TATA", "GGG", "CCC", "AATAAA", "GAATTC", "GGATCC"]
  motifs.filter(fn(motif) { has_motif(sequence, motif) })
}

///|
pub fn motif_density(sequence : String, motif : String) -> Double {
  if sequence.length() == 0 {
    0.0
  } else {
    motif_count(sequence, motif).to_double() *
    100.0 /
    sequence.length().to_double()
  }
}

///|
pub fn sequence_complexity(sequence : String) -> Double {
  composition_report(sequence).entropy
}

///|
pub fn sequence_is_low_complexity(
  sequence : String,
  threshold : Double,
) -> Bool {
  sequence_complexity(sequence) < threshold
}

///|
pub fn sequence_quality_flags(sequence : String) -> Array[String] {
  let flags : Array[String] = []
  let normalized = normalize_sequence(sequence)
  if normalized == "" {
    flags.push("empty")
  }
  if !is_valid_dna(normalized) {
    flags.push("invalid_base")
  }
  if normalized.length() > 0 && composition_report(normalized).ambiguous > 0 {
    flags.push("ambiguous_base")
  }
  if sequence_is_low_complexity(normalized, 1.0) {
    flags.push("low_complexity")
  }
  if longest_homopolymer(normalized) >= 8 {
    flags.push("long_homopolymer")
  }
  flags
}

///|
pub fn sequence_quality_ok(sequence : String) -> Bool {
  sequence_quality_flags(sequence).length() == 0
}

///|
pub fn translate_all_frames(
  sequence : String,
  table : CodonTable,
) -> Array[String] {
  [
    translate_frame(sequence, 1, table),
    translate_frame(sequence, 2, table),
    translate_frame(sequence, 3, table),
  ]
}

///|
pub fn reverse_complement_frames(
  sequence : String,
  table : CodonTable,
) -> Array[String] {
  translate_all_frames(reverse_complement(sequence), table)
}

///|
pub fn sequence_signature(sequence : String) -> String {
  let report = composition_report(sequence)
  report.length.to_string() +
  ":" +
  report.a.to_string() +
  ":" +
  report.c.to_string() +
  ":" +
  report.g.to_string() +
  ":" +
  report.t.to_string() +
  ":" +
  report.ambiguous.to_string()
}