-----------------------------------------------------------------------------
-- Copyright 2020, Ideas project team. This file is distributed under the
-- terms of the Apache License 2.0. For more information, see the files
-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution.
-----------------------------------------------------------------------------
module Domain.Hypothesis.Rules
   ( addAlphaRule, addDfRule, addDfBetweenWithin
   , addH0FromHARule, addH0FromHAEqualSignRule, addHARule
   , addHypothesesRule, addHypothesesChiSquaredRule
   , addObservedTotals, addExpectedFrequencies
   , addConclusionPValueRule, addRejectionRule
   , addTestFormulaRule, addTestValueRule
   , chooseTTestPairedRule, chooseTTestRule, chooseTTestTwoRule
   , chooseZTestRule, chooseRPearsonRule, chooseAnovaRule, chooseChiSquaredRule
   , computePValueTTest, computePValueZTest
   , criticalConclusionRule
   , hypothesesConclusionCriticalRule, hypothesesConclusionPValueRule
   , addStandardErrorSigma, addStandardErrorSD
   , determineSided
   , lookupTValueRule, lookupZValueRule, lookupRValueRule, lookupFValueRule
   , lookupChiValueRule
     -------------
   , inferSidedness, inferTestChoice, inferTestChoices, inferRejectionCritical
   , inferConclusionCritical, inferConclusionPValue, inferDf, inferVar
   , inferTestFormula, inferCriticalZWith, inferCriticalTWith
   , inferCriticalRWith, inferCriticalFWith, inferCriticalChiWith
   , inferTestValue, inferDfBetweenWithin, chiSquaredDf, chiSquaredTestValue
   , computeTotals, computeExpectedFrequencies, getTable
   ) where

import Control.Monad
import Data.List
import Data.Maybe
import Domain.Hypothesis.Common
import Domain.Math.Data.Relation
import Domain.Math.Expr hiding ((.*.), (./.), (^))
import Domain.Math.Numeric.Views
import Domain.Statistics.ComponentSet
import Domain.Statistics.Views
import Ideas.Common.Library

----------------------------------------------------------
-- Rules for determining the confidence level

addAlphaRule :: Rule ComponentSet
addAlphaRule =
  describe "Rule for adding the alpha component" .
  makeRule "component.alpha" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (derived cs `doesNotContain` SignificanceLevel)
      let alpha = case get SignificanceLevel (initials cs) of
                     Just (CExpr a) -> a
                     _ -> toExpr pickAlpha
      return $ append SignificanceLevel (CExpr alpha) cs

----------------------------------------------------------
-- Rules for constructing the hypotheses

determineSided :: Rule ComponentSet
determineSided =
  describe "Rule for determine one-/two-sided testing" .
  makeRule "component.sided" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
       guard (cs `doesNotContain` Sidedness)
       sided <- 
          case getTestType TestChoice (initials cs) of
             Just Anova -> return RightSided
             _ -> do
                ha <- getRelation AlternativeHypothesis cs
                return (sidedFromHA ha)
       return $ append Sidedness (CChoice $ Sided sided) cs

addHypothesesRule :: Rule ComponentSet
addHypothesesRule = 
   describe "Add null hypothesis and alternative hypothesis, in one step" $
   makeRule "component.hypotheses" f
 where
   f :: ComponentSet -> [ComponentSet]
   f = applyAll $ (addH0FromHARule ./. addH0FromHAEqualSignRule) .*. addHARule

addH0FromHARule :: Rule ComponentSet
addH0FromHARule =
  describe "Rule for adding the H0 component based on HA" .
  makeRule "component.h0-from-ha" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (cs `doesNotContain` NullHypothesis)
      ha <- getRelation AlternativeHypothesis cs
      return . append NullHypothesis
        (CRelation $ h0FromHA ha) $ cs

addH0FromHAEqualSignRule :: Rule ComponentSet
addH0FromHAEqualSignRule =
  describe "Rule for adding the H0 component based on HA; use equal sign (by convention)" .
  makeRule "component.h0-from-ha-eq" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (cs `doesNotContain` NullHypothesis)
      ha <- getRelation AlternativeHypothesis cs   
      return . append NullHypothesis
        (CRelation $ h0FromHAEqualSign ha) $ cs

addHARule :: Rule ComponentSet
addHARule =
  describe "Rule for adding the HA component" .
  makeRule "component.ha" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
       guard (derived cs `doesNotContain` AlternativeHypothesis)
       ha <- get AlternativeHypothesis (initials cs) 
       return $ append AlternativeHypothesis ha cs

----------------------------------------------------------
-- Rules for determining the properties of the data

chooseTTestRule :: Rule ComponentSet
chooseTTestRule =
  describe "Rule for choosing the type of test" .
  makeRule "component.test.t-test" $
     addTestChoice TTestOne

chooseTTestTwoRule :: Rule ComponentSet
chooseTTestTwoRule =
  describe "Rule for choosing the type of test" .
  makeRule "component.test.t-test-two" $
     addTestChoice TTestTwo

chooseTTestPairedRule :: Rule ComponentSet
chooseTTestPairedRule =
  describe "Rule for choosing the type of test" .
  makeRule "component.test.t-test-paired" $
     addTestChoice TTestPaired

-- | If the standard deviation of the population is known then the z-test is
-- /always/ chosen, thus a t-test is never chosen in that case. When this
-- information is unknown, then an estimation of the standard deviation of the
-- population has to be made based on the sample, in that case a t-test is
-- chosen.
--
-- When the sample size becomes large enough* it is also possible to choose a
-- z-test because the z-distribution looks like the t-distribution for large
-- sample sizes.
--
-- NOTE*: The threshold for what `large' means may vary. For now this threshold
--        is fixed at 100.
--
-- Consisely: If PopulationSdev known: ZTest
--            else if SampleSize big enough: ZTest or TTest
--            else TTest
chooseZTestRule :: Rule ComponentSet
chooseZTestRule =
  describe "Rule for choosing the type of test" .
  makeRule "component.test.z-test" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs =
      let largeThreshold = 100
      in do
        guard (derived cs `doesNotContain` TestChoice)
        n <- match naturalView $ fromMaybe 0 (getExpr SampleSize cs)
        if cs `contains` PopulationSdev || n >= largeThreshold 
           then return (append TestChoice (CChoice $ TestType ZTest) cs)
           else do
              -- This branch is similar to the previous rule body in revision 10549
              let tests = validTests cs
              guard (ZTest `elem` tests)
              return $ append TestChoice (CChoice $ TestType ZTest) cs

addTestFormulaRule :: Rule ComponentSet
addTestFormulaRule =
  describe "Rule for adding the test formula" .
  makeRule "component.test-formula" $ f
  where
    f :: ComponentSet -> [ComponentSet]
    f cs = do
      guard (derived cs `doesNotContain` TestFormula)
      rel <- inferTestFormula cs
      return $ append TestFormula (CRelation rel) cs

addTestValueRule :: Rule ComponentSet
addTestValueRule =
  describe "Rule for adding the test value (from the formula)" .
  makeRule "component.test-value" $ f
  where
    f :: ComponentSet -> [ComponentSet]
    f cs = do
      guard (derived cs `doesNotContain` TestValue)
      tv <- inferTestValue cs
      return $ append TestValue (CRelation tv) cs

inferTestValue :: MonadPlus m => ComponentSet -> m (Relation Expr)
inferTestValue cs =
   case getRelation TestValue (initials cs) of
      Just initialTestValue -> 
         return initialTestValue
      _ | inferTestChoice cs == Just ChiSquared -> 
         chiSquaredTestValue cs
      _ -> do
         new <- msum (map return $ applyAll addTestFormulaRule cs)
         let cs' = substitute new
         var  <- leftHandSide  <$> getRelation TestFormula new
         expr <- rightHandSide <$> getRelation TestFormula cs'
         val  <- matchM doubleView expr
         return $ var .==. toExpr val

----------------------------------------------------------
-- Rules for performing a T-Test

addDfRule :: Rule ComponentSet
addDfRule =
  describe "Rule for adding the degrees of freedom component" .
  makeRule "component.df" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (derived cs `doesNotContain` Df)
      case get Df (initials cs) of
         Just initialDf -> 
            return $ append Df initialDf cs
         _ | inferTestChoice cs == Just ChiSquared -> do
            df <- chiSquaredDf cs
            return $ append Df (CExpr (toExpr df)) cs
         _ -> do 
            test <- inferTestChoice cs
            df   <- degreesOfFreedomFromTest test (chooseTypeOfTest cs)
            val  <- matchM doubleView (getSubstitution cs |-> df)
            return $ append Df (CExpr (toExpr val)) cs

lookupTValueRule :: Rule ComponentSet
lookupTValueRule =
  describe "Rule for looking up a t-value" .
  makeRule "component.critical.t-value" $ 
     inferCriticalTWith $ \sided alpha df -> 
        [computeCriticalT sided alpha df]

-- shared function for computing t-value and r-value
inferCriticalGenericWith :: (TestType -> Bool) -> Expr -> (Sided -> Double -> Double -> [Double]) 
                         -> ComponentSet -> [ComponentSet]
inferCriticalGenericWith forTestType var compute cs = do
   guard (cs `doesNotContain` Critical)
   guard (derived cs `contains` AlternativeHypothesis)
   guard (any forTestType (inferTestChoices cs))
   cs'   <- matchM substitutedView cs
   alpha <- matchM doubleView <=< getExpr SignificanceLevel $ cs   
   let cs'' = case inferDf cs' of 
                 Just df -> substitute (append Df (CExpr df) cs')
                 Nothing -> cs'
   df <- matchM doubleView =<< getExpr Df cs''
   sided <- inferSidedness cs
   value <- compute sided alpha df
   return $ append Critical (CRelation $ var .==. fromDouble value) cs

inferCriticalTWith :: (Sided -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet]
inferCriticalTWith = inferCriticalGenericWith isTTest (Var "tcrit")

computePValueTTest :: Rule ComponentSet
computePValueTTest =
  describe "Rule for computing the p-value for a t-test" .
  makeRule "component.p-value.t-test" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (cs `doesNotContain` PValue)
      guard (derived cs `contains` AlternativeHypothesis)
      guard (maybe False isTTest (inferTestChoice cs))
      let cs'   = substitute cs
      testStatistic <- match doubleView =<< fmap rightHandSide (getRelation TestValue cs')

      let cs'' = case inferDf cs' of 
                    Just df -> substitute (append Df (CExpr df) cs')
                    Nothing -> cs'
      df <- matchM doubleView =<< getExpr Df cs''
      sided         <- inferSidedness cs
      let value = computePValueT sided testStatistic df
      return $ append PValue (CExpr $ fromDouble value) cs

----------------------------------------------------------
-- Rules for performing a Z-Test

lookupZValueRule :: Rule ComponentSet
lookupZValueRule =
  describe "Rule for looking up a z-value" .
  makeRule "component.critical.z-value" $ 
     inferCriticalZWith $ \sided alpha -> 
        [computeCriticalZ sided alpha]

inferCriticalZWith :: (Sided -> Double -> [Double]) -> ComponentSet -> [ComponentSet]
inferCriticalZWith compute cs = do
  guard (cs `doesNotContain` Critical)
  guard (derived cs `contains` AlternativeHypothesis)
  guard (maybe False (ZTest ==) (inferTestChoice cs))
  alpha <- matchM doubleView <=< getExpr SignificanceLevel $ cs
  sided <- inferSidedness cs
  value <- compute sided alpha
  return $ append Critical (CRelation $ Var "zcrit" .==. fromDouble value) cs

computePValueZTest :: Rule ComponentSet
computePValueZTest =
  describe "Rule for computing the p-value for a z-test" .
  makeRule "component.p-value.z-test" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (cs `doesNotContain` PValue)
      guard (maybe False (ZTest ==) (inferTestChoice cs))
      guard (derived cs `contains` AlternativeHypothesis)
      let cs'   = substitute cs
      -- TestStatistic was renamed to TestFormula, which is now a CRelation
      testStatistic <- match doubleView =<< fmap rightHandSide (getRelation TestValue cs')
      sided <- inferSidedness cs
      let value = computePValueZ sided testStatistic
      return $ append PValue (CExpr $ fromDouble value) cs

----------------------------------------------------------
-- Rules for performing R-Pearson

lookupRValueRule :: Rule ComponentSet
lookupRValueRule =
  describe "Rule for looking up a r-value" .
  makeRule "component.critical.r-value" $ 
     inferCriticalRWith $ \sided alpha df -> 
        [computeCriticalR sided alpha df]

chooseRPearsonRule :: Rule ComponentSet
chooseRPearsonRule =
  describe "Rule for choosing the type of test" .
  makeRule "component.test.r-pearson" $
     addTestChoice RPearson

addTestChoice :: TestType -> ComponentSet -> Maybe ComponentSet
addTestChoice testType cs = do
  guard (derived cs `doesNotContain` TestChoice)
  let tests = validTests cs
  guard (testType `elem` tests)
  return $ append TestChoice (CChoice $ TestType testType) cs

inferCriticalRWith :: (Sided -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet]
inferCriticalRWith = inferCriticalGenericWith (== RPearson) (Var "rcrit")

----------------------------------------------------------
-- Rules for performing Anova

lookupFValueRule :: Rule ComponentSet
lookupFValueRule =
  describe "Rule for looking up a F-value" .
  makeRule "component.critical.f-value" $ inferCriticalFWith computeCriticalF

chooseAnovaRule :: Rule ComponentSet
chooseAnovaRule =
  describe "Rule for choosing the type of test" .
  makeRule "component.test.anova" $
     addTestChoice Anova

addDfBetweenWithin :: Rule ComponentSet
addDfBetweenWithin = describe "Add df between and within (for Anova)" $
  makeRule "component.df-anova" f
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (derived cs `doesNotContain` DfBetween)
      guard (derived cs `doesNotContain` DfWithin)
      (between, within) <- inferDfBetweenWithin cs
      return $ append DfBetween (CExpr $ toExpr between) $ 
               append DfWithin  (CExpr $ toExpr within)  cs

inferDfBetweenWithin :: Monad m => ComponentSet -> m (Double, Double)
inferDfBetweenWithin cs = do
   n     <- matchM doubleView <=< getExpr SampleSize $ cs
   let nrOfGroups = 2 -- always 2, for now
       dfBetween  = nrOfGroups - 1
       dfWithin   = n - nrOfGroups
   return (dfBetween, dfWithin)

inferCriticalFWith :: (Double -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet]
inferCriticalFWith compute cs = do
  guard (cs `doesNotContain` Critical)
  guard (derived cs `contains` AlternativeHypothesis)
  guard (maybe False (Anova ==) (inferTestChoice cs))
  alpha <- matchM doubleView <=< getExpr SignificanceLevel $ cs
  (dfBetween, dfWithin) <- inferDfBetweenWithin cs
  value <- compute dfBetween dfWithin alpha
  return $ append Critical (CRelation $ Var "Fcrit" .==. fromDouble value) cs

----------------------------------------------------------
-- Rules for performing Chi-Squared

chooseChiSquaredRule :: Rule ComponentSet
chooseChiSquaredRule = describe "Rule for choosing the type of test" .
  makeRule "component.test.chi-squared" $
     addTestChoice ChiSquared

addHypothesesChiSquaredRule :: Rule ComponentSet
addHypothesesChiSquaredRule = describe "Add hypotheses (null and alternative) for chi-squared" $
   makeRule "component.hypotheses-chi-squared" f
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (derived cs `doesNotContain` NullHypothesis)
      guard (derived cs `doesNotContain` AlternativeHypothesis)
      guard (maybe False (ChiSquared ==) (inferTestChoice cs))
      return $ append NullHypothesis (CExpr independent) 
             $ append AlternativeHypothesis (CExpr dependent) cs

addObservedTotals :: Rule ComponentSet
addObservedTotals = describe "Add totals (rows and columns) for observed frequencies" $ 
   makeRule "component.observed-totals" f
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do  
      guard (derived cs `doesNotContain` ObservedColumnTotals)
      guard (derived cs `doesNotContain` ObservedRowTotals)
      guard (derived cs `doesNotContain` ObservedTotal)
      table <- getTable ObservedFrequencies cs
      let (rowTotals, columnTotals, total) = computeTotals table
      return $ append ObservedRowTotals (CExpr $ toExpr rowTotals) 
             $ append ObservedColumnTotals (CExpr $ toExpr columnTotals) 
             $ append ObservedTotal (CExpr $ toExpr total) cs

type ChiSquaredTotals = ([Int], [Int], Int)

computeTotals :: [[Int]] -> ChiSquaredTotals
computeTotals table = 
   let rowTotals    = map sum table
       columnTotals = map sum (transpose table)
       total        = sum rowTotals
   in (rowTotals, columnTotals, total)

computeExpectedFrequencies :: ChiSquaredTotals -> [[Double]]
computeExpectedFrequencies (rowTotals, columnTotals, total) = table
 where
   xss = map (replicate (length columnTotals)) rowTotals
   yss = replicate (length rowTotals) columnTotals
   table = zipWith (zipWith f) xss yss
   f x y = fromIntegral (x*y) / fromIntegral total

addExpectedFrequencies :: Rule ComponentSet
addExpectedFrequencies = describe "Add expected frequencies" $ 
   makeRule "component.expected-frequencies" f
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (derived cs `doesNotContain` ExpectedFrequencies)
      observed <- getTable ObservedFrequencies cs
      let totals   = computeTotals observed
          expected = computeExpectedFrequencies totals
      return $ append ExpectedFrequencies (CExpr $ toExpr expected) cs
   
lookupChiValueRule :: Rule ComponentSet
lookupChiValueRule = describe "Rule for looking up a chi^2-value" .
  makeRule "component.critical.chi-value" $ 
     inferCriticalChiWith $ \sided alpha df -> 
        computeCriticalChi sided alpha df

inferCriticalChiWith :: (Sided -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet]
inferCriticalChiWith = inferCriticalGenericWith (== ChiSquared) (Var "chicrit")

chiSquaredTestValue :: MonadPlus m => ComponentSet -> m (Relation Expr)
chiSquaredTestValue cs = do 
   observed <- getTable ObservedFrequencies cs
   expected <- getDoubleTable ExpectedFrequencies cs
   let table = zipWith (zipWith f) observed expected
       f o e = (fromIntegral o-e)^(2 :: Int) / e
       value = sum (map sum table)
   return (chiSquared .==. toExpr value)
   
chiSquaredDf :: MonadPlus m => ComponentSet -> m Int
chiSquaredDf cs = do
   observed <- getTable ObservedFrequencies cs
   guard (not $ null observed)
   let r = length observed
       c = length $ head observed
   return ((r-1)*(c-1))
   
getTable :: MonadPlus m => ComponentId -> ComponentSet -> m [[Int]]
getTable n cs = getExpr n cs >>= fromExpr

getDoubleTable :: MonadPlus m => ComponentId -> ComponentSet -> m [[Double]]
getDoubleTable n cs = getExpr n cs >>= fromExpr

----------------------------------------------------------
-- Rules for constructing the rejection region

-- lookup the Sidedness component; if not present, try to infer (using the rule)
inferSidedness :: MonadPlus m => ComponentSet -> m Sided
inferSidedness cs = 
   getSided Sidedness cs `mplus` do
      testType <- inferTestChoice cs
      guard (testType `elem` [Anova, ChiSquared])
      return RightSided
    `mplus` do
      cs' <- applyM determineSided cs
      getSided Sidedness cs'

inferVar :: MonadPlus m => ComponentSet -> m Expr
inferVar cs = do 
   testType <- inferTestChoice cs
   return (varForTestType testType)

varForTestType :: TestType -> Expr
varForTestType testType = 
   case testType of 
      ZTest      -> Var "z"
      RPearson   -> Var "r"
      Anova      -> Var "F"
      ChiSquared -> chiSquared
      _          -> Var "t"

inferTestChoice :: MonadPlus m => ComponentSet -> m TestType
inferTestChoice cs =
   case validTests cs of
      hd:_ -> return hd
      _    -> fail "no valid test choice"

inferTestChoices :: ComponentSet -> [TestType]
inferTestChoices cs = 
   case getTestType TestChoice cs of
      Just test -> [test]
      Nothing   -> validTests cs

inferTestFormula :: ComponentSet -> [Relation Expr]
inferTestFormula cs = 
   case getRelation TestFormula cs of
      Just formula -> [formula]
      Nothing -> do
         testType <- inferTestChoices cs
         let var  = varForTestType testType
         let test = chooseTypeOfTest cs
         t <- testFormulaFromTest testType test
         return (var .==. t)

inferRejectionCritical :: MonadPlus m => ComponentSet -> m (Relation Expr)
inferRejectionCritical cs = 
   getRelation RejectionCritical cs `mplus` do
      cs' <- applyM addRejectionRule cs
      getRelation RejectionCritical cs'

inferConclusionCritical :: MonadPlus m => ComponentSet -> m Bool
inferConclusionCritical cs = 
   getConclusion ConclusionCritical cs `mplus` do
      cs' <- applyM criticalConclusionRule cs
      getConclusion ConclusionCritical cs'

inferConclusionPValue :: MonadPlus m => ComponentSet -> m (Relation Expr)
inferConclusionPValue cs = 
   getRelation ConclusionPValue cs `mplus` do
      cs' <- applyM addConclusionPValueRule cs
      getRelation ConclusionPValue cs'

inferDf :: MonadPlus m => ComponentSet -> m Expr
inferDf cs =
   getExpr Df cs `mplus` do
      cs' <- applyM addDfRule cs
      getExpr Df cs'

addRejectionRule :: Rule ComponentSet
addRejectionRule =
  describe "Rule for constructing the rejection critical component" .
  makeRule "component.rejection.critical" $ f
  where
    f :: ComponentSet -> [ComponentSet]
    f cs = do
      guard (cs `doesNotContain` RejectionCritical)
      guard (derived cs `contains` AlternativeHypothesis)
      let cs' = substitute cs
      sided <- inferSidedness cs'
      testType <- inferTestChoices cs
      let rel = case testType of 
                   ZTest      -> sidedRelation sided (Var "z") (Var "zcrit")
                   RPearson   -> sidedRelation sided (Var "r") (Var "rcrit")
                   Anova      -> sidedRelation sided (Var "F") (Var "Fcrit")
                   ChiSquared -> sidedRelation sided chiSquared (Var "chicrit")
                   _          -> sidedRelation sided (Var "t") (Var "tcrit")
      return . append RejectionCritical
        (CRelation rel) $ cs

addConclusionPValueRule :: Rule ComponentSet
addConclusionPValueRule =
  describe "Rule for constructing the conclusion p-value component" .
  makeRule "component.conclusion.p-value" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      guard (cs `doesNotContain` ConclusionPValue)
      pv    <- matchM doubleView =<< getExpr PValue cs 
      alpha <- matchM doubleView =<< getExpr SignificanceLevel cs
      let relType = if pv > alpha then GreaterThan else LessThanOrEqualTo
      return . append ConclusionPValue
       (CRelation $ makeType relType (Var "p") (Var "alpha")) $ cs

----------------------------------------------------------
-- Rules for making a conclusion

criticalConclusionRule :: Rule ComponentSet
criticalConclusionRule =
  describe ("Rule for creating a conclusion based on a critical value " ++
  "and test statistic") .
  makeRule "component.critical-conclusion" $ f
  where
    f :: ComponentSet -> Maybe ComponentSet
    f cs = do
      rej <- inferRejectionCritical cs
      -- to do: rejection critical is added to component set only to get the substituted relation
      let cs' = substitute (append RejectionCritical (CRelation rej) cs)
      rejection <- getRelation RejectionCritical cs'
      guard (cs `doesNotContain` ConclusionCritical)
      lhs       <- match doubleView $ leftHandSide  rejection
      rhs       <- match doubleView $ rightHandSide rejection
      let result = eval (relationType rejection) lhs rhs
      return $ append ConclusionCritical (CChoice $ Conclusion result) cs

hypothesesConclusionCriticalRule :: Rule ComponentSet
hypothesesConclusionCriticalRule = 
   describe "derive the hypotheses conclusion from the critical conclusion" $
   makeRule "component.hypotheses-conclusion-critical" f 
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (cs `doesNotContain` ConclusionHypotheses)
      concl <- inferConclusionCritical cs
      let rejhyp = if concl then RejectH0 else DontRejectH0
      return $ append ConclusionHypotheses (CChoice (RejectionHypotheses rejhyp)) cs

hypothesesConclusionPValueRule :: Rule ComponentSet
hypothesesConclusionPValueRule = 
   describe "derive the hypotheses conclusion from the p-value" $
   makeRule "component.hypotheses-conclusion-pvalue" f
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (cs `doesNotContain` ConclusionHypotheses)
      rel <- inferConclusionPValue cs
      let rejhyp = if relationType rel == LessThanOrEqualTo then RejectH0 else DontRejectH0
      return $ append ConclusionHypotheses (CChoice (RejectionHypotheses rejhyp)) cs

addStandardErrorSigma :: Rule ComponentSet
addStandardErrorSigma = 
   describe "derive standard error from population standard deviation and sample size" $ 
   makeRule "component.standard-error-sigma" f 
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (cs `doesNotContain` StandardError)
      n   <- matchM doubleView =<< getExpr SampleSize cs
      psdev <- matchM doubleView =<< getRhsExpr PopulationSdev cs
      let se = Var "sigmaM" .==. toExpr (psdev / sqrt n)
      return $ append StandardError (CRelation se) cs
      
addStandardErrorSD :: Rule ComponentSet
addStandardErrorSD = 
   describe "derive standard error from sample standard deviation and sample size" $ 
   makeRule "component.standard-error-sd" f 
 where
   f :: ComponentSet -> Maybe ComponentSet
   f cs = do
      guard (cs `doesNotContain` StandardError)
      n   <- matchM doubleView =<< getExpr SampleSize cs
      sdev <- matchM doubleView =<< getRhsExpr SampleSdev cs
      let se = Var "SEM" .==. toExpr (sdev / sqrt n)
      return $ append StandardError (CRelation se) cs