Commit ff87a474 authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

Update inverse problems exercise

parent 417df820
%% Cell type:markdown id: tags:
# Exercises: inverse problems
The *Volterra operator* $V$ is defined by
$$V f(x) = \int_{x_0}^x f(x) \,dx.$$
A simple discrete approximation is
$$(V f)_i = h \sum_{j=0}^i f_j,$$
where $f_j = f(x_j)$ are the values of $f$ at equidistant sampling points $(x_j)_{j=1}^{N-1}$ with distance $h = x_1 - x_0$. As a matrix, this can simply be written as
$$V = h \begin{pmatrix}
1 & 0 & 0 & 0 & \ldots & 0 \\
1 & 1 & 0 & 0 & \ldots & 0 \\
1 & 1 & 1 & 0 & \ldots & 0 \\
\vdots & & & \ddots & & \\
1 & 1 & 1 & 1 & \ldots & 1
\end{pmatrix}$$
- Implement the Volterra operator with a parameter `h :: Double` as an instance of the `Operator` type class from the lecture, but with the `Arr` type alias replaced by
```haskell
type Arr r = Array r DIM1 Double
```
Though it is possible to do this using a matrix formulation and e.g. `mmultS`, this unneccessarily inefficient. A better possibility is to note the connection between the Volterra operator and Haskell's `scanl` function. Moreover, there is an implementation of `scanl` for unboxed `Vector`s from `Data.Vector.Unboxed` (*not* HMatrix's `Numeric.Container`) as well as conversion functions `toUnboxed` and `fromUnboxed` for Repa `U` arrays to and from these `Vector`s. Use these to implement the operator.
- Implement the Volterra operator with parameter `h :: Double` and its adjoint.
Though it is possible to do this using a matrix formulation and e.g. `mmultS`, this unneccessarily inefficient. A better possibility is to note the connection between the Volterra operator and Haskell's `scanl` function. Moreover, there is an implementation of `scanl` for unboxed `Vector`s from `Data.Vector.Unboxed` as well as O(1) conversion functions `toUnboxed` and `fromUnboxed` for Repa `U` arrays to and from these `Vector`s. Use these to implement the operator.
*Note*: The adjoint (i.e. transpose) of the operator can be found most easily by noting the similarity between the Volterra matrix above and its transpose: the adjoint simply proceeds in the opposite direction.
- Use the regularized CG implementation from the lecture to compute regularized *derivatives* of the functions
$$\begin{align}
f_1(x) &= e^{-x^2} \\
f_2(x) &= \text{sgn}(x) = \begin{cases} -1 & x < 0 \\ 1 & x \geq 0 \end{cases}
\end{align}$$
sampled at 200 points between -2 and 2 and distorted by random noise between -0.1 and 0.1. The CG tolerance should be set suffieciently low (e.g. 1e-6).
- Plot the results together with the noisy functions and, in the case of $f_1$, the analytical derivative.
%% Cell type:code id: tags:
``` haskell
```
......
%% Cell type:markdown id: tags:
# Exercises: inverse problems
The *Volterra operator* $V$ is defined by
$$V f(x) = \int_{x_0}^x f(x) \,dx.$$
A simple discrete approximation is
$$(V f)_i = h \sum_{j=0}^i f_j,$$
where $f_j = f(x_j)$ are the values of $f$ at equidistant sampling points $(x_j)_{j=1}^{N-1}$ with distance $h = x_1 - x_0$. As a matrix, this can simply be written as
$$V = h \begin{pmatrix}
1 & 0 & 0 & 0 & \ldots & 0 \\
1 & 1 & 0 & 0 & \ldots & 0 \\
1 & 1 & 1 & 0 & \ldots & 0 \\
\vdots & & & \ddots & & \\
1 & 1 & 1 & 1 & \ldots & 1
\end{pmatrix}$$
- Implement the Volterra operator with a parameter `h :: Double` as an instance of the `Operator` type class from the lecture, but with the `Arr` type alias replaced by
```haskell
type Arr r = Array r DIM1 Double
```
Though it is possible to do this using a matrix formulation and e.g. `mmultS`, this unneccessarily inefficient. A better possibility is to note the connection between the Volterra operator and Haskell's `scanl` function. Moreover, there is an implementation of `scanl` for unboxed `Vector`s from `Data.Vector.Unboxed` (*not* HMatrix's `Numeric.Container`) as well as conversion functions `toUnboxed` and `fromUnboxed` for Repa `U` arrays to and from these `Vector`s. Use these to implement the operator.
- Implement the Volterra operator with parameter `h :: Double` and its adjoint.
Though it is possible to do this using a matrix formulation and e.g. `mmultS`, this unneccessarily inefficient. A better possibility is to note the connection between the Volterra operator and Haskell's `scanl` function. Moreover, there is an implementation of `scanl` for unboxed `Vector`s from `Data.Vector.Unboxed` as well as O(1) conversion functions `toUnboxed` and `fromUnboxed` for Repa `U` arrays to and from these `Vector`s. Use these to implement the operator.
*Note*: The adjoint (i.e. transpose) of the operator can be found most easily by noting the similarity between the Volterra matrix above and its transpose: the adjoint simply proceeds in the opposite direction.
- Use the regularized CG implementation from the lecture to compute regularized *derivatives* of the functions
$$\begin{align}
f_1(x) &= e^{-x^2} \\
f_2(x) &= \text{sgn}(x) = \begin{cases} -1 & x < 0 \\ 1 & x \geq 0 \end{cases}
\end{align}$$
sampled at 200 points between -2 and 2 and distorted by random noise between -0.1 and 0.1. The CG tolerance should be set suffieciently low (e.g. 1e-6).
- Plot the results together with the noisy functions and, in the case of $f_1$, the analytical derivative.
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
import Control.Monad.Identity
import Data.Array.Repa as R hiding ((++))
import Data.Array.Repa as R
import qualified Data.Vector.Unboxed as V
import Data.Array.Repa.Helpers (runCGreg)
import Data.Random.Normal (normalsIO')
import Numeric.LinearAlgebra.Helpers (linspace)
```
%% Cell type:code id: tags:
``` haskell
-- Copied from the lecture, except for changing DIM2 to DIM1 in the type Arr
type Arr r = Array r DIM1 Double
class Operator a where
evalOp :: a -> Arr U -> Arr D
adjointOp :: a -> Arr U -> Arr D
data CGState = CGState { cgx :: Arr U
, cgp :: Arr U
, cgr :: Arr U
, cgr2 :: Double
}
cgreg :: Operator a => a -> Double -> Arr U -> Arr U -> [CGState]
cgreg op reg rhs initial =
runIdentity $ do
(res :: Arr U) <- computeP $ rhs -^ evalOp op initial
rInit <- computeP $ adjointOp op res
r2Init <- normSquaredP rInit
return $ iterate cgStep (CGState initial rInit rInit r2Init)
where normSquaredP = sumAllP . R.map (^(2::Int))
scale a = R.map (* a)
cgStep :: CGState -> CGState
cgStep (CGState x p r r2) =
runIdentity $ do
!(q :: Arr U) <- computeP $ evalOp op p
!p2 <- normSquaredP p
!q2 <- normSquaredP q
let alpha = r2 / (q2 + reg*p2)
!x' <- computeP $ x +^ scale alpha p
!(s :: Arr U) <- computeP $ adjointOp op q
!r' <- computeP $ r -^ scale alpha (s +^ scale reg p)
!r2' <- normSquaredP r'
let beta = r2' / r2
!p' <- computeP $ r' +^ scale beta p
return $ CGState x' p' r' r2'
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil _ [] = []
takeUntil predicate (x:xs)
| predicate x = [x]
| otherwise = x : takeUntil predicate xs
process :: Monad m => [a] -> (a -> m ()) -> m a
process xs f = foldM (\_ x -> f x >> return x) undefined xs
runCG :: Operator a => Double -> a -> Double -> Arr U -> IO (Arr U)
runCG tol op reg rhs = do
let initial = computeS $ fromFunction (extent rhs) (const 0)
let steps' = cgreg op reg rhs initial
let r20 = cgr2 $ head steps'
let steps = takeUntil (\x -> sqrt (cgr2 x / r20) < tol) steps'
result <- process (zip [(1::Int)..] steps) $ \(n, cgs) ->
putStrLn $ show n ++ " " ++ show (sqrt $ cgr2 cgs / r20)
return $ cgx . snd $ result
```
vec2arr :: Shape sh => (V.Vector Double -> V.Vector Double) -> Array U sh Double -> Array U sh Double
vec2arr f arr = fromUnboxed (extent arr) . f . toUnboxed $ arr
%% Cell type:code id: tags:
volterra :: Double -> Array U DIM1 Double -> Array U DIM1 Double
volterra h = vec2arr $ V.postscanl' (\s x -> s+h*x) 0
``` haskell
-- Volterra operator for a given stepsize, defined on V.Vectors
volterraOp :: Double -> V.Vector Double -> V.Vector Double
volterraOp h = V.map (* h) . V.tail . V.scanl (+) 0
-- Function that takes an operation on V.Vectors and returns the
-- corresponding operation on Arrays
wrapVectorOp :: (V.Vector Double -> V.Vector Double) -> Arr U -> Arr D
wrapVectorOp op array = delay . fromUnboxed (extent array) . op . toUnboxed $ array
-- Volterra operator type and instance
data Volterra = Volterra Double
instance Operator Volterra where
evalOp (Volterra h) = wrapVectorOp (volterraOp h)
-- The adjoint is this case is equivalent to reversing the matrix both up-down and left-right,
-- which is the same as reversing the input and output vectors
adjointOp (Volterra h) = wrapVectorOp (V.reverse . (volterraOp h) . V.reverse)
volterraAdjoint :: Double -> Array U DIM1 Double -> Array U DIM1 Double
volterraAdjoint h = vec2arr $ V.postscanr' (\x s -> s+h*x) 0
```
%% Cell type:code id: tags:
``` haskell
-- This was not required by the exercise, but is nice to have: test the
-- implementation of the adjoint by checking whether y^T A x = (A^T y)^T x
-- is fulfilled numerically for 50 random vectors of a given dimension.
-- Note that randomishDoubleArray is yet another way to generate random numbers.
import Data.Array.Repa.Algorithms.Randomish
testAdjoint :: Operator a => a -> DIM1 -> Double
testAdjoint op shape = maximum $ Prelude.map testAdjoint' [1..50]
testAdjoint :: (Array U DIM1 Double -> Array U DIM1 Double) -> (Array U DIM1 Double -> Array U DIM1 Double) -> DIM1 -> Double
testAdjoint operator adjoint shape = maximum $ Prelude.map testAdjoint' [1..50]
where testAdjoint' seed = abs $ sumAllS (y *^ fx) - sumAllS (fty *^ x)
where x = randomishDoubleArray shape (-1) 1 seed
fx = evalOp op x
fx = operator x
y = randomishDoubleArray (extent fx) (-1) 1 (seed+1)
fty = adjointOp op y
fty = adjoint y
testAdjoint (Volterra 1) (Z:.100)
testAdjoint (volterra 1) (volterraAdjoint 1) (Z:.100)
```
%% Cell type:code id: tags:
``` haskell
-- The function performing the main computations; returns the noisy function and
-- the regularizied solution for given function, sample points, noise level and
-- regularization parameter.
doDeriv :: (Double -> Double) -> Arr U -> Double -> Double -> IO (Arr U, Arr U)
doDeriv :: (Double -> Double) -> Array U DIM1 Double -> Double -> Double -> IO (Array U DIM1 Double, Array U DIM1 Double)
doDeriv f xs noiselevel reg = do
let func = computeS $ R.map f xs :: Arr U
let noise = randomishDoubleArray (extent func) (-noiselevel) noiselevel 0 :: Arr U
let rhs = computeS $ func +^ noise :: Arr U
let func = computeUnboxedS $ R.map f xs
let noise = randomishDoubleArray (extent func) (-noiselevel) noiselevel 0
let rhs = computeUnboxedS $ func +^ noise
let h = (xs ! (Z:.1)) - (xs ! (Z:.0))
print h
deriv <- runCG 1e-6 (Volterra h) reg rhs
deriv <- runCGreg 1e-6 (volterra h) (volterraAdjoint h) reg rhs
return (rhs, deriv)
```
%% Cell type:code id: tags:
``` haskell
linspace :: Int -> Double -> Double -> Arr U
linspace n a b = computeS $ fromFunction (Z:.n) $ \(Z:.i) ->
a + fromIntegral i / (fromIntegral n - 1) * (b - a)
xs = linspace 200 (-2) 2
```
%% Cell type:code id: tags:
``` haskell
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo
xs :: Array U DIM1 Double
xs = computeS $ linspace 200 (-2, 2)
exact = computeUnboxedS $ R.map (\x -> -2*x*exp(-x^2)) xs
(rhs, deriv) <- doDeriv (\x -> exp (-x^2)) xs 0.1 0.07
exact = computeS $ R.map (\x -> -2*x*exp(-x^2)) xs :: Arr U
toRenderable $ do
plot (line "noisy" [zip (toList xs) (toList rhs)])
plot (line "exact" [zip (toList xs) (toList exact)])
plot (line "deriv" [zip (toList xs) (toList deriv)])
```
%% Cell type:code id: tags:
``` haskell
(rhs, deriv) <- doDeriv signum xs 0.1 0.05
toRenderable $ do
plot (line "noisy" [zip (toList xs) (toList rhs)])
plot (line "derivative" [zip (toList xs) (toList deriv)])
```
%% Cell type:code id: tags:
``` haskell
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment