Commit 589fd8cd authored by sof's avatar sof
Browse files

[project @ 1997-06-05 23:08:40 by sof]

new regression test
parent 46cc1b32
\begin{code}
module Arr (
module Array,
safezipWith, safezip,
row,
sum1, map2, map3,
mapat, mapat2, mapat3,
mapindexed, mapindexed2, mapindexed3,
-- zipArr, sumArr, scaleArr,
arraySize,
matvec, inner,
outerVector,
Vector (Vector), toVector, fromVector, listVector, vectorList, vector,
zipVector, scaleVector, sumVector, vectorNorm2, vectorSize,
Matrix (Matrix), toMatrix, fromMatrix, listMatrix, matrixList, matrix,
zipMatrix, scaleMatrix, sumMatrix,
augment,
trMatrix,
-- showsVector,
-- showsMatrix,
-- showsVecList, showsMatList
-- spy,
) where
import Array
import Numeric
--import Trace
--import IOExtensions(unsafePerformIO)
\end{code}
@Vector@ and @Matrix@ are 1-based arrays with read/show in form of Lists.
\begin{code}
data Vector a = Vector (Array Int a) deriving (Eq) --, Show)
toVector :: Array Int a -> Vector a
toVector x = Vector x
fromVector :: Vector a -> Array Int a
fromVector (Vector x) = x
instance Functor (Vector) where
map fn x = toVector (map fn (fromVector x))
{-instance Eq a => Eq (Vector a) where
-- (Vector x) == (Vector y) = x == y
-}
instance Show a => Show (Vector a) where
showsPrec p x = showsPrec p (elems (fromVector x))
instance Read a => Read (Vector a) where
readsPrec p = readParen False
(\r -> [(listVector s, t) | (s, t) <- reads r])
instance Num b => Num (Vector b) where
(+) = zipVector "+" (+)
(-) = zipVector "-" (-)
negate = map negate
abs = map abs
signum = map signum
-- (*) = matMult -- works only for matrices!
-- fromInteger = map fromInteger
\end{code}
Convert a list to 1-based vector.
\begin{code}
listVector :: [a] -> Vector a
listVector x = toVector (listArray (1,length x) x)
vectorList :: Vector a -> [a]
vectorList = elems . fromVector
vector (l,u) x | l == 1 = toVector (array (l,u) x)
| otherwise = error "vector: l != 1"
zipVector :: String -> (b -> c -> d) -> Vector b -> Vector c -> Vector d
zipVector s f (Vector a) (Vector b)
| bounds a == bounds b = vector (bounds a) [(i, f (a!i) (b!i)) | i <- indices a]
| otherwise = error ("zipVector: " ++ s ++ ": unconformable arrays")
scaleVector :: Num a => a -> Vector a -> Vector a
scaleVector a = map (* a)
sumVector :: Num a => Vector a -> a
sumVector = sum . elems . fromVector
vectorNorm2 :: Num a => Vector a -> a
vectorNorm2 x = inner x x
vectorSize :: Vector a -> Int
vectorSize (Vector x) = rangeSize (bounds x)
\end{code}
==============
\begin{code}
data Matrix a = Matrix (Array (Int, Int) a) deriving Eq
toMatrix :: Array (Int, Int) a -> Matrix a
toMatrix x = Matrix x
fromMatrix :: Matrix a -> Array (Int, Int) a
fromMatrix (Matrix x) = x
instance Functor (Matrix) where
map fn x = toMatrix (map fn (fromMatrix x))
--instance Eq a => Eq (Matrix a) where
-- (Matrix x) == (Matrix y) = x == y
instance Show a => Show (Matrix a) where
showsPrec p x = vertl (matrixList x)
vertl [] = showString "[]"
vertl (x:xs) = showChar '[' . shows x . showl xs
where showl [] = showChar ']'
showl (x:xs) = showString ",\n" . shows x . showl xs
instance Read a => Read (Matrix a) where
readsPrec p = readParen False
(\r -> [(listMatrix s, t) | (s, t) <- reads r])
instance Num b => Num (Matrix b) where
(+) = zipMatrix "+" (+)
(-) = zipMatrix "-" (-)
negate = map negate
abs = map abs
signum = map signum
x * y = toMatrix (matMult (fromMatrix x) (fromMatrix y)) -- works only for matrices!
-- fromInteger = map fromInteger
\end{code}
Convert a nested list to a matrix.
\begin{code}
listMatrix :: [[a]] -> Matrix a
listMatrix x = Matrix (listArray ((1, 1), (length x, length (x!!0))) (concat x))
matrixList :: Matrix a -> [[a]]
matrixList (Matrix x) = [ [x!(i,j) | j <- range (l',u')] | i <- range (l,u)]
where ((l,l'),(u,u')) = bounds x
matrix ((l,l'),(u,u')) x | l == 1 && l' == 1 = toMatrix (array ((l,l'),(u,u')) x)
| otherwise = error "matrix: l != 1"
zipMatrix :: String -> (b -> c -> d) -> Matrix b -> Matrix c -> Matrix d
zipMatrix s f (Matrix a) (Matrix b)
| bounds a == bounds b = matrix (bounds a) [(i, f (a!i) (b!i)) | i <- indices a]
| otherwise = error ("zipMatrix: " ++ s ++ ": unconformable arrays")
scaleMatrix :: Num a => a -> Matrix a -> Matrix a
scaleMatrix a = map (* a)
sumMatrix :: Num a => Matrix a -> a
sumMatrix = sum . elems . fromMatrix
\end{code}
============
\begin{code}
safezipWith :: String -> (a -> b -> c) -> [a] -> [b] -> [c]
safezipWith _ _ [] [] = []
safezipWith s f (x:xs) (y:ys) = f x y : safezipWith s f xs ys
safezipWith s _ _ _ = error ("safezipWith: " ++ s ++ ": unconformable vectors")
safezip :: [a] -> [b] -> [(a,b)]
safezip = safezipWith "(,)" (,)
trMatrix :: Matrix a -> Matrix a
trMatrix (Matrix x) = matrix ((l,l'),(u',u)) [((j,i), x!(i,j)) | j <- range (l',u'), i <- range (l,u)]
where ((l,l'),(u,u')) = bounds x
row :: (Ix a, Ix b) => a -> Array (a,b) c -> Array b c
row i x = ixmap (l',u') (\j->(i,j)) x where ((l,l'),(u,u')) = bounds x
zipArr :: (Ix a) => String -> (b -> c -> d) -> Array a b -> Array a c -> Array a d
zipArr s f a b | bounds a == bounds b = array (bounds a) [(i, f (a!i) (b!i)) | i <- indices a]
| otherwise = error ("zipArr: " ++ s ++ ": unconformable arrays")
\end{code}
Valid only for b -> c -> b functions.
\begin{code}
zipArr' :: (Ix a) => String -> (b -> c -> b) -> Array a b -> Array a c -> Array a b
zipArr' s f a b | bounds a == bounds b = accum f a (assocs b)
| otherwise = error ("zipArr': " ++ s ++ ": unconformable arrays")
\end{code}
Overload arithmetical operators to work on lists.
\begin{code}
instance Num a => Num [a] where
(+) = safezipWith "+" (+)
(-) = safezipWith "-" (-)
negate = map negate
abs = map abs
signum = map signum
-- (*) = undefined
-- fromInteger = undefined
\end{code}
\begin{code}
sum1 :: (Num a) => [a] -> a
sum1 = foldl1 (+)
--main = print (sum1 [[4,1,1], [5,1,2], [6,1,3,4]])
\end{code}
\begin{code}
map2 f = map (map f)
map3 f = map (map2 f)
\end{code}
Map function f at position n only. Out of range indices are silently
ignored.
\begin{code}
mapat n f x = mapat1 0 f x where
mapat1 _ _ [] = []
mapat1 i f (x:xs) = (if i == n then f x else x) : mapat1 (i + 1) f xs
mapat2 (i,j) = mapat i . mapat j
mapat3 (i,j,k) = mapat i . mapat j . mapat k
-- main = print (mapat 2 (10+) [1,2,3,4])
-- main = print (mapat2 (1,0) (1000+) ginp)
-- main = print (mapat3 (1,0,1) (1000+) gw)
\end{code}
\begin{code}
mapindexed f x = mapindexed1 f 0 x where
mapindexed1 _ _ [] = []
mapindexed1 f n (x:xs) = f n x : mapindexed1 f (n + 1) xs
mapindexed2 f = mapindexed (\i -> mapindexed (\j -> f (i, j)))
mapindexed3 f = mapindexed (\i -> mapindexed (\j -> mapindexed (\k -> f (i, j, k))))
-- main = print (mapindexed (\x y -> mapat (10+) [1,2,3,4] y) [1,2,3,4])
-- main = print (mapindexed2 (\(i,j) x -> 100*i + 10*j + x) ginp)
-- main = print (mapindexed3 (\(i,j,k) x -> 1000*i + 100*j + 10*k + x) gw)
\end{code}
Overload arithmetical operators to work on arrays.
\begin{code}
instance (Ix a, Show a, Num b) => Num (Array a b) where
(+) = zipArr "+" (+)
(-) = zipArr "-" (-)
negate = map negate
abs = map abs
signum = map signum
-- (*) = matMult -- works only for matrices!
-- fromInteger = map fromInteger
\end{code}
\begin{xcode}
scaleArr :: (Ix i, Num a) => a -> Array i a -> Array i a
scaleArr a = map (*a)
sumArr :: (Ix i, Num a) => Array i a -> a
sumArr = sum . elems
\end{xcode}
\begin{code}
arraySize :: (Ix i) => Array i a -> Int
arraySize = rangeSize . bounds
\end{code}
\begin{code}
matMult :: (Ix a, Ix b, Ix c, Num d) =>
Array (a,b) d -> Array (b,c) d -> Array (a,c) d
matMult x y = array resultBounds
[((i,j), sum [x!(i,k) * y!(k,j) | k <- range (lj,uj)])
| i <- range (li,ui),
j <- range (lj',uj') ]
where ((li,lj),(ui,uj)) = bounds x
((li',lj'),(ui',uj')) = bounds y
resultBounds
| (lj,uj)==(li',ui') = ((li,lj'),(ui,uj'))
| otherwise = error "matMult: incompatible bounds"
\end{code}
Inner product of two vectors.
\begin{code}
inner :: Num a => Vector a -> Vector a -> a
inner (Vector v) (Vector w) = if b == bounds w
then sum [v!i * w!i | i <- range b]
else error "nn.inner: inconformable vectors"
where b = bounds v
\end{code}
Outer product of two vectors $v \dot w^\mathrm{T}$.
\begin{code}
outerVector :: Num b => Vector b -> Vector b -> Matrix b
outerVector (Vector v) (Vector w) = if (l,u) == (l',u')
then matrix ((l,l'),(u,u')) [((i,j), v!i * w!j) | i <- range (l,u), j <- range (l',u')]
else error "nn.outer: inconformable vectors"
where ((l,u),(l',u')) = (bounds v, bounds w)
\end{code}
\begin{code}
outerArr :: (Ix a, Num b) => Array a b -> Array a b -> Array (a,a) b
outerArr v w = if (l,u) == (l',u')
then array ((l,l'),(u,u')) [((i,j), v!i * w!j) | i <- range (l,u), j <- range (l',u')]
else error "nn.outer: inconformable vectors"
where ((l,u),(l',u')) = (bounds v, bounds w)
\end{code}
Inner product of a matrix and a vector.
\begin{code}
matvec :: (Ix a, Num b) => Array (a,a) b -> Array a b -> Array a b
matvec w x | bounds x == (l',u') =
array (l,u) [(i, sum [w!(i,j) * x!j | j <- range (l',u')])
| i <- range (l,u)]
| otherwise = error "nn.matvec: inconformable arrays"
where ((l,l'),(u,u')) = bounds w
\end{code}
Append to a vector.
\begin{code}
augment :: (Num a) => Vector a -> a -> Vector a
augment (Vector x) y = Vector (array (a,b') ((b',y) : assocs x))
where (a,b) = bounds x
b' = b + 1
\end{code}
Older approach (x!!i!!j fails in ghc-2.03).
\begin{code}
toMatrix' :: [[a]] -> Matrix a
toMatrix' x = Matrix (array ((1,1),(u,u')) [((i,j), (x!!(i-1))!!(j-1))
| i <- range (1,u), j <- range (1,u')])
where (u,u') = (length x,length (x!!0))
\end{code}
Matrix 2D printout.
\begin{code}
padleft :: Int -> String -> String
padleft n x | n <= length x = x
| otherwise = replicate (n - length x) ' ' ++ x
\end{code}
\begin{code}
padMatrix :: RealFloat a => Int -> Matrix a -> Matrix String
padMatrix n x = let ss = map (\a -> showFFloat (Just n) a "") x
maxw = maximum (map length (elems (fromMatrix ss)))
in map (padleft maxw) ss
\end{code}
\begin{xcode}
showsVector :: (RealFloat a) => Int -> Vector a -> ShowS
showsVector n x z1 = let x' = padArr n x
(l,u) = bounds x' in
concat (map (\ (i, s) -> if i == u then s ++ "\n" else s ++ " ") (assocs x')) ++ z1
\end{xcode}
\begin{xcode}
showsMatrix :: RealFloat a => Int -> Matrix a -> ShowS
showsMatrix n x z1 = let x' = padMatrix n x
((l,l'),(u,u')) = bounds x' in
concat (map (\ ((i,j), s) -> if j == u' then s ++ "\n" else s ++ " ") (assocs x')) ++ z1
\end{xcode}
{-
showsVecList n x s = foldr (showsVector n) s x
showsMatList n x s = foldr (showsMatrix n) s x
-}
\begin{code}
--spy :: Show a => String -> a -> a
--spy msg x = trace ('<' : msg ++ ": " ++ shows x ">\n") x
--spy x = seq (unsafePerformIO (putStr ('<' : shows x ">\n"))) x
--spy x = traceShow "z" x
\end{code}
\begin{code}
module Main(main) where
import Arr
\end{code}
\begin{code}
type F a = Vector a -> a
type DF a = Vector a -> Vector a
\end{code}
\begin{code}
data (Eval a) => ScgData a = ScgData {k :: !Int, err :: !a,
w, p, r :: !(Vector a),
delta, pnorm2, lambda, lambdabar :: !a,
success :: !Bool}
\end{code}
\begin{code}
calculate2order :: Floating a => ScgData a -> a -> DF a -> ScgData a
calculate2order d sigma1 df =
let pnorm2' = vectorNorm2 (p d)
sigma = sigma1 / (sqrt pnorm2')
s = scaleVector (1/sigma) (df ((w d) + (scaleVector sigma (p d))) - df (w d))
in d {pnorm2 = pnorm2', delta = inner (p d) s}
\end{code}
\begin{code}
hessPosDef :: (Floating a, Ord a) => ScgData a -> ScgData a
hessPosDef d =
let delta' = delta d + (lambda d - lambdabar d) * pnorm2 d {- step 3 -}
in if delta' <= 0 {- step 4 -}
then let lambdabar' = 2.0 * (lambda d - delta' / pnorm2 d)
in d {delta = -delta' + lambda d * pnorm2 d, lambda = lambdabar', lambdabar = lambdabar'}
else d {delta = delta'}
\end{code}
\begin{code}
reduceError :: (Floating a, Ord a) => ScgData a -> DF a -> Bool -> a -> a -> ScgData a
reduceError d df restart bdelta mu =
let r' = negate (df (w d))
p' = if restart
then r'
else let beta = (vectorNorm2 r' - inner r' (r d)) / mu
in r' + scaleVector beta (p d)
in d {p = p', r = r', lambda = if bdelta >= 0.75 then lambda d / 4 else lambda d
}
\end{code}
\begin{code}
data ScgInput a = ScgInput (F a) (DF a) (Vector a)
\end{code}
\begin{code}
scgIter :: (Floating a, Ord a) => ScgInput a -> [ScgData a]
scgIter (ScgInput f df w1) =
let p1 = negate (df w1) {- step 1 -}
r1 = p1
pnorm21 = vectorNorm2 p1
n = vectorSize w1
sigma1 = 1.0e-4
lambda1 = 1.0e-6
err1 = f w1
in iterate (\d ->
let d1 = if success d {- step 2 -}
then calculate2order d sigma1 df
else d
d2 = hessPosDef d1
mu = inner (p d2) (r d2) {- step 5 -}
alpha = mu / (delta d2)
w' = (w d2) + scaleVector alpha (p d2)
err' = f w'
bdelta = 2 * (delta d2) * ((err d2) - err') / (mu^2) {- step 6 -}
success' = (bdelta >= 0) {- step 7 -}
restart = ((k d) `mod` n == 0)
d3 = if success'
then (reduceError (d2 {w = w'}) df restart bdelta mu)
{err = err', lambdabar = 0}
else d2 {lambdabar = lambda d2}
d4 = if bdelta < 0.25 {- step 8 -}
then d3 {lambda = (lambda d3) + (delta d3) * (1 - bdelta) / (pnorm2 d3)}
else d3
in d4 {k = k d4 + 1, success = success'}
)
(ScgData 1 err1 w1 p1 r1 0.0 pnorm21 lambda1 0.0 True)
\end{code}
\begin{code}
rosenbrock = ScgInput
(\ (Vector x) -> 100 * (x!2 - x!1^2)^2 + (1 - x!1)^2)
(\ (Vector x) -> listVector [-2 * (1 - x!1) - 400 * x!1 * (x!2 - x!1^2),
200 * (x!2 -x!1^2)])
(listVector [-1,-1.0])
\end{code}
\begin{code}
main = print (w ((scgIter rosenbrock)!!1))
\end{code}
{-
Date: Tue, 20 May 1997 05:10:04 GMT
From: Tomasz Cholewo <tjchol01@mecca.spd.louisville.edu>
ghc-2.03 cannot compile the following code, which I think is correct
according to the Report
data X = A {a :: Int} | B {a :: Int}
The error message is:
Conflicting definitions for: a
Defined at bug4.lhs:2
Defined at bug4.lhs:2
In addition the following snippet
data X = A {a :: Int}
y = let A {a} = x
in a
fails with:
bug4.lhs:4:5: Not a valid LHS on input: "in"
-}
--module Main(main) where
data X = A {a :: Int} | B {a :: Int}
f x = let A {a} = x
in a
main = print (f (A {a = 3}))
-- Test newtype derived instances
newtype Age = MkAge Int deriving (Eq, Show)
instance Num Age where
(+) (MkAge a) (MkAge b) = MkAge (a+b)
main = print (MkAge 3 + MkAge 5)
{- Tests let-expressions in do-statments -}
module Main( main ) where
foo = do
putStr "a"
let x = "b" in putStr x
putStr "c"
main = do
putStr "a"
foo
let x = "b" in putStrLn x
--!!! Empty comments terminating a file..
main = print "Hello" --
\ No newline at end of file
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