Reimplement calcscore.js in C++

[ihh]

Motivated by #7, the idea is to accelerate score calculations inside the AWS Lambda call.

Essentially you'd be reimplementing calcscore.js but reading the history in as a CigarTree directly from a JSON file. The idea will be that the Lambda script will then pass data to this program (probably by saving the CigarTree and sequences as a temporary local file).

C++ is probably easiest. I'd probably use Nick Lohmann's JSON parser. I think you could get away without using a matrix library since the only matrix operations are multiplications and it's probably more efficient to just code row-column multiplication by hand.

Originally posted by @ihh in #7

It is probably worth just confirming that you can reproduce the timing behavior locally by running calcscore.js.

A quicker workaround might be to do some performance profiling of the JavaScript and use that to drive some data-guided optimizations. However in the longer term the C++ route will be faster.

The fiddliest algorithm involved is Felsenstein's pruning algorithm. I do have various C++ implementations of Felsenstein's algorithm kicking around (e.g. here) but I suspect they'll have too much baggage to save you real time; you might as well implement it yourself, it's a fun coding exercise.

Originally posted by @ihh in #7

[ihh]

I've added a stub, calcscore.cpp, for this. Currently it just parses in some JSON from stdin, and regurgitates it in expanded form to stdout.

In calcscore.cpp's version of subLogLike, instead of using precomputed internalsByColumn, leavesByColumn, and branchesByColumn lists, we'll need to compute them on the fly for each column. We do (in calcscore.cpp) know the root node for each column, so we can just iterate upwards from the leaves, ignoring any nodes that have a gap in the current column.

[ihh]

Here's the prompt I used to start this off (I needed to troubleshoot/fix it a bunch)

Here’s an example of a JSON dialect I call “sequence-decorated CIGAR tree”.

(insert data/align.seqcigartree.json here)

It represents a binary phylogenetic tree with pairiwse alignments on the branches and sequences at the leaves. The JSON syntax tree reflects the phylogeny, via "child" properties at internal nodes (always arrays of length 2). A CIGAR string (run-length encoded pairwise alignment) on each branch, indicating parent->child alignment (for internal nodes; at the root node the CIGAR string must be a run of insertions the same length as the root sequence).

Here is a formal schema for sequence-decorated CIGAR tree JSON. You may not need to reference it directly in generated code, since I am not asking for code that validates the input just yet. I want something as lightweight as possible, and for now we will trust the input.

{
"$schema": "https://json-schema.org/draft/2020-12/schema",
"$id": "https://example.org/cigartree.schema.json",
"title": "CIGAR Tree Root Node",
"type": "object",
"required": ["child"],
"properties": {
"id": { "type": "string" },
"child": {
"type": "array",
"minItems": 1,
"items": { "$ref": "#/$defs/internalOrLeaf" }
}
},
"additionalProperties": false,
"$defs": {
"internalOrLeaf": {
"type": "object",
"required": ["cigar", "distance"],
"properties": {
"id": { "type": "string" },
"cigar": {
"type": "string",
"pattern": "^([0-9]+[MID])+$"
},
"distance": {
"type": "number",
"minimum": 0
},
"child": {
"type": "array",
"minItems": 1,
"items": { "$ref": "#/$defs/internalOrLeaf" }
},
"seq": { "type": "string" }
},
"additionalProperties": false,
"allOf": [
{
"if": {
"required": ["child"]
},
"then": {
"not": { "required": ["seq"] }
}
},
{
"if": {
"required": ["seq"]
},
"then": {
"not": { "required": ["child"] }
}
}
]
}
}
}

Here is a schema for parameter files:
{
"$schema": "https://json-schema.org/draft/2020-12/schema",
"title": "TKF92Model",
"type": "object",
"required": ["alphabet", "seqlen", "inslen", "reslife", "root", "evecs_l", "evals", "evecs_r"],
"properties": {
"alphabet": {
"type": "string",
"description": "A string of single-character symbols, e.g. amino acids or nucleotides"
},
"seqlen": {
"type": "number",
"description": "Expected sequence length at equilibrium",
"minimum": 0
},
"inslen": {
"type": "number",
"description": "Expected insertion length",
"minimum": 0
},
"reslife": {
"type": "number",
"minimum": 0
},
"root": {
"type": "array",
"description": "Stationary distribution over alphabet symbols",
"minItems": 1,
"items": {
"type": "number"
}
},
"evecs_l": {
"type": "array",
"description": "Left eigenvectors of substitution rate matrix (each is a vector of floats)",
"items": {
"type": "array",
"items": {
"type": "number"
}
}
},
"evals": {
"type": "array",
"description": "Eigenvalues of substitution rate matrix",
"items": {
"type": "number"
}
},
"evecs_r": {
"type": "array",
"description": "Right eigenvectors of substitution rate matrix (each is a vector of floats)",
"items": {
"type": "array",
"items": {
"type": "number"
}
}
}
},
"additionalProperties": false
}

Here is an example of a parameter file:

(insert data/lg08evol.json here)

Here is some JavaScript code to expand a CIGAR tree to a multiple alignment and a tree:

(insert expandCigarTree from js/cigartree.js)

Using all the above as an example/reference, here is what I want:

A C++ program with as few well-trusted dependencies as possible (e.g. gason for JSON parsing + STL; try to avoid anything as heavy as Boost or GSL) that reads its entire standard input into a string, parses that as JSON, extracts "cigartree" and "params" properties from top-level object, parses "cigartree" as CIGAR tree JSON and creates multiple alignment & tree data structures (as returned by expandCigarTree, above), parses the "params" object as parameter JSON, pulls out the evecs_l and evecs_r matrices (e.g. as vectors of vectors) and the root & evals vectors, along with the other params ("seqlen", "inslen", "reslife"), pulls out the alphabet as a string, and then prints the entire data structure back out as a JSON object (combining a CIGAR tree and parameters) without referring back to the originally parsed JSON �(i.e. using only the internal variables & structures that the JSON-specified info has been parsed into), as a test.

[ihh]

In my last commit (6fe3a76) I updated calcscore.js so that it can be used in exactly the same way as the planned C++ program, reading one combined JSON parameter+data file from standard input. This is an optional mode to make it easier to test the C++ code (you need specify the -p option on the command line, otherwise it will look for separate tree, alignment, and parameter files)

[ihh]

Also added some Rust code - might be better to do this in Rust (better language to build on, easier webassembly integration) but probably doesn't matter too much for this task

[ShreyaKakhandiki]

Here is a table showing some basic benchmarking I did on calcscore.js. Based on this, it doesn't seem like scoring should be a huge bottleneck, as the largest file does not take longer than 2 seconds.

tree_file	align_file	tree_size	align_size	total_ms	expand_ms	leaf_check_ms	count_gaps_ms	subll_compute_ms	subll_sum_ms	ggi_params_ms	gapll_compute_ms	gapll_sum_ms	null_counts_ms	null_sub_ms	null_gap_ms	rss_mb	alignment_cols	branch_count
PF00003.mini_dummy.nh	PF00003.mini.fa	117	3604	31.06475	1.333125	0.047916	0.316833	28.722584	0.00075	0.068334	0.421083	0.000292	0.1115	0.011958	0.006083	142	6	6
PF00003.small_dummy.nh	PF00003.small.fa	806	25233	129.159041	4.883375	0.079042	1.834083	120.89675	0.000875	0.070583	0.8715	0.000333	0.479	0.013333	0.002792	145	36	36
PF00003.med_dummy.nh	PF00003.med.fa	3498	108195	514.471042	17.901417	0.103958	6.352834	485.395708	0.001542	0.079792	3.057291	0.001	1.511	0.01925	0.00425	147	151	151
PF00003.28_dummy.nh	PF00003.28.fa	12927	387373	1947.914292	69.317333	0.225166	25.109416	1837.12075	0.001417	0.086917	10.84175	0.0015	5.098625	0.046292	0.005208	167	539	539