srtree 3.0.0.6 → 3.0.0.7
raw patch · 5 files changed
+65/−15 lines, 5 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- ChangeLog.md +18/−0
- src/Algorithm/SRTree/AD/Unboxed.hs +1/−1
- src/Algorithm/SRTree/ConfidenceIntervals.hs +22/−13
- srtree.cabal +1/−1
- test/Spec.hs +23/−0
ChangeLog.md view
@@ -1,5 +1,23 @@ # Changelog for srtree +## 3.0.0.7++- **Bound profile-likelihood walk restarts** (`ConfidenceIntervals`): the+ per-parameter profile walk restarted from every marginally better optimum it+ found and `getAllProfiles`' restart limit was a no-op (both branches+ identical), so profiling weakly-identified parameters could run effectively+ forever. `getProfile` now caps internal restarts (5) and `getAllProfiles`+ gives up after its restart limit instead of re-profiling indefinitely.+- **`PowerAbs` static-base gradient fix** (`Algorithm.SRTree.AD.Unboxed`): in+ the multi-threaded `evalGradVec` backward pass, the static-left/dynamic-right+ case of opcode 5 (`PowerAbs`, i.e. `|l|^r`) passed the static left child's+ buffer offset to `bwdST` instead of the dynamic right child's. The computed+ exponent gradient was written to the static node's unused adjoint slot, so+ parameters in the exponent (e.g. `|x0|^(t0*x0)`) received no gradient. This+ silently degraded NLopt fits for such expressions (wrong eggp fitness and+ refits); the gradient is now accumulated on the correct right child.+ Added a finite-difference regression test for this case.+ ## 3.0.0.6 - **AD theta-length guards** (`Algorithm.SRTree.AD.CompiledAD` +
src/Algorithm/SRTree/AD/Unboxed.hs view
@@ -684,7 +684,7 @@ (2, 3) -> bwdST nb stSrc fwd adj (\dx _ gy fg -> (dx / gy, dx * (negate fg / gy))) kb rb (2, 4) -> bwdST nb stSrc fwd adj (\dx fx gy fg -> (fixNaN (dx * gy * fg / fx), fixNaN (dx * fg * log fx))) kb rb (2, 5) -> bwdST nb stSrc fwd adj (\dx fx gy fg ->- let v2 = abs fx in (fixNaN (dx * (fx * gy) * fg / (v2 * v2)), fixNaN (dx * fg * log (abs fx)))) kb lb+ let v2 = abs fx in (fixNaN (dx * (fx * gy) * fg / (v2 * v2)), fixNaN (dx * fg * log (abs fx)))) kb rb (2, 6) -> bwdST nb stSrc fwd adj (\dx fx gy _ -> let dxl = dx * (recip . sqrt . (+1) . (^(2::Int))) gy dxy = fx * gy * dxl ^ (3::Int)
src/Algorithm/SRTree/ConfidenceIntervals.hs view
@@ -228,8 +228,11 @@ where restartLimit = 5 :: Int + -- Give up after too many restarts instead of re-profiling forever (a+ -- profile walk that keeps finding marginally better optima would not+ -- terminate otherwise). go restarts et' theta' stdErr' estCIs'- | restarts >= restartLimit = profileAll restarts et' theta' stdErr' estCIs'+ | restarts >= restartLimit = pure [] | otherwise = profileAll restarts et' theta' stdErr' estCIs' profileAll restarts et' theta' stdErr' estCIs' = go' 0 []@@ -289,9 +292,9 @@ nllOptRef <- newIORef nll_opt0 thetaOptRef <- newIORef theta_opt0 - negDelta <- go kmax (-stdErr_i / 8) 0 1 mempty nllOptRef thetaOptRef+ negDelta <- go restartLimit kmax (-stdErr_i / 8) 0 1 mempty nllOptRef thetaOptRef thetaOpt1 <- readIORef thetaOptRef- posDelta <- go kmax (stdErr_i / 8) 0 1 ([0], [thetaOpt1]) nllOptRef thetaOptRef+ posDelta <- go restartLimit kmax (stdErr_i / 8) 0 1 ([0], [thetaOpt1]) nllOptRef thetaOptRef thetaOpt2 <- readIORef thetaOptRef let optTh' = thetaOpt2 U.! ix@@ -301,13 +304,14 @@ (tau2theta, theta2tau) = createSplines taus thetas stdErr_i tau_max ix optTh' pure $ ProfileT taus thetas optTh' tau2theta theta2tau where- kmax = 500+ kmax = 500+ restartLimit = 5 :: Int nll_opt0 = ctNLL et theta_opt0 theta_opt0 = ctOptimizer et theta minimizer = ctOptimizerFixed et ix - go 0 _delta _t _inv_slope acc _nllRef _thetaRef = pure acc- go k delta t inv_slope acc@(taus, thetas) nllOptRef thetaOptRef = do+ go _ 0 _delta _t _inv_slope acc _nllRef _thetaRef = pure acc+ go restarts k delta t inv_slope acc@(taus, thetas) nllOptRef thetaOptRef = do nllOpt <- readIORef nllOptRef thetaOpt <- readIORef thetaOptRef let t_delta = (thetaOpt U.! ix) + delta * (t + inv_slope)@@ -330,14 +334,19 @@ then acc else (tau:taus, theta_t:thetas) if | not validDelta || isNaN inv_slope -> pure acc- | nll_cond < nllOpt - 1e-6 * abs nllOpt -> do- -- Better optimum found: update references and restart walk- -- from the new MLE, discarding previously collected points.- writeIORef nllOptRef nll_cond- writeIORef thetaOptRef theta_t- go kmax delta 0 1 mempty nllOptRef thetaOptRef+ | nll_cond < nllOpt - 1e-6 * abs nllOpt ->+ if restarts <= 0+ -- Restart budget exhausted: keep the point found so far+ -- instead of walking/restarting indefinitely.+ then pure acc'+ else do+ -- Better optimum found: update references and restart walk+ -- from the new MLE, discarding previously collected points.+ writeIORef nllOptRef nll_cond+ writeIORef thetaOptRef theta_t+ go (restarts - 1) kmax delta 0 1 mempty nllOptRef thetaOptRef | abs tau > tau_max -> pure acc'- | otherwise -> go (k-1) delta (t + inv_slope) inv_slope' acc' nllOptRef thetaOptRef+ | otherwise -> go restarts (k-1) delta (t + inv_slope) inv_slope' acc' nllOptRef thetaOptRef -- Based on https://insysbio.github.io/LikelihoodProfiler.jl/latest/ -- Borisov, Ivan, and Evgeny Metelkin. "Confidence intervals by constrained optimization—An algorithm and software package for practical identifiability analysis in systems biology." PLOS Computational Biology 16.12 (2020): e1008495.
srtree.cabal view
@@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack name: srtree-version: 3.0.0.6+version: 3.0.0.7 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
test/Spec.hs view
@@ -62,6 +62,28 @@ assertBool "evalGrad gradient != finite diff" (and (zipWith (\a b -> abs (a - b) < 1e-4) (VS.toList g0) gfd)) +-- Regression for the PowerAbs (opcode 5) static-left/dynamic-right backward+-- case: the exponent parameter's gradient used to be written to the static+-- left child's unused adjoint slot, silently dropping it from evalGradVec.+test_grad_powerabs_static_base :: Test+test_grad_powerabs_static_base = TestCase $ do+ let xss = [VU.fromList [1.5, 2.0, 3.0, 4.0]]+ y = VU.fromList [1.0, 2.0, 3.0, 4.0]+ tree = Fix (Bin PowerAbs (var 0) (param 0 * var 0)) -- |x0| ** (t0*x0)+ theta = VS.fromList [0.5]+ ct = compileTree xss y Nothing tree+ (f0, g0) = evalGrad ct theta+ (f1, g1) = evalGradVec ct theta+ h = 1e-6+ (fp, _) = evalGradVec ct (VS.fromList [0.5 + h])+ (fm, _) = evalGradVec ct (VS.fromList [0.5 - h])+ gfd = (fp - fm) / (2 * h)+ assertBool "PowerAbs ST objective mismatch" (abs (f1 - f0) < 1e-9 * max 1 (abs f0))+ assertBool ("PowerAbs ST evalGradVec gradient mismatch: " ++ show (VS.toList g1, gfd))+ (abs (g1 VS.! 0 - gfd) < 1e-4 * max 1 (abs gfd))+ assertBool ("PowerAbs ST evalGrad gradient mismatch: " ++ show (VS.toList g0, gfd))+ (abs (g0 VS.! 0 - gfd) < 1e-4 * max 1 (abs gfd))+ test_benchgrad :: Test test_benchgrad = TestCase $ do let genTerm = do coin <- tossBiased 0.4@@ -107,6 +129,7 @@ counts <- runTestTT $ TestList [ t1 , t2+ , TestLabel "grad-powerabs-st" test_grad_powerabs_static_base , TestLabel "benchgrad" test_benchgrad , TestLabel "eqsat" EqSatTests.tests , TestLabel "store" StoreTests.tests