I have some chemical models, which are based on sparse matrices of some coefficients. So, given the model parameters, I generate F# code based on only non-zero elements of these coefficients. The generated model is then fed to ALGLIB (http://www.alglib.net/) ODE solver. The matrices of coefficients are about 99.9% to 99.99% sparse, which means that only 0.01% to 0.1% of coefficients are not exact zeros. Below is a very simplified example of how a generated F# model file looks like. It is the function update (x : array<double>) : array<double> that is fed to ALGLIB ODE solver using 64-bit FSI.
Now, the ALGLIB ODE solver is perfectly capable of handling at least 1M variables for a simple input function. I've tested that and it works with no issues. I have under 10K variables for a typical model.
However, when I increase the model size, I start getting stack overflow exception at run time: the model with about 100K LOC works fine, but the model with about 150K LOC fails with stack overflow exception.
I am guessing that this is related to how initialization / processing of large "hard-coded" arrays is handled and I wonder how should I tweak the generated code OR how can I increase the stack size for FSI and/or F# 64-bit program, let's say to 1 GB.
I shall stress that this is not a typical recursive function stack overflow issue, but just the overall size of the model, which causes the problem.
If you look at update function, you will see that it has a generated array, each element of which is produced by taking another array and applying |> Array.sum. This becomes huge for a large model and I am guessing that that could be causing a stack overflow.
Thanks a lot!
PS Below is a very simplified example of the model. It does have all the necessary structures that appear in the real model.
namespace Model
open Clm.Substances
open Clm.Model
open Clm.ReactionTypes
module ModelData =
let seedValue = 123456
let numberOfAminoAcids = NumberOfAminoAcids.OneAminoAcid
let maxPeptideLength = MaxPeptideLength.TwoMax
let numberOfSubstances = 7
let aminoAcids = AminoAcid.getAminoAcids numberOfAminoAcids
let chiralAminoAcids = ChiralAminoAcid.getAminoAcids numberOfAminoAcids
let peptides = Peptide.getPeptides maxPeptideLength numberOfAminoAcids
let allSubst =
[ Substance.food ]
@
(chiralAminoAcids |> List.map (fun a -> Chiral a))
@
(peptides |> List.map (fun p -> PeptideChain p))
let allInd = allSubst |> List.mapi (fun i s -> (s, i)) |> Map.ofList
let getTotalSubst (x : array<double>) =
[|
x.[0] // Y
x.[1] // A
x.[2] // a
2.0 * x.[3] // AA
2.0 * x.[4] // Aa
2.0 * x.[5] // aA
2.0 * x.[6] // aa
|]
|> Array.sum
let getTotals (x : array<double>) =
[|
// A
(
[|
x.[1] // A
2.0 * x.[3] // AA
x.[4] // Aa
x.[5] // aA
|]
|> Array.sum
,
[|
x.[2] // a
x.[4] // Aa
x.[5] // aA
2.0 * x.[6] // aa
|]
|> Array.sum
)
|]
let update (x : array<double>) : array<double> =
let xSum = (x |> Array.sum) - x.[0]
let xSumN =
[|
1.0 * x.[1] // A
1.0 * x.[2] // a
2.0 * x.[3] // AA
2.0 * x.[4] // Aa
2.0 * x.[5] // aA
2.0 * x.[6] // aa
|]
|> Array.sum
let xSumSquaredN =
[|
1.0 * x.[1] * x.[1] // A
1.0 * x.[2] * x.[2] // a
2.0 * x.[3] * x.[3] // AA
2.0 * x.[4] * x.[4] // Aa
2.0 * x.[5] * x.[5] // aA
2.0 * x.[6] * x.[6] // aa
|]
|> Array.sum
[|
// 0 - Y
[|
0.0001 * x.[2] // a | SynthesisName: Y <-> a
-0.001 * x.[0] // Y | SynthesisName: Y <-> a
0.0001 * x.[1] // A | SynthesisName: Y <-> A
-0.001 * x.[0] // Y | SynthesisName: Y <-> A
|]
|> Array.sum
// 1 - A
[|
0.0001 * x.[5] // aA | LigationName: a + A <-> aA
-0.001 * x.[2] * x.[1] // a + A | LigationName: a + A <-> aA
0.0001 * x.[4] // Aa | LigationName: A + a <-> Aa
-0.001 * x.[1] * x.[2] // A + a | LigationName: A + a <-> Aa
0.0001 * x.[3] // AA | LigationName: A + A <-> AA
0.0001 * x.[3] // AA | LigationName: A + A <-> AA
-0.001 * x.[1] * x.[1] // A + A | LigationName: A + A <-> AA
-0.001 * x.[1] * x.[1] // A + A | LigationName: A + A <-> AA
-0.0001 * x.[1] // A | SynthesisName: Y <-> A
0.001 * x.[0] // Y | SynthesisName: Y <-> A
|]
|> Array.sum
// 2 - a
[|
0.0001 * x.[5] // aA | LigationName: a + A <-> aA
-0.001 * x.[2] * x.[1] // a + A | LigationName: a + A <-> aA
0.0001 * x.[4] // Aa | LigationName: A + a <-> Aa
-0.001 * x.[1] * x.[2] // A + a | LigationName: A + a <-> Aa
0.0001 * x.[6] // aa | LigationName: a + a <-> aa
0.0001 * x.[6] // aa | LigationName: a + a <-> aa
-0.001 * x.[2] * x.[2] // a + a | LigationName: a + a <-> aa
-0.001 * x.[2] * x.[2] // a + a | LigationName: a + a <-> aa
-0.0001 * x.[2] // a | SynthesisName: Y <-> a
0.001 * x.[0] // Y | SynthesisName: Y <-> a
|]
|> Array.sum
// 3 - AA
[|
-0.0001 * x.[3] // AA | LigationName: A + A <-> AA
0.001 * x.[1] * x.[1] // A + A | LigationName: A + A <-> AA
|]
|> Array.sum
// 4 - Aa
[|
-0.0001 * x.[4] // Aa | LigationName: A + a <-> Aa
0.001 * x.[1] * x.[2] // A + a | LigationName: A + a <-> Aa
|]
|> Array.sum
// 5 - aA
[|
-0.0001 * x.[5] // aA | LigationName: a + A <-> aA
0.001 * x.[2] * x.[1] // a + A | LigationName: a + A <-> aA
|]
|> Array.sum
// 6 - aa
[|
-0.0001 * x.[6] // aa | LigationName: a + a <-> aa
0.001 * x.[2] * x.[2] // a + a | LigationName: a + a <-> aa
|]
|> Array.sum
|]
let modelDataParams =
{
numberOfSubstances = 7
numberOfAminoAcids = OneAminoAcid
maxPeptideLength = TwoMax
getTotals = getTotals
getTotalSubst = getTotalSubst
allSubst = allSubst
allInd = allInd
allRawReactions =
[
]
allReactions =
[
(SynthesisName, 2)
(LigationName, 4)
]
}
To summarize our findings.
After a series of trial an errors testing, tracing and debugging different hypothesis @Konstantin was able to discover that the issue was due to the JIT compiler. Apparently it was trying to compile the
updatefunction prior to its first execution. This function was too large and that was causing the Stack Overflow.Splitting the function into smaller ones was the solution.
Bravo Konstantin!