diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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` + 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 @@ -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) 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 @@ -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. 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.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 diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -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