+ PsiVector::Iterator** iters = this->getIterators();
+
+ usint runs = this->countRuns();
+ usint samples = 0, minimal_sum = 0, nonstrict_minimal_sum = 0, minimal_samples = 0;
+ usint max_samples = runs + (this->data_size - runs) / sample_rate;
+ sampled_values = new pair_type[max_samples];
+
+ std::list<Triple> matches;
+ for(usint j = 0; j < this->number_of_sequences; j++)
+ {
+ pair_type seq_range = this->getSequenceRange(j);
+ usint first_sample = samples; // First minimal sample of the current sequence.
+ usint minimal = seq_range.second + 1; // Last minimal position.
+ usint i, x = this->sa_samples->inverseSA(seq_range.first), next_x, prev_x = 0;
+
+ // Invariant: x == inverseSA(i)
+ for(i = seq_range.first; i <= seq_range.second; i++, x = next_x)
+ {
+ usint c = this->getCharacter(x);
+
+ // T[i,n] is lexicographically the first suffix beginning with c.
+ if(x == this->index_ranges[c].first)
+ {
+ if(!matches.empty())
+ {
+ for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end(); ++iter)
+ {
+ nonstrict_minimal_sum += i - (*iter).first;
+ }
+ matches.clear();
+ }
+
+ // Minimal sample: SA[x] = i, LCP[x] = 0
+ sampled_values[samples] = pair_type(x, 0); minimal = i;
+ samples++; minimal_samples++;
+ next_x = this->psiUnsafe(x, c, *(iters[c])) - this->number_of_sequences; prev_x = x;
+ continue;
+ }
+
+ // Process the previous left matches we are still following.
+ usint low = this->index_ranges[c].first + this->number_of_sequences;
+ pair_type potential_sample(this->data_size, this->data_size);
+ for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end();)
+ {
+ Triple match = *iter;
+ if(match.second < low) // These no longer match the current character.
+ {
+ if(potential_sample.first != this->data_size) { nonstrict_minimal_sum += potential_sample.second; }
+
+ // Potential minimal sample: SA[match.third] = match.first, LCP[match.third] = i - match.first
+ potential_sample = pair_type(match.third, i - match.first);
+ iter = matches.erase(iter);
+ }
+ else
+ {
+ (*iter).second = this->psiUnsafe(match.second - this->number_of_sequences, c, *(iters[c]));
+ ++iter;
+ }
+ }
+ if(potential_sample.first != this->data_size)
+ {
+ // Last potential sample is minimal.
+ sampled_values[samples] = potential_sample;
+ samples++; minimal_sum += potential_sample.second; minimal_samples++;
+ }
+
+ // If PLCP[i] is minimal, we add the left match to the list.
+ // We add triples of type (j, inverseSA(j) + n_of_s, x) as inverseSA is < 0 for end markers.
+ // x is the original minimal position.
+ usint next_y = this->psiUnsafe(x - 1, c, *(iters[c]));
+ next_x = this->psiUnsafeNext(c, *(iters[c]));
+ if(next_y != next_x - 1 || next_x < this->number_of_sequences)
+ {
+ matches.push_back(Triple(i, next_y, x)); minimal = i;
+ }
+ next_x -= this->number_of_sequences; prev_x = x;
+ }
+
+ // Are we still following something?
+ if(!matches.empty())
+ {
+ pair_type potential_sample(this->data_size, this->data_size);
+ for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end(); ++iter)
+ {
+ Triple match = *iter;
+ if(potential_sample.first != this->data_size) { nonstrict_minimal_sum += potential_sample.second; }
+
+ // Potential minimal sample: SA[match.third] = match.first, LCP[match.third] = i - match.first
+ potential_sample = pair_type(match.third, i - match.first);
+ }
+ matches.clear();
+ if(potential_sample.first != this->data_size)
+ {
+ // Last potential sample is minimal.
+ sampled_values[samples] = potential_sample;
+ samples++; minimal_sum += potential_sample.second; minimal_samples++;
+ }
+ }
+
+ // Add the non-minimal samples.
+ if(sample_rate <= this->data_size)
+ {
+ usint last_sample = samples - 1;
+ i = seq_range.first; x = this->sa_samples->inverseSA(seq_range.first);
+ for(usint current_sample = first_sample; current_sample <= last_sample; current_sample++)
+ {
+ // Find the next minimal sample and add nonminimal samples if needed.
+ pair_type first_nonminimal(i + sample_rate - 1, samples);
+ pair_type next_nonminimal = first_nonminimal;
+ while(x != sampled_values[current_sample].first)
+ {
+ if(i == next_nonminimal.first)
+ {
+ sampled_values[samples] = pair_type(x, 0); samples++;
+ next_nonminimal.first += sample_rate; next_nonminimal.second++;
+ }
+ i++; x = this->psi(x, iters) - this->number_of_sequences;
+ }
+
+ // Reduce the nonminimal samples to the current minimal sample.
+ for(next_nonminimal = first_nonminimal; next_nonminimal.second < samples; next_nonminimal.first += sample_rate, next_nonminimal.second++)
+ {
+ sampled_values[next_nonminimal.second].second = sampled_values[current_sample].second + i - next_nonminimal.first;
+ }
+ i++; x = this->psi(x, iters) - this->number_of_sequences;
+ }
+ }
+ }
+
+ std::sort(sampled_values, sampled_values + samples);
+
+ if(report)
+ {
+ std::cout << "Samples: " << samples << " (total) / " << minimal_samples << " (minimal)" << std::endl;
+ std::cout << "Upper bounds: " << max_samples << " (total) / " << runs << " (minimal)" << std::endl;
+ std::cout << "Sum of minimal samples: " << (minimal_sum + nonstrict_minimal_sum) << " (total) / " << minimal_sum << " (strict)" << std::endl;
+ std::cout << std::endl;
+ }
+
+ this->deleteIterators(iters);
+
+ return samples;