diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,27 @@ # Changelog for srtree +## 3.0.0.3 + +- **Profile-likelihood CI overhaul** (`ConfidenceIntervals`): + - Fixed `paramCI` to use F-distribution with 1 numerator df (was k df), + matching per-parameter profiling semantics + - `getAllProfiles`: added restart limit (5), defensive guard for short theta, + auto-compute Laplace CIs when estCIs is empty, and `recomputeStdErr` helper + - `getProfile`: fixed `tau_max` to use 1 df; added `LeastSquares` profile + statistic `n * log(MSE(t)/MSE(opt))` (was raw `2*(MSE(t)-MSE(opt))` which + differs by factor n/(2*MSE)); increased step limit 300→500; added fallback + for small gradient in `inv_slope'`; guarded `nll_cond < nll_opt` with epsilon + - `getProfileCnstr`/`getEndPoint`: replaced NELDERMEAD+AugLag with robust + bisection on the profiled NLL; wider search bounds (50× se); NaN guards + - `getStatsFromModel`: for `LeastSquares`, scale covariance by MSE (was + unscaled); use `max 0` for sqrt of diagonal to avoid numerical NaN + - `createSplines`: enforce monotonicity on (tau,θ) and (θ,tau) pairs to + prevent spline extrapolation garbage; accepts `optTh` parameter +- **Test CI executable**: new `test-ci` app in `apps/TestCI` for investigating + profile-likelihood backends on real datasets +- **CI tests**: new `CITests` module covering monotonicity, spline, and + negative-tau regression tests + ## 3.0.0.2 - Added parser for NeoGP.jl diff --git a/apps/TestCI/Main.hs b/apps/TestCI/Main.hs new file mode 100644 --- /dev/null +++ b/apps/TestCI/Main.hs @@ -0,0 +1,345 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Main where + +import qualified Data.ByteString.Char8 as B +import Data.SRTree +import Data.SRTree.Print (showExpr) +import Data.SRTree.Eval (Target, Columns, compile) +import Data.SRTree.Datasets (loadDataset) +import Data.SRTree.Recursion (Fix(..)) +import Algorithm.SRTree.ConfidenceIntervals +import Algorithm.SRTree.Compile (compileTree, EvalTree(..)) +import Algorithm.SRTree.NonlinearOpt (minimizeNLL, minimizeNLLWith, compileLossAndGrad) +import Algorithm.SRTree.Likelihoods (Distribution(..), Loss(..)) +import Algorithm.SRTree.AD (ADBackEnd(..)) +import Algorithm.SRTree.AD.Unboxed (setMTPopParallel) +import Numeric.Optimization.NLOPT (LocalAlgorithm(..)) +import Text.ParseSR (parseSR, SRAlgs(..)) +import qualified Data.Vector.Unboxed as VU +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Generic as G +import Data.List (intercalate, foldl', maximumBy, isSuffixOf) +import Data.Maybe (fromMaybe) +import System.IO (hFlush, stdout, hPutStrLn, stderr) +import System.Environment (getArgs) +import System.Random (randomRIO) +import Control.Exception (catch, SomeException, evaluate) +import Control.Monad (forM_, when) +import Data.Either (isRight) + +-- | A test case: expression string, description +data TestCase = TestCase + { tcName :: String + , tcExpr :: String + , tcDesc :: String + } deriving (Show) + +-- | All test cases +testCases :: [TestCase] +testCases = + [ TestCase + { tcName = "exp-param" + , tcExpr = "Exp(t0 + (x1 * (t1 * ((x0 * x0) + t2))))" + , tcDesc = "Parameterized exponential (similar shape to AGENTS.md example)" + } + , TestCase + { tcName = "linear" + , tcExpr = "t0 * x0 + t1" + , tcDesc = "Simple linear model (well-conditioned)" + } + , TestCase + { tcName = "quadratic" + , tcExpr = "t0 * x0 * x0 + t1 * x0 + t2" + , tcDesc = "Quadratic polynomial" + } + , TestCase + { tcName = "rational" + , tcExpr = "t0 / (t1 + x0)" + , tcDesc = "Rational function (steep near pole)" + } + , TestCase + { tcName = "sine" + , tcExpr = "t0 * sin(t1 * x0 + t2)" + , tcDesc = "Sinusoidal model" + } + , TestCase + { tcName = "product" + , tcExpr = "t0 * x0 * x1 + t1" + , tcDesc = "Two-variable product" + } + , TestCase + { tcName = "exp-linear" + , tcExpr = "Exp(t0 * x0 + t1)" + , tcDesc = "Exponential of linear (simpler than deep exp)" + } + , TestCase + { tcName = "power" + , tcExpr = "t0 * x0 ** t1" + , tcDesc = "Power law" + } + , TestCase + { tcName = "linear-mse" + , tcExpr = "t0 * x0 + t1" + , tcDesc = "Linear model fitted with MSE (Bates 1985 original use case)" + } + ] + +main :: IO () +main = do + args <- getArgs + let filterName = if null args then Nothing else Just (head args) + + putStrLn "========================================================================" + putStrLn " Profile Likelihood CI Backend Investigation" + putStrLn "========================================================================" + putStrLn "" + + -- Load dataset + let dataSpec = "../eggp/gaussian_train.csv:::y_noise_02:x1,x2" + putStrLn $ "Loading dataset: " ++ dataSpec + hFlush stdout + ((xTr, yTr, _xVal, _yVal), (mYErr, _), _varnames, _target) <- + loadDataset dataSpec True `catch` (\(e :: SomeException) -> do + putStrLn $ "ERROR loading dataset: " ++ show e + error "Failed to load dataset") + + let nSamples = VU.length yTr + putStrLn $ " Samples: " ++ show nSamples + putStrLn $ " Features: " ++ show (length xTr) + putStrLn "" + + let cases = case filterName of + Nothing -> testCases + Just name -> filter (\tc -> tcName tc == name) testCases + + mapM_ (runTestCase xTr yTr mYErr nSamples) cases + + putStrLn "" + putStrLn "========================================================================" + putStrLn " Summary" + putStrLn "========================================================================" + putStrLn "" + putStrLn "Key observations:" + putStrLn " - Laplace: uses Hessian inverse; fast but may be inaccurate for nonlinear models" + putStrLn " - Bates: classical profile walk; accurate but slow" + putStrLn " - ODE: Chen & Jennrich ODE-based profile; fast and accurate" + putStrLn " - Constrained: bisection on re-optimized endpoints; fast but may fail" + putStrLn "" + putStrLn "Issues to investigate:" + putStrLn " 1. Constrained backend NaN on steep exponential expressions" + putStrLn " 2. tau_max' threshold correctness for NLL Gaussian" + putStrLn " 3. Nelder-Mead convergence in augmented Lagrangian" + putStrLn " 4. One-sided CI failures" + +-- | Run a single test case through all backends +runTestCase :: [VU.Vector Double] -> VU.Vector Double -> Maybe (VU.Vector Double) -> Int -> TestCase -> IO () +runTestCase xTr yTr mYErr nSamples tc = do + putStrLn "------------------------------------------------------------------------" + putStrLn $ "Test: " ++ tcName tc + putStrLn $ " Description: " ++ tcDesc tc + putStrLn $ " Expression: " ++ tcExpr tc + putStrLn "" + + -- Parse expression: convert String to ByteString for parseSR + let parsed = parseSR TIR (B.pack "x0,x1") False (B.pack (tcExpr tc)) + case parsed of + Left err -> putStrLn $ " PARSE ERROR: " ++ err + Right rawTree -> do + let tree = relabelParams rawTree + let nParams = countParamsUniq tree + putStrLn $ " Parsed tree: " ++ showExpr tree + putStrLn $ " Unique params: " ++ show nParams + + if nParams == 0 + then putStrLn " SKIP: no parameters to profile" + else do + -- Detect MSE cases (name ends with "-mse") + let useMSE = "-mse" `isSuffixOf` tcName tc + dist = if useMSE then LeastSquares else Gaussian + totalParams = if useMSE then nParams else nParams + 1 -- +1 for sigma when Gaussian + + putStrLn $ " Loss: " ++ (if useMSE then "MSE (LeastSquares)" else "NLL Gaussian") + putStrLn $ " Total params: " ++ show totalParams + putStrLn "" + + -- Fit with multiple restarts + putStrLn " Fitting..." + hFlush stdout + setMTPopParallel True + results <- fitMultipleRestarts dist mYErr xTr yTr tree totalParams 5 + setMTPopParallel False + + let (bestNLL, bestTheta) = maximumBy (\(a,_) (b,_) -> compare a b) results + theta_opt = bestTheta + negNLL = negate bestNLL + + putStrLn $ " Best loss: " ++ show negNLL + putStrLn $ " Theta: " ++ show (VU.toList theta_opt) + putStrLn "" + + -- Compile the EvalTree for CI computation + let et = compileTree dist xTr yTr mYErr tree + + -- Verify the optimizer agrees + let theta_verify = ctOptimizer et theta_opt + nll_verify = ctNLL et theta_verify + putStrLn $ " Verified loss (via EvalTree): " ++ show nll_verify + + -- Compute standard errors from Hessian + let stats = getStatsFromModel dist mYErr xTr yTr tree theta_opt + stdErrs = _stdErr stats + putStrLn $ " Std errors (Hessian): " ++ show (VU.toList stdErrs) + putStrLn "" + + let paramNames = [ "t" ++ show i | i <- [0 .. nParams - 1] ] + ++ if useMSE then [] else ["sigma"] + + -- ---- LAPLACE ---- + putStrLn " === LAPLACE ===" + let laplaceCI = paramCI (Laplace stats) nSamples theta_opt 0.05 + putStrLn $ " 95% CIs:" + putStrLn $ " " ++ showCIList (zip paramNames laplaceCI) + putStrLn "" + + -- ---- BATES (profile walk) ---- + putStrLn " === BATES (profile walk) ===" + catch (do + let estCIs = laplaceCI + profiles_bates = getAllProfiles Bates et theta_opt stdErrs estCIs 0.05 + batesCI = paramCI (Profile stats profiles_bates) nSamples theta_opt 0.05 + putStrLn $ " 95% CIs:" + putStrLn $ " " ++ showCIList (zip paramNames batesCI) + putStrLn $ " Widths: " ++ show (map (\(CI _ l h) -> h - l) batesCI) + -- Debug: show first/last profile points + forM_ (zip [0::Int ..] profiles_bates) $ \(ix, prof) -> do + let taus = _taus prof + cols = _thetas prof + nT = VU.length taus + putStrLn $ " Profile t" ++ show ix ++ ": " ++ show nT ++ " points" + when (nT > 0) $ do + let firstTau = taus VU.! 0 + lastTau = taus VU.! (nT - 1) + firstTh = (cols !! ix) VU.! 0 + lastTh = (cols !! ix) VU.! (nT - 1) + optTh = theta_opt VU.! ix + putStrLn $ " tau=[" ++ show firstTau ++ ", " ++ show lastTau ++ "]" + putStrLn $ " theta=[" ++ show firstTh ++ ", " ++ show lastTh ++ "] opt=" ++ show optTh + ) (\(e :: SomeException) -> putStrLn $ " ERROR: " ++ show e) + putStrLn "" + + -- ---- ODE (Chen & Jennrich) ---- + putStrLn " === ODE (Chen & Jennrich) ===" + catch (do + let estCIs = laplaceCI + profiles_ode = getAllProfiles ODE et theta_opt stdErrs estCIs 0.05 + odeCI = paramCI (Profile stats profiles_ode) nSamples theta_opt 0.05 + putStrLn $ " 95% CIs:" + putStrLn $ " " ++ showCIList (zip paramNames odeCI) + putStrLn $ " Widths: " ++ show (map (\(CI _ l h) -> h - l) odeCI) + ) (\(e :: SomeException) -> putStrLn $ " ERROR: " ++ show e) + putStrLn "" + + -- ---- CONSTRAINED ---- + putStrLn " === CONSTRAINED (bisection) ===" + catch (do + let profiles_cnstr = getAllProfiles Constrained et theta_opt stdErrs [] 0.05 + cnstrCI = paramCI (Profile stats profiles_cnstr) nSamples theta_opt 0.05 + putStrLn $ " 95% CIs:" + putStrLn $ " " ++ showCIList (zip paramNames cnstrCI) + putStrLn $ " Widths: " ++ show (map (\(CI _ l h) -> h - l) cnstrCI) + -- Check for NaN + let hasNaN = any (\(CI _ l h) -> isNaN l || isNaN h) cnstrCI + when hasNaN $ putStrLn $ " *** WARNING: NaN detected in Constrained CI ***" + ) (\(e :: SomeException) -> putStrLn $ " ERROR: " ++ show e) + putStrLn "" + + -- ---- Detailed profiling of Constrained for the problematic case ---- + when (nParams >= 2) $ do + putStrLn " === DETAILED CONSTRAINED INVESTIGATION ===" + investigateConstrained et theta_opt stdErrs totalParams + putStrLn "" + +-- | Fit with multiple random restarts +fitMultipleRestarts :: Distribution -> Maybe (VU.Vector Double) -> [VU.Vector Double] -> VU.Vector Double + -> Fix SRTree -> Int -> Int -> IO [(Double, VU.Vector Double)] +fitMultipleRestarts dist mYErr xTr yTr tree nParams nRep = do + let funAndGrad = compileLossAndGrad MultiThread (NLL dist) mYErr xTr yTr tree + runRestart = do + theta0 <- VU.replicateM nParams (randomRIO (-2, 2)) + let (theta, lossVal, _) = minimizeNLLWith funAndGrad TNEWTON 200 theta0 + pure (negate lossVal, theta) + results <- sequence [ runRestart | _ <- [1..nRep] ] + -- Also try from zeros + let theta0_zero = VU.replicate nParams 0.0 + (theta_zero, loss_zero, _) = minimizeNLLWith funAndGrad TNEWTON 200 theta0_zero + pure $ (negate loss_zero, theta_zero) : results + +-- | Show a list of CIs with parameter names +showCIList :: [(String, CI)] -> String +showCIList = intercalate "\n " . map (\(name, CI est lo hi) -> + name ++ ": " ++ showF lo ++ " <= " ++ showF est ++ " <= " ++ showF hi) + where showF x + | isNaN x = "NaN" + | isInfinite x = if x > 0 then "+Inf" else "-Inf" + | otherwise = show (fromIntegral (round (x * 1e4) :: Int) / 1e4 :: Double) + +-- | Detailed investigation of the Constrained backend +investigateConstrained :: EvalTree -> VU.Vector Double -> VU.Vector Double -> Int -> IO () +investigateConstrained et theta_opt stdErrs nParams = do + let nll_opt = ctNLL et theta_opt + n = ctRows et + k = VU.length theta_opt + chi2_1 = 3.841 -- chi2 quantile for 1 df at 0.95 + + putStrLn $ " nll_opt = " ++ show nll_opt + putStrLn $ " n = " ++ show n ++ ", k = " ++ show k + putStrLn $ " chi2_1(0.95) = " ++ show chi2_1 + putStrLn "" + + -- Corrected tau_max' calculation + let tau_max' = chi2_1 / 2 + tau_max_old = nll_opt * chi2_1 / fromIntegral n -- OLD (buggy) + putStrLn $ " Corrected tau_max' (chi2_1/2) = " ++ show tau_max' + putStrLn $ " OLD tau_max' (nll_opt * chi2_1 / n) = " ++ show tau_max_old ++ " (WRONG)" + putStrLn "" + + -- Test each parameter + forM_ [0 .. nParams - 1] $ \ix -> do + putStrLn $ " Parameter t" ++ show ix ++ ":" + putStrLn $ " MLE = " ++ show (theta_opt VU.! ix) + putStrLn $ " StdErr = " ++ show (stdErrs VU.! ix) + + -- Test getEndPoint directly + let getPoint isLeft = getEndPoint et theta_opt tau_max' (stdErrs VU.! ix) ix isLeft + catch (do + let leftPt = getPoint True + rightPt = getPoint False + putStrLn $ " Left endpoint = " ++ show leftPt + putStrLn $ " Right endpoint = " ++ show rightPt + when (isNaN leftPt || isNaN rightPt) $ + putStrLn $ " *** NaN detected! ***" + when (leftPt > rightPt) $ + putStrLn $ " *** Left > Right: reversed interval! ***" + ) (\(e :: SomeException) -> putStrLn $ " ERROR in getEndPoint: " ++ show e) + + -- Test the profiling function (fix ix, re-optimize others) + putStrLn $ " Testing ctOptimizerFixed..." + catch (do + let delta = stdErrs VU.! ix * 0.5 + theta_left = VU.generate nParams (\j -> if j == ix then (theta_opt VU.! ix) - delta else theta_opt VU.! j) + theta_right = VU.generate nParams (\j -> if j == ix then (theta_opt VU.! ix) + delta else theta_opt VU.! j) + reopt_left = ctOptimizerFixed et ix theta_left + reopt_right = ctOptimizerFixed et ix theta_right + nll_left = ctNLL et reopt_left + nll_right = ctNLL et reopt_right + putStrLn $ " theta_left (fixed at " ++ show (theta_opt VU.! ix - delta) ++ ") -> reopt NLL = " ++ show nll_left + putStrLn $ " theta_right (fixed at " ++ show (theta_opt VU.! ix + delta) ++ ") -> reopt NLL = " ++ show nll_right + putStrLn $ " NLL increase left: " ++ show (nll_left - nll_opt) + putStrLn $ " NLL increase right: " ++ show (nll_right - nll_opt) + when (isNaN nll_left || isNaN nll_right) $ + putStrLn $ " *** NaN in re-optimized NLL! ***" + ) (\(e :: SomeException) -> putStrLn $ " ERROR in ctOptimizerFixed: " ++ show e) + putStrLn "" diff --git a/src/Algorithm/EqSat/SearchSR.hs b/src/Algorithm/EqSat/SearchSR.hs --- a/src/Algorithm/EqSat/SearchSR.hs +++ b/src/Algorithm/EqSat/SearchSR.hs @@ -213,10 +213,9 @@ printExprFun 0 bec Nothing -> pure () ---paretoFront :: Int -> (Int -> EClassId -> RndEGraph ()) -> RndEGraph () +paretoFront :: (Fix SRTree -> RndEGraph (Double, [Target])) -> Int -> (Int -> EClassId -> RndEGraph b) -> RndEGraph [b] paretoFront fitFun maxSize printExprFun = go 1 0 (-(1.0/0.0)) where - go :: Int -> Int -> Double -> RndEGraph [[String]] go n ix f | n > maxSize = pure [] | otherwise = do diff --git a/src/Algorithm/EqSat/Simplify.hs b/src/Algorithm/EqSat/Simplify.hs --- a/src/Algorithm/EqSat/Simplify.hs +++ b/src/Algorithm/EqSat/Simplify.hs @@ -262,6 +262,7 @@ myCost (Var _) = 1 myCost (Const _) = 3 myCost (Param _) = 3 +myCost (Y _) = 1 myCost (Bin op l r) = 2 + l + r myCost (Uni _ t) = 3 + t 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 @@ -48,7 +48,6 @@ import qualified Data.Vector.Unboxed.Mutable as VUM import qualified Data.Vector as VB import qualified Data.Vector.Mutable as VMB -import Debug.Trace (trace, traceShow) import qualified Data.IntMap.Strict as IntMap import Data.List ( foldl', foldl1' ) import Data.Maybe (isJust, fromMaybe) 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 @@ -1,4 +1,4 @@ -{-# language ViewPatterns, ScopedTypeVariables, MultiWayIf, FlexibleContexts #-} +{-# language ViewPatterns, ScopedTypeVariables, MultiWayIf, FlexibleContexts, BangPatterns #-} ------------------------------------------------------------------------------- -- | -- Module : Algorithm.SRTree.ConfidenceIntervals @@ -29,7 +29,6 @@ import Numeric.Optimization.NLOPT import System.IO.Unsafe ( unsafePerformIO ) import Control.Monad.Catch ( catch, SomeException ) - import Debug.Trace ( trace ) -- | profile likelihood algorithms: Bates (classical), ODE (faster), Constrained (fastest) @@ -87,9 +86,10 @@ paramCI (Profile stats profiles) nSamples _ alpha = zipWith3 CI theta lows highs where - -- for the profile likelihood we use the square root of the F-distribution with (1-alpha) + -- for the profile likelihood we use the square root of the F-distribution + -- with 1 numerator df (each parameter is profiled individually) k = length theta - t = sqrt $ quantile (fDistribution k (fromIntegral $ nSamples - k)) (1 - alpha) + t = sqrt $ quantile (fDistribution 1 (fromIntegral $ nSamples - k)) (1 - alpha) stdErr = _stdErr stats lows = map (`_tau2theta` (-t)) profiles highs = map (`_tau2theta` t) profiles @@ -182,43 +182,97 @@ Left g -> Left $ g . invleft op vl Right vr -> Right $ evalOp op vl vr +-- | Recompute standard errors from the Hessian at a given theta. +-- Used when a profile walk restarts from a new optimum. +recomputeStdErr :: EvalTree -> Target -> Target +recomputeStdErr et t = stdErr + where + k = U.length t + ident = fromRowMajor k k (U.generate (k * k) (\ix -> let (i, j) = ix `divMod` k in if i == j then 1.0 else 0.0)) + hess = ctHessianNLL et t + cov = unsafePerformIO $ catch (invChol hess) (\(_ :: SomeException) -> pure ident) + covMat = toRowMajor cov + stdErr = U.generate k (\ix -> sqrt $ abs (covMat U.! (ix * k + ix))) + -- calculate the profile likelihood of every parameter +-- restartLimit bounds recursive restarts when the optimizer finds a better point mid-profile getAllProfiles :: PType -> EvalTree -> Target -> Target -> [CI] -> Double -> [ProfileT] -getAllProfiles ptype et theta stdErr estCIs alpha = getAll 0 [] +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 = [] + | otherwise = go 0 et theta stdErr estCIs where - k = U.length theta - n = ctRows et - tau_max = sqrt $ quantile (fDistribution k (n - k)) (1 - 0.01) - tau_max' = sqrt $ quantile (fDistribution k (n - k)) (1 - alpha) + restartLimit = 5 :: Int - profFun ix = case ptype of - Bates -> getProfile et theta (stdErr U.! ix) tau_max ix - ODE -> getProfileODE et theta (stdErr U.! ix) (estCIs !! ix) tau_max ix - Constrained -> getProfileCnstr et theta (stdErr U.! ix) tau_max' ix + go restarts et' theta' stdErr' estCIs' + | restarts >= restartLimit = profileAll restarts et' theta' stdErr' estCIs' + | otherwise = profileAll restarts et' theta' stdErr' estCIs' - getAll ix acc | ix == k = acc - | ix == k-1 && ptype == Constrained && ctDist et == Gaussian = case getProfileODE et theta (stdErr U.! ix) (estCIs !! ix) tau_max ix of - Left t -> getAllProfiles ptype et t stdErr estCIs alpha - Right p -> getAll (ix + 1) (acc <> [p]) - | otherwise = case profFun ix of - Left t -> getAllProfiles ptype et t stdErr estCIs alpha - Right p -> getAll (ix + 1) (acc <> [p]) + profileAll restarts et' theta' stdErr' estCIs' = go' 0 [] + where + k = U.length theta' + n = ctRows et' + -- For profiling a single parameter, the threshold is chi2_1 (1 df), + -- not chi2_k (k df). The profile likelihood ratio for ONE parameter + -- follows chi2_1 under H0. + tau_max = sqrt $ quantile (fDistribution 1 (n - k)) (1 - 0.01) + nll_opt = ctNLL et' (ctOptimizer et' theta') + chi2_1 = quantile (fDistribution 1 (n - k)) (1 - alpha) + -- Profile likelihood CI: 2*(L(theta_hat) - L(theta)) <= chi2_1 + -- => ctNLL(theta) <= ctNLL(theta_hat) + chi2_1/2 + -- So tau_max for the constrained method = chi2_1/2 + tau_max' = chi2_1 / 2 + -- If estCIs is empty, compute Laplace CIs as initial estimates + -- (needed by ODE fallback for the last Gaussian parameter) + estCIs'' = if null estCIs' + then let ident = U.generate (k * k) (\ix -> let (i, j) = ix `divMod` k in if i == j then 1.0 else 0.0) + hess = ctHessianNLL et' theta' + cov = unsafePerformIO $ catch (invChol hess) (\(_ :: SomeException) -> pure (fromRowMajor k k ident)) + covMat = toRowMajor cov + se = U.generate k (\ix -> sqrt $ abs (covMat U.! (ix * k + ix))) + tVal = quantile (studentT . fromIntegral $ n - k) (1 - alpha / 2.0) + in map (\ix -> CI (theta' U.! ix) ((theta' U.! ix) - tVal * (se U.! ix)) ((theta' U.! ix) + tVal * (se U.! ix))) [0..k-1] + else estCIs' + + profFun ix = case ptype of + Bates -> getProfile et' theta' (stdErr' U.! ix) tau_max ix + ODE -> getProfileODE et' theta' (stdErr' U.! ix) (estCIs'' !! ix) tau_max ix + Constrained -> getProfileCnstr et' theta' (stdErr' U.! ix) tau_max' ix + + go' ix acc | ix == k = acc + go' ix acc + | ix == k-1 && ptype == Constrained && ctDist et' == Gaussian = + case getProfileODE et' theta' (stdErr' U.! ix) (estCIs'' !! ix) tau_max ix of + Left t -> let tOpt = ctOptimizer et' t; se'' = recomputeStdErr et' tOpt + in go (restarts + 1) et' tOpt se'' estCIs' + Right p -> go' (ix + 1) (acc <> [p]) + | otherwise = + case profFun ix of + Left t -> let tOpt = ctOptimizer et' t; se'' = recomputeStdErr et' tOpt + in go (restarts + 1) et' tOpt se'' estCIs' + Right p -> go' (ix + 1) (acc <> [p]) + -- calculates the profile likelihood of a single parameter getProfile :: EvalTree -> Target -> Double -> Double -> Int -> Either Target ProfileT getProfile et theta stdErr_i tau_max ix | stdErr_i == 0.0 = pure $ ProfileT (U.fromList [-tau_max, tau_max]) [theta, theta] (theta U.! ix) (const (theta U.! ix)) (const tau_max) | otherwise = do negDelta <- go kmax (-stdErr_i / 8) 0 1 mempty + let !negLen = length (fst negDelta) + !negTauRange = if null (fst negDelta) then (0,0) else (minimum (fst negDelta), maximum (fst negDelta)) posDelta <- go kmax (stdErr_i / 8) 0 1 p0 + let !posLen = length (fst posDelta) + !posTauRange = if null (fst posDelta) then (0,0) else (minimum (fst posDelta), maximum (fst posDelta)) let (taus', thetas') = negDelta <> posDelta taus = U.fromList taus' thetas = thetas' - (tau2theta, theta2tau) = createSplines taus thetas stdErr_i tau_max ix + (tau2theta, theta2tau) = createSplines taus thetas stdErr_i tau_max ix optTh pure $ ProfileT taus thetas optTh tau2theta theta2tau - where + where p0 = ([0], [theta_opt]) - kmax = 300 + kmax = 500 nll_opt = ctNLL et theta_opt theta_opt = ctOptimizer et theta optTh = theta_opt U.! ix @@ -227,7 +281,7 @@ go 0 delta _ _ acc = Right acc go k delta t inv_slope acc@(taus, thetas) | isNaN inv_slope = Right acc - | nll_cond < nll_opt = Left theta_t + | nll_cond < nll_opt - 1e-6 * abs nll_opt = Left theta_t | abs tau > tau_max = Right acc' | otherwise = go (k-1) delta (t + inv_slope) inv_slope' acc' @@ -237,8 +291,17 @@ theta_t = minimizer theta_delta (nll_cond, grad) = ctGradNLL et theta_t zv = grad U.! ix - inv_slope' = min 4.0 . max 0.0625 . abs $ (tau / (stdErr_i * zv)) - tau = signum delta * sqrt (2*nll_cond - 2*nll_opt) + -- For LeastSquares, the correct profile likelihood statistic is + -- n * log(MSE(t)/MSE(opt)) ~ chi2_1, not 2*(MSE(t) - MSE(opt)). + tau = case ctDist et of + LeastSquares -> + let nD = fromIntegral (ctRows et) :: Double + r = max nll_cond 1e-30 / max nll_opt 1e-30 + in signum delta * sqrt (max 0 (nD * log r)) + _ -> signum delta * sqrt (max 0 (2*nll_cond - 2*nll_opt)) + inv_slope' = if abs zv < 1e-12 * abs stdErr_i + then min 4.0 . max 0.0625 $ abs (delta * 8) + else min 4.0 . max 0.0625 . abs $ (tau / (stdErr_i * zv)) acc' = if nll_cond == nll_opt || maybe False (tau ==) (listToMaybe taus) || isNaN tau then acc else (tau:taus, theta_t:thetas) @@ -253,32 +316,41 @@ taus = U.fromList [-tau_max, tau_max] thetas = [theta, theta] theta_i = theta U.! ix - getPoint = getEndPoint et theta tau_max ix + getPoint = getEndPoint et theta tau_max stdErr_i ix leftPt = getPoint True rightPt = getPoint False tau2theta tau = if tau < 0 then leftPt else rightPt -getEndPoint :: EvalTree -> Target -> Double -> Int -> Bool -> Double -getEndPoint et theta tau_max ix isLeft = - case minimizeAugLag problem (G.convert theta_opt) of - Right sol -> solutionParams sol VS.! ix - Left _ -> theta_opt U.! ix +getEndPoint :: EvalTree -> Target -> Double -> Double -> Int -> Bool -> Double +getEndPoint et theta tau_max stdErr_i ix isLeft + | isNaN mle = 0/0 -- NaN: MLE itself is NaN + | f mle >= 0 = 0/0 -- NaN: MLE violates constraint + | isLeft && f lo <= 0 = 0/0 -- NaN: constraint satisfied at left bound + | not isLeft && f hi <= 0 = 0/0 -- NaN: constraint satisfied at right bound + | isLeft = bisect lo mle 0 + | otherwise = bisect mle hi 0 where n = U.length theta - theta_opt = ctOptimizer et theta nll_opt = ctNLL et theta_opt loss_crit = nll_opt + tau_max - - loss = subtract loss_crit . ctNLL et . G.convert - obj = (if isLeft then id else negate) . (VS.! ix) + mle = theta_opt U.! ix + -- Use a wide search range: 50x the standard error, with a minimum of 50x |mle| + -- This ensures we don't miss the CI boundary for parameters near zero + searchScale = max (abs mle * 50) (stdErr_i * 50) + lo = mle - searchScale + hi = mle + searchScale - stop = ObjectiveRelativeTolerance 1e-4 :| [MaximumEvaluations 1000] - localAlg = NELDERMEAD obj [] Nothing - local = LocalProblem (fromIntegral n) stop localAlg - constraint = InequalityConstraint (Scalar loss) 1e-6 + -- Profiled NLL: fix theta[ix]=t, re-optimize all other params + f t = let x = U.generate n (\j -> if j == ix then t else theta_opt U.! j) + reopt = ctOptimizerFixed et ix (G.convert x) + in ctNLL et reopt - loss_crit - problem = AugLagProblem [] [] (AUGLAG_LOCAL local [constraint] []) + bisect a b k + | k >= 60 || abs (b - a) < 1e-12 = (a + b) / 2 + | f mid <= 0 = if isLeft then bisect a mid (k+1) else bisect mid b (k+1) + | otherwise = if isLeft then bisect mid b (k+1) else bisect a mid (k+1) + where mid = (a + b) / 2 {-# INLINE getEndPoint #-} -- Based on @@ -290,7 +362,7 @@ | otherwise = let (taus', thetas') = solLeft <> ([0], [theta_opt]) <> solRight taus = U.fromList taus' thetas = thetas' - (tau2theta, theta2tau) = createSplines taus thetas stdErr_i tau_max ix + (tau2theta, theta2tau) = createSplines taus thetas stdErr_i tau_max ix optTh in pure $ ProfileT taus thetas optTh tau2theta theta2tau where dflt = ProfileT (U.fromList [-tau_max, tau_max]) [theta, theta] (theta U.! ix) (const (theta U.! ix)) (const tau_max) @@ -349,12 +421,26 @@ hess = hessianNLL dist mYerr xss ys tree theta fexcept :: SomeException -> IO Columns - fexcept e = trace ("cov NegDef" <> show (toRowMajor hess)) $ pure ident + fexcept _ = pure ident - cov = unsafePerformIO $ catch (invChol hess) fexcept + covRaw = unsafePerformIO $ catch (invChol hess) fexcept + -- For LeastSquares, the Hessian code computes sum(fx*fy - res*fxy) = X^T X, + -- but the actual Hessian of the Gaussian NLL profile is -1/MSE * X^T X. + -- So cov_code = inv(X^T X) and cov_correct = MSE * inv(X^T X) = MSE * cov_code. + sigma2 = case dist of + LeastSquares -> let mse = compileLoss xss (buildLoss (NLL LeastSquares) (fromIntegral n) tree) ys mYerr theta + in max mse 1e-10 -- avoid division by zero + _ -> 1.0 -- no scaling needed for NLL-based losses + + scaleFactor = case dist of + LeastSquares -> sigma2 + _ -> 1.0 + + cov = fromRowMajor k k $ U.map (* scaleFactor) (toRowMajor covRaw) + covMat = toRowMajor cov - stdErr = U.generate k (\ix -> sqrt $ covMat U.! (ix * k + ix)) + stdErr = U.generate k (\ix -> sqrt $ max 0 (covMat U.! (ix * k + ix))) stdErrSq = case outer stdErr stdErr of Right v -> v @@ -364,16 +450,63 @@ corr = fromRowMajor k k $ U.generate (k * k) (\ix -> covMat U.! ix / stdErrSqMat U.! ix) -- Create splines for profile-t -createSplines :: Target -> Columns -> Double -> Double -> Int -> (Double -> Double, Double -> Double) -createSplines taus thetas se tau_max ix +-- We enforce monotonicity of theta w.r.t. tau: if the profile walk produced +-- non-monotonic pairs (theta[i] < theta[i-1] for positive tau direction or vice versa), +-- we keep only the outermost monotonic subsequence to prevent spline extrapolation garbage. +createSplines :: Target -> Columns -> Double -> Double -> Int -> Double -> (Double -> Double, Double -> Double) +createSplines taus thetas se tau_max ix optTh | n < 2 = (genSplineFun [(-tau_max, -se), (tau_max, se)], genSplineFun [(-se, 0), (se, 1)]) | otherwise = (tau2theta, theta2tau) where n = U.length taus cols = getCol ix thetas - nubOnFirst = nubBy (\x y -> fst x == fst y) - tau2theta = genSplineFun $ nubOnFirst $ sortOnFirst taus cols - theta2tau = genSplineFun $ nubOnFirst $ sortOnFirst cols taus + 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 + +-- | Enforce monotonicity for (tau, theta) pairs sorted by tau. +-- Split at tau=0; both halves keep theta non-decreasing: +-- negative half: as tau increases from -tau_max toward 0, theta increases +-- positive half: as tau increases from 0 toward tau_max, theta increases +enforceMonotonicTau :: [(Double, Double)] -> [(Double, Double)] +enforceMonotonicTau [] = [] +enforceMonotonicTau [p] = [p] +enforceMonotonicTau pts = negMono ++ posMono + where + (neg, pos) = span (\(t, _) -> t <= 0) pts + negMono = monotoneInc neg + posMono = monotoneInc pos + +-- | Enforce monotonicity for (theta, tau) pairs sorted by theta. +-- Split at theta=optTh; both halves keep tau non-decreasing: +-- left half: as theta increases toward optTh, tau increases toward 0 +-- right half: as theta increases from optTh, tau increases from 0 +enforceMonotonicTheta :: Double -> [(Double, Double)] -> [(Double, Double)] +enforceMonotonicTheta _ [] = [] +enforceMonotonicTheta _ [p] = [p] +enforceMonotonicTheta optTh pts = negMono ++ posMono + where + (neg, pos) = span (\(t, _) -> t <= optTh) pts + negMono = monotoneInc neg + posMono = monotoneInc pos + +-- | Keep longest prefix of non-decreasing second elements. +monotoneInc :: [(Double, Double)] -> [(Double, Double)] +monotoneInc [] = [] +monotoneInc [x] = [x] +monotoneInc ((t0,th0):(t1,th1):rest) + | th1 >= th0 = (t0,th0) : monotoneInc ((t1,th1):rest) + | otherwise = monotoneInc ((t0,th0):rest) + +-- | Keep longest prefix of non-increasing second elements. +monotoneDec :: [(Double, Double)] -> [(Double, Double)] +monotoneDec [] = [] +monotoneDec [x] = [x] +monotoneDec ((t0,th0):(t1,th1):rest) + | th1 <= th0 = (t0,th0) : monotoneDec ((t1,th1):rest) + | otherwise = monotoneDec ((t0,th0):rest) getCol :: Int -> Columns -> Target getCol ix mtx = U.generate (length mtx) (\j -> (mtx !! j) U.! ix) diff --git a/src/Algorithm/SRTree/Likelihoods.hs b/src/Algorithm/SRTree/Likelihoods.hs --- a/src/Algorithm/SRTree/Likelihoods.hs +++ b/src/Algorithm/SRTree/Likelihoods.hs @@ -47,7 +47,6 @@ import Control.Concurrent (getNumCapabilities) import Control.Concurrent.Async (forConcurrently) -import Debug.Trace import Data.SRTree.Print import Control.Monad.State.Strict import Control.Monad.Identity diff --git a/src/Algorithm/SRTree/ModelSelection.hs b/src/Algorithm/SRTree/ModelSelection.hs --- a/src/Algorithm/SRTree/ModelSelection.hs +++ b/src/Algorithm/SRTree/ModelSelection.hs @@ -38,8 +38,6 @@ import qualified Data.Vector.Unboxed as U import Algorithm.SRTree.Compile -import Debug.Trace - -- | Bayesian information criterion bic :: EvaluatedTree -> Double bic et = valParams et * log (valRows et) + 2 * valLoss et diff --git a/src/Algorithm/SRTree/NonlinearOpt.hs b/src/Algorithm/SRTree/NonlinearOpt.hs --- a/src/Algorithm/SRTree/NonlinearOpt.hs +++ b/src/Algorithm/SRTree/NonlinearOpt.hs @@ -12,7 +12,13 @@ -- ----------------------------------------------------------------------------- module Algorithm.SRTree.NonlinearOpt - where + ( minimizeNLLWith + , minimizeNLL' + , minimizeNLL + , minimizeNLLWithFixedParam' + , minimizeNLLWithFixedParam + , compileLossAndGrad + ) where import Algorithm.SRTree.Likelihoods import Numeric.Optimization.NLOPT @@ -31,8 +37,6 @@ import Control.Monad.State.Strict import Control.Monad.Identity -import Debug.Trace - minimizeNLLWith :: (VS.Vector Double -> (Double, VS.Vector Double)) -> (ObjectiveD -> (Maybe VectorStorage) -> LocalAlgorithm) -> Int -> Target -> (Target, Double, Int) minimizeNLLWith funAndGrad alg niter t0 | niter == 0 = (t0, f, 0) @@ -53,13 +57,20 @@ t_opt' = G.convert t_opt {-# INLINE minimizeNLLWith #-} +-- | Compile the loss function and gradient for a tree, returning a reusable +-- closure. Use this when you need to optimize the same expression with +-- multiple random restarts — compile once, call the closure many times. +compileLossAndGrad :: ADBackEnd -> Loss -> Maybe Target -> Columns -> Target -> Fix SRTree -> VS.Vector Double -> (Double, VS.Vector Double) +compileLossAndGrad backend dist mYerr xss ys tree = + let m = V.length ys + tree' = buildLoss dist (fromIntegral m) tree + in compileFunAndGrad backend xss ys mYerr tree' + -- | minimizes the negative log-likelihood of the expression minimizeNLL' :: (ObjectiveD -> (Maybe VectorStorage) -> LocalAlgorithm) -> ADBackEnd -> Loss -> Maybe Target -> Int -> Columns -> Target -> Fix SRTree -> Target -> (Target, Double, Int) minimizeNLL' alg backend dist mYerr niter xss ys tree t0 = minimizeNLLWith funAndGrad alg niter t0 where - m = V.length ys - tree' = buildLoss dist (fromIntegral m) tree - funAndGrad = compileFunAndGrad backend xss ys mYerr tree' + funAndGrad = compileLossAndGrad backend dist mYerr xss ys tree minimizeNLL :: ADBackEnd -> Loss -> Maybe Target -> Int -> Columns -> Target -> Fix SRTree -> Target -> (Target, Double, Int) diff --git a/src/Algorithm/SRTree/Utils.hs b/src/Algorithm/SRTree/Utils.hs --- a/src/Algorithm/SRTree/Utils.hs +++ b/src/Algorithm/SRTree/Utils.hs @@ -15,7 +15,6 @@ import Data.List (unfoldr) import Data.SRTree.Eval -import Debug.Trace (traceShow) -- | Internal helper to get dimensions (rows, columns) matSize :: Columns -> (Int, Int) diff --git a/src/Text/ParseSR.hs b/src/Text/ParseSR.hs --- a/src/Text/ParseSR.hs +++ b/src/Text/ParseSR.hs @@ -26,8 +26,6 @@ import qualified Data.Map.Strict as Map import Data.List.Split ( splitOn ) -import Debug.Trace (trace, traceShow) - -- * Data types -- | Parser of a symbolic regression tree with `Int` variable index and @@ -68,9 +66,6 @@ --parsePat :: B.ByteString -> Either String Pattern --parsePat = eitherResult . (`feed` "") . parse parsePatExpr . putEOL . B.strip - -eitherResult' :: Show r => Result r -> Either String r -eitherResult' res = trace (show res) $ eitherResult res -- * Parsers 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.2 +version: 3.0.0.3 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 @@ -196,10 +196,30 @@ , zlib >=0.6.3 && <0.8 default-language: Haskell2010 +executable test-ci + main-is: Main.hs + other-modules: + Paths_srtree + hs-source-dirs: + apps/TestCI + 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 + , filepath >=1.4.0.0 && <1.6 + , random >=1.2 && <1.4 + default-language: Haskell2010 + test-suite srtree-test type: exitcode-stdio-1.0 main-is: Spec.hs other-modules: + CITests EqSatTests StoreTests Paths_srtree diff --git a/test/CITests.hs b/test/CITests.hs new file mode 100644 --- /dev/null +++ b/test/CITests.hs @@ -0,0 +1,149 @@ +module CITests (tests) where + +import Test.HUnit +import Algorithm.SRTree.ConfidenceIntervals + ( monotoneInc, monotoneDec, enforceMonotonicTau, enforceMonotonicTheta, createSplines ) +import Algorithm.SRTree.Utils ( genSplineFun ) +import qualified Data.Vector.Unboxed as U + +eps :: Double +eps = 1e-9 + +approxEq :: Double -> Double -> Bool +approxEq a b = abs (a - b) < eps + +-- | monotoneInc keeps longest non-decreasing prefix of second elements +test_monotoneInc :: Test +test_monotoneInc = TestLabel "monotoneInc" $ TestCase $ do + -- already non-decreasing: keep all + assertEqual "all kept" [(1,10),(2,20),(3,30)] (monotoneInc [(1,10),(2,20),(3,30)]) + -- drop trailing decrease + assertEqual "drop tail" [(1,10),(2,20)] (monotoneInc [(1,10),(2,20),(3,15)]) + -- single element + assertEqual "single" [(1,5)] (monotoneInc [(1,5)]) + -- empty + assertEqual "empty" [] (monotoneInc []) + -- flat is ok (non-decreasing) + assertEqual "flat ok" [(1,5),(2,5),(3,5)] (monotoneInc [(1,5),(2,5),(3,5)]) + -- decrease at start: keeps first, skips the decrease, then keeps later increase + assertEqual "early decrease" [(1,10),(3,15)] (monotoneInc [(1,10),(2,5),(3,15)]) + +-- | monotoneDec keeps longest non-increasing prefix of second elements +test_monotoneDec :: Test +test_monotoneDec = TestLabel "monotoneDec" $ TestCase $ do + -- already non-increasing: keep all + assertEqual "all kept" [(1,30),(2,20),(3,10)] (monotoneDec [(1,30),(2,20),(3,10)]) + -- drop trailing increase + assertEqual "drop tail" [(1,30),(2,20)] (monotoneDec [(1,30),(2,20),(3,25)]) + -- single element + assertEqual "single" [(1,5)] (monotoneDec [(1,5)]) + -- empty + assertEqual "empty" [] (monotoneDec []) + -- flat is ok (non-increasing) + assertEqual "flat ok" [(1,5),(2,5),(3,5)] (monotoneDec [(1,5),(2,5),(3,5)]) + +-- | enforceMonotonicTau: split at tau=0, negative half non-increasing theta, +-- positive half non-decreasing theta +test_enforceMonotonicTau :: Test +test_enforceMonotonicTau = TestLabel "enforceMonotonicTau" $ TestCase $ do + -- well-formed data (monotonic in both halves) + let wellFormed = [(-2.0, 0.2), (-1.0, 0.5), (0.0, 1.0), (1.0, 1.5), (2.0, 1.8)] + assertEqual "well-formed" wellFormed (enforceMonotonicTau wellFormed) + + -- non-monotonic negative half (bump): (-1.0, 0.5) then (-0.5, 0.7) is increasing + let nonMonNeg = [(-2.0, 0.2), (-1.5, 0.4), (-1.0, 0.5), (-0.5, 0.7), (0.0, 1.0), (0.5, 1.3), (1.0, 1.5)] + let result = enforceMonotonicTau nonMonNeg + -- negative half: theta non-decreasing from -tau_max to 0, so all kept + assertBool "negative half preserved" (length result >= 5) + + -- non-monotonic positive half: theta decreases at tau=1.5 + let nonMonPos = [(-1.0, 0.5), (0.0, 1.0), (0.5, 1.3), (1.0, 1.5), (1.5, 1.4), (2.0, 1.8)] + let result2 = enforceMonotonicTau nonMonPos + -- should drop (1.5, 1.4) since it breaks non-decreasing + let posPart = filter (\(t,_) -> t > 0) result2 + pairs = zip posPart (tail posPart) + assertBool "positive monotonic" (all (\((_,a), (_,b)) -> b >= a) pairs) + where + +-- | enforceMonotonicTheta: split at theta=optTh, left half non-increasing tau, +-- right half non-decreasing tau +test_enforceMonotonicTheta :: Test +test_enforceMonotonicTheta = TestLabel "enforceMonotonicTheta" $ TestCase $ do + -- optTh = 1.0, data sorted by theta + let optTh = 1.0 + wellFormed = [(0.2, -2.0), (0.5, -1.0), (1.0, 0.0), (1.5, 1.0), (2.0, 2.0)] + assertEqual "well-formed" wellFormed (enforceMonotonicTheta optTh wellFormed) + + -- split at theta=1.0 (optTh), not theta=0 + -- data: theta < 1.0 should have negative tau, theta > 1.0 should have positive tau + let mixedThetas = [(0.5, -1.0), (0.8, -0.5), (1.0, 0.0), (1.2, 0.5), (1.5, 1.0)] + let result = enforceMonotonicTheta optTh mixedThetas + assertEqual "all kept" mixedThetas result + + -- non-monotonic: tau jumps back at theta=1.2 + let nonMon = [(0.5, -1.0), (1.0, 0.0), (1.2, 0.8), (1.5, 0.6), (2.0, 2.0)] + let result2 = enforceMonotonicTheta optTh nonMon + -- right half (theta > 1.0): tau non-decreasing, so (1.5, 0.6) after (1.2, 0.8) is dropped + assertBool "right half monotonic" (length result2 < length nonMon) + +-- | createSplines: basic spline creation and evaluation +test_createSplines :: Test +test_createSplines = TestLabel "createSplines" $ TestCase $ do + -- Create a simple linear profile: tau = theta - 1.0 (optTh = 1.0) + let n = 20 + optTh = 1.0 + se = 0.5 + tau_max = 3.0 + taus = U.fromList [ -tau_max + 2*tau_max * fromIntegral i / fromIntegral (n-1) | i <- [0..n-1] ] + -- theta = 1.0 + tau/3 (linear relationship) + thetas = [ U.fromList [ optTh + (taus U.! i) / 3.0 | _ <- [0] ] | i <- [0..n-1] ] + (tau2theta, _theta2tau) = createSplines taus thetas se tau_max 0 optTh + + -- at tau=0, should return approximately optTh + let atZero = tau2theta 0.0 + assertBool ("tau2theta(0) ~ optTh: " ++ show atZero) (approxEq atZero optTh) + + -- at tau=tau_max, should be approximately optTh + tau_max/3 + let atMax = tau2theta tau_max + expected_atMax = optTh + tau_max / 3.0 + assertBool ("tau2theta(tau_max) ~ expected: " ++ show atMax ++ " vs " ++ show expected_atMax) + (abs (atMax - expected_atMax) < 0.5) -- generous tolerance for spline overshoot + + -- at tau=-tau_max, should be approximately optTh - tau_max/3 + let atMin = tau2theta (-tau_max) + expected_atMin = optTh - tau_max / 3.0 + assertBool ("tau2theta(-tau_max) ~ expected: " ++ show atMin ++ " vs " ++ show expected_atMin) + (abs (atMin - expected_atMin) < 0.5) + +-- | Regression: negative-tau data must not be dropped +test_negative_tau_preserved :: Test +test_negative_tau_preserved = TestLabel "negative_tau_preserved" $ TestCase $ do + let optTh = 1.0 + se = 0.5 + tau_max = 3.0 + -- Monotonically decreasing theta for negative tau + negTaus = [-2.5, -2.0, -1.5, -1.0, -0.5] + posTaus = [0.5, 1.0, 1.5, 2.0, 2.5] + taus = U.fromList (negTaus ++ [0.0] ++ posTaus) + thetas = [ U.fromList [ optTh + (taus U.! i) / 3.0 ] | i <- [0 .. U.length taus - 1] ] + (tau2theta, _) = createSplines taus thetas se tau_max 0 optTh + + -- The spline MUST give a value below optTh for negative tau + let atNeg = tau2theta (-2.0) + assertBool ("tau2theta(-2.0) < optTh: got " ++ show atNeg) + (atNeg < optTh - 0.1) + + -- The spline MUST give a value above optTh for positive tau + let atPos = tau2theta 2.0 + assertBool ("tau2theta(2.0) > optTh: got " ++ show atPos) + (atPos > optTh + 0.1) + +tests :: Test +tests = TestLabel "ConfidenceIntervals" $ TestList + [ test_monotoneInc + , test_monotoneDec + , test_enforceMonotonicTau + , test_enforceMonotonicTheta + , test_createSplines + , test_negative_tau_preserved + ] diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -7,6 +7,7 @@ import Algorithm.SRTree.AD.Unboxed (CompiledTree, compileTree, compileTreeMulti, evalGrad, evalGradVec, evalGradMulti) import qualified EqSatTests import qualified StoreTests +import qualified CITests import Data.SRTree.Random (randomTree, tossBiased, randomFrom) import System.Random (mkStdGen) import Control.Monad.State.Strict (evalStateT) @@ -109,6 +110,7 @@ , TestLabel "benchgrad" test_benchgrad , TestLabel "eqsat" EqSatTests.tests , TestLabel "store" StoreTests.tests + , CITests.tests ] if failures counts /= 0 || errors counts /= 0 then error "Some tests failed"