diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,20 @@ # Changelog for srtree +## 3.0.0.6 + +- **AD theta-length guards** (`Algorithm.SRTree.AD.CompiledAD` + + `Algorithm.SRTree.AD.Unboxed`): `CompiledTree` gains `ctMaxParam` and + `evalGrad`/`evalGradVec`/`evalLossVec` fail soft with NaN when the theta + vector is shorter than the tree's largest `Param` index, instead of writing + out of bounds (which corrupted the heap). +- **Param bounds guards** (`Data.SRTree.Eval`): `compile`, `compileLoss` and + `evalTree` return NaN for a `Param` index beyond the theta length. +- **Confidence intervals** (`Algorithm.SRTree.ConfidenceIntervals`): added + `maxParamOf` and `distExtraParams`; `getAllProfiles` now validates theta + against the tree's max param plus the distribution's extra parameters, and + supports 1-parameter models (previously required >= 2). Removed a leftover + `Debug.Trace`. + ## 3.0.0.5 - **Fix NLopt re-entered thousands of times per fit** (`NonlinearOpt` + `SearchSR`): diff --git a/src/Algorithm/SRTree/AD/CompiledAD.hs b/src/Algorithm/SRTree/AD/CompiledAD.hs --- a/src/Algorithm/SRTree/AD/CompiledAD.hs +++ b/src/Algorithm/SRTree/AD/CompiledAD.hs @@ -36,4 +36,5 @@ , ctFcode :: !(VU.Vector Int) -- id -> Uni: fromEnum Function , ctOcode :: !(VU.Vector Int) -- id -> Bin: fromEnum Op , ctVars :: !(VB.Vector (VU.Vector Double)) -- leaf source columns xss ++ [y, yErr] (referenced, not copied) + , ctMaxParam :: !Int -- largest Param index in the tree (-1 if none) } diff --git a/src/Algorithm/SRTree/AD/Unboxed.hs b/src/Algorithm/SRTree/AD/Unboxed.hs --- a/src/Algorithm/SRTree/AD/Unboxed.hs +++ b/src/Algorithm/SRTree/AD/Unboxed.hs @@ -69,8 +69,21 @@ compileTree :: [VU.Vector Double] -> VU.Vector Double -> Maybe (VU.Vector Double) -> Fix SRTree -> CompiledTree compileTree xss ys mYErr tree = CompiledTree { ctNodes = nodes, ctRoot = root, ctDyn = dynArr, ctStatic = staticArr, ctStaticBase = staticBaseArr, ctM = m, ctNPred = root + 1 - , ctKind = kindArr, ctArg = argArr, ctArg2 = arg2Arr, ctFcode = fcodeArr, ctOcode = ocodeArr, ctVars = vars } + , ctKind = kindArr, ctArg = argArr, ctArg2 = arg2Arr, ctFcode = fcodeArr, ctOcode = ocodeArr, ctVars = vars + , ctMaxParam = maxParam } where + -- Largest Param index referenced by the tree (-1 when there is none). + -- The evaluation kernels index theta by Param index, so a theta shorter + -- than maxParam + 1 would read/write out of bounds; the eval entry points + -- use this value to fail soft (NaN) instead of corrupting memory. + maxParam = maximum (negate 1 : cata algP tree) + where algP :: SRTree [Int] -> [Int] + algP (Param ix) = [ix] + algP (Var _) = [] + algP (Const _) = [] + algP (Y _) = [] + algP (Uni _ t) = t + algP (Bin _ l r) = l <> r -- yErr is only defined when mYErr is present (a tree referencing Var -2 -- always pairs with mYErr = Just e, see the likelihood loss wrappers). The -- ctVars list must stay well-defined for every column even when mYErr is @@ -304,7 +317,11 @@ -- arrays) while staying in plain ST: allocation drops from O(nodes * m) -- to O(nodes + params), and the tight inner loops are all unboxed. evalGrad :: CompiledTree -> V.Vector Double -> (Double, V.Vector Double) -evalGrad ct theta = runST $ do +evalGrad ct theta + -- theta shorter than the largest Param index would make the kernel + -- read/write out of bounds; fail soft with NaNs (NLopt rejects NaN). + | ctMaxParam ct >= V.length theta = (0/0, V.replicate (V.length theta) (0/0)) + | otherwise = runST $ do fwd <- VUM.new (root + 1) -- node id -> forward value, current row adj <- VUM.new (root + 1) -- node id -> adjoint (dL/dnode), current row gradM <- VUM.replicate p 0 -- accumulated per-parameter gradient @@ -396,7 +413,11 @@ -- trades the O(nodes + params) scratch of the row-fused version for the -- O(nodes * m) fwd/adj columns of the massiv-style whole-column kernel. evalGradVec :: CompiledTree -> V.Vector Double -> (Double, V.Vector Double) -evalGradVec ct theta = runST $ do +evalGradVec ct theta + -- see the note in `evalGrad`: a theta shorter than the tree's largest + -- Param index must never reach the unsafe kernel accesses. + | ctMaxParam ct >= V.length theta = (0/0, V.replicate (V.length theta) (0/0)) + | otherwise = runST $ do -- Per-chunk buffers of O(stride * chunk) instead of one O(stride * m) -- allocation per call: the fwd/adj matrices are streamed one chunk of -- `chunk` rows at a time, so the per-call allocation drops ~m/chunk x @@ -505,7 +526,9 @@ -- but it cuts the per-call allocation ~30x so this is cheap enough for the -- val-eval hot path that runs once per explored expression. evalLossVec :: CompiledTree -> V.Vector Double -> Double -evalLossVec ct theta = runST $ do +evalLossVec ct theta + | ctMaxParam ct >= V.length theta = 0/0 + | otherwise = runST $ do buf <- VUM.new (stride * chunk) go buf 0 0 where diff --git a/src/Algorithm/SRTree/ConfidenceIntervals.hs b/src/Algorithm/SRTree/ConfidenceIntervals.hs --- a/src/Algorithm/SRTree/ConfidenceIntervals.hs +++ b/src/Algorithm/SRTree/ConfidenceIntervals.hs @@ -31,7 +31,6 @@ import System.IO.Unsafe ( unsafePerformIO ) import Control.Monad.Catch ( catch, SomeException ) import Control.Exception ( evaluate ) -import Debug.Trace ( trace ) -- | profile likelihood algorithms: Bates (classical), ODE (faster), Constrained (fastest) -- The Constrained approach returns only the endpoints. @@ -196,6 +195,26 @@ covMat = toRowMajor cov stdErr = U.generate k (\ix -> sqrt $ abs (covMat U.! (ix * k + ix))) +-- | Largest `Param` index of a tree (-1 if none). Theta vectors are indexed +-- by `Param` index, so profiling requires `U.length theta > maxParamOf tree`, +-- plus the extra parameters appended by the loss (e.g. Gaussian's log-noise). +maxParamOf :: Fix SRTree -> Int +maxParamOf = maximum . (negate 1 :) . cata alg + where alg :: SRTree [Int] -> [Int] + alg (Param ix) = [ix] + alg (Var _) = [] + alg (Const _) = [] + alg (Y _) = [] + alg (Uni _ t) = t + alg (Bin _ l r) = l <> r + +-- | Extra parameters appended by 'buildDistLoss' for a distribution +-- (Gaussian: log-noise; ROXY: 3 auxiliary parameters). +distExtraParams :: Distribution -> Int +distExtraParams Gaussian = 1 +distExtraParams ROXY = 3 +distExtraParams _ = 0 + -- | Calculate the profile likelihood of every parameter. -- restartLimit bounds recursive restarts when getProfileODE finds a better point mid-profile. -- For Bates, getProfile handles restarts internally and never returns Left. @@ -203,7 +222,8 @@ getAllProfiles ptype et theta stdErr estCIs alpha -- Defensive: if theta is too short for the EvalTree's distribution, -- return empty profiles instead of crashing (e.g. MSE loss with Gaussian dist) - | U.length theta < 2 = pure [] + | U.null theta = pure [] + | U.length theta < maxParamOf (ctTree et) + 1 + distExtraParams (ctDist et) = pure [] | otherwise = go 0 et theta stdErr estCIs where restartLimit = 5 :: Int @@ -475,7 +495,6 @@ cols = getCol ix thetas rawPairs = sortOnFirst taus cols monoPairs = enforceMonotonicTau rawPairs - _ = trace ("createSplines: raw=" ++ show (length rawPairs) ++ " mono=" ++ show (length monoPairs) ++ " head=" ++ show (take 3 monoPairs) ++ " last=" ++ show (reverse $ take 3 $ reverse monoPairs)) () tau2theta = genSplineFun monoPairs theta2tau = genSplineFun $ enforceMonotonicTheta optTh $ sortOnFirst cols taus diff --git a/src/Data/SRTree/Eval.hs b/src/Data/SRTree/Eval.hs --- a/src/Data/SRTree/Eval.hs +++ b/src/Data/SRTree/Eval.hs @@ -163,7 +163,7 @@ alg (Var (-1)) = Static y alg (Var (-2)) = Static yErr alg (Var i) = Static (dataset !! i) - alg (Param i) = Dynamic (\th !idx -> th `V.unsafeIndex` i) + alg (Param i) = Dynamic (\th !idx -> if i < V.length th then th `V.unsafeIndex` i else 0/0) -- 2. Univariate Functions alg (Uni f (Scl c)) = Scl (evalFun f c) @@ -215,7 +215,7 @@ alg (Const c) = Scl c alg (Var i) = Static (dataset !! i) -- Look at this! No more V.replicate. It just fetches the scalar directly. - alg (Param i) = Dynamic (\th !idx -> th `V.unsafeIndex` i) + alg (Param i) = Dynamic (\th !idx -> if i < V.length th then th `V.unsafeIndex` i else 0/0) alg (Y i) = undefined -- this shouldn't be called -- 2. Univariate Functions @@ -261,7 +261,7 @@ evalTree xss params = cata $ \case Var ix -> xss !! ix - Param ix -> replicateAs xss $ params V.! ix + Param ix -> replicateAs xss $ if ix < V.length params then params V.! ix else 0/0 Const c -> replicateAs xss c Y _ -> undefined Uni g t -> evalFun g t diff --git a/srtree.cabal b/srtree.cabal --- a/srtree.cabal +++ b/srtree.cabal @@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack name: srtree -version: 3.0.0.5 +version: 3.0.0.6 synopsis: A general library to work with Symbolic Regression expression trees. description: A Symbolic Regression Tree data structure to work with mathematical expressions with support to first order derivative and simplification; license: BSD3 @@ -60,28 +60,28 @@ other-modules: Paths_srtree build-depends: - async >=2.2 && <2.3 + async ==2.2.* , attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.19 && <5 - , binary >=0.8 && <0.9 + , binary ==0.8.* , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.9 , deepseq >=1.4 && <1.6 - , directory >=1.3 && <1.4 - , exceptions >=0.10 && <0.11 + , directory ==1.3.* + , exceptions ==0.10.* , filepath >=1.4.0.0 && <1.6 , hashable >=1.4 && <1.6 - , ieee754 >=0.8 && <0.9 + , ieee754 ==0.8.* , lens >=5.0 && <6 , mtl >=2.2 && <2.4 - , parallel >=3.2 && <3.4 + , parallel >=3.2 && <3.4 , primitive >=0.8 && <0.10 , random >=1.2 && <1.4 , split >=0.2.5 && <0.3 , statistics >=0.15 && <0.17 , time >=1.9 && <1.15 - , unordered-containers >=0.2 && <0.3 + , unordered-containers ==0.2.* , vector >=0.12 && <0.14 , zlib >=0.6.3 && <0.8 hs-source-dirs: @@ -99,29 +99,29 @@ apps/Bench ghc-options: -threaded -rtsopts -with-rtsopts=-N -O2 -fllvm -pgmlo opt-20 -pgmlc llc-20 -optlo-O3 -optlc-mcpu=native -mavx2 -mfma -fspec-constr -fmax-simplifier-iterations=20 -fexpose-all-unfoldings build-depends: - async >=2.2 && <2.3 + async ==2.2.* , attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.19 && <5 - , binary >=0.8 && <0.9 + , binary ==0.8.* , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.9 , criterion >=1.5 && <2 , deepseq >=1.4 && <1.6 - , directory >=1.3 && <1.4 - , exceptions >=0.10 && <0.11 + , directory ==1.3.* + , exceptions ==0.10.* , filepath >=1.4.0.0 && <1.6 , hashable >=1.4 && <1.6 - , ieee754 >=0.8 && <0.9 + , ieee754 ==0.8.* , lens >=5.0 && <6 , mtl >=2.2 && <2.4 - , parallel >=3.2 && <3.4 + , parallel >=3.2 && <3.4 , primitive >=0.8 && <0.10 , random >=1.2 && <1.4 , split >=0.2.5 && <0.3 , srtree , statistics >=0.15 && <0.17 - , unordered-containers >=0.2 && <0.3 + , unordered-containers ==0.2.* , vector >=0.12 && <0.14 , zlib >=0.6.3 && <0.8 default-language: Haskell2010 @@ -134,29 +134,29 @@ apps/BenchEqSat ghc-options: -threaded -rtsopts -with-rtsopts=-N -O2 build-depends: - async >=2.2 && <2.3 + async ==2.2.* , attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.19 && <5 - , binary >=0.8 && <0.9 + , binary ==0.8.* , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.9 , criterion >=1.5 && <2 , deepseq >=1.4 && <1.6 - , directory >=1.3 && <1.4 - , exceptions >=0.10 && <0.11 + , directory ==1.3.* + , exceptions ==0.10.* , filepath >=1.4.0.0 && <1.6 , hashable >=1.4 && <1.6 - , ieee754 >=0.8 && <0.9 + , ieee754 ==0.8.* , lens >=5.0 && <6 , mtl >=2.2 && <2.4 - , parallel >=3.2 && <3.4 + , parallel >=3.2 && <3.4 , primitive >=0.8 && <0.10 , random >=1.2 && <1.4 , split >=0.2.5 && <0.3 , srtree , statistics >=0.15 && <0.17 - , unordered-containers >=0.2 && <0.3 + , unordered-containers ==0.2.* , vector >=0.12 && <0.14 , zlib >=0.6.3 && <0.8 default-language: Haskell2010 @@ -169,29 +169,29 @@ apps/Report ghc-options: -threaded -rtsopts -with-rtsopts=-N -O2 build-depends: - async >=2.2 && <2.3 + async ==2.2.* , attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.19 && <5 - , binary >=0.8 && <0.9 + , binary ==0.8.* , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.9 , deepseq >=1.4 && <1.6 - , directory >=1.3 && <1.4 - , exceptions >=0.10 && <0.11 + , directory ==1.3.* + , exceptions ==0.10.* , filepath >=1.4.0.0 && <1.6 , hashable >=1.4 && <1.6 - , ieee754 >=0.8 && <0.9 + , ieee754 ==0.8.* , lens >=5.0 && <6 , mtl >=2.2 && <2.4 , optparse-applicative >=0.16 && <0.20 - , parallel >=3.2 && <3.4 + , parallel >=3.2 && <3.4 , primitive >=0.8 && <0.10 , random >=1.2 && <1.4 , split >=0.2.5 && <0.3 , srtree , statistics >=0.15 && <0.17 - , unordered-containers >=0.2 && <0.3 + , unordered-containers ==0.2.* , vector >=0.12 && <0.14 , zlib >=0.6.3 && <0.8 default-language: Haskell2010 @@ -205,14 +205,14 @@ ghc-options: -threaded -rtsopts -with-rtsopts=-N -O2 build-depends: base >=4.19 && <5 - , srtree , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.9 - , vector >=0.12 && <0.14 - , exceptions >=0.10 && <0.11 - , directory >=1.3 && <1.4 + , directory ==1.3.* + , exceptions ==0.10.* , filepath >=1.4.0.0 && <1.6 , random >=1.2 && <1.4 + , srtree + , vector >=0.12 && <0.14 default-language: Haskell2010 test-suite srtree-test @@ -227,29 +227,29 @@ test ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends: - HUnit >=1.6 && <1.7 - , async >=2.2 && <2.3 + HUnit ==1.6.* + , async ==2.2.* , attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.19 && <5 - , binary >=0.8 && <0.9 + , binary ==0.8.* , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.9 , deepseq >=1.4 && <1.6 - , directory >=1.3 && <1.4 - , exceptions >=0.10 && <0.11 + , directory ==1.3.* + , exceptions ==0.10.* , filepath >=1.4.0.0 && <1.6 , hashable >=1.4 && <1.6 - , ieee754 >=0.8 && <0.9 + , ieee754 ==0.8.* , lens >=5.0 && <6 , mtl >=2.2 && <2.4 - , parallel >=3.2 && <3.4 + , parallel >=3.2 && <3.4 , primitive >=0.8 && <0.10 , random >=1.2 && <1.4 , split >=0.2.5 && <0.3 , srtree , statistics >=0.15 && <0.17 - , unordered-containers >=0.2 && <0.3 + , unordered-containers ==0.2.* , vector >=0.12 && <0.14 , zlib >=0.6.3 && <0.8 default-language: Haskell2010