srtree 3.0.0.5 → 3.0.0.6

raw patch · 6 files changed

+109/−51 lines, 6 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• ChangeLog.md +15/−0
• src/Algorithm/SRTree/AD/CompiledAD.hs +1/−0
• src/Algorithm/SRTree/AD/Unboxed.hs +27/−4
• src/Algorithm/SRTree/ConfidenceIntervals.hs +22/−3
• src/Data/SRTree/Eval.hs +3/−3
• srtree.cabal +41/−41

@@ -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`):

@@ -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)   }

@@ -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

@@ -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 

@@ -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

@@ -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