{-# LANGUAGE BangPatterns #-}
module Numeric.Integration.SimplexCubature.Internal
where
import Data.Array.IO (IOUArray, getBounds, getElems
, mapIndices, newArray, newArray_
, newListArray, readArray
, writeArray)
import Data.List (foldl', foldl1', sort)
import Control.Monad ((<$!>), when)
import qualified Data.Foldable as DF
import Data.Maybe (fromJust)
import Data.Sequence (Seq, index, update, (><))
import qualified Data.Sequence as S
import Data.Vector.Unboxed (Vector, fromList, toList
, unsafeFreeze)
import qualified Data.Vector.Unboxed as UV
import Data.Vector.Unboxed.Mutable (IOVector, new, unsafeRead
, unsafeWrite)
import qualified Data.Vector.Unboxed.Mutable as UMV
import Numeric.Integration.Simplex.Simplex (Simplex, simplexVolume)
type IOMatrix = IOUArray (Int, Int) Double
type IO1dArray = IOUArray Int Double
type IO3dArray = IOUArray (Int, Int, Int) Double
type IOVectorD = IOVector Double
type IOVectorI = IOVector Int
type VectorD = Vector Double
toDbl :: Int -> Double
toDbl :: Int -> Double
toDbl = forall a b. (Integral a, Num b) => a -> b
fromIntegral
pow :: Double -> Int -> Double
pow :: Double -> Int -> Double
pow Double
x Int
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a. Int -> a -> [a]
replicate Int
n Double
x)
permuteMultiset :: (Eq a, Ord a) => [a] -> [[a]]
permuteMultiset :: forall a. (Eq a, Ord a) => [a] -> [[a]]
permuteMultiset [a]
list = forall a. (a -> Maybe a) -> a -> [a]
unfold1 forall {a}. Ord a => [a] -> Maybe [a]
next (forall a. Ord a => [a] -> [a]
sort [a]
list) where
unfold1 :: (a -> Maybe a) -> a -> [a]
unfold1 :: forall a. (a -> Maybe a) -> a -> [a]
unfold1 a -> Maybe a
f a
x = case a -> Maybe a
f a
x of
Maybe a
Nothing -> [a
x]
Just a
y -> a
x forall a. a -> [a] -> [a]
: forall a. (a -> Maybe a) -> a -> [a]
unfold1 a -> Maybe a
f a
y
next :: [a] -> Maybe [a]
next [a]
xs = case forall {a}. Ord a => ([a], [a]) -> Maybe ([a], [a])
findj (forall a. [a] -> [a]
reverse [a]
xs, []) of
Maybe ([a], [a])
Nothing -> forall a. Maybe a
Nothing
Just ( a
l:[a]
ls, [a]
rs) -> forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall {a}. Ord a => a -> [a] -> ([a], [a]) -> [a]
inc a
l [a]
ls (forall a. [a] -> [a]
reverse [a]
rs, [])
Just ( [], [a]
_ ) -> forall a. HasCallStack => [Char] -> a
error [Char]
"permute: should not happen"
findj :: ([a], [a]) -> Maybe ([a], [a])
findj ( xxs :: [a]
xxs@(a
x:[a]
xs), yys :: [a]
yys@(a
y:[a]
_) ) = if a
x forall a. Ord a => a -> a -> Bool
>= a
y
then ([a], [a]) -> Maybe ([a], [a])
findj ( [a]
xs, a
x forall a. a -> [a] -> [a]
: [a]
yys )
else forall a. a -> Maybe a
Just ( [a]
xxs, [a]
yys )
findj ( a
x:[a]
xs , [] ) = ([a], [a]) -> Maybe ([a], [a])
findj ( [a]
xs , [a
x] )
findj ( [] , [a]
_ ) = forall a. Maybe a
Nothing
inc :: a -> [a] -> ([a], [a]) -> [a]
inc !a
u [a]
us ( a
x:[a]
xs, [a]
yys ) = if a
u forall a. Ord a => a -> a -> Bool
>= a
x
then a -> [a] -> ([a], [a]) -> [a]
inc a
u [a]
us ( [a]
xs, a
x forall a. a -> [a] -> [a]
: [a]
yys )
else forall a. [a] -> [a]
reverse (a
xforall a. a -> [a] -> [a]
:[a]
us) forall a. [a] -> [a] -> [a]
++ forall a. [a] -> [a]
reverse (a
uforall a. a -> [a] -> [a]
:[a]
yys) forall a. [a] -> [a] -> [a]
++ [a]
xs
inc a
_ [a]
_ ( [], [a]
_ ) = forall a. HasCallStack => [Char] -> a
error [Char]
"permute: should not happen"
square :: Double -> Double
square :: Double -> Double
square Double
x = Double
xforall a. Num a => a -> a -> a
*Double
x
smprms :: Int -> Int -> IO (IOMatrix, Seq VectorD, [Int])
smprms :: Int -> Int -> IO (IOMatrix, Seq VectorD, [Int])
smprms Int
n Int
key = do
let (Int
rls, Int
gms, Int
wts) | Int
key forall a. Eq a => a -> a -> Bool
== Int
1 = (Int
3, Int
2, Int
3) :: (Int, Int, Int)
| Int
key forall a. Eq a => a -> a -> Bool
== Int
2 = (Int
5, Int
4, Int
6) :: (Int, Int, Int)
| Int
key forall a. Eq a => a -> a -> Bool
== Int
3 = (Int
7, Int
7, Int
11) :: (Int, Int, Int)
| Int
key forall a. Eq a => a -> a -> Bool
== Int
4 = (if Int
n forall a. Eq a => a -> a -> Bool
== Int
2 then (Int
7, Int
11, Int
20) else (Int
7, Int
12, Int
21))
:: (Int, Int, Int)
| Bool
otherwise = forall a. HasCallStack => [Char] -> a
error [Char]
"this should not happen"
IOMatrix
w <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
wts, Int
rls)) Double
0 :: IO IOMatrix
IOVectorI
pts <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
UMV.replicate Int
wts Int
0 :: IO IOVectorI
IOMatrix
g <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
nforall a. Num a => a -> a -> a
+Int
1, Int
wts)) Double
0 :: IO IOMatrix
let np :: Int
np = Int
nforall a. Num a => a -> a -> a
+Int
1
n2 :: Int
n2 = Int
np forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
2)
n4 :: Int
n4 = Int
n2 forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
3) forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
4)
n6 :: Int
n6 = Int
n4 forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
5) forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
6)
n8 :: Int
n8 = Int
n6 forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
7) forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
8)
o :: Int
o = forall a. Integral a => a -> a -> a
div (Int
nforall a. Num a => a -> a -> a
*Int
np) Int
2
ndbl :: Double
ndbl = Int -> Double
toDbl Int
n
sqrt15 :: Double
sqrt15 = Double
3.872983346207416885179265399782399611
r1 :: Double
r1 = (Double
ndbl forall a. Num a => a -> a -> a
+ Double
4 forall a. Num a => a -> a -> a
- Double
sqrt15) forall a. Fractional a => a -> a -> a
/ (Double
ndblforall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
+ Double
8forall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
+ Double
1)
s1 :: Double
s1 = Double
1 forall a. Num a => a -> a -> a
- Double
ndblforall a. Num a => a -> a -> a
*Double
r1
l1 :: Double
l1 = Double
s1 forall a. Num a => a -> a -> a
- Double
r1
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
1) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
1))) [Int
1 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
0 Int
1
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
1) Double
s1
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
1) Double
r1) [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
gms Int
np
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Ord a => a -> a -> Bool
< Int
4) forall a b. (a -> b) -> a -> b
$ do
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
1, Int
rls) Double
1
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
rlsforall a. Num a => a -> a -> a
-Int
1) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
1))
let iw :: Int
iw = if Int
key forall a. Ord a => a -> a -> Bool
< Int
4 then Int
rlsforall a. Num a => a -> a -> a
-Int
2 else Int
rls
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
2) (Double
3forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3))
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
2) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3))) [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
1 Int
np
let n2double :: Double
n2double = Int -> Double
toDbl Int
n2
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iw) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
2 forall a. Fractional a => a -> a -> a
/ (Double
4forall a. Num a => a -> a -> a
*Double
n2double))
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Ord a => a -> a -> Bool
> Int
1) forall a b. (a -> b) -> a -> b
$ do
if Int
n forall a. Eq a => a -> a -> Bool
== Int
2
then do
let l2 :: Double
l2 = Double
0.62054648267200632589046034361711
r1' :: Double
r1' = (Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt(Double
0.5 forall a. Num a => a -> a -> a
- Double
l2forall a. Num a => a -> a -> a
*Double
l2))forall a. Fractional a => a -> a -> a
/Double
3
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
1) (Double
1 forall a. Num a => a -> a -> a
- Double
2forall a. Num a => a -> a -> a
*Double
r1')
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
1) Double
r1') [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
gms Int
3
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
1) (Double
1forall a. Fractional a => a -> a -> a
/Double
6)
let r2 :: Double
r2 = (Double
1forall a. Num a => a -> a -> a
-Double
l2)forall a. Fractional a => a -> a -> a
/Double
3
s2 :: Double
s2 = Double
1 forall a. Num a => a -> a -> a
- Double
2forall a. Num a => a -> a -> a
*Double
r2
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
s2
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
r2) [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
1) Int
3
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
1) (Double
1forall a. Fractional a => a -> a -> a
/Double
6)
else do
let r2 :: Double
r2 = (Double
ndblforall a. Num a => a -> a -> a
+Double
4forall a. Num a => a -> a -> a
+Double
sqrt15) forall a. Fractional a => a -> a -> a
/ (Double
ndblforall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
8forall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
1)
s2 :: Double
s2 = Double
1 forall a. Num a => a -> a -> a
- Double
ndblforall a. Num a => a -> a -> a
*Double
r2
l2 :: Double
l2 = Double
s2 forall a. Num a => a -> a -> a
- Double
r2
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
s2
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
r2) [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
1) Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
1) ((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3)forall a. Num a => a -> a -> a
-Double
l1)forall a. Fractional a => a -> a -> a
/(Double
n2doubleforall a. Num a => a -> a -> a
*(Double
l2forall a. Num a => a -> a -> a
-Double
l1)forall a. Num a => a -> a -> a
*Double
l2forall a. Num a => a -> a -> a
*Double
l2))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
1) ((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3)forall a. Num a => a -> a -> a
-Double
l2)forall a. Fractional a => a -> a -> a
/(Double
n2doubleforall a. Num a => a -> a -> a
*(Double
l1forall a. Num a => a -> a -> a
-Double
l2)forall a. Num a => a -> a -> a
*Double
l1forall a. Num a => a -> a -> a
*Double
l1))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
3) (Double
5forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
3) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))) [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
2 Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
4) (Double
3forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
4) (Double
3forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
4) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))) [Int
3 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
3 Int
o
let tmp :: Double
tmp = Int -> Double
toDbl (Int
16forall a. Num a => a -> a -> a
*Int
n4)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
2) (- Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
5 forall a. Fractional a => a -> a -> a
/ Double
tmp)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
2) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
4 forall a. Fractional a => a -> a -> a
/ Double
tmp)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
2) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
4 forall a. Fractional a => a -> a -> a
/ Double
tmp)
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Ord a => a -> a -> Bool
> Int
2) forall a b. (a -> b) -> a -> b
$ do
let tmp' :: Double
tmp' = Double
ndblforall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
+ Double
14forall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
- Double
11
u1 :: Double
u1 = (Double
ndblforall a. Num a => a -> a -> a
+Double
7forall a. Num a => a -> a -> a
+Double
2forall a. Num a => a -> a -> a
*Double
sqrt15) forall a. Fractional a => a -> a -> a
/ Double
tmp'
v1 :: Double
v1 = Double
0.5forall a. Num a => a -> a -> a
*(Double
1forall a. Num a => a -> a -> a
-(Double
ndblforall a. Num a => a -> a -> a
-Double
1)forall a. Num a => a -> a -> a
*Double
u1)
d1 :: Double
d1 = Double
v1 forall a. Num a => a -> a -> a
- Double
u1
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
3) Double
v1
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
3) Double
v1
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
3) Double
u1) [Int
3 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
2) Int
o
let u2 :: Double
u2 = (Double
ndbl forall a. Num a => a -> a -> a
+ Double
7 forall a. Num a => a -> a -> a
- Double
2forall a. Num a => a -> a -> a
*Double
sqrt15) forall a. Fractional a => a -> a -> a
/ Double
tmp'
v2 :: Double
v2 = (Double
1forall a. Num a => a -> a -> a
-(Double
ndblforall a. Num a => a -> a -> a
-Double
1)forall a. Num a => a -> a -> a
*Double
u2)forall a. Fractional a => a -> a -> a
/Double
2
d2 :: Double
d2 = Double
v2 forall a. Num a => a -> a -> a
- Double
u2
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
4) Double
v2
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
4) Double
v2
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
4) Double
u2) [Int
3 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
3) Int
o
if Int
n forall a. Eq a => a -> a -> Bool
== Int
2
then do
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
155forall a. Num a => a -> a -> a
-Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
1200)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
155forall a. Num a => a -> a -> a
+Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
1200)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
3) (Double
9forall a. Fractional a => a -> a -> a
/Double
40)
else
if Int
n forall a. Eq a => a -> a -> Bool
== Int
3
then do
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
2665forall a. Num a => a -> a -> a
+Double
14forall a. Num a => a -> a -> a
*Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
37800)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
2665forall a. Num a => a -> a -> a
-Double
14forall a. Num a => a -> a -> a
*Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
37800)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
3) (Double
10forall a. Fractional a => a -> a -> a
/Double
189)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
3) Int
0
else do
let r2 :: Double
r2 = (Double
ndblforall a. Num a => a -> a -> a
+Double
4forall a. Num a => a -> a -> a
+Double
sqrt15) forall a. Fractional a => a -> a -> a
/ (Double
ndblforall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
8forall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
1)
l2 :: Double
l2 = Double
1 forall a. Num a => a -> a -> a
- (Double
ndblforall a. Num a => a -> a -> a
+Double
1)forall a. Num a => a -> a -> a
*Double
r2
n4dbl :: Double
n4dbl = Int -> Double
toDbl Int
n4
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
3)
((Double
2forall a. Num a => a -> a -> a
*(Double
27forall a. Num a => a -> a -> a
-Double
ndbl)forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
l2forall a. Num a => a -> a -> a
*(Double
13forall a. Num a => a -> a -> a
-Double
ndbl)) forall a. Fractional a => a -> a -> a
/
(Double -> Int -> Double
pow Double
l1 Int
4 forall a. Num a => a -> a -> a
* (Double
l1forall a. Num a => a -> a -> a
-Double
l2) forall a. Num a => a -> a -> a
* Double
n4dbl))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
3)
((Double
2forall a. Num a => a -> a -> a
*(Double
27forall a. Num a => a -> a -> a
-Double
ndbl)forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
l1forall a. Num a => a -> a -> a
*(Double
13forall a. Num a => a -> a -> a
-Double
ndbl)) forall a. Fractional a => a -> a -> a
/
(Double -> Int -> Double
pow Double
l2 Int
4 forall a. Num a => a -> a -> a
*(Double
l2forall a. Num a => a -> a -> a
-Double
l1)forall a. Num a => a -> a -> a
*Double
n4dbl))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
3)
((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
d2) forall a. Fractional a => a -> a -> a
/ (Double
n4dbl forall a. Num a => a -> a -> a
* (Double
d1forall a. Num a => a -> a -> a
-Double
d2) forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow Double
d1 Int
4))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
3)
((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
d1) forall a. Fractional a => a -> a -> a
/ (Double
n4dbl forall a. Num a => a -> a -> a
* (Double
d2forall a. Num a => a -> a -> a
-Double
d1) forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow Double
d2 Int
4))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
5) (Double
7forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7))
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
5) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7))) [Int
2 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
4 Int
np
let invnp7 :: Double
invnp7 = Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
6) (Double
5forall a. Num a => a -> a -> a
*Double
invnp7)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
6) (Double
3forall a. Num a => a -> a -> a
*Double
invnp7)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
6) Double
invnp7) [Int
3 .. Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
5 (Int
npforall a. Num a => a -> a -> a
*Int
n)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
7) (Double
3forall a. Num a => a -> a -> a
*Double
invnp7)) [Int
1, Int
2, Int
3]
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
np forall a. Ord a => a -> a -> Bool
> Int
3) forall a b. (a -> b) -> a -> b
$
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
7) Double
invnp7) [Int
4..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
6 (forall a. Integral a => a -> a -> a
div ((Int
nforall a. Num a => a -> a -> a
-Int
1)forall a. Num a => a -> a -> a
*Int
nforall a. Num a => a -> a -> a
*Int
np) Int
6)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
4) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
7 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
128forall a. Num a => a -> a -> a
*Int
n4forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
+Int
5)))
let tmp'' :: Double
tmp'' = Int -> Double
toDbl (Int
64forall a. Num a => a -> a -> a
*Int
n6)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
4) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
7 forall a. Fractional a => a -> a -> a
/ Double
tmp'')) [Int
3, Int
4]
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
4) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Double
tmp'')) [Int
5, Int
6, Int
7]
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Eq a => a -> a -> Bool
== Int
4) forall a b. (a -> b) -> a -> b
$ do
let sg :: Double
sg = Double
1 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
23328forall a. Num a => a -> a -> a
*Int
n6)
u5 :: Double
u5 = -Double
216 forall a. Num a => a -> a -> a
* Double
sg forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
52212 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*(Int
6353 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
1934forall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
*Int
27)))
u6 :: Double
u6 = Double
1296 forall a. Num a => a -> a -> a
* Double
sg forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
7884 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*(Int
1541 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
9))
u7 :: Double
u7 = -Double
7776 forall a. Num a => a -> a -> a
* Double
sg forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
8292 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*(Int
1139 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
3))forall a. Fractional a => a -> a -> a
/(Double
ndbl forall a. Num a => a -> a -> a
+ Double
7)
p0 :: Double
p0 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -Int
144 forall a. Num a => a -> a -> a
* (Int
142528 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
23073 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
115))
p1 :: Double
p1 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -Int
12 forall a. Num a => a -> a -> a
* (Int
6690556 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
2641189 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
245378 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
1495)))
p2 :: Double
p2 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -Int
16 forall a. Num a => a -> a -> a
* (Int
6503401 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
4020794forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
787281forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
47323forall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
*Int
385))))
p3 :: Double
p3 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -(Int
6386660 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
4411997forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
951821forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
61659forall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
*Int
665))))forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
+Int
7)
a :: Double
a = Double
p2forall a. Fractional a => a -> a -> a
/(Double
3forall a. Num a => a -> a -> a
*Double
p3)
p :: Double
p = Double
aforall a. Num a => a -> a -> a
*(Double
p1forall a. Fractional a => a -> a -> a
/Double
p2 forall a. Num a => a -> a -> a
- Double
a)
q :: Double
q = Double
aforall a. Num a => a -> a -> a
*(Double
2forall a. Num a => a -> a -> a
*Double
aforall a. Num a => a -> a -> a
*Double
a forall a. Num a => a -> a -> a
- Double
p1forall a. Fractional a => a -> a -> a
/Double
p3) forall a. Num a => a -> a -> a
+ Double
p0forall a. Fractional a => a -> a -> a
/Double
p3
th :: Double
th = forall a. Floating a => a -> a
acos(-Double
qforall a. Fractional a => a -> a -> a
/(Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt(-Double
pforall a. Num a => a -> a -> a
*Double
pforall a. Num a => a -> a -> a
*Double
p))) forall a. Fractional a => a -> a -> a
/ Double
3
r :: Double
r = Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
sqrt(-Double
p)
tp :: Double
tp = Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
3
a1 :: Double
a1 = -Double
a forall a. Num a => a -> a -> a
+ Double
rforall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
th
a3 :: Double
a3 = -Double
a forall a. Num a => a -> a -> a
+ Double
rforall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(Double
thforall a. Num a => a -> a -> a
+Double
tp)
a2 :: Double
a2 = -Double
3forall a. Num a => a -> a -> a
*Double
a forall a. Num a => a -> a -> a
- Double
a1 forall a. Num a => a -> a -> a
- Double
a3
npdbl :: Double
npdbl = Int -> Double
toDbl Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
5) ((Double
1forall a. Num a => a -> a -> a
-Double
ndblforall a. Num a => a -> a -> a
*Double
a1)forall a. Fractional a => a -> a -> a
/Double
npdbl)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
5) ((Double
1forall a. Num a => a -> a -> a
+Double
a1)forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
4) Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
6) ((Double
1forall a. Num a => a -> a -> a
-Double
ndblforall a. Num a => a -> a -> a
*Double
a2)forall a. Fractional a => a -> a -> a
/Double
npdbl)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
6) ((Double
1forall a. Num a => a -> a -> a
+Double
a2)forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
5) Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
7) ((Double
1forall a. Num a => a -> a -> a
-Double
ndblforall a. Num a => a -> a -> a
*Double
a3)forall a. Fractional a => a -> a -> a
/Double
npdbl)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
7) ((Double
1forall a. Num a => a -> a -> a
+Double
a3)forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
6) Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
5, Int
iwforall a. Num a => a -> a -> a
-Int
5)
((Double
u7forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
u6forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
*Double
u5)forall a. Fractional a => a -> a -> a
/(Double
a1forall a. Num a => a -> a -> a
*Double
a1forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
a1forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a3) forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a1 Int
5)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
6, Int
iwforall a. Num a => a -> a -> a
-Int
5)
((Double
u7forall a. Num a => a -> a -> a
-(Double
a1forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
u6forall a. Num a => a -> a -> a
+Double
a1forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
*Double
u5)forall a. Fractional a => a -> a -> a
/(Double
a2forall a. Num a => a -> a -> a
*Double
a2forall a. Num a => a -> a -> a
-(Double
a1forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
a2forall a. Num a => a -> a -> a
+Double
a1forall a. Num a => a -> a -> a
*Double
a3) forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a2 Int
5)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
7, Int
iwforall a. Num a => a -> a -> a
-Int
5)
((Double
u7forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a1)forall a. Num a => a -> a -> a
*Double
u6forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a1forall a. Num a => a -> a -> a
*Double
u5)forall a. Fractional a => a -> a -> a
/(Double
a3forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a1)forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a1) forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a3 Int
5)
let invnp7' :: Double
invnp7' = Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
8) (Double
4forall a. Num a => a -> a -> a
*Double
invnp7')
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
8) (Double
4forall a. Num a => a -> a -> a
*Double
invnp7')
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
8) Double
invnp7') [Int
3..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
7) Int
o
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
8, Int
iwforall a. Num a => a -> a -> a
-Int
5) (Double
10 forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
729forall a. Num a => a -> a -> a
*Int
n6))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
9) (Double
5.5forall a. Num a => a -> a -> a
*Double
invnp7')
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
9) (Double
2.5forall a. Num a => a -> a -> a
*Double
invnp7')
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
9) Double
invnp7') [Int
3..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
8) (Int
npforall a. Num a => a -> a -> a
*Int
n)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
9, Int
iwforall a. Num a => a -> a -> a
-Int
5) (Double
64 forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
6561forall a. Num a => a -> a -> a
*Int
n6))
let tmp''' :: Double
tmp''' = Int -> Double
toDbl (Int
64forall a. Num a => a -> a -> a
*Int
n6)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
5) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
7 forall a. Fractional a => a -> a -> a
/ Double
tmp''')
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
7, Int
iwforall a. Num a => a -> a -> a
-Int
5) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Double
tmp''')
let invnp9 :: Double
invnp9 = Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
9)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
8) (Double
9forall a. Num a => a -> a -> a
*Double
invnp9)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
8) Double
invnp9) [Int
2..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
7 Int
np
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
9) (Double
7forall a. Num a => a -> a -> a
*Double
invnp9)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
9) (Double
3forall a. Num a => a -> a -> a
*Double
invnp9)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
9) Double
invnp9) [Int
3..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
8 (Int
npforall a. Num a => a -> a -> a
*Int
n)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
10) (Double
5forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
1, Int
2]
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
10) Double
invnp9) [Int
3..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
9 Int
o
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
11) (Double
5forall a. Num a => a -> a -> a
*Double
invnp9)
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
11) (Double
3forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
2, Int
3]
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
np forall a. Ord a => a -> a -> Bool
> Int
3) forall a b. (a -> b) -> a -> b
$
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
11) Double
invnp9) [Int
4..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
10 (Int
oforall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
-Int
1))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
6) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
9 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
1536forall a. Num a => a -> a -> a
*Int
n6))
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
6)
(Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
9 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
512forall a. Num a => a -> a -> a
*Int
n6forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
+Int
7)))) [Int
3, Int
4]
let tmp'''' :: Double
tmp'''' = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ Int
256forall a. Num a => a -> a -> a
*Int
n8
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
6) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
9 forall a. Fractional a => a -> a -> a
/ Double
tmp'''')) [Int
5, Int
6, Int
7]
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
6) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
9) Int
8 forall a. Fractional a => a -> a -> a
/ Double
tmp'''')) [Int
8..Int
11]
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
n forall a. Ord a => a -> a -> Bool
> Int
2) forall a b. (a -> b) -> a -> b
$ do
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
12) (Double
3forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
1..Int
4]
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
12) Double
invnp9) [Int
5..Int
np]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
11 (forall a. Integral a => a -> a -> a
div (Int
npforall a. Num a => a -> a -> a
*Int
nforall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
-Int
1)forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
-Int
2)) Int
24)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
12, Int
iwforall a. Num a => a -> a -> a
-Int
6) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
9) Int
8 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
256forall a. Num a => a -> a -> a
*Int
n8))
Seq IO1dArray
rowsIO <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IOMatrix -> Int -> IO IO1dArray
extractRow IOMatrix
w) (forall a. [a] -> Seq a
S.fromList [Int
2..Int
wts])
Seq VectorD
rows <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM IO1dArray -> IO VectorD
array1dToVectorD Seq IO1dArray
rowsIO
let cols :: Seq VectorD
cols = Seq VectorD -> Seq VectorD
transpose Seq VectorD
rows
Vector Int
pts_out <- forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
UV.unsafeFreeze IOVectorI
pts
let ptsU :: VectorD
ptsU = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Int -> Double
toDbl Vector Int
pts_out
row1 :: Seq Double
row1 = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\VectorD
col -> Double
1 forall a. Num a => a -> a -> a
- forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0
(forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) (forall a. Unbox a => Vector a -> Vector a
UV.tail VectorD
ptsU) VectorD
col))
(forall a. Int -> Seq a -> Seq a
S.take Int
rls Seq VectorD
cols)
wcols :: Seq VectorD
wcols = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Int
j -> forall a. Unbox a => a -> Vector a -> Vector a
UV.cons (forall a. Seq a -> Int -> a
index Seq Double
row1 Int
j) (forall a. Seq a -> Int -> a
index Seq VectorD
cols Int
j))
(forall a. [a] -> Seq a
S.fromList [Int
0..(Int
rlsforall a. Num a => a -> a -> a
-Int
1)])
col1 :: VectorD
col1 = forall a. Seq a -> Int -> a
index Seq VectorD
wcols Int
0
wcols2 :: Seq VectorD
wcols2 = forall a. a -> Seq a -> Seq a
(S.<|) VectorD
col1
(forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\VectorD
col -> forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
col VectorD
col1) (forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
wcols))
col2 :: VectorD
col2 = forall a. Seq a -> Int -> a
index Seq VectorD
wcols2 Int
1
nb :: Double
nb = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
col1))
ratio :: Double
ratio = Double
nb forall a. Fractional a => a -> a -> a
/ forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
col2))
wcol2 :: VectorD
wcol2 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
ratio) VectorD
col2
wcols3 :: Seq VectorD
wcols3 = forall a. a -> Seq a -> Seq a
(S.<|) VectorD
col1 (forall a. a -> Seq a -> Seq a
(S.<|) VectorD
wcol2 (forall a. Int -> Seq a -> Seq a
S.drop Int
2 Seq VectorD
wcols2))
let updateW :: Seq VectorD -> Int -> Seq VectorD
updateW :: Seq VectorD -> Int -> Seq VectorD
updateW Seq VectorD
cols' Int
k = forall a. Int -> a -> Seq a -> Seq a
update (Int
kforall a. Num a => a -> a -> a
-Int
1) VectorD
wknew Seq VectorD
cols'
where
ptsW :: VectorD
ptsW = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Fractional a => a -> a -> a
/Double
nb) (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a. Seq a -> Int -> a
index Seq VectorD
cols' (Int
kforall a. Num a => a -> a -> a
-Int
1)))
slice :: Seq VectorD
slice = forall a. Int -> Seq a -> Seq a
S.drop Int
1 (forall a. Int -> Seq a -> Seq a
S.take (Int
kforall a. Num a => a -> a -> a
-Int
1) Seq VectorD
cols')
prod1 :: VectorD
prod1 = (forall a. Unbox a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) forall a b. (a -> b) -> a -> b
$
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsW) Seq VectorD
slice
rows' :: Seq VectorD
rows' = Seq VectorD -> Seq VectorD
transpose Seq VectorD
slice
prod2 :: VectorD
prod2 = (forall a. Unbox a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) forall a b. (a -> b) -> a -> b
$
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
prod1) Seq VectorD
rows'
wk :: VectorD
wk = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) (forall a. Seq a -> Int -> a
index Seq VectorD
cols' (Int
kforall a. Num a => a -> a -> a
-Int
1)) VectorD
prod2
ratio' :: Double
ratio' = Double
nb forall a. Fractional a => a -> a -> a
/
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
wk))
wknew :: VectorD
wknew = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
ratio') VectorD
wk
wcolsnew :: Seq VectorD
wcolsnew = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Seq VectorD -> Int -> Seq VectorD
updateW Seq VectorD
wcols3 [Int
3..Int
rls]
forall (m :: * -> *) a. Monad m => a -> m a
return (IOMatrix
g, Seq VectorD -> Seq VectorD
transpose Seq VectorD
wcolsnew, forall a. Unbox a => Vector a -> [a]
toList forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => (a -> Bool) -> Vector a -> Vector Int
UV.findIndices (forall a. Eq a => a -> a -> Bool
/= Int
0) Vector Int
pts_out)
transpose :: Seq VectorD -> Seq VectorD
transpose :: Seq VectorD -> Seq VectorD
transpose Seq VectorD
cols =
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Int
i -> (forall a. Unbox a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Unbox a => Vector a -> Int -> a
UV.!Int
i) Seq VectorD
cols)
(forall a. [a] -> Seq a
S.fromList [Int
0..(forall a. Unbox a => Vector a -> Int
UV.length (forall a. Seq a -> Int -> a
index Seq VectorD
cols Int
0) forall a. Num a => a -> a -> a
- Int
1)])
matprod :: IOMatrix -> VectorD -> IO VectorD
matprod :: IOMatrix -> VectorD -> IO VectorD
matprod IOMatrix
mat VectorD
x = do
((Int, Int)
_, (Int
m, Int
n)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
mat
IOVectorD
out <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
UMV.new Int
m :: IO IOVectorD
let step :: Int -> IO VectorD
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
mforall a. Num a => a -> a -> a
+Int
1 = forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze IOVectorD
out
| Bool
otherwise = do
!Double
coef <- Int -> Double -> IO Double
innerstep Int
1 Double
0
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
out (Int
iforall a. Num a => a -> a -> a
-Int
1) Double
coef
Int -> IO VectorD
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
where
innerstep :: Int -> Double -> IO Double
innerstep :: Int -> Double -> IO Double
innerstep Int
j !Double
s | Int
j forall a. Eq a => a -> a -> Bool
== Int
nforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return Double
s
| Bool
otherwise = do
Double
mat_ij <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
mat (Int
i, Int
j)
Int -> Double -> IO Double
innerstep (Int
jforall a. Num a => a -> a -> a
+Int
1) (Double
s forall a. Num a => a -> a -> a
+ Double
mat_ij forall a. Num a => a -> a -> a
* (VectorD
x forall a. Unbox a => Vector a -> Int -> a
UV.! (Int
jforall a. Num a => a -> a -> a
-Int
1)))
Int -> IO VectorD
step Int
1
smpsms :: IOMatrix -> Int -> (VectorD -> VectorD) -> IO1dArray
-> Double -> IO IO1dArray
smpsms :: IOMatrix
-> Int
-> (VectorD -> VectorD)
-> IO1dArray
-> Double
-> IO IO1dArray
smpsms IOMatrix
vertex Int
nf VectorD -> VectorD
f IO1dArray
g Double
scalar = do
[Double]
gAsList <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems IO1dArray
g
[VectorD]
f_gPermuts <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap VectorD -> VectorD
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. IOMatrix -> VectorD -> IO VectorD
matprod IOMatrix
vertex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Unbox a => [a] -> Vector a
fromList)
(forall a. (Eq a, Ord a) => [a] -> [[a]]
permuteMultiset [Double]
gAsList)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> [e] -> m (a i e)
newListArray (Int
1, Int
nf)
(forall a. Unbox a => Vector a -> [a]
toList (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*Double
scalar) (forall a. (a -> a -> a) -> [a] -> a
foldl1' (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) [VectorD]
f_gPermuts)))
extractColumn :: IOMatrix -> Int -> IO IO1dArray
extractColumn :: IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
m Int
j = do
((Int, Int)
_, (Int
nrow, Int
_)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
nrow) (\Int
i -> (Int
i, Int
j)) IOMatrix
m
extractRow :: IOMatrix -> Int -> IO IO1dArray
IOMatrix
m Int
i = do
((Int, Int)
_, (Int
_, Int
ncol)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
ncol) (\Int
j -> (Int
i, Int
j)) IOMatrix
m
outerProduct :: IO1dArray -> VectorD -> IO IOMatrix
outerProduct :: IO1dArray -> VectorD -> IO IOMatrix
outerProduct IO1dArray
x1 VectorD
x2 = do
(Int
_, Int
n1) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO1dArray
x1
let n2 :: Int
n2 = forall a. Unbox a => Vector a -> Int
UV.length VectorD
x2
IOMatrix
out <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ ((Int
1, Int
1), (Int
n1, Int
n2)) :: IO IOMatrix
let step :: Int -> IO IOMatrix
step :: Int -> IO IOMatrix
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
n1forall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return IOMatrix
out
| Bool
otherwise = do
Double
x1_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
x1 Int
i
Double -> Int -> IO IOMatrix
inner Double
x1_i Int
0
where
inner :: Double -> Int -> IO IOMatrix
inner !Double
x Int
j | Int
j forall a. Eq a => a -> a -> Bool
== Int
n2 = Int -> IO IOMatrix
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
| Bool
otherwise = do
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
out (Int
i, Int
jforall a. Num a => a -> a -> a
+Int
1) (Double
x forall a. Num a => a -> a -> a
* (VectorD
x2 forall a. Unbox a => Vector a -> Int -> a
UV.! Int
j))
Double -> Int -> IO IOMatrix
inner Double
x (Int
jforall a. Num a => a -> a -> a
+Int
1)
Int -> IO IOMatrix
step Int
1
sumMatrices :: [IOMatrix] -> IO IOMatrix
sumMatrices :: [IOMatrix] -> IO IOMatrix
sumMatrices [IOMatrix]
matrices = do
((Int, Int)
_, (Int
n1, Int
n2)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds (forall a. [a] -> a
head [IOMatrix]
matrices)
IOMatrix
out <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ ((Int
1, Int
1), (Int
n1, Int
n2)) :: IO IOMatrix
let step :: Int -> IO IOMatrix
step :: Int -> IO IOMatrix
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
n1forall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return IOMatrix
out
| Bool
otherwise = Int -> IO IOMatrix
inner Int
1
where
inner :: Int -> IO IOMatrix
inner :: Int -> IO IOMatrix
inner Int
j | Int
j forall a. Eq a => a -> a -> Bool
== Int
n2forall a. Num a => a -> a -> a
+Int
1 = Int -> IO IOMatrix
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
| Bool
otherwise = do
[Double]
coefs <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\IOMatrix
m -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
m (Int
i, Int
j)) [IOMatrix]
matrices
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
out (Int
i, Int
j) (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
coefs)
Int -> IO IOMatrix
inner (Int
jforall a. Num a => a -> a -> a
+Int
1)
Int -> IO IOMatrix
step Int
1
smprul :: IOMatrix -> Int -> (VectorD -> VectorD) -> Double -> IOMatrix
-> Seq VectorD -> [Int] -> IO (IO1dArray, IO1dArray)
smprul :: IOMatrix
-> Int
-> (VectorD -> VectorD)
-> Double
-> IOMatrix
-> Seq VectorD
-> [Int]
-> IO (IO1dArray, IO1dArray)
smprul IOMatrix
vrts Int
nf VectorD -> VectorD
f Double
vol IOMatrix
g Seq VectorD
w [Int]
pospts = do
let rtmn :: Double
rtmn = Double
0.1
small :: Double
small = Double
1e-12
errcof :: Double
errcof = Double
8
rls :: Int
rls = forall a. Unbox a => Vector a -> Int
UV.length (forall a. Seq a -> Int -> a
index Seq VectorD
w Int
0)
[IOMatrix]
toSum <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\Int
k -> do
IO1dArray
g_colk <- IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
g (Int
kforall a. Num a => a -> a -> a
+Int
1)
IO1dArray
sms <- IOMatrix
-> Int
-> (VectorD -> VectorD)
-> IO1dArray
-> Double
-> IO IO1dArray
smpsms IOMatrix
vrts Int
nf VectorD -> VectorD
f IO1dArray
g_colk Double
vol
IO1dArray -> VectorD -> IO IOMatrix
outerProduct IO1dArray
sms (forall a. Seq a -> Int -> a
index Seq VectorD
w Int
k))
[Int]
pospts
IOMatrix
rule <- [IOMatrix] -> IO IOMatrix
sumMatrices [IOMatrix]
toSum
IO1dArray
basval <- IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
rule Int
1
IO1dArray
rgnerr <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray (Int
1, Int
nf) Double
0 :: IO IO1dArray
let step :: Int -> IO ()
step :: Int -> IO ()
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
nfforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
Double
basval_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
basval Int
i
let nmbs :: Double
nmbs = forall a. Num a => a -> a
abs Double
basval_i
(Double
rt, Double
nmcp) <- Int -> Double -> Double -> Double -> IO (Double, Double)
inner Int
rls Double
rtmn Double
0 Double
nmbs
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
rt forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Int
rls forall a. Ord a => a -> a -> Bool
> Int
3) forall a b. (a -> b) -> a -> b
$
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (Double
rtforall a. Num a => a -> a -> a
*Double
nmcp)
Double
rgnerr_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
rgnerr Int
i
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (forall a. Ord a => a -> a -> a
max (Double
errcofforall a. Num a => a -> a -> a
*Double
rgnerr_i) (Double
smallforall a. Num a => a -> a -> a
*Double
nmbs))
Int -> IO ()
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
where
inner :: Int -> Double -> Double -> Double -> IO (Double, Double)
inner :: Int -> Double -> Double -> Double -> IO (Double, Double)
inner Int
k !Double
x !Double
y !Double
z | Int
k forall a. Eq a => a -> a -> Bool
== Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return (Double
x, Double
y)
| Bool
otherwise = do
Double
rule_ik <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
rule (Int
i, Int
k)
Double
rule_ikm1 <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
rule (Int
i, Int
kforall a. Num a => a -> a -> a
-Int
1)
let nmrl :: Double
nmrl = forall a. Ord a => a -> a -> a
max (forall a. Num a => a -> a
abs Double
rule_ik) (forall a. Num a => a -> a
abs Double
rule_ikm1)
Double
rgnerr_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
rgnerr Int
i
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (forall a. Ord a => a -> a -> a
max Double
nmrl Double
rgnerr_i)
if Double
nmrl forall a. Ord a => a -> a -> Bool
> Double
smallforall a. Num a => a -> a -> a
*Double
z Bool -> Bool -> Bool
&& Int
k forall a. Ord a => a -> a -> Bool
< Int
rls
then Int -> Double -> Double -> Double -> IO (Double, Double)
inner (Int
kforall a. Num a => a -> a -> a
-Int
2) (forall a. Ord a => a -> a -> a
max (Double
nmrlforall a. Fractional a => a -> a -> a
/Double
y) Double
x) Double
nmrl Double
z
else Int -> Double -> Double -> Double -> IO (Double, Double)
inner (Int
kforall a. Num a => a -> a -> a
-Int
2) Double
x Double
nmrl Double
z
Int -> IO ()
step Int
1
forall (m :: * -> *) a. Monad m => a -> m a
return (IO1dArray
basval, IO1dArray
rgnerr)
rowMeans :: IOMatrix -> IO VectorD
rowMeans :: IOMatrix -> IO VectorD
rowMeans IOMatrix
m = do
((Int, Int)
_, (Int
nrow, Int
ncol)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
IOVectorD
outIO <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
new Int
nrow :: IO IOVectorD
let step :: Int -> IO ()
step :: Int -> IO ()
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
nrowforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
!Double
sum_i <- Int -> Double -> IO Double
inner Int
1 Double
0
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
outIO (Int
iforall a. Num a => a -> a -> a
-Int
1) Double
sum_i
Int -> IO ()
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
where
inner :: Int -> Double -> IO Double
inner :: Int -> Double -> IO Double
inner Int
j !Double
x | Int
j forall a. Eq a => a -> a -> Bool
== Int
ncolforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return (Double
x forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl Int
ncol)
| Bool
otherwise = do
Double
coef <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
m (Int
i, Int
j)
Int -> Double -> IO Double
inner (Int
jforall a. Num a => a -> a -> a
+Int
1) (Double
x forall a. Num a => a -> a -> a
+ Double
coef)
Int -> IO ()
step Int
1
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze IOVectorD
outIO
array1dToVectorD :: IO1dArray -> IO VectorD
array1dToVectorD :: IO1dArray -> IO VectorD
array1dToVectorD IO1dArray
array = forall (m :: * -> *) a b. Monad m => (a -> b) -> m a -> m b
(<$!>) forall a. Unbox a => [a] -> Vector a
fromList (forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems IO1dArray
array)
getVectors :: Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors :: Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
n IO3dArray
m Int
k Int
j1 Int
j2 = do
VectorD
out1U <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
n) (\Int
i -> (Int
i, Int
j1, Int
k)) IO3dArray
m)
VectorD
out2U <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
n) (\Int
i -> (Int
i, Int
j2, Int
k)) IO3dArray
m)
forall (m :: * -> *) a. Monad m => a -> m a
return (VectorD
out1U, VectorD
out2U)
smpdfs :: Int -> (VectorD -> VectorD) -> Int -> Int
-> IO3dArray -> IO (Int, IO3dArray)
smpdfs :: Int
-> (VectorD -> VectorD)
-> Int
-> Int
-> IO3dArray
-> IO (Int, IO3dArray)
smpdfs Int
nd VectorD -> VectorD
f Int
top Int
sbs IO3dArray
vrts = do
let cuttf :: Double
cuttf = Double
2.0
cuttb :: Double
cuttb = Double
8.0
IOMatrix
v <- forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1)) (\(Int
i, Int
j) -> (Int
i, Int
j, Int
top)) IO3dArray
vrts
VectorD
cn <- IOMatrix -> IO VectorD
rowMeans IOMatrix
v
let fc :: VectorD
fc = VectorD -> VectorD
f VectorD
cn
dfmd :: Double
dfmd = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map forall a. Num a => a -> a
abs VectorD
fc)
IOMatrix
frthdf <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1)) Double
0 :: IO IOMatrix
IOVectorI
iejeitjtisjsls <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
new Int
7 :: IO IOVectorI
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
1
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
2
IOVectorD
dfmxdfnx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
UMV.replicate Int
2 Double
0 :: IO IOVectorD
let step :: Int -> Double -> IO ()
step :: Int -> Double -> IO ()
step Int
i Double
x | Int
i forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
!Double
emx <- Int -> Double -> IO Double
inner (Int
iforall a. Num a => a -> a -> a
+Int
1) Double
x
Int -> Double -> IO ()
step (Int
iforall a. Num a => a -> a -> a
+Int
1) Double
emx
where
inner :: Int -> Double -> IO Double
inner :: Int -> Double -> IO Double
inner Int
j !Double
y | Int
j forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
2 = forall (m :: * -> *) a. Monad m => a -> m a
return Double
y
| Bool
otherwise = do
VectorD
vi <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
i)
VectorD
vj <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
j)
let h :: VectorD
h = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*(Double
2forall a. Fractional a => a -> a -> a
/(Double
5forall a. Num a => a -> a -> a
*(Int -> Double
toDbl Int
nd forall a. Num a => a -> a -> a
+Double
1))))
(forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
vi VectorD
vj)
ewd :: Double
ewd = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map forall a. Num a => a -> a
abs VectorD
h)
twoh :: VectorD
twoh = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*Double
2) VectorD
h
t1 :: VectorD
t1 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
cn VectorD
twoh)
t2 :: VectorD
t2 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
cn VectorD
twoh)
t3 :: VectorD
t3 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*Double
6) VectorD
fc
t4 :: VectorD
t4 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
cn VectorD
h)
t5 :: VectorD
t5 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
cn VectorD
h)
t6 :: VectorD
t6 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*(-Double
4))forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
t4 VectorD
t5
tsum :: VectorD
tsum = forall a. (a -> a -> a) -> [a] -> a
foldl1' (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) [VectorD
t1, VectorD
t2, VectorD
t3, VectorD
t6]
dfr1 :: Double
dfr1 = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map forall a. Num a => a -> a
abs VectorD
tsum)
dfr2 :: Double
dfr2 = if Double
dfmdforall a. Num a => a -> a -> a
+Double
dfr1forall a. Fractional a => a -> a -> a
/Double
8 forall a. Eq a => a -> a -> Bool
== Double
dfmd then Double
0 else Double
dfr1
dfr3 :: Double
dfr3 = Double
dfr2forall a. Num a => a -> a -> a
*Double
ewd
Double
dfmx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
0
if Double
dfr3 forall a. Ord a => a -> a -> Bool
>= Double
dfmx
then do
Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
is
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
3 Int
js
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
i
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
j
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
1 Double
dfmx
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
0 Double
dfr3
else do
Double
dfnx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
1
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
dfr3 forall a. Ord a => a -> a -> Bool
>= Double
dfnx) forall a b. (a -> b) -> a -> b
$ do
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
i
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
3 Int
j
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
1 Double
dfr3
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
frthdf (Int
i, Int
j) Double
dfr3
if Double
ewd forall a. Ord a => a -> a -> Bool
>= Double
y
then do
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
0 Int
i
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
1 Int
j
Int -> Double -> IO Double
inner (Int
jforall a. Num a => a -> a -> a
+Int
1) Double
ewd
else Int -> Double -> IO Double
inner (Int
jforall a. Num a => a -> a -> a
+Int
1) Double
y
Int -> Double -> IO ()
step Int
1 Double
0
Double
dfmx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
0
Double
dfnx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
1
let nregions :: Int
nregions = if Double
dfnx forall a. Ord a => a -> a -> Bool
> Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttf then Int
4 else Int
3
if Double
dfnx forall a. Ord a => a -> a -> Bool
> Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttf
then forall (m :: * -> *) a. Monad m => a -> m a
return ()
else
if Double
dfmx forall a. Eq a => a -> a -> Bool
== Double
0
then do
Int
ie <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
0
Int
je <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
1
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
ie
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
je
else do
let loop :: Int -> Double -> Int -> IO Int
loop :: Int -> Double -> Int -> IO Int
loop Int
l !Double
x !Int
ls | Int
l forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
2 = forall (m :: * -> *) a. Monad m => a -> m a
return Int
ls
| Bool
otherwise = do
Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
if (Int
l forall a. Eq a => a -> a -> Bool
/= Int
is) Bool -> Bool -> Bool
&& (Int
l forall a. Eq a => a -> a -> Bool
/= Int
js)
then do
let it :: Int
it = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Int
l, Int
is, Int
js]
jt :: Int
jt = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Int
l, Int
is, Int
js]
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
it
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
3 Int
jt
let lt :: Int
lt = Int
isforall a. Num a => a -> a -> a
+Int
jsforall a. Num a => a -> a -> a
+Int
lforall a. Num a => a -> a -> a
-Int
itforall a. Num a => a -> a -> a
-Int
jt
Double
dfr1 <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (Int
it, Int
lt)
Double
dfr2 <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (Int
lt, Int
jt)
let dfr :: Double
dfr = Double
dfr1 forall a. Num a => a -> a -> a
+ Double
dfr2
if Double
dfr forall a. Ord a => a -> a -> Bool
>= Double
x
then Int -> Double -> Int -> IO Int
loop (Int
lforall a. Num a => a -> a -> a
+Int
1) Double
dfr Int
l
else Int -> Double -> Int -> IO Int
loop (Int
lforall a. Num a => a -> a -> a
+Int
1) Double
x Int
ls
else Int -> Double -> Int -> IO Int
loop (Int
lforall a. Num a => a -> a -> a
+Int
1) Double
x Int
ls
!Int
ls <- Int -> Double -> Int -> IO Int
loop Int
1 Double
0 Int
0
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
6 Int
ls
Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
Double
difil <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (forall a. Ord a => a -> a -> a
min Int
is Int
ls, forall a. Ord a => a -> a -> a
max Int
is Int
ls)
Double
diflj <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (forall a. Ord a => a -> a -> a
min Int
js Int
ls, forall a. Ord a => a -> a -> a
max Int
js Int
ls)
let dfnx' :: Double
dfnx' = forall a. Ord a => a -> a -> a
max Double
difil Double
diflj
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
1 Double
dfnx'
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttb forall a. Ord a => a -> a -> Bool
< Double
dfnx' Bool -> Bool -> Bool
&& Double
difil forall a. Ord a => a -> a -> Bool
> Double
diflj) forall a b. (a -> b) -> a -> b
$ do
Int
is' <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
Int
js' <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
Int
it <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
2
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
is'
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
js'
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
it
IO3dArray
vrts2 <- forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1, Int
sbsforall a. Num a => a -> a -> a
+Int
nregionsforall a. Num a => a -> a -> a
-Int
1))
(\(Int
i, Int
j, Int
k) -> (Int
i, Int
j, if Int
k forall a. Ord a => a -> a -> Bool
<= Int
sbs then Int
k else Int
top)) IO3dArray
vrts
Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
VectorD
vti <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
is)
VectorD
vtj <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
js)
if Int
nregions forall a. Eq a => a -> a -> Bool
== Int
4
then do
let vt :: VectorD
vt = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vti VectorD
vtj
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
top) (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vt
IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vti
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
2) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vt
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vtj
Int
it <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
2
Int
jt <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
3
(VectorD
vtit, VectorD
vtjt) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 Int
top Int
it Int
jt
let vtt :: VectorD
vtt = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vtit VectorD
vtjt
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
jt, Int
top) (Int
it, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vtt
IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
jt, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vtjt
(VectorD
vti2, VectorD
vtj2) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 (Int
sbsforall a. Num a => a -> a -> a
+Int
2) Int
it Int
jt
let vt2 :: VectorD
vt2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vti2 VectorD
vtj2
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
jt, Int
sbsforall a. Num a => a -> a -> a
+Int
2) (Int
it, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vt2
IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
jt, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vtj2
else do
let vt :: VectorD
vt = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (\Double
x Double
y -> (Double
2forall a. Num a => a -> a -> a
*Double
xforall a. Num a => a -> a -> a
+Double
y)forall a. Fractional a => a -> a -> a
/Double
3) VectorD
vti VectorD
vtj
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
top) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vt
if Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttf forall a. Ord a => a -> a -> Bool
< Double
dfnx
then do
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
1) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
Int
ls <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
6
(VectorD
vtj1, VectorD
vtl1) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 (Int
sbsforall a. Num a => a -> a -> a
+Int
1) Int
js Int
ls
let vt1 :: VectorD
vt1 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vtj1 VectorD
vtl1
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
ls, Int
sbsforall a. Num a => a -> a -> a
+Int
1) (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vt1
IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
ls, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtl1
else do
let vv :: VectorD
vv = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (\Double
x Double
y -> (Double
xforall a. Num a => a -> a -> a
+Double
2forall a. Num a => a -> a -> a
*Double
y)forall a. Fractional a => a -> a -> a
/Double
3) VectorD
vti VectorD
vtj
IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
1) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vv
IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
nregions, IO3dArray
vrts2)
where
replaceDimension :: IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension :: IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
m (Int
j, Int
k) VectorD
v = do
let loop :: Int -> IO ()
loop :: Int -> IO ()
loop Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO3dArray
m (Int
i, Int
j, Int
k) (forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iforall a. Num a => a -> a -> a
-Int
1))
Int -> IO ()
loop (Int
iforall a. Num a => a -> a -> a
+Int
1)
Int -> IO ()
loop Int
1
replaceDimensions :: IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions :: IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
m (Int
j1, Int
k1) (Int
j2, Int
k2) VectorD
v = do
let loop :: Int -> IO ()
loop :: Int -> IO ()
loop Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO3dArray
m (Int
i, Int
j1, Int
k1) (forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iforall a. Num a => a -> a -> a
-Int
1))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO3dArray
m (Int
i, Int
j2, Int
k2) (forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iforall a. Num a => a -> a -> a
-Int
1))
Int -> IO ()
loop (Int
iforall a. Num a => a -> a -> a
+Int
1)
Int -> IO ()
loop Int
1
smpchc :: Int -> Int -> Int
smpchc :: Int -> Int -> Int
smpchc Int
nd Int
key
| Int
key forall a. Eq a => a -> a -> Bool
== Int
3
= forall a. Integral a => a -> a -> a
div ((Int
ndforall a. Num a => a -> a -> a
+Int
4)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
3)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
2)) Int
6 forall a. Num a => a -> a -> a
+ (Int
ndforall a. Num a => a -> a -> a
+Int
2)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
1)
| Int
key forall a. Eq a => a -> a -> Bool
== Int
1
= Int
2forall a. Num a => a -> a -> a
*Int
nd forall a. Num a => a -> a -> a
+ Int
3
| Int
key forall a. Eq a => a -> a -> Bool
== Int
2
= forall a. Integral a => a -> a -> a
div ((Int
ndforall a. Num a => a -> a -> a
+Int
3)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
2)) Int
2 forall a. Num a => a -> a -> a
+ Int
2forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
1)
| Bool
otherwise
= forall a. Integral a => a -> a -> a
div ((Int
ndforall a. Num a => a -> a -> a
+Int
5)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
4)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
3)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
2)) Int
24 forall a. Num a => a -> a -> a
+ Int
5 forall a. Num a => a -> a -> a
* forall a. Integral a => a -> a -> a
div ((Int
ndforall a. Num a => a -> a -> a
+Int
2)forall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
1)) Int
2
check :: Int -> Int -> Int -> Double -> Double -> Int -> Int -> Int
check :: Int -> Int -> Int -> Double -> Double -> Int -> Int -> Int
check Int
nd Int
nf Int
mxfs Double
ea Double
er Int
sbs Int
key
| Double
ea forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
er forall a. Ord a => a -> a -> Bool
< Double
0
= Int
5
| Int
nf forall a. Ord a => a -> a -> Bool
< Int
1
= Int
4
| Int
nd forall a. Ord a => a -> a -> Bool
< Int
2
= Int
3
| Int
key forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Int
key forall a. Ord a => a -> a -> Bool
> Int
4
= Int
2
| Int
mxfs forall a. Ord a => a -> a -> Bool
< Int
sbs forall a. Num a => a -> a -> a
* Int -> Int -> Int
smpchc Int
nd Int
key
= Int
1
| Bool
otherwise
= Int
0
adsimp :: Int -> Int -> Int -> (VectorD -> VectorD) -> Double -> Double
-> Int -> IO3dArray -> IO (VectorD, VectorD, Int, Bool)
adsimp :: Int
-> Int
-> Int
-> (VectorD -> VectorD)
-> Double
-> Double
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
adsimp Int
nd Int
nf Int
mxfs VectorD -> VectorD
f Double
ea Double
er Int
key IO3dArray
vrts = do
((Int, Int, Int)
_, (Int
_, Int
_, Int
sbs)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO3dArray
vrts
case Int -> Int -> Int -> Double -> Double -> Int -> Int -> Int
check Int
nd Int
nf Int
mxfs Double
ea Double
er Int
sbs Int
key of
Int
0 -> Int
-> Int
-> (VectorD -> VectorD)
-> Int
-> Double
-> Double
-> Int
-> Int
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
smpsad Int
nd Int
nf VectorD -> VectorD
f Int
mxfs Double
ea Double
er Int
key (Int -> Int -> Int
smpchc Int
nd Int
key) Int
sbs IO3dArray
vrts
Int
1 -> Int
-> Int
-> (VectorD -> VectorD)
-> Int
-> Double
-> Double
-> Int
-> Int
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
smpsad Int
nd Int
nf VectorD -> VectorD
f (Int
sbs forall a. Num a => a -> a -> a
* Int -> Int -> Int
smpchc Int
nd Int
key) Double
ea Double
er Int
key (Int -> Int -> Int
smpchc Int
nd Int
key) Int
sbs IO3dArray
vrts
Int
2 -> forall a. HasCallStack => [Char] -> a
error [Char]
"integration rule must be between 1 and 4"
Int
3 -> forall a. HasCallStack => [Char] -> a
error [Char]
"dimension must be at least 2"
Int
4 -> forall a. HasCallStack => [Char] -> a
error [Char]
"number of components must be at least 1"
Int
5 -> forall a. HasCallStack => [Char] -> a
error [Char]
"requested errors must be positive"
Int
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"this should not happen"
type Params = (Bool, Int, Int, Seq VectorD, Seq VectorD,
VectorD, VectorD, IO3dArray, Seq Double)
smpsad :: Int -> Int -> (VectorD -> VectorD) -> Int -> Double -> Double -> Int
-> Int -> Int -> IO3dArray -> IO (VectorD, VectorD, Int, Bool)
smpsad :: Int
-> Int
-> (VectorD -> VectorD)
-> Int
-> Double
-> Double
-> Int
-> Int
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
smpsad Int
nd Int
nf VectorD -> VectorD
f Int
mxfs Double
ea Double
er Int
key Int
rcls Int
sbs IO3dArray
vrts = do
let dfcost :: Int
dfcost = Int
1 forall a. Num a => a -> a -> a
+ Int
2forall a. Num a => a -> a -> a
*Int
ndforall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
1)
(IOMatrix
g, Seq VectorD
w, [Int]
pospts) <- Int -> Int -> IO (IOMatrix, Seq VectorD, [Int])
smprms Int
nd Int
key
[Simplex]
simplices <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO3dArray -> Int -> Int -> IO Simplex
toSimplex IO3dArray
vrts (Int
ndforall a. Num a => a -> a -> a
+Int
1)) [Int
1..Int
sbs]
let vol :: Seq Double
vol = forall a. [a] -> Seq a
S.fromList forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map Simplex -> Double
simplexVolume [Simplex]
simplices
nv :: Int
nv = Int
sbsforall a. Num a => a -> a -> a
*Int
rcls
Seq IOMatrix
matrices <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM
(\Int
k -> forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1)) (\(Int
i, Int
j) -> (Int
i, Int
j, Int
k)) IO3dArray
vrts)
(forall a. [a] -> Seq a
S.fromList [Int
1..Int
sbs])
Seq (IO1dArray, IO1dArray)
br <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\(IOMatrix
m, Double
v) -> IOMatrix
-> Int
-> (VectorD -> VectorD)
-> Double
-> IOMatrix
-> Seq VectorD
-> [Int]
-> IO (IO1dArray, IO1dArray)
smprul IOMatrix
m Int
nf VectorD -> VectorD
f Double
v IOMatrix
g Seq VectorD
w [Int]
pospts) (forall a b. Seq a -> Seq b -> Seq (a, b)
S.zip Seq IOMatrix
matrices Seq Double
vol)
Seq VectorD
aes <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> b
snd) Seq (IO1dArray, IO1dArray)
br
Seq VectorD
vls <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> a
fst) Seq (IO1dArray, IO1dArray)
br
let vl :: VectorD
vl = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
vls
ae :: VectorD
ae = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
aes
fl :: Bool
fl = VectorD -> VectorD -> Bool
getFL VectorD
ae VectorD
vl
let loop :: Params -> IO (VectorD, VectorD, Int, Bool)
loop :: Params -> IO (VectorD, VectorD, Int, Bool)
loop !Params
params | Bool -> Bool
not Bool
fl' Bool -> Bool -> Bool
|| Int
nv'forall a. Num a => a -> a -> a
+Int
dfcostforall a. Num a => a -> a -> a
+Int
4forall a. Num a => a -> a -> a
*Int
rcls forall a. Ord a => a -> a -> Bool
> Int
mxfs =
forall (m :: * -> *) a. Monad m => a -> m a
return (VectorD
vl', VectorD
ae', Int
nv', Bool
fl')
| Bool
otherwise = do
let maxs :: Seq Double
maxs = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. (Unbox a, Ord a) => Vector a -> a
UV.maximum Seq VectorD
aes'
imax :: Int
imax = forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> Seq a -> Maybe Int
S.findIndexL (forall a. Eq a => a -> a -> Bool
== forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum Seq Double
maxs) Seq Double
maxs
vl0 :: VectorD
vl0 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
vl' (forall a. Seq a -> Int -> a
index Seq VectorD
vls' Int
imax)
ae0 :: VectorD
ae0 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
ae' (forall a. Seq a -> Int -> a
index Seq VectorD
aes' Int
imax)
(Int
nregions, IO3dArray
vrts2) <- Int
-> (VectorD -> VectorD)
-> Int
-> Int
-> IO3dArray
-> IO (Int, IO3dArray)
smpdfs Int
nd VectorD -> VectorD
f (Int
imaxforall a. Num a => a -> a -> a
+Int
1) Int
sbs' IO3dArray
vrts'
let vi :: Double
vi = forall a. Seq a -> Int -> a
index Seq Double
vol' Int
imax forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl Int
nregions
nv2 :: Int
nv2 = Int
nv' forall a. Num a => a -> a -> a
+ Int
nregionsforall a. Num a => a -> a -> a
*Int
rcls forall a. Num a => a -> a -> a
+ Int
dfcost
sbs2 :: Int
sbs2 = Int
sbs' forall a. Num a => a -> a -> a
+ Int
nregionsforall a. Num a => a -> a -> a
-Int
1
Seq IOMatrix
matrices2 <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM
(\Int
k -> forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1))
(\(Int
i, Int
j) -> (Int
i, Int
j, Int
k)) IO3dArray
vrts2)
(forall a. [a] -> Seq a
S.fromList ((Int
imaxforall a. Num a => a -> a -> a
+Int
1)forall a. a -> [a] -> [a]
:[(Int
sbs'forall a. Num a => a -> a -> a
+Int
1)..Int
sbs2]))
Seq (IO1dArray, IO1dArray)
br2 <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\IOMatrix
m -> IOMatrix
-> Int
-> (VectorD -> VectorD)
-> Double
-> IOMatrix
-> Seq VectorD
-> [Int]
-> IO (IO1dArray, IO1dArray)
smprul IOMatrix
m Int
nf VectorD -> VectorD
f Double
vi IOMatrix
g Seq VectorD
w [Int]
pospts) Seq IOMatrix
matrices2
Seq VectorD
rgnerrs <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> b
snd) Seq (IO1dArray, IO1dArray)
br2
Seq VectorD
basvals <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> a
fst) Seq (IO1dArray, IO1dArray)
br2
let vl2 :: VectorD
vl2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
vl0
(forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
basvals)
ae2 :: VectorD
ae2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
ae0
(forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
rgnerrs)
aes2 :: Seq VectorD
aes2 = forall a. Int -> a -> Seq a -> Seq a
update Int
imax (forall a. Seq a -> Int -> a
index Seq VectorD
rgnerrs Int
0) Seq VectorD
aes'
forall a. Seq a -> Seq a -> Seq a
>< forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
rgnerrs
vls2 :: Seq VectorD
vls2 = forall a. Int -> a -> Seq a -> Seq a
update Int
imax (forall a. Seq a -> Int -> a
index Seq VectorD
basvals Int
0) Seq VectorD
vls'
forall a. Seq a -> Seq a -> Seq a
>< forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
basvals
fl2 :: Bool
fl2 = VectorD -> VectorD -> Bool
getFL VectorD
ae2 VectorD
vl2
vol2 :: Seq Double
vol2 = forall a. Int -> a -> Seq a -> Seq a
update Int
imax Double
vi Seq Double
vol'
forall a. Seq a -> Seq a -> Seq a
>< forall a. Int -> a -> Seq a
S.replicate (Int
nregionsforall a. Num a => a -> a -> a
-Int
1) Double
vi
Params -> IO (VectorD, VectorD, Int, Bool)
loop (Bool
fl2, Int
nv2, Int
sbs2, Seq VectorD
aes2, Seq VectorD
vls2, VectorD
ae2, VectorD
vl2, IO3dArray
vrts2, Seq Double
vol2)
where
(Bool
fl', Int
nv', Int
sbs', Seq VectorD
aes', Seq VectorD
vls', VectorD
ae', VectorD
vl', IO3dArray
vrts', Seq Double
vol') = Params
params
Params -> IO (VectorD, VectorD, Int, Bool)
loop (Bool
fl, Int
nv, Int
sbs, Seq VectorD
aes, Seq VectorD
vls, VectorD
ae, VectorD
vl, IO3dArray
vrts, Seq Double
vol)
where
getFL :: VectorD -> VectorD -> Bool
getFL VectorD
a VectorD
v = forall a. Unbox a => (a -> Bool) -> Vector a -> Bool
UV.any (forall a. Ord a => a -> a -> Bool
> forall a. Ord a => a -> a -> a
max Double
ea (forall a. (Unbox a, Ord a) => Vector a -> a
UV.maximum (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map ((forall a. Num a => a -> a -> a
*Double
er)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a
abs) VectorD
v))) VectorD
a
toSimplex :: IO3dArray -> Int -> Int -> IO Simplex
toSimplex :: IO3dArray -> Int -> Int -> IO Simplex
toSimplex IO3dArray
m Int
n Int
k = do
((Int, Int, Int)
_, (Int
nrow, Int
_, Int
_)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO3dArray
m
let getColumn :: Int -> IO IO1dArray
getColumn :: Int -> IO IO1dArray
getColumn Int
col = forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
nrow) (\Int
i -> (Int
i, Int
col, Int
k)) IO3dArray
m
[IO1dArray]
columns <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM Int -> IO IO1dArray
getColumn [Int
1..Int
n]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems [IO1dArray]
columns