Compare commits

...
Author SHA1 Message Date
Eric Coissac a8841c203d Release 4.4.46 2026-07-02 08:50:00 +02:00
Eric Coissac f1a424343c Update ignore patterns, configure Serena, and refactor FASTSEQ parser
Adds `.gitignore` rules to exclude benchmark outputs and local build artifacts. Initializes Serena AI configuration for the Go project with UTF-8 encoding and gitignore-based file exclusion. Refactors the FASTSEQ header parser to fix regex syntax, optimize alternation ordering, improve control flow readability, and introduce explicit float-to-int coercion that preserves fractional values.
2026-07-02 08:49:26 +02:00
coissacandGitHub 4dae7f708d Merge pull request #119 from metabarcoding/push-zrmnkuwsztxx
Release 4.4.45
2026-06-02 15:02:42 +02:00
Eric Coissac 578445f38a Release 4.4.45 2026-06-02 15:01:32 +02:00
Eric Coissac 2edb33ad08 fix: correct duplicated typos in GVal function names
Corrects duplicated typos in the registered GVal function names, changing "which_maxwhichmax" and "which_minwhichmin" to "which_max" and "which_min". The underlying obiutils.WhichMax/WhichMin logic and its int-to-float64 index conversion remain unchanged.
2026-06-02 15:00:48 +02:00
coissacandGitHub b108d4978e Merge pull request #118 from metabarcoding/push-twztopkvxyop
Release 4.4.44
2026-06-02 14:42:00 +02:00
Eric Coissac e9210e28a3 Release 4.4.44 2026-06-02 14:40:09 +02:00
Eric Coissac 13a93fce11 feat: add which_max and which_min to retrieve extreme element indices
Implement reflection-based WhichMax and WhichMin to dynamically find the index or key of the maximum/minimum element in slices, arrays, or maps. Functions validate orderability, handle empty collections, and dispatch via reflect.Kind. Expose as which_max and which_min GVal functions, with float64 type assertions for compatibility and preserved error handling.
2026-06-02 14:39:31 +02:00
Eric Coissac 14064c919e fix(obiutils): correctly unwrap interface values in min/max
Introduces an `unwrapInterface` reflection helper to dereference `interface{}`-wrapped values before type validation. Updates slice and map iteration loops in min/max functions to apply this helper, ensuring `isOrderedKind` accurately identifies underlying concrete types instead of incorrectly rejecting `reflect.Interface` elements.
2026-06-02 14:34:00 +02:00
coissacandGitHub 1dfd68aa6d Merge pull request #117 from metabarcoding/push-wvlmzvomslzv
Release 4.4.43
2026-06-01 14:14:40 +02:00
Eric Coissac 930fe5f1ba Release 4.4.43 2026-06-01 13:22:58 +02:00
Eric Coissac dcdaf9e372 feat: support map and slice types in OBI attributes
Extends OBI header parsing to recognize and deserialize JSON-like arrays and objects. Introduces safe conversion utilities in `obiutils` to cast generic interface values into typed maps, and exposes them via new `BioSequence` methods. Header values are now marshaled, quote-normalized, and formatted for map and slice types.
2026-06-01 13:21:11 +02:00
Eric Coissac af7ae3d60c Correct Shannon entropy bias for canonical k-mers
Multiple raw k-mers collapsing into identical circular canonical forms introduce bias into complexity estimates. This change pre-computes `log(class_size)` tables and per-word-size maximum entropy bounds. The `KmerEntropy` function and `KmerEntropyFilter` are updated to apply the corrected formula `(log(N) + Σf·log(s) - Σf·log(f))/N / emax`, ensuring accurate sequence complexity estimation.
2026-05-17 14:54:57 +08:00
Eric Coissac cecf90fa40 feat: add min/max filtering and saturating subtraction utilities
Introduce generic and reflection-based utilities for filtering slices and maps by minimum/maximum thresholds, along with saturating subtraction. The `obiutils` package provides type-safe generic implementations alongside dynamic reflection dispatchers to handle arbitrary ordered and numeric types. These are exposed as GVAL expression functions in `obiseq`, extending the language's built-in filtering and numeric capabilities.
2026-05-14 20:58:24 +08:00
Eric Coissac a186bd1c92 fix: validate non-empty sequence IDs in FASTA and FASTQ writers
Adds a pre-processing guard that checks for empty sequence identifiers before formatting. This prevents malformed FASTA output and stops downstream processing of invalid FASTQ data by terminating early. The check is placed before existing sequence-length validations to enforce non-empty IDs during batch processing.
2026-05-05 18:07:58 +02:00
coissacandGitHub 46d60c1a44 Merge pull request #115 from metabarcoding/push-lkzqoskvyqtr
[4.4.2] Enhanced taxonomy handling, input robustness & PCR tag validation
2026-04-30 16:59:49 +02:00
14 changed files with 790 additions and 91 deletions
+2
View File
@@ -21,6 +21,8 @@ xx
.rhistory .rhistory
/.vscode /.vscode
/benchmarck
/build /build
/bugs /bugs
autodoc autodoc
+2
View File
@@ -0,0 +1,2 @@
/cache
/project.local.yml
+133
View File
@@ -0,0 +1,133 @@
# the name by which the project can be referenced within Serena
project_name: "obitools4"
# list of languages for which language servers are started; choose from:
# al angular ansible bash clojure
# cpp cpp_ccls crystal csharp csharp_omnisharp
# dart elixir elm erlang fortran
# fsharp go groovy haskell haxe
# hlsl html java json julia
# kotlin lean4 lua luau markdown
# matlab msl nix ocaml pascal
# perl php php_phpactor powershell python
# python_jedi python_ty r rego ruby
# ruby_solargraph rust scala scss solidity
# svelte swift systemverilog terraform toml
# typescript typescript_vts vue yaml zig
# (This list may be outdated. For the current list, see values of Language enum here:
# https://github.com/oraios/serena/blob/main/src/solidlsp/ls_config.py
# For some languages, there are alternative language servers, e.g. csharp_omnisharp, ruby_solargraph.)
# Note:
# - For C, use cpp
# - For JavaScript, use typescript
# - For Angular projects, use angular (subsumes typescript+html; requires `npm install` in the project root)
# - For Svelte projects, use svelte (subsumes typescript/javascript for .svelte projects; requires npm)
# - For SCSS / Sass / plain CSS, use scss (some-sass-language-server handles all three)
# - For Free Pascal/Lazarus, use pascal
# Special requirements:
# Some languages require additional setup/installations.
# See here for details: https://oraios.github.io/serena/01-about/020_programming-languages.html#language-servers
# When using multiple languages, the first language server that supports a given file will be used for that file.
# The first language is the default language and the respective language server will be used as a fallback.
# Note that when using the JetBrains backend, language servers are not used and this list is correspondingly ignored.
languages:
- go
# the encoding used by text files in the project
# For a list of possible encodings, see https://docs.python.org/3.11/library/codecs.html#standard-encodings
encoding: "utf-8"
# line ending convention to use when writing source files.
# Possible values: unset (use global setting), "lf", "crlf", or "native" (platform default)
# This does not affect Serena's own files (e.g. memories and configuration files), which always use native line endings.
line_ending:
# The language backend to use for this project.
# If not set, the global setting from serena_config.yml is used.
# Valid values: LSP, JetBrains
# Note: the backend is fixed at startup. If a project with a different backend
# is activated post-init, an error will be returned.
language_backend:
# whether to use project's .gitignore files to ignore files
ignore_all_files_in_gitignore: true
# advanced configuration option allowing to configure language server-specific options.
# Maps the language key to the options.
# Have a look at the docstring of the constructors of the LS implementations within solidlsp (e.g., for C# or PHP) to see which options are available.
# No documentation on options means no options are available.
ls_specific_settings: {}
# list of additional workspace folder paths for cross-package reference support (e.g. in monorepos).
# Paths can be absolute or relative to the project root.
# Each folder is registered as an LSP workspace folder, enabling language servers to discover
# symbols and references across package boundaries.
# Currently supported for: TypeScript.
# Example:
# additional_workspace_folders:
# - ../sibling-package
# - ../shared-lib
additional_workspace_folders: []
# list of additional paths to ignore in this project.
# Same syntax as gitignore, so you can use * and **.
# Note: global ignored_paths from serena_config.yml are also applied additively.
ignored_paths: []
# whether the project is in read-only mode
# If set to true, all editing tools will be disabled and attempts to use them will result in an error
# Added on 2025-04-18
read_only: false
# list of tool names to exclude.
# This extends the existing exclusions (e.g. from the global configuration)
# Find the list of tools here: https://oraios.github.io/serena/01-about/035_tools.html
excluded_tools: []
# list of tools to include that would otherwise be disabled (particularly optional tools that are disabled by default).
# This extends the existing inclusions (e.g. from the global configuration).
# Find the list of tools here: https://oraios.github.io/serena/01-about/035_tools.html
included_optional_tools: []
# fixed set of tools to use as the base tool set (if non-empty), replacing Serena's default set of tools.
# This cannot be combined with non-empty excluded_tools or included_optional_tools.
# Find the list of tools here: https://oraios.github.io/serena/01-about/035_tools.html
fixed_tools: []
# list of mode names that are to be activated by default, overriding the setting in the global configuration.
# The full set of modes to be activated is base_modes (from global config) + default_modes + added_modes.
# If the setting is undefined/empty, the default_modes from the global configuration (serena_config.yml) apply.
# Otherwise, this overrides the setting from the global configuration (serena_config.yml).
# Therefore, you can set this to [] if you do not want the default modes defined in the global config to apply
# for this project.
# This setting can, in turn, be overridden by CLI parameters (--mode).
# See https://oraios.github.io/serena/02-usage/050_configuration.html#modes
default_modes:
# list of mode names to be activated additionally for this project, e.g. ["query-projects"]
# The full set of modes to be activated is base_modes (from global config) + default_modes + added_modes.
# See https://oraios.github.io/serena/02-usage/050_configuration.html#modes
added_modes:
# initial prompt for the project. It will always be given to the LLM upon activating the project
# (contrary to the memories, which are loaded on demand).
initial_prompt: ""
# time budget (seconds) per tool call for the retrieval of additional symbol information
# such as docstrings or parameter information.
# This overrides the corresponding setting in the global configuration; see the documentation there.
# If null or missing, use the setting from the global configuration.
symbol_info_budget:
# list of regex patterns which, when matched, mark a memory entry as readonly.
# Extends the list from the global configuration, merging the two lists.
read_only_memory_patterns: []
# list of regex patterns for memories to completely ignore.
# Matching memories will not appear in list_memories or activate_project output
# and cannot be accessed via read_memory or write_memory.
# To access ignored memory files, use the read_file tool on the raw file path.
# Extends the list from the global configuration, merging the two lists.
# Example: ["_archive/.*", "_episodes/.*"]
ignored_memory_patterns: []
+82 -5
View File
@@ -17,7 +17,7 @@ import (
) )
var __obi_header_value_string_pattern__ = regexp.MustCompile(`^'\s*([^']*'|"[^"]*")\s*;`) var __obi_header_value_string_pattern__ = regexp.MustCompile(`^'\s*([^']*'|"[^"]*")\s*;`)
var __obi_header_value_numeric_pattern__ = regexp.MustCompile(`^\s*([+-]?\.\d+|[+-]?\d+(\.\d*)?([eE][+-]?\d+)?)\s*;`) var __obi_header_value_numeric_pattern__ = regexp.MustCompile(`^\s*[+-]?(\d+(\.\d*)?|\.\d+)([eE][+-]?\d+)?\s*;`)
var __obi_header_map_int_key__ = regexp.MustCompile("([{,])([0-9]+):") var __obi_header_map_int_key__ = regexp.MustCompile("([{,])([0-9]+):")
func __match__dict__(text []byte) []int { func __match__dict__(text []byte) []int {
@@ -146,6 +146,65 @@ func __match__key__(text []byte) []int {
return []int{} // Not a key return []int{} // Not a key
} }
func __match__array__(text []byte) []int {
state := 0
level := 0
start := 0
instring := byte(0)
for i, r := range text {
if state == 2 {
if r == ';' {
return []int{start, i + 1}
}
if r != ' ' && r != '\t' {
return []int{}
}
}
if state == 0 {
if r == '[' {
level++
state++
start = i
continue
}
if r != ' ' && r != '\t' {
return []int{}
}
continue
}
// state == 1: inside the array
if instring != 0 {
if r == instring {
instring = 0
}
continue
}
if r == '"' || r == '\'' {
instring = r
continue
}
if r == '[' || r == '{' {
level++
continue
}
if r == ']' || r == '}' {
level--
if level == 0 {
state++
}
}
}
return []int{}
}
func __match__general__(text []byte) []int { func __match__general__(text []byte) []int {
for i, r := range text { for i, r := range text {
@@ -242,6 +301,21 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
stop = m[1] + 1 stop = m[1] + 1
} else { } else {
// array value
m = __match__array__(part)
if len(m) > 0 {
bvalue = bytes.TrimSpace(part[m[0]:(m[1] - 1)])
j := bytes.ReplaceAll(bvalue, []byte("'"), []byte(`"`))
j = __obi_header_map_int_key__.ReplaceAll(j, []byte(`$1"$2":`))
arr, err := _parse_json_array_interface(j)
if err != nil {
value = string(bvalue)
} else {
value = arr
}
stop = m[1] + 1
} else {
// Generic value // Generic value
// m = __obi_header_value_general_pattern__.FindIndex(part) // m = __obi_header_value_general_pattern__.FindIndex(part)
@@ -264,6 +338,7 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
// no value // no value
break break
} // End of No value } // End of No value
} // End of not array
} // End of not dict } // End of not dict
} // End of not string } // End of not string
} // End of not numeric } // End of not numeric
@@ -272,6 +347,8 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
case float64: case float64:
if vt == math.Floor(vt) { if vt == math.Floor(vt) {
annotations[key] = int(vt) annotations[key] = int(vt)
} else {
annotations[key] = vt
} }
default: default:
annotations[key] = value annotations[key] = value
@@ -327,9 +404,8 @@ func WriteFastSeqOBIHeade(buffer *bytes.Buffer, sequence *obiseq.BioSequence) {
buffer.WriteString(fmt.Sprintf("%s=", key)) buffer.WriteString(fmt.Sprintf("%s=", key))
buffer.Write(tv) buffer.Write(tv)
buffer.WriteString("; ") buffer.WriteString("; ")
case map[string]int, default:
map[string]string, if obiutils.IsAMap(value) || obiutils.IsASlice(value) || obiutils.IsAnArray(value) {
map[string]interface{}:
tv, err := obiutils.JsonMarshal(t) tv, err := obiutils.JsonMarshal(t)
if err != nil { if err != nil {
log.Fatalf("Cannot convert %v value", value) log.Fatalf("Cannot convert %v value", value)
@@ -338,11 +414,12 @@ func WriteFastSeqOBIHeade(buffer *bytes.Buffer, sequence *obiseq.BioSequence) {
buffer.WriteString(fmt.Sprintf("%s=", key)) buffer.WriteString(fmt.Sprintf("%s=", key))
buffer.Write(tv) buffer.Write(tv)
buffer.WriteString("; ") buffer.WriteString("; ")
default: } else {
buffer.WriteString(fmt.Sprintf("%s=%v; ", key, value)) buffer.WriteString(fmt.Sprintf("%s=%v; ", key, value))
} }
} }
} }
}
if sequence.HasDefinition() { if sequence.HasDefinition() {
buffer.WriteByte(' ') buffer.WriteByte(' ')
+3
View File
@@ -90,6 +90,9 @@ func FormatFastaBatch(batch obiiter.BioSequenceBatch, formater FormatHeader, ski
log.Debugf("FormatFastaBatch: #%d : %d seqs", batch.Order(), batch.Len()) log.Debugf("FormatFastaBatch: #%d : %d seqs", batch.Order(), batch.Len())
for _, seq := range batch.Slice() { for _, seq := range batch.Slice() {
if len(seq.Id()) == 0 {
log.Fatalf("Sequence identifier is empty")
}
if seq.Len() > 0 { if seq.Len() > 0 {
// Write header directly into bs — no intermediate string // Write header directly into bs — no intermediate string
bs.WriteByte('>') bs.WriteByte('>')
+3
View File
@@ -64,6 +64,9 @@ func FormatFastqBatch(batch obiiter.BioSequenceBatch,
first := true first := true
for _, seq := range batch.Slice() { for _, seq := range batch.Slice() {
if len(seq.Id()) == 0 {
log.Fatalf("Sequence identifier is empty")
}
if seq.Len() > 0 { if seq.Len() > 0 {
_formatFastq(&bs, seq, formater) _formatFastq(&bs, seq, formater)
+75 -40
View File
@@ -4,22 +4,21 @@ import "math"
// KmerEntropy computes the entropy of a single encoded k-mer. // KmerEntropy computes the entropy of a single encoded k-mer.
// //
// The algorithm mirrors the lowmask entropy calculation: it decodes the k-mer // The algorithm mirrors the Rust obiskbuilder entropy: it decodes the k-mer
// to a DNA sequence, extracts all sub-words of each size from 1 to levelMax, // to a DNA sequence, extracts all sub-words of each size from 1 to levelMax,
// normalizes them by circular canonical form, counts their frequencies, and // normalizes them by circular canonical form, counts their frequencies, and
// computes Shannon entropy normalized by the maximum possible entropy. // computes Shannon entropy corrected for class sizes, normalized by the
// maximum possible entropy over 4^ws raw bins.
// The returned value is the minimum entropy across all word sizes. // The returned value is the minimum entropy across all word sizes.
// //
// Correction for small sequences: the raw entropy H = log(N) - Σ f·log(f)/N
// under-estimates the true complexity when many raw words collapse to the same
// canonical form. Adding Σ f·log(class_size)/N recovers the entropy of the
// underlying uncollapsed distribution (assuming uniform mixing within each
// equivalence class).
//
// A value close to 0 indicates very low complexity (e.g. "AAAA..."), // A value close to 0 indicates very low complexity (e.g. "AAAA..."),
// while a value close to 1 indicates high complexity. // while a value close to 1 indicates high complexity.
//
// Parameters:
// - kmer: the encoded k-mer (2 bits per base)
// - k: the k-mer size
// - levelMax: maximum sub-word size for entropy (typically 6)
//
// Returns:
// - minimum normalized entropy across all word sizes 1..levelMax
func KmerEntropy(kmer uint64, k int, levelMax int) float64 { func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
if k < 1 || levelMax < 1 { if k < 1 || levelMax < 1 {
return 1.0 return 1.0
@@ -35,7 +34,7 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
var seqBuf [32]byte var seqBuf [32]byte
seq := DecodeKmer(kmer, k, seqBuf[:]) seq := DecodeKmer(kmer, k, seqBuf[:])
// Pre-compute nLogN lookup (same as lowmask) // Pre-compute nLogN lookup
nLogN := make([]float64, k+1) nLogN := make([]float64, k+1)
for i := 1; i <= k; i++ { for i := 1; i <= k; i++ {
nLogN[i] = float64(i) * math.Log(float64(i)) nLogN[i] = float64(i) * math.Log(float64(i))
@@ -51,6 +50,23 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
} }
} }
// Build ln(class_size) tables: for each canonical form, how many raw
// words map to it under circular normalization.
classLogSizeTables := make([][]float64, levelMax+1)
for ws := 1; ws <= levelMax; ws++ {
tableSize := 1 << (ws * 2)
classSize := make([]int, tableSize)
for code := 0; code < tableSize; code++ {
classSize[normTables[ws][code]]++
}
classLogSizeTables[ws] = make([]float64, tableSize)
for j := 0; j < tableSize; j++ {
if classSize[j] > 0 {
classLogSizeTables[ws][j] = math.Log(float64(classSize[j]))
}
}
}
minEntropy := math.MaxFloat64 minEntropy := math.MaxFloat64
for ws := 1; ws <= levelMax; ws++ { for ws := 1; ws <= levelMax; ws++ {
@@ -75,23 +91,13 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
table[normWord]++ table[normWord]++
} }
// Compute Shannon entropy // Compute emax over 4^ws raw bins (uncollapsed distribution).
floatNwords := float64(nwords) floatNwords := float64(nwords)
logNwords := math.Log(floatNwords) logNwords := math.Log(floatNwords)
na := tableSize // 4^ws
var sumNLogN float64
for j := 0; j < tableSize; j++ {
n := table[j]
if n > 0 {
sumNLogN += nLogN[n]
}
}
// Compute emax (maximum possible entropy for this word size)
na := CanonicalCircularKmerCount(ws)
var emax float64 var emax float64
if nwords < na { if nwords < na {
emax = math.Log(float64(nwords)) emax = logNwords
} else { } else {
cov := nwords / na cov := nwords / na
remains := nwords - (na * cov) remains := nwords - (na * cov)
@@ -105,7 +111,19 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
continue continue
} }
entropy := (logNwords - sumNLogN/floatNwords) / emax // Accumulate Σ f·log(f) and Σ f·log(class_size) over canonical forms.
classLogSize := classLogSizeTables[ws]
var sumNLogN, sumClassLogN float64
for j := 0; j < tableSize; j++ {
n := table[j]
if n > 0 {
sumNLogN += nLogN[n]
sumClassLogN += float64(n) * classLogSize[j]
}
}
// Corrected entropy: H_raw ≈ log(N) + (Σf·log(s) - Σf·log(f)) / N
entropy := (logNwords + sumClassLogN/floatNwords - sumNLogN/floatNwords) / emax
if entropy < 0 { if entropy < 0 {
entropy = 0 entropy = 0
} }
@@ -134,6 +152,7 @@ type KmerEntropyFilter struct {
threshold float64 threshold float64
nLogN []float64 nLogN []float64
normTables [][]int normTables [][]int
classLogSizeTables [][]float64
emaxValues []float64 emaxValues []float64
logNwords []float64 logNwords []float64
// Pre-allocated frequency tables reused across Entropy() calls. // Pre-allocated frequency tables reused across Entropy() calls.
@@ -142,11 +161,6 @@ type KmerEntropyFilter struct {
} }
// NewKmerEntropyFilter creates an entropy filter with pre-computed tables. // NewKmerEntropyFilter creates an entropy filter with pre-computed tables.
//
// Parameters:
// - k: the k-mer size
// - levelMax: maximum sub-word size for entropy (typically 6)
// - threshold: entropy threshold (k-mers with entropy <= threshold are rejected)
func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter { func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter {
if levelMax >= k { if levelMax >= k {
levelMax = k - 1 levelMax = k - 1
@@ -169,20 +183,38 @@ func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter
} }
} }
// ln(class_size) for each canonical form under circular normalization.
classLogSizeTables := make([][]float64, levelMax+1)
for ws := 1; ws <= levelMax; ws++ {
tableSize := 1 << (ws * 2)
classSize := make([]int, tableSize)
for code := 0; code < tableSize; code++ {
classSize[normTables[ws][code]]++
}
classLogSizeTables[ws] = make([]float64, tableSize)
for j := 0; j < tableSize; j++ {
if classSize[j] > 0 {
classLogSizeTables[ws][j] = math.Log(float64(classSize[j]))
}
}
}
// Pre-compute emax and logNwords per word size.
// emax uses 4^ws raw bins to match the corrected entropy.
emaxValues := make([]float64, levelMax+1) emaxValues := make([]float64, levelMax+1)
logNwords := make([]float64, levelMax+1) logNwords := make([]float64, levelMax+1)
for ws := 1; ws <= levelMax; ws++ { for ws := 1; ws <= levelMax; ws++ {
nw := k - ws + 1 nw := k - ws + 1
na := CanonicalCircularKmerCount(ws) na := 1 << (ws * 2) // 4^ws raw bins
floatNw := float64(nw)
logNwords[ws] = math.Log(floatNw)
if nw < na { if nw < na {
logNwords[ws] = math.Log(float64(nw)) emaxValues[ws] = logNwords[ws]
emaxValues[ws] = math.Log(float64(nw))
} else { } else {
cov := nw / na cov := nw / na
remains := nw - (na * cov) remains := nw - (na * cov)
f1 := float64(cov) / float64(nw) f1 := float64(cov) / floatNw
f2 := float64(cov+1) / float64(nw) f2 := float64(cov+1) / floatNw
logNwords[ws] = math.Log(float64(nw))
emaxValues[ws] = -(float64(na-remains)*f1*math.Log(f1) + emaxValues[ws] = -(float64(na-remains)*f1*math.Log(f1) +
float64(remains)*f2*math.Log(f2)) float64(remains)*f2*math.Log(f2))
} }
@@ -200,6 +232,7 @@ func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter
threshold: threshold, threshold: threshold,
nLogN: nLogN, nLogN: nLogN,
normTables: normTables, normTables: normTables,
classLogSizeTables: classLogSizeTables,
emaxValues: emaxValues, emaxValues: emaxValues,
logNwords: logNwords, logNwords: logNwords,
freqTables: freqTables, freqTables: freqTables,
@@ -236,7 +269,7 @@ func (ef *KmerEntropyFilter) Entropy(kmer uint64) float64 {
// Count circular-canonical sub-word frequencies // Count circular-canonical sub-word frequencies
tableSize := 1 << (ws * 2) tableSize := 1 << (ws * 2)
table := ef.freqTables[ws] table := ef.freqTables[ws]
clear(table) // reset to zero clear(table)
mask := (1 << (ws * 2)) - 1 mask := (1 << (ws * 2)) - 1
normTable := ef.normTables[ws] normTable := ef.normTables[ws]
@@ -251,19 +284,21 @@ func (ef *KmerEntropyFilter) Entropy(kmer uint64) float64 {
table[normWord]++ table[normWord]++
} }
// Compute Shannon entropy
floatNwords := float64(nwords) floatNwords := float64(nwords)
logNwords := ef.logNwords[ws] logNwords := ef.logNwords[ws]
classLogSize := ef.classLogSizeTables[ws]
var sumNLogN float64 var sumNLogN, sumClassLogN float64
for j := 0; j < tableSize; j++ { for j := 0; j < tableSize; j++ {
n := table[j] n := table[j]
if n > 0 { if n > 0 {
sumNLogN += ef.nLogN[n] sumNLogN += ef.nLogN[n]
sumClassLogN += float64(n) * classLogSize[j]
} }
} }
entropy := (logNwords - sumNLogN/floatNwords) / emax // Corrected entropy: H_raw ≈ log(N) + (Σf·log(s) - Σf·log(f)) / N
entropy := (logNwords + sumClassLogN/floatNwords - sumNLogN/floatNwords) / emax
if entropy < 0 { if entropy < 0 {
entropy = 0 entropy = 0
} }
+1 -1
View File
@@ -3,7 +3,7 @@ package obioptions
// Version is automatically updated by the Makefile from version.txt // Version is automatically updated by the Makefile from version.txt
// The patch number (third digit) is incremented on each push to the repository // The patch number (third digit) is incremented on each push to the repository
var _Version = "Release 4.4.42" var _Version = "Release 4.4.46"
// Version returns the version of the obitools package. // Version returns the version of the obitools package.
// //
+18
View File
@@ -364,6 +364,24 @@ func (s *BioSequence) GetIntSlice(key string) ([]int, bool) {
return val, ok return val, ok
} }
func (s *BioSequence) GetMapOfIntSlice(key string) (map[string][]int, bool) {
v, ok := s.GetAttribute(key)
if !ok {
return nil, false
}
val, err := obiutils.InterfaceToMapOfIntSlice(v)
return val, err == nil
}
func (s *BioSequence) GetMapOfStringSlice(key string) (map[string][]string, bool) {
v, ok := s.GetAttribute(key)
if !ok {
return nil, false
}
val, err := obiutils.InterfaceToMapOfStringSlice(v)
return val, err == nil
}
// Count returns the value of the "count" attribute of the BioSequence. // Count returns the value of the "count" attribute of the BioSequence.
// //
// The count of a sequence is the number of times it has been observed in the dataset. // The count of a sequence is the number of times it has been observed in the dataset.
+1 -1
View File
@@ -103,7 +103,7 @@ func TestNewBioSequence(t *testing.T) {
// Return type: None. // Return type: None.
func TestNewBioSequenceWithQualities(t *testing.T) { func TestNewBioSequenceWithQualities(t *testing.T) {
id := "123" id := "123"
sequence := []byte("ATGC") sequence := []byte("atgc")
definition := "DNA sequence" definition := "DNA sequence"
qualities := []byte("1234") qualities := []byte("1234")
+27
View File
@@ -141,6 +141,33 @@ var OBILang = gval.NewLanguage(
gval.Function("max", func(args ...interface{}) (interface{}, error) { gval.Function("max", func(args ...interface{}) (interface{}, error) {
return obiutils.Max(args[0]) return obiutils.Max(args[0])
}), }),
gval.Function("whichmax", func(args ...interface{}) (interface{}, error) {
result, err := obiutils.WhichMax(args[0])
if idx, ok := result.(int); ok {
return float64(idx), nil
}
return result, err
}),
gval.Function("whichmin", func(args ...interface{}) (interface{}, error) {
result, err := obiutils.WhichMin(args[0])
if idx, ok := result.(int); ok {
return float64(idx), nil
}
return result, err
}),
gval.Function("filtermin", func(args ...interface{}) (interface{}, error) {
return obiutils.FilterMin(args[0], args[1])
}),
gval.Function("filtermax", func(args ...interface{}) (interface{}, error) {
return obiutils.FilterMax(args[0], args[1])
}),
gval.Function("saturatingsub", func(args ...interface{}) (interface{}, error) {
return obiutils.SaturatingSub(args[0], args[1])
}),
gval.Function("contains", func(args ...interface{}) (interface{}, error) { gval.Function("contains", func(args ...interface{}) (interface{}, error) {
if obiutils.IsAMap(args[0]) { if obiutils.IsAMap(args[0]) {
val := reflect.ValueOf(args[0]).MapIndex(reflect.ValueOf(args[1])) val := reflect.ValueOf(args[0]).MapIndex(reflect.ValueOf(args[1]))
+38
View File
@@ -276,6 +276,44 @@ func InterfaceToStringMap(i interface{}) (val map[string]string, err error) {
return return
} }
func InterfaceToMapOfIntSlice(i interface{}) (val map[string][]int, err error) {
err = nil
switch m := i.(type) {
case map[string][]int:
val = m
case map[string]interface{}:
val = make(map[string][]int, len(m))
for k, v := range m {
val[k], err = InterfaceToIntSlice(v)
if err != nil {
return
}
}
default:
err = &NotAMapInt{"value attribute cannot be casted to a map[string][]int"}
}
return
}
func InterfaceToMapOfStringSlice(i interface{}) (val map[string][]string, err error) {
err = nil
switch m := i.(type) {
case map[string][]string:
val = m
case map[string]interface{}:
val = make(map[string][]string, len(m))
for k, v := range m {
val[k], err = InterfaceToStringSlice(v)
if err != nil {
return
}
}
default:
err = &NotAMapInt{"value attribute cannot be casted to a map[string][]string"}
}
return
}
func InterfaceToStringSlice(i interface{}) (val []string, err error) { func InterfaceToStringSlice(i interface{}) (val []string, err error) {
err = nil err = nil
+365 -4
View File
@@ -34,6 +34,26 @@ func MinMaxSlice[T constraints.Ordered](vec []T) (min, max T) {
return return
} }
func FilterMinSlice[T constraints.Ordered](vec []T, minimum T) []T {
result := make([]T, 0, len(vec))
for _, v := range vec {
if v >= minimum {
result = append(result, v)
}
}
return result
}
func FilterMaxSlice[T constraints.Ordered](vec []T, maximum T) []T {
result := make([]T, 0, len(vec))
for _, v := range vec {
if v <= maximum {
result = append(result, v)
}
}
return result
}
func MaxMap[K comparable, T constraints.Ordered](values map[K]T) (K, T, error) { func MaxMap[K comparable, T constraints.Ordered](values map[K]T) (K, T, error) {
var maxKey K var maxKey K
var maxValue T var maxValue T
@@ -73,6 +93,46 @@ func MinMap[K comparable, T constraints.Ordered](values map[K]T) (K, T, error) {
return minKey, minValue, nil return minKey, minValue, nil
} }
func FilterMinMap[K comparable, T constraints.Ordered](values map[K]T, minimum T) map[K]T {
result := make(map[K]T)
for k, v := range values {
if v >= minimum {
result[k] = v
}
}
return result
}
func FilterMaxMap[K comparable, T constraints.Ordered](values map[K]T, maximum T) map[K]T {
result := make(map[K]T)
for k, v := range values {
if v <= maximum {
result[k] = v
}
}
return result
}
func SaturatingSubSlice[T Numeric](vec []T, sub T) []T {
result := make([]T, len(vec))
for i, v := range vec {
if v > sub {
result[i] = v - sub
}
}
return result
}
func SaturatingSubMap[K comparable, T Numeric](values map[K]T, sub T) map[K]T {
result := make(map[K]T)
for k, v := range values {
if v > sub {
result[k] = v - sub
}
}
return result
}
// Min returns the smallest element in a slice/array or map, // Min returns the smallest element in a slice/array or map,
// or the value itself if data is a single comparable value. // or the value itself if data is a single comparable value.
// Returns an error if the container is empty or the type is unsupported. // Returns an error if the container is empty or the type is unsupported.
@@ -135,11 +195,121 @@ func Max(data interface{}) (interface{}, error) {
} }
} }
func FilterMin(data interface{}, minimum interface{}) (interface{}, error) {
v := reflect.ValueOf(data)
switch v.Kind() {
case reflect.Slice, reflect.Array:
if v.Len() == 0 {
return nil, errors.New("empty slice or array")
}
return filterMinFromIterable(v, minimum)
case reflect.Map:
if v.Len() == 0 {
return nil, errors.New("empty map")
}
return filterMinFromMap(v, minimum)
default:
if !isOrderedKind(v.Kind()) {
return nil, fmt.Errorf("unsupported type: %s", v.Kind())
}
return data, nil
}
}
func FilterMax(data interface{}, maximum interface{}) (interface{}, error) {
v := reflect.ValueOf(data)
switch v.Kind() {
case reflect.Slice, reflect.Array:
if v.Len() == 0 {
return nil, errors.New("empty slice or array")
}
return filterMaxFromIterable(v, maximum)
case reflect.Map:
if v.Len() == 0 {
return nil, errors.New("empty map")
}
return filterMaxFromMap(v, maximum)
default:
if !isOrderedKind(v.Kind()) {
return nil, fmt.Errorf("unsupported type: %s", v.Kind())
}
return data, nil
}
}
func SaturatingSub(data interface{}, sub interface{}) (interface{}, error) {
v := reflect.ValueOf(data)
switch v.Kind() {
case reflect.Slice, reflect.Array:
return saturatingSubFromIterable(v, sub)
case reflect.Map:
return saturatingSubFromMap(v, sub)
default:
if !isNumericKind(v.Kind()) {
return nil, fmt.Errorf("unsupported type: %s", v.Kind())
}
r, err := saturatingSubValues(v, reflect.ValueOf(sub))
if err != nil {
return nil, err
}
return r.Interface(), nil
}
}
func saturatingSubFromIterable(v reflect.Value, sub interface{}) (interface{}, error) {
subVal := reflect.ValueOf(sub)
result := reflect.MakeSlice(v.Type(), v.Len(), v.Len())
for i := 0; i < v.Len(); i++ {
r, err := saturatingSubValues(v.Index(i), subVal)
if err != nil {
return nil, err
}
result.Index(i).Set(r)
}
return result.Interface(), nil
}
func saturatingSubFromMap(v reflect.Value, sub interface{}) (interface{}, error) {
subVal := reflect.ValueOf(sub)
result := reflect.MakeMap(v.Type())
for _, key := range v.MapKeys() {
r, err := saturatingSubValues(v.MapIndex(key), subVal)
if err != nil {
return nil, err
}
if !r.IsZero() {
result.SetMapIndex(key, r)
}
}
return result.Interface(), nil
}
func saturatingSubValues(a, b reflect.Value) (reflect.Value, error) {
result := reflect.New(a.Type()).Elem()
switch a.Kind() {
case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
if av, bv := a.Int(), b.Int(); av > bv {
result.SetInt(av - bv)
}
case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64:
if av, bv := a.Uint(), b.Uint(); av > bv {
result.SetUint(av - bv)
}
case reflect.Float32, reflect.Float64:
if av, bv := a.Float(), b.Float(); av > bv {
result.SetFloat(av - bv)
}
default:
return reflect.Value{}, fmt.Errorf("unsupported type for saturating subtraction: %s", a.Kind())
}
return result, nil
}
// maxFromIterable scans a slice/array to find the maximum. // maxFromIterable scans a slice/array to find the maximum.
func maxFromIterable(v reflect.Value) (interface{}, error) { func maxFromIterable(v reflect.Value) (interface{}, error) {
var best reflect.Value var best reflect.Value
for i := 0; i < v.Len(); i++ { for i := 0; i < v.Len(); i++ {
elem := v.Index(i) elem := unwrapInterface(v.Index(i))
if !isOrderedKind(elem.Kind()) { if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind()) return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
} }
@@ -154,7 +324,7 @@ func maxFromIterable(v reflect.Value) (interface{}, error) {
func minFromIterable(v reflect.Value) (interface{}, error) { func minFromIterable(v reflect.Value) (interface{}, error) {
var minVal reflect.Value var minVal reflect.Value
for i := 0; i < v.Len(); i++ { for i := 0; i < v.Len(); i++ {
elem := v.Index(i) elem := unwrapInterface(v.Index(i))
if !isOrderedKind(elem.Kind()) { if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind()) return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
} }
@@ -165,12 +335,182 @@ func minFromIterable(v reflect.Value) (interface{}, error) {
return minVal.Interface(), nil return minVal.Interface(), nil
} }
func filterMinFromIterable(v reflect.Value, minimum interface{}) (interface{}, error) {
minVal := reflect.ValueOf(minimum)
result := reflect.MakeSlice(v.Type(), 0, v.Len())
for i := 0; i < v.Len(); i++ {
elem := unwrapInterface(v.Index(i))
if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if !less(elem, minVal) { // elem >= minimum
result = reflect.Append(result, elem)
}
}
return result.Interface(), nil
}
func filterMaxFromIterable(v reflect.Value, maximum interface{}) (interface{}, error) {
maxVal := reflect.ValueOf(maximum)
result := reflect.MakeSlice(v.Type(), 0, v.Len())
for i := 0; i < v.Len(); i++ {
elem := unwrapInterface(v.Index(i))
if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if !greater(elem, maxVal) { // elem <= maximum
result = reflect.Append(result, elem)
}
}
return result.Interface(), nil
}
// whichMaxFromIterable returns the index of the maximum element in a slice/array.
func whichMaxFromIterable(v reflect.Value) (int, error) {
var best reflect.Value
bestIdx := 0
for i := 0; i < v.Len(); i++ {
elem := unwrapInterface(v.Index(i))
if !isOrderedKind(elem.Kind()) {
return 0, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if i == 0 || greater(elem, best) {
best = elem
bestIdx = i
}
}
return bestIdx, nil
}
// whichMinFromIterable returns the index of the minimum element in a slice/array.
func whichMinFromIterable(v reflect.Value) (int, error) {
var minVal reflect.Value
minIdx := 0
for i := 0; i < v.Len(); i++ {
elem := unwrapInterface(v.Index(i))
if !isOrderedKind(elem.Kind()) {
return 0, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if i == 0 || less(elem, minVal) {
minVal = elem
minIdx = i
}
}
return minIdx, nil
}
// whichMaxFromMap returns the key associated with the maximum value in a map.
func whichMaxFromMap(v reflect.Value) (interface{}, error) {
var best reflect.Value
var bestKey reflect.Value
first := true
for _, key := range v.MapKeys() {
elem := unwrapInterface(v.MapIndex(key))
if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if first || greater(elem, best) {
best = elem
bestKey = key
first = false
}
}
return bestKey.Interface(), nil
}
// whichMinFromMap returns the key associated with the minimum value in a map.
func whichMinFromMap(v reflect.Value) (interface{}, error) {
var minVal reflect.Value
var minKey reflect.Value
first := true
for _, key := range v.MapKeys() {
elem := unwrapInterface(v.MapIndex(key))
if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if first || less(elem, minVal) {
minVal = elem
minKey = key
first = false
}
}
return minKey.Interface(), nil
}
// WhichMax returns the key (for a map) or index (for a slice/array) of the maximum value.
func WhichMax(data interface{}) (interface{}, error) {
v := reflect.ValueOf(data)
switch v.Kind() {
case reflect.Slice, reflect.Array:
if v.Len() == 0 {
return nil, errors.New("empty slice or array")
}
return whichMaxFromIterable(v)
case reflect.Map:
if v.Len() == 0 {
return nil, errors.New("empty map")
}
return whichMaxFromMap(v)
default:
return nil, fmt.Errorf("unsupported type: %s", v.Kind())
}
}
// WhichMin returns the key (for a map) or index (for a slice/array) of the minimum value.
func WhichMin(data interface{}) (interface{}, error) {
v := reflect.ValueOf(data)
switch v.Kind() {
case reflect.Slice, reflect.Array:
if v.Len() == 0 {
return nil, errors.New("empty slice or array")
}
return whichMinFromIterable(v)
case reflect.Map:
if v.Len() == 0 {
return nil, errors.New("empty map")
}
return whichMinFromMap(v)
default:
return nil, fmt.Errorf("unsupported type: %s", v.Kind())
}
}
func filterMinFromMap(v reflect.Value, minimum interface{}) (interface{}, error) {
minVal := reflect.ValueOf(minimum)
result := reflect.MakeMap(v.Type())
for _, key := range v.MapKeys() {
elem := unwrapInterface(v.MapIndex(key))
if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if !less(elem, minVal) { // elem >= minimum
result.SetMapIndex(key, elem)
}
}
return result.Interface(), nil
}
func filterMaxFromMap(v reflect.Value, maximum interface{}) (interface{}, error) {
maxVal := reflect.ValueOf(maximum)
result := reflect.MakeMap(v.Type())
for _, key := range v.MapKeys() {
elem := unwrapInterface(v.MapIndex(key))
if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
}
if !greater(elem, maxVal) { // elem <= maximum
result.SetMapIndex(key, elem)
}
}
return result.Interface(), nil
}
// maxFromMap scans map values to find the maximum. // maxFromMap scans map values to find the maximum.
func maxFromMap(v reflect.Value) (interface{}, error) { func maxFromMap(v reflect.Value) (interface{}, error) {
var best reflect.Value var best reflect.Value
first := true first := true
for _, key := range v.MapKeys() { for _, key := range v.MapKeys() {
elem := v.MapIndex(key) elem := unwrapInterface(v.MapIndex(key))
if !isOrderedKind(elem.Kind()) { if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind()) return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
} }
@@ -187,7 +527,7 @@ func minFromMap(v reflect.Value) (interface{}, error) {
var minVal reflect.Value var minVal reflect.Value
first := true first := true
for _, key := range v.MapKeys() { for _, key := range v.MapKeys() {
elem := v.MapIndex(key) elem := unwrapInterface(v.MapIndex(key))
if !isOrderedKind(elem.Kind()) { if !isOrderedKind(elem.Kind()) {
return nil, fmt.Errorf("unsupported element type: %s", elem.Kind()) return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
} }
@@ -199,6 +539,27 @@ func minFromMap(v reflect.Value) (interface{}, error) {
return minVal.Interface(), nil return minVal.Interface(), nil
} }
func isNumericKind(k reflect.Kind) bool {
switch k {
case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64,
reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64,
reflect.Float32, reflect.Float64:
return true
default:
return false
}
}
// unwrapInterface returns v.Elem() when v holds an interface value, otherwise v unchanged.
// This is necessary when iterating map[string]interface{} or []interface{} via reflection:
// the element Kind is reflect.Interface, not the underlying concrete type.
func unwrapInterface(v reflect.Value) reflect.Value {
if v.Kind() == reflect.Interface {
return v.Elem()
}
return v
}
// isOrderedKind reports whether k supports comparison ordering. // isOrderedKind reports whether k supports comparison ordering.
func isOrderedKind(k reflect.Kind) bool { func isOrderedKind(k reflect.Kind) bool {
switch k { switch k {
+1 -1
View File
@@ -1 +1 @@
4.4.42 4.4.46