The coiterative comonad generated by a comonad

Coiterative comonads represent non-terminating, productive computations.

They are the dual notion of iterative monads. While iterative computations produce no values or eventually terminate with one, coiterative computations constantly produce values and they never terminate.

It's simpler form, Coiter, is an infinite stream of data. CoiterT extends this so that each step of the computation can be performed in a comonadic context.

This is literate Haskell! To run the example, open the source and copy this comment block into a new file with '.lhs' extension.

Many numerical approximation methods compute infinite sequences of results; each, hopefully, more accurate than the previous one.

Newton's method to find zeroes of a function is one such algorithm.

 {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-}

 import Control.Comonad.Trans.Coiter
 import Control.Comonad.Env
 import Control.Applicative
 import Data.Foldable (toList, find)

 data Function = Function {
   -- Function to find zeroes of
   function   :: Double -> Double,
   -- Derivative of the function
   derivative :: Double -> Double
 }
 
 data Result = Result {
   -- Estimated zero of the function
   value  :: Double,
   -- Estimated distance to the actual zero
   xerror :: Double,
   -- How far is value from being an actual zero; that is,
   -- the difference between @0@ and @f value@
   ferror :: Double
 } deriving (Show)
 
 data Outlook = Outlook { result :: Result,
                          -- Whether the result improves in future steps
                          progress :: Bool } deriving (Show)

To make our lives easier, we will store the problem at hand using the Env environment comonad.

 type Solution a = CoiterT (Env Function) a

Problems consist of a function and its derivative as the environment, and an initial value.

 type Problem = Env Function Double

We can express an iterative algorithm using unfold over an initial environment.

 newton :: Problem -> Solution Double
 newton = unfold (\wd ->
                     let  f  = asks function wd in
                     let df  = asks derivative wd in
                     let  x  = extract wd in
                     x - f x / df x)
 
 

To estimate the error, we look forward one position in the stream. The next value will be much more precise than the current one, so we can consider it as the actual result.

We know that the exact value of a function at one of it's zeroes is 0. So, ferror can be computed exactly as abs (f a - f 0) == abs (f a)

 estimateError :: Solution Double -> Result
 estimateError s =
   let a:a':_ = toList s in
   let f = asks function s in
   Result { value = a,
            xerror = abs $ a - a',
            ferror = abs $ f a
          }

To get a sense of when the algorithm is making any progress, we can sample the future and check if the result improves at all.

 estimateOutlook :: Int -> Solution Result -> Outlook
 estimateOutlook sampleSize solution =
   let sample = map ferror $ take sampleSize $ tail $ toList solution in
   let result = extract solution in
   Outlook { result = result,
             progress = ferror result > minimum sample } 

To compute the square root of c, we solve the equation x*x - c = 0. We will stop whenever the accuracy of the result doesn't improve in the next 5 steps.

The starting value for our algorithm is c itself. One could compute a better estimate, but the algorithm converges fast enough that it's not really worth it.

 squareRoot :: Double -> Maybe Result
 squareRoot c = let problem = flip env c (Function { function = (\x -> x*x - c),
                                                     derivative = (\x -> 2*x) })
                in 
                fmap result $ find (not . progress) $ 
                  newton problem =>> estimateError =>> estimateOutlook 5

This program will output the result together with the error.

 main :: IO ()
 main = putStrLn $ show $ squareRoot 4