Skip to content

Commit

Permalink
split out new package diagrams-solve
Browse files Browse the repository at this point in the history
Everything formerly in `Diagrams.Solve` is now in
`Diagrams.Solve.Polynomial` in the `diagrams-solve` package; the
tridiagonal solving code in `Diagrams.CubicSpline.Internal` was moved to
`Diagrams.Solve.Tridiagonal`.

Closes #235.
  • Loading branch information
byorgey committed Feb 26, 2015
1 parent 1877e88 commit fa0834b
Show file tree
Hide file tree
Showing 7 changed files with 44 additions and 287 deletions.
2 changes: 1 addition & 1 deletion diagrams-lib.cabal
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@ Library
Diagrams.Query,
Diagrams.Segment,
Diagrams.Size,
Diagrams.Solve,
Diagrams.Tangent,
Diagrams.ThreeD,
Diagrams.ThreeD.Align,
Expand Down Expand Up @@ -101,6 +100,7 @@ Library
monoid-extras >= 0.3 && < 0.4,
dual-tree >= 0.2 && < 0.3,
diagrams-core >= 1.2 && < 1.3,
diagrams-solve >= 0.1 && < 0.2,
active >= 0.1 && < 0.2,
colour >= 2.3.2 && < 2.4,
data-default-class < 0.1,
Expand Down
66 changes: 3 additions & 63 deletions src/Diagrams/CubicSpline/Internal.hs
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
-----------------------------------------------------------------------------
-- |
-- Module : Diagrams.CubicSpline
Expand All @@ -17,71 +14,14 @@
module Diagrams.CubicSpline.Internal
(
-- * Solving for spline coefficents
solveTriDiagonal
, solveCyclicTriDiagonal
, solveCubicSplineDerivatives
solveCubicSplineDerivatives
, solveCubicSplineDerivativesClosed
, solveCubicSplineCoefficients
) where

import Data.List

-- | Solves a system of the form 'A*X=D' for 'x' where 'A' is an
-- 'n' by 'n' matrix with 'bs' as the main diagonal and
-- 'as' the diagonal below and 'cs' the diagonal above.
-- See: <http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm>
solveTriDiagonal :: Fractional a => [a] -> [a] -> [a] -> [a] -> [a]
solveTriDiagonal as (b0:bs) (c0:cs) (d0:ds) = h cs' ds'
where
cs' = c0 / b0 : f cs' as bs cs
f _ [_] _ _ = []
f (c':cs') (a:as) (b:bs) (c:cs) = c / (b - c' * a) : f cs' as bs cs
f _ _ _ _ = error "solveTriDiagonal.f: impossible!"

ds' = d0 / b0 : g ds' as bs cs' ds
g _ [] _ _ _ = []
g (d':ds') (a:as) (b:bs) (c':cs') (d:ds) = (d - d' * a)/(b - c' * a) : g ds' as bs cs' ds
g _ _ _ _ _ = error "solveTriDiagonal.g: impossible!"
import Diagrams.Solve.Tridiagonal

h _ [d] = [d]
h (c:cs) (d:ds) = let xs@(x:_) = h cs ds in d - c * x : xs
h _ _ = error "solveTriDiagonal.h: impossible!"

solveTriDiagonal _ _ _ _ = error "arguments 2,3,4 to solveTriDiagonal must be nonempty"

-- Helper that applies the passed function only to the last element of a list
modifyLast :: (a -> a) -> [a] -> [a]
modifyLast _ [] = []
modifyLast f [a] = [f a]
modifyLast f (a:as) = a : modifyLast f as

-- Helper that builds a list of length n of the form: '[s,m,m,...,m,m,e]'
sparseVector :: Int -> a -> a -> a -> [a]
sparseVector n s m e
| n < 1 = []
| otherwise = s : h (n - 1)
where
h 1 = [e]
h n = m : h (n - 1)

-- | Solves a system similar to the tri-diagonal system using a special case
-- of the Sherman-Morrison formula <http://en.wikipedia.org/wiki/Sherman-Morrison_formula>.
-- This code is based on /Numerical Recpies in C/'s @cyclic@ function in section 2.7.
solveCyclicTriDiagonal :: Fractional a => [a] -> [a] -> [a] -> [a] -> a -> a -> [a]
solveCyclicTriDiagonal as (b0:bs) cs ds alpha beta = zipWith ((+) . (fact *)) zs xs
where
l = length ds
gamma = -b0
us = sparseVector l gamma 0 alpha

bs' = (b0 - gamma) : modifyLast (subtract (alpha*beta/gamma)) bs

xs@(x:_) = solveTriDiagonal as bs' cs ds
zs@(z:_) = solveTriDiagonal as bs' cs us

fact = -(x + beta * last xs / gamma) / (1.0 + z + beta * last zs / gamma)

solveCyclicTriDiagonal _ _ _ _ _ _ = error "second argument to solveCyclicTriDiagonal must be nonempty"
import Data.List

-- | Use the tri-diagonal solver with the appropriate parameters for an open cubic spline.
solveCubicSplineDerivatives :: Fractional a => [a] -> [a]
Expand Down
12 changes: 6 additions & 6 deletions src/Diagrams/Segment.hs
Original file line number Diff line number Diff line change
Expand Up @@ -63,23 +63,23 @@ module Diagrams.Segment

) where

import Control.Lens (Rewrapped, Wrapped (..), iso, makeLenses, op,
over, Each (..))
import Control.Lens (Each (..), Rewrapped, Wrapped (..),
iso, makeLenses, op, over)
import Data.FingerTree
import Data.Monoid.MList
import Data.Semigroup
import Numeric.Interval.Kaucher (Interval (..))
import qualified Numeric.Interval.Kaucher as I
import Numeric.Interval.Kaucher (Interval (..))
import qualified Numeric.Interval.Kaucher as I

import Linear.Affine
import Linear.Metric
import Linear.Vector

import Control.Applicative
import Diagrams.Core hiding (Measured)
import Diagrams.Core hiding (Measured)
import Diagrams.Located
import Diagrams.Parametric
import Diagrams.Solve
import Diagrams.Solve.Polynomial


------------------------------------------------------------
Expand Down
187 changes: 0 additions & 187 deletions src/Diagrams/Solve.hs

This file was deleted.

6 changes: 3 additions & 3 deletions src/Diagrams/ThreeD/Shapes.hs
Original file line number Diff line number Diff line change
Expand Up @@ -23,16 +23,16 @@ module Diagrams.ThreeD.Shapes
) where

import Control.Applicative
import Control.Lens (review, (^.), _1)
import Control.Lens (review, (^.), _1)
import Data.Typeable

import Data.Semigroup
import Diagrams.Angle
import Diagrams.Core
import Diagrams.Solve
import Diagrams.Points
import Diagrams.Solve.Polynomial
import Diagrams.ThreeD.Types
import Diagrams.ThreeD.Vector
import Diagrams.Points

import Linear.Affine
import Linear.Metric
Expand Down
Loading

0 comments on commit fa0834b

Please sign in to comment.