{-# 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 = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral

pow :: Double -> Int -> Double
pow :: Double -> Int -> Double
pow Double
x Int
n = [Double] -> Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
x)

permuteMultiset :: (Eq a, Ord a) => [a] -> [[a]] 
permuteMultiset :: [a] -> [[a]]
permuteMultiset [a]
list = ([a] -> Maybe [a]) -> [a] -> [[a]]
forall a. (a -> Maybe a) -> a -> [a]
unfold1 [a] -> Maybe [a]
forall a. Ord a => [a] -> Maybe [a]
next ([a] -> [a]
forall a. Ord a => [a] -> [a]
sort [a]
list) where
  unfold1 :: (a -> Maybe a) -> a -> [a]
  unfold1 :: (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 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> Maybe a) -> a -> [a]
forall a. (a -> Maybe a) -> a -> [a]
unfold1 a -> Maybe a
f a
y
  -- next :: [a] -> Maybe [a]
  next :: [a] -> Maybe [a]
next [a]
xs = case ([a], [a]) -> Maybe ([a], [a])
forall a. Ord a => ([a], [a]) -> Maybe ([a], [a])
findj ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
xs, []) of 
    Maybe ([a], [a])
Nothing -> Maybe [a]
forall a. Maybe a
Nothing
    Just ( a
l:[a]
ls, [a]
rs) -> [a] -> Maybe [a]
forall a. a -> Maybe a
Just ([a] -> Maybe [a]) -> [a] -> Maybe [a]
forall a b. (a -> b) -> a -> b
$ a -> [a] -> ([a], [a]) -> [a]
forall a. Ord a => a -> [a] -> ([a], [a]) -> [a]
inc a
l [a]
ls ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
rs, []) 
    Just ( [], [a]
_ ) -> [Char] -> Maybe [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"permute: should not happen"
  -- we use simple list zippers: (left,right)
  -- findj :: ([a],[a]) -> Maybe ([a],[a])   
  findj :: ([a], [a]) -> Maybe ([a], [a])
findj ( xxs :: [a]
xxs@(a
x:[a]
xs), yys :: [a]
yys@(a
y:[a]
_) ) = if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y 
    then ([a], [a]) -> Maybe ([a], [a])
findj ( [a]
xs, a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
yys )
    else ([a], [a]) -> Maybe ([a], [a])
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]
_ ) = Maybe ([a], [a])
forall a. Maybe a
Nothing
  -- inc :: a -> [a] -> ([a],[a]) -> [a]
  inc :: a -> [a] -> ([a], [a]) -> [a]
inc !a
u [a]
us ( a
x:[a]
xs, [a]
yys ) = if a
u a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
x
    then a -> [a] -> ([a], [a]) -> [a]
inc a
u [a]
us ( [a]
xs, a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
yys ) 
    else [a] -> [a]
forall a. [a] -> [a]
reverse (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
us)  [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a] -> [a]
forall a. [a] -> [a]
reverse (a
ua -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
yys) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
xs
  inc a
_ [a]
_ ( [], [a]
_ ) = [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"permute: should not happen"

square :: Double -> Double
square :: Double -> Double
square Double
x = Double
xDouble -> Double -> Double
forall 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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Int
3, Int
2, Int
3) :: (Int, Int, Int)
                      | Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 = (Int
5, Int
4, Int
6) :: (Int, Int, Int)
                      | Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 = (Int
7, Int
7, Int
11) :: (Int, Int, Int)
                      | Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
4 = (if Int
n Int -> Int -> Bool
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 = [Char] -> (Int, Int, Int)
forall a. HasCallStack => [Char] -> a
error [Char]
"this should not happen"
  IOMatrix
w <- ((Int, Int), (Int, Int)) -> Double -> IO IOMatrix
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 <- Int -> Int -> IO (MVector (PrimState IO) Int)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
UMV.replicate Int
wts Int
0 :: IO IOVectorI
  IOMatrix
g <- ((Int, Int), (Int, Int)) -> Double -> IO IOMatrix
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
wts)) Double
0 :: IO IOMatrix
  let np :: Int
np = Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
      n2 :: Int
n2 = Int
np Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)
      n4 :: Int
n4 = Int
n2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4)
      n6 :: Int
n6 = Int
n4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
6)
      n8 :: Int
n8 = Int
n6 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
8)
      o :: Int
o = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div (Int
nInt -> Int -> Int
forall 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 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sqrt15) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndbl Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
8Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndbl Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
      s1 :: Double
s1 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
r1
      l1 :: Double
l1 = Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
r1
  (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1))) [Int
1 .. Int
np]
  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
0 Int
1
  IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
s1
  (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
r1) [Int
2 .. Int
np]
  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
gms Int
np
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
4) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
rlsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1))
  let iw :: Int
iw = if Int
key Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
4 then Int
rlsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2 else Int
rls
  IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3))
  (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3))) [Int
2 .. Int
np]
  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
1 Int
np
  let n2double :: Double
n2double = Int -> Double
toDbl Int
n2
  IOMatrix -> (Int, Int) -> Double -> IO ()
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
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3) Int
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
n2double))
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
1) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
      then do
        let l2 :: Double
l2 = Double
0.62054648267200632589046034361711
            r1' :: Double
r1' = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double -> Double
forall a. Floating a => a -> a
sqrt(Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l2))Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
r1')
        (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
r1') [Int
2 .. Int
np]
        MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
gms Int
3
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6)
        let r2 :: Double
r2 = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3
            s2 :: Double
s2 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
r2
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Double
s2
        (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Double
r2) [Int
2 .. Int
np]
        MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
3
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6)
      else do
        let r2 :: Double
r2 = (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt15) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
8Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)
            s2 :: Double
s2 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
r2
            l2 :: Double
l2 = Double
s2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
r2
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Double
s2
        (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Double
r2) [Int
2 .. Int
np]
        MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
np
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) ((Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3)Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
n2doubleDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l2))
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) ((Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3)Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
n2doubleDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
l1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l1))
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
5Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5))
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5))) [Int
2 .. Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
2 Int
np
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5))
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5))
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5))) [Int
3 .. Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
3 Int
o
    let tmp :: Double
tmp = Int -> Double
toDbl (Int
16Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n4)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) (- Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3) Int
5 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
3, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5) Int
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
4, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5) Int
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp)
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
2) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
    let tmp' :: Double
tmp' = Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndbl Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
14Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndbl Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
11
        u1 :: Double
u1 = (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt15) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp'
        v1 :: Double
v1 = Double
0.5Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u1)
        d1 :: Double
d1 = Double
v1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u1
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) Double
v1
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) Double
v1
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) Double
u1) [Int
3 .. Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Int
o
    let u2 :: Double
u2 = (Double
ndbl Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
7 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt15) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp'
        v2 :: Double
v2 = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2
        d2 :: Double
d2 = Double
v2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u2
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4) Double
v2
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4) Double
v2
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4) Double
u2) [Int
3 .. Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) Int
o
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
      then do
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) ((Double
155Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
sqrt15)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
1200)
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) ((Double
155Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt15)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
1200)
        IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
1, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) (Double
9Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
40)
      else
        if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3
          then do
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) ((Double
2665Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
14Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt15)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
37800)
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) ((Double
2665Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
14Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt15)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
37800)
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) (Double
10Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
189)
            MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) Int
0
          else do
            let r2 :: Double
r2 = (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt15) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
8Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)
                l2 :: Double
l2 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
r2
                n4dbl :: Double
n4dbl = Int -> Double
toDbl Int
n4
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
27Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndbl)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5)Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
13Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndbl)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                         (Double -> Int -> Double
pow Double
l1 Int
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
l1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
n4dbl))
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
27Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndbl)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5)Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l1Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
13Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndbl)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                         (Double -> Int -> Double
pow Double
l2 Int
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
n4dbl))
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5)Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
d2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n4dbl Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
d1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
d2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow Double
d1 Int
4))
            IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5)Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
d1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n4dbl Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
d2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
d1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow Double
d2 Int
4))
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
7Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7))
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7))) [Int
2 .. Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
4 Int
np
    let invnp7 :: Double
invnp7 = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7)
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7)
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
5 (Int
npInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7)) [Int
1, Int
2, Int
3]
    Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
np Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
3) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
      (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
6 (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div ((Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
np) Int
6)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
4) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3) Int
7 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
128Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n4Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5)))
    let tmp'' :: Double
tmp'' = Int -> Double
toDbl (Int
64Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
4) (-Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5) Int
7 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp'')) [Int
3, Int
4]
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
4) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7) Int
6 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp'')) [Int
5, Int
6, Int
7]
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
4) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
    let sg :: Double
sg = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
23328Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6)
        u5 :: Double
u5 = -Double
216 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sg Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
52212 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
6353 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
1934Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
27)))
        u6 :: Double
u6 = Double
1296 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sg Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
7884 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
1541 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
9))
        u7 :: Double
u7 = -Double
7776 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sg Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
8292 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
1139 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
3))Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndbl Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
7)
        p0 :: Double
p0 = Int -> Double
toDbl (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ -Int
144 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
142528 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
23073 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
115))
        p1 :: Double
p1 = Int -> Double
toDbl (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ -Int
12 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
6690556 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
2641189 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
245378 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
1495)))
        p2 :: Double
p2 = Int -> Double
toDbl (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ -Int
16 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
6503401 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
4020794Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
787281Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
47323Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
385))))
        p3 :: Double
p3 = Int -> Double
toDbl (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ -(Int
6386660 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
4411997Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
951821Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
61659Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
665))))Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7)
        a :: Double
a = Double
p2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
p3)
        p :: Double
p = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
p1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
p2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)
        q :: Double
q = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
p3) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
p3
        th :: Double
th = Double -> Double
forall a. Floating a => a -> a
acos(-Double
qDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt(-Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
p))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
        r :: Double
r = Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
sqrt(-Double
p)
        tp :: Double
tp = Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3
        a1 :: Double
a1 = -Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
th
        a3 :: Double
a3 = -Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos(Double
thDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
tp)
        a2 :: Double
a2 = -Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a3 -- a3 = -a + r*(cos(th+tp))
        npdbl :: Double
npdbl = Int -> Double
toDbl Int
np
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5) ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
npdbl)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5) ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4) Int
np
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
6) ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
npdbl)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
6) ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5) Int
np
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7) ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
npdbl)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7) ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a3)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
6) Int
np
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5)
               ((Double
u7Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u6Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u5)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a1 Int
5)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
6, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5)
               ((Double
u7Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u6Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u5)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a2 Int
5)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5)
               ((Double
u7Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u6Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
u5)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
a3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a3Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a3 Int
5)
    let invnp7' :: Double
invnp7' = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
8) (Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7')
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
8) (Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7')
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
8) Double
invnp7') [Int
3..Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7) Int
o
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
8, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5) (Double
10 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7) Int
6 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
729Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6))
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
9) (Double
5.5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7')
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
9) (Double
2.5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp7')
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
9) Double
invnp7') [Int
3..Int
np]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
8) (Int
npInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
9, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5) (Double
64 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7) Int
6 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
6561Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6))
    let tmp''' :: Double
tmp''' = Int -> Double
toDbl (Int
64Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6)
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
4, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5) (-Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5) Int
7 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp''')
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
7, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7) Int
6 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp''')
    let invnp9 :: Double
invnp9 = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
9)
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
9Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
7 Int
np
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
7Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
8 (Int
npInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
1, Int
2]
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
9 Int
o
    IOMatrix -> (Int, Int) -> Double -> IO ()
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
5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
2, Int
3]
    Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
np Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
3) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
      (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
    MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
10 (Int
oInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
    IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
6) (-Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
3) Int
9 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
1536Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6))
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
6)
                            (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
5) Int
9 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
512Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n6Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
7)))) [Int
3, Int
4]
    let tmp'''' :: Double
tmp'''' = Int -> Double
toDbl (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ Int
256Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n8
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
6) (-Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
7) Int
9 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp'''')) [Int
5, Int
6, Int
7]
    (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
6) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
9) Int
8 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tmp'''')) [Int
8..Int
11]
    Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
2) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
1..Int
4]
      (Int -> IO ()) -> [Int] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> IOMatrix -> (Int, Int) -> Double -> IO ()
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]
      MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
pts Int
11 (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div (Int
npInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2)) Int
24)
      IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
12, Int
iwInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
6) (Double -> Int -> Double
pow (Double
ndblDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
9) Int
8 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
256Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n8))
  Seq IO1dArray
rowsIO <- (Int -> IO IO1dArray) -> Seq Int -> IO (Seq IO1dArray)
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) ([Int] -> Seq Int
forall a. [a] -> Seq a
S.fromList [Int
2..Int
wts])
  Seq VectorD
rows <- (IO1dArray -> IO VectorD) -> Seq IO1dArray -> IO (Seq VectorD)
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 <- MVector (PrimState IO) Int -> IO (Vector Int)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
UV.unsafeFreeze IOVectorI
MVector (PrimState IO) Int
pts
  let ptsU :: VectorD
ptsU = (Int -> Double) -> Vector Int -> VectorD
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 = (VectorD -> Double) -> Seq VectorD -> Seq Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\VectorD
col -> Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0
                           ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) (VectorD -> VectorD
forall a. Unbox a => Vector a -> Vector a
UV.tail VectorD
ptsU) VectorD
col))
                  (Int -> Seq VectorD -> Seq VectorD
forall a. Int -> Seq a -> Seq a
S.take Int
rls Seq VectorD
cols)
      wcols :: Seq VectorD
wcols = (Int -> VectorD) -> Seq Int -> Seq VectorD
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Int
j -> Double -> VectorD -> VectorD
forall a. Unbox a => a -> Vector a -> Vector a
UV.cons (Seq Double -> Int -> Double
forall a. Seq a -> Int -> a
index Seq Double
row1 Int
j) (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
cols Int
j))
                   ([Int] -> Seq Int
forall a. [a] -> Seq a
S.fromList [Int
0..(Int
rlsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)])
      col1 :: VectorD
col1 = Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
wcols Int
0
      wcols2 :: Seq VectorD
wcols2 = VectorD -> Seq VectorD -> Seq VectorD
forall a. a -> Seq a -> Seq a
(S.<|) VectorD
col1
                      ((VectorD -> VectorD) -> Seq VectorD -> Seq VectorD
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\VectorD
col -> (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
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) (Int -> Seq VectorD -> Seq VectorD
forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
wcols))
      col2 :: VectorD
col2 = Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
wcols2 Int
1
      nb :: Double
nb = (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) VectorD
ptsU ((Double -> Double) -> VectorD -> VectorD
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 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) VectorD
ptsU ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
col2))
      wcol2 :: VectorD
wcol2 = (Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
ratio) VectorD
col2
      wcols3 :: Seq VectorD
wcols3 = VectorD -> Seq VectorD -> Seq VectorD
forall a. a -> Seq a -> Seq a
(S.<|) VectorD
col1 (VectorD -> Seq VectorD -> Seq VectorD
forall a. a -> Seq a -> Seq a
(S.<|) VectorD
wcol2 (Int -> Seq VectorD -> Seq VectorD
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 = Int -> VectorD -> Seq VectorD -> Seq VectorD
forall a. Int -> a -> Seq a -> Seq a
update (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) VectorD
wknew Seq VectorD
cols'
         where
          ptsW :: VectorD
ptsW = (Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
nb) ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) VectorD
ptsU (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
cols' (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)))
          slice :: Seq VectorD
slice = Int -> Seq VectorD -> Seq VectorD
forall a. Int -> Seq a -> Seq a
S.drop Int
1 (Int -> Seq VectorD -> Seq VectorD
forall a. Int -> Seq a -> Seq a
S.take (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Seq VectorD
cols')
          prod1 :: VectorD
prod1 = ([Double] -> VectorD
forall a. Unbox a => [a] -> Vector a
fromList ([Double] -> VectorD)
-> (Seq Double -> [Double]) -> Seq Double -> VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Seq Double -> [Double]
forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) (Seq Double -> VectorD) -> Seq Double -> VectorD
forall a b. (a -> b) -> a -> b
$ -- ou alors pas de vecteurs
                  (VectorD -> Double) -> Seq VectorD -> Seq Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 (VectorD -> Double) -> (VectorD -> VectorD) -> VectorD -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
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 = ([Double] -> VectorD
forall a. Unbox a => [a] -> Vector a
fromList ([Double] -> VectorD)
-> (Seq Double -> [Double]) -> Seq Double -> VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Seq Double -> [Double]
forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) (Seq Double -> VectorD) -> Seq Double -> VectorD
forall a b. (a -> b) -> a -> b
$
                  (VectorD -> Double) -> Seq VectorD -> Seq Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 (VectorD -> Double) -> (VectorD -> VectorD) -> VectorD -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) VectorD
prod1) Seq VectorD
rows'
          wk :: VectorD
wk = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
cols' (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) VectorD
prod2
          ratio' :: Double
ratio' = Double
nb Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                   (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) VectorD
ptsU ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
wk))
          wknew :: VectorD
wknew = (Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
ratio') VectorD
wk
      wcolsnew :: Seq VectorD
wcolsnew = (Seq VectorD -> Int -> Seq VectorD)
-> Seq VectorD -> [Int] -> Seq VectorD
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]
  (IOMatrix, Seq VectorD, [Int]) -> IO (IOMatrix, Seq VectorD, [Int])
forall (m :: * -> *) a. Monad m => a -> m a
return (IOMatrix
g, Seq VectorD -> Seq VectorD
transpose Seq VectorD
wcolsnew, Vector Int -> [Int]
forall a. Unbox a => Vector a -> [a]
toList (Vector Int -> [Int]) -> Vector Int -> [Int]
forall a b. (a -> b) -> a -> b
$ (Int -> Bool) -> Vector Int -> Vector Int
forall a. Unbox a => (a -> Bool) -> Vector a -> Vector Int
UV.findIndices (Int -> Int -> Bool
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 =
  (Int -> VectorD) -> Seq Int -> Seq VectorD
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Int
i -> ([Double] -> VectorD
forall a. Unbox a => [a] -> Vector a
fromList ([Double] -> VectorD)
-> (Seq Double -> [Double]) -> Seq Double -> VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Seq Double -> [Double]
forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) (Seq Double -> VectorD) -> Seq Double -> VectorD
forall a b. (a -> b) -> a -> b
$ (VectorD -> Double) -> Seq VectorD -> Seq Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (VectorD -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
UV.!Int
i) Seq VectorD
cols)
       ([Int] -> Seq Int
forall a. [a] -> Seq a
S.fromList [Int
0..(VectorD -> Int
forall a. Unbox a => Vector a -> Int
UV.length (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
cols Int
0) Int -> Int -> Int
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)) <- IOMatrix -> IO ((Int, Int), (Int, Int))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
mat
  IOVectorD
out <- Int -> IO (MVector (PrimState IO) Double)
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = MVector (PrimState IO) Double -> IO VectorD
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze IOVectorD
MVector (PrimState IO) Double
out
             | Bool
otherwise = do
              !Double
coef <- Int -> Double -> IO Double
innerstep Int
1 Double
0
              MVector (PrimState IO) Double -> Int -> Double -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
MVector (PrimState IO) Double
out (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Double
coef
              Int -> IO VectorD
step (Int
iInt -> Int -> Int
forall 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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
s
                               | Bool
otherwise = do
                                 Double
mat_ij <- IOMatrix -> (Int, Int) -> IO Double
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
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mat_ij Double -> Double -> Double
forall a. Num a => a -> a -> a
* (VectorD
x VectorD -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
UV.! (Int
jInt -> Int -> Int
forall 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 <- IO1dArray -> IO [Double]
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems IO1dArray
g
  [VectorD]
f_gPermuts <- ([Double] -> IO VectorD) -> [[Double]] -> IO [VectorD]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM ((VectorD -> VectorD) -> IO VectorD -> IO VectorD
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap VectorD -> VectorD
f (IO VectorD -> IO VectorD)
-> ([Double] -> IO VectorD) -> [Double] -> IO VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IOMatrix -> VectorD -> IO VectorD
matprod IOMatrix
vertex (VectorD -> IO VectorD)
-> ([Double] -> VectorD) -> [Double] -> IO VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> VectorD
forall a. Unbox a => [a] -> Vector a
fromList)
                     ([Double] -> [[Double]]
forall a. (Eq a, Ord a) => [a] -> [[a]]
permuteMultiset [Double]
gAsList)
  (Int, Int) -> [Double] -> IO IO1dArray
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> [e] -> m (a i e)
newListArray (Int
1, Int
nf)
               (VectorD -> [Double]
forall a. Unbox a => Vector a -> [a]
toList ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
scalar) ((VectorD -> VectorD -> VectorD) -> [VectorD] -> VectorD
forall a. (a -> a -> a) -> [a] -> a
foldl1' ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
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
_)) <- IOMatrix -> IO ((Int, Int), (Int, Int))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
  (Int, Int) -> (Int -> (Int, Int)) -> IOMatrix -> IO IO1dArray
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
extractRow :: IOMatrix -> Int -> IO IO1dArray
extractRow IOMatrix
m Int
i = do
  ((Int, Int)
_, (Int
_, Int
ncol)) <- IOMatrix -> IO ((Int, Int), (Int, Int))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
  (Int, Int) -> (Int -> (Int, Int)) -> IOMatrix -> IO IO1dArray
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) <- IO1dArray -> IO (Int, Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO1dArray
x1
  let n2 :: Int
n2 = VectorD -> Int
forall a. Unbox a => Vector a -> Int
UV.length VectorD
x2
  IOMatrix
out <- ((Int, Int), (Int, Int)) -> IO IOMatrix
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = IOMatrix -> IO IOMatrix
forall (m :: * -> *) a. Monad m => a -> m a
return IOMatrix
out
             | Bool
otherwise = do
                Double
x1_i <- IO1dArray -> Int -> IO Double
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n2 = Int -> IO IOMatrix
step (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                           | Bool
otherwise = do
                              IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
out (Int
i, Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (VectorD
x2 VectorD -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
UV.! Int
j))
                              Double -> Int -> IO IOMatrix
inner Double
x (Int
jInt -> Int -> Int
forall 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)) <- IOMatrix -> IO ((Int, Int), (Int, Int))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds ([IOMatrix] -> IOMatrix
forall a. [a] -> a
head [IOMatrix]
matrices)
  IOMatrix
out <- ((Int, Int), (Int, Int)) -> IO IOMatrix
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = IOMatrix -> IO IOMatrix
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n2Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = Int -> IO IOMatrix
step (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                         | Bool
otherwise = do
                           [Double]
coefs <- (IOMatrix -> IO Double) -> [IOMatrix] -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\IOMatrix
m -> IOMatrix -> (Int, Int) -> IO Double
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
                           IOMatrix -> (Int, Int) -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
out (Int
i, Int
j) ([Double] -> Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
coefs)
                           Int -> IO IOMatrix
inner (Int
jInt -> Int -> Int
forall 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 = VectorD -> Int
forall a. Unbox a => Vector a -> Int
UV.length (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
w Int
0)
  [IOMatrix]
toSum <- (Int -> IO IOMatrix) -> [Int] -> IO [IOMatrix]
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
kInt -> Int -> Int
forall 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 (Seq VectorD -> Int -> VectorD
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 <- (Int, Int) -> Double -> IO IO1dArray
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nfInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
             | Bool
otherwise = do
                Double
basval_i <- IO1dArray -> Int -> IO Double
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 = Double -> Double
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
                Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
rt Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Int
rls Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
3) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
                  IO1dArray -> Int -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (Double
rtDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
nmcp)
                Double
rgnerr_i <- IO1dArray -> Int -> IO Double
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
rgnerr Int
i
                IO1dArray -> Int -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double
errcofDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rgnerr_i) (Double
smallDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
nmbs))
                Int -> IO ()
step (Int
iInt -> Int -> Int
forall 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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Double, Double) -> IO (Double, Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Double
x, Double
y)
                                | Bool
otherwise = do
                                  Double
rule_ik <- IOMatrix -> (Int, Int) -> IO Double
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 <- IOMatrix -> (Int, Int) -> IO Double
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
rule (Int
i, Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
                                  let nmrl :: Double
nmrl = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
rule_ik) (Double -> Double
forall a. Num a => a -> a
abs Double
rule_ikm1)
                                  Double
rgnerr_i <- IO1dArray -> Int -> IO Double
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
rgnerr Int
i
                                  IO1dArray -> Int -> Double -> IO ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
nmrl Double
rgnerr_i)
                                  if Double
nmrl Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
smallDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
z Bool -> Bool -> Bool
&& Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
rls
                                    then Int -> Double -> Double -> Double -> IO (Double, Double)
inner (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double
nmrlDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
y) Double
x) Double
nmrl Double
z
                                    else Int -> Double -> Double -> Double -> IO (Double, Double)
inner (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) Double
x Double
nmrl Double
z
  Int -> IO ()
step Int
1
  (IO1dArray, IO1dArray) -> IO (IO1dArray, IO1dArray)
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)) <- IOMatrix -> IO ((Int, Int), (Int, Int))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
  IOVectorD
outIO <- Int -> IO (MVector (PrimState IO) Double)
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nrowInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
             | Bool
otherwise = do
                !Double
sum_i <- Int -> Double -> IO Double
inner Int
1 Double
0
                MVector (PrimState IO) Double -> Int -> Double -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
MVector (PrimState IO) Double
outIO (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Double
sum_i
                Int -> IO ()
step (Int
iInt -> Int -> Int
forall 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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
ncolInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl Int
ncol)
                           | Bool
otherwise = do
                             Double
coef <- IOMatrix -> (Int, Int) -> IO Double
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
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
coef)
  Int -> IO ()
step Int
1
  MVector (PrimState IO) Double -> IO VectorD
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze IOVectorD
MVector (PrimState IO) Double
outIO

array1dToVectorD :: IO1dArray -> IO VectorD
array1dToVectorD :: IO1dArray -> IO VectorD
array1dToVectorD IO1dArray
array = ([Double] -> VectorD) -> IO [Double] -> IO VectorD
forall (m :: * -> *) a b. Monad m => (a -> b) -> m a -> m b
(<$!>) [Double] -> VectorD
forall a. Unbox a => [a] -> Vector a
fromList (IO1dArray -> IO [Double]
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) -- pas de gain
getVectors :: Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
n IO3dArray
m Int
k Int
j1 Int
j2 = do
  VectorD
out1U <- (IO1dArray -> IO VectorD) -> IO IO1dArray -> IO VectorD
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD ((Int, Int) -> (Int -> (Int, Int, Int)) -> IO3dArray -> IO IO1dArray
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 <- (IO1dArray -> IO VectorD) -> IO IO1dArray -> IO VectorD
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD ((Int, Int) -> (Int -> (Int, Int, Int)) -> IO3dArray -> IO IO1dArray
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)
  -- out1 <- new n :: IO IOVectorD
  -- out2 <- new n :: IO IOVectorD
  -- let loop :: Int -> IO ()
  --     loop i | i == n+1 = return ()
  --            | otherwise = do
  --              coef1 <- readArray m (i, j1, k)
  --              unsafeWrite out1 (i-1) coef1
  --              coef2 <- readArray m (i, j2, k)
  --              unsafeWrite out2 (i-1) coef2
  --              loop (i+1)
  -- loop 1
  -- out1U <- UV.unsafeFreeze out1
  -- out2U <- UV.unsafeFreeze out2
  (VectorD, VectorD) -> IO (VectorD, VectorD)
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 <- ((Int, Int), (Int, Int))
-> ((Int, Int) -> (Int, Int, Int)) -> IO3dArray -> IO IOMatrix
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
ndInt -> Int -> Int
forall 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 = (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
forall a. Num a => a -> a
abs VectorD
fc)
  IOMatrix
frthdf <- ((Int, Int), (Int, Int)) -> Double -> IO IOMatrix
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
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Double
0 :: IO IOMatrix
  IOVectorI
iejeitjtisjsls <- Int -> IO (MVector (PrimState IO) Int)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
new Int
7 :: IO IOVectorI
  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4 Int
1
  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5 Int
2
  IOVectorD
dfmxdfnx <- Int -> Double -> IO (MVector (PrimState IO) Double)
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
               | Bool
otherwise = do
                  !Double
emx <- Int -> Double -> IO Double
inner (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
x
                  Int -> Double -> IO ()
step (Int
iInt -> Int -> Int
forall 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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2 = Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
y
                             | Bool
otherwise = do
                              VectorD
vi <- (IO1dArray -> IO VectorD) -> IO IO1dArray -> IO VectorD
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 <- (IO1dArray -> IO VectorD) -> IO IO1dArray -> IO VectorD
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 = (Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
5Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Int -> Double
toDbl Int
nd Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1))))
                                             ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
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 = (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
forall a. Num a => a -> a
abs VectorD
h)
                                  twoh :: VectorD
twoh = (Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
2) VectorD
h
                                  t1 :: VectorD
t1 = VectorD -> VectorD
f ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
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 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) VectorD
cn VectorD
twoh)
                                  t3 :: VectorD
t3 = (Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
6) VectorD
fc
                                  t4 :: VectorD
t4 = VectorD -> VectorD
f ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
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 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) VectorD
cn VectorD
h)
                                  t6 :: VectorD
t6 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((Double -> Double -> Double
forall a. Num a => a -> a -> a
*(-Double
4))(Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.)((Double -> Double) -> Double -> Double)
-> (Double -> Double -> Double) -> Double -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) VectorD
t4 VectorD
t5
                                  tsum :: VectorD
tsum = (VectorD -> VectorD -> VectorD) -> [VectorD] -> VectorD
forall a. (a -> a -> a) -> [a] -> a
foldl1' ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) [VectorD
t1, VectorD
t2, VectorD
t3, VectorD
t6]
                                  dfr1 :: Double
dfr1 = (Double -> Double -> Double) -> Double -> VectorD -> Double
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
forall a. Num a => a -> a
abs VectorD
tsum)
                                  dfr2 :: Double
dfr2 = if Double
dfmdDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
dfr1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
8 Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
dfmd then Double
0 else Double
dfr1
                                  dfr3 :: Double
dfr3 = Double
dfr2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ewd
                              Double
dfmx <- MVector (PrimState IO) Double -> Int -> IO Double
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
0
                              if Double
dfr3 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
dfmx
                                then do
                                 Int
is <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4
                                 Int
js <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5
                                 MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
2 Int
is
                                 MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
3 Int
js
                                 MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4 Int
i
                                 MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5 Int
j
                                 MVector (PrimState IO) Double -> Int -> Double -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
1 Double
dfmx
                                 MVector (PrimState IO) Double -> Int -> Double -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
0 Double
dfr3
                                else do
                                 Double
dfnx <- MVector (PrimState IO) Double -> Int -> IO Double
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
1
                                 Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
dfr3 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
dfnx) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
                                  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
2 Int
i
                                  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
3 Int
j
                                  MVector (PrimState IO) Double -> Int -> Double -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
1 Double
dfr3
                              IOMatrix -> (Int, Int) -> Double -> IO ()
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 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
y
                                then do
                                  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
0 Int
i
                                  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
1 Int
j
                                  Int -> Double -> IO Double
inner (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
ewd
                                else Int -> Double -> IO Double
inner (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
y
  Int -> Double -> IO ()
step Int
1 Double
0
  Double
dfmx <- MVector (PrimState IO) Double -> Int -> IO Double
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
0
  Double
dfnx <- MVector (PrimState IO) Double -> Int -> IO Double
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
1
  let nregions :: Int
nregions = if Double
dfnx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
dfmxDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cuttf then Int
4 else Int
3
  if Double
dfnx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
dfmxDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cuttf
    then () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
    else
      if Double
dfmx Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
        then do
          Int
ie <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
0
          Int
je <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
1
          MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4 Int
ie
          MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2 = Int -> IO Int
forall (m :: * -> *) a. Monad m => a -> m a
return Int
ls
                            | Bool
otherwise = do
                              Int
is <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4
                              Int
js <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5
                              if (Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
is) Bool -> Bool -> Bool
&& (Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
js)
                                then do
                                  let it :: Int
it = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Int
l, Int
is, Int
js]
                                      jt :: Int
jt = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Int
l, Int
is, Int
js]
                                  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
2 Int
it
                                  MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
3 Int
jt
                                  let lt :: Int
lt = Int
isInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
jsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
lInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
itInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
jt
                                  Double
dfr1 <- IOMatrix -> (Int, Int) -> IO Double
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 <- IOMatrix -> (Int, Int) -> IO Double
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 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dfr2
                                  if Double
dfr Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
x
                                    then Int -> Double -> Int -> IO Int
loop (Int
lInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
dfr Int
l
                                    else Int -> Double -> Int -> IO Int
loop (Int
lInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Double
x Int
ls
                                else Int -> Double -> Int -> IO Int
loop (Int
lInt -> Int -> Int
forall 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
          MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
6 Int
ls
          Int
is <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4
          Int
js <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5
          Double
difil <- IOMatrix -> (Int, Int) -> IO Double
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
is Int
ls, Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
is Int
ls)
          Double
diflj <- IOMatrix -> (Int, Int) -> IO Double
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
js Int
ls, Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
js Int
ls)
          let dfnx' :: Double
dfnx' = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
difil Double
diflj
          MVector (PrimState IO) Double -> Int -> Double -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
MVector (PrimState IO) Double
dfmxdfnx Int
1 Double
dfnx'
          Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
dfmxDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cuttb Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
dfnx' Bool -> Bool -> Bool
&& Double
difil Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
diflj) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
            Int
is' <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4
            Int
js' <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5
            Int
it <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
2
            MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
2 Int
is'
            MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4 Int
js'
            MVector (PrimState IO) Int -> Int -> Int -> IO ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5 Int
it
  IO3dArray
vrts2 <- ((Int, Int, Int), (Int, Int, Int))
-> ((Int, Int, Int) -> (Int, Int, Int))
-> IO3dArray
-> IO IO3dArray
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
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nregionsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
                      (\(Int
i, Int
j, Int
k) -> (Int
i, Int
j, if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
sbs then Int
k else Int
top)) IO3dArray
vrts
  -- vrts2 <- newArray_ ((1, 1, 1), (nd, nd+1, sbs+nregions-1)) :: IO IO3dArray
  -- let go1 :: Int -> IO ()
  --     go1 i | i == nd+1 = return ()
  --           | otherwise = do
  --             let go2 j | j == nd+2 = go1 (i+1)
  --                       | otherwise = do
  --                         let go3 k | k == sbs+nregions = go2 (j+1)
  --                                   | otherwise = do
  --                                     case k <= sbs of
  --                                       True -> do
  --                                         coef <- readArray vrts (i, j, k)
  --                                         writeArray vrts2 (i, j, k) coef
  --                                       False -> do
  --                                         coef <- readArray v (i, j)
  --                                         writeArray vrts2 (i, j, k) coef
  --                                     go3 (k+1)
  --                         go3 1
  --             go2 1
  -- go1 1
  Int
is <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
4
  Int
js <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
5
  VectorD
vti <- (IO1dArray -> IO VectorD) -> IO IO1dArray -> IO VectorD
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 <- (IO1dArray -> IO VectorD) -> IO IO1dArray -> IO VectorD
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
4
    then do
      let vt :: VectorD
vt = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
0.5)(Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.)((Double -> Double) -> Double -> Double)
-> (Double -> Double -> Double) -> Double -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Double -> Double
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
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) VectorD
vt
      IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
is, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) VectorD
vti
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
is, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) (Int
is, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) VectorD
vt
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) (Int
js, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) VectorD
vtj
      Int
it <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
2
      Int
jt <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
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
      -- vtit <- (=<<) array1dToVectorD (mapIndices (1, nd) (\i -> (i, it, top)) vrts2)
      -- vtjt <- (=<<) array1dToVectorD (mapIndices (1, nd) (\i -> (i, jt, top)) vrts2)
      -- vtit <- array1dToVectorD vtitIO -- (=<<) (return . fromList) (getElems vtitIO)
      -- vtjt <- array1dToVectorD vtjtIO -- (=<<) (return . fromList) (getElems vtjtIO)
      let vtt :: VectorD
vtt = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
0.5)(Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.)((Double -> Double) -> Double -> Double)
-> (Double -> Double -> Double) -> Double -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Double -> Double
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
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) VectorD
vtt
      IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
jt, Int
sbsInt -> Int -> Int
forall 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
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Int
it Int
jt
      -- vti2 <- (=<<) array1dToVectorD (mapIndices (1, nd) (\i -> (i, it, sbs+2)) vrts2)
      -- vtj2 <- (=<<) array1dToVectorD (mapIndices (1, nd) (\i -> (i, jt, sbs+2)) vrts2) -- :: IO IO1dArray
      -- vti2 <- (=<<) (return . fromList) (getElems vti2IO)
      -- vtj2 <- (=<<) (return . fromList) (getElems vtj2IO)
      let vt2 :: VectorD
vt2 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
0.5)(Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.)((Double -> Double) -> Double -> Double)
-> (Double -> Double -> Double) -> Double -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) VectorD
vti2 VectorD
vtj2
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
jt, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) (Int
it, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) VectorD
vt2
      IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
jt, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) VectorD
vtj2
    else do
      let vt :: VectorD
vt = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
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
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
y)Double -> Double -> Double
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
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) VectorD
vt
      if Double
dfmxDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cuttf Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
dfnx
        then do
          IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
is, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
          IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
js, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
          Int
ls <- MVector (PrimState IO) Int -> Int -> IO Int
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
MVector (PrimState IO) Int
iejeitjtisjsls Int
6
          (VectorD
vtj1, VectorD
vtl1) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 (Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
js Int
ls
          -- vtj1 <- (=<<) array1dToVectorD (mapIndices (1, nd) (\i -> (i, js, sbs+1)) vrts2)
          -- vtl1 <- (=<<) array1dToVectorD (mapIndices (1, nd) (\i -> (i, ls, sbs+1)) vrts2)
          -- vtj1 <- (=<<) (return . fromList) (getElems vtj1IO)
          -- vtl1 <- (=<<) (return . fromList) (getElems vtl1IO)
          let vt1 :: VectorD
vt1 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
0.5)(Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.)((Double -> Double) -> Double -> Double)
-> (Double -> Double -> Double) -> Double -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) VectorD
vtj1 VectorD
vtl1
          IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
ls, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
js, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) VectorD
vt1
          IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
ls, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) VectorD
vtl1
        else do
          let vv :: VectorD
vv = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
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
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y)Double -> Double -> Double
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
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
is, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) VectorD
vv
          IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
js, Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
  (Int, IO3dArray) -> IO (Int, IO3dArray)
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
                 | Bool
otherwise = do
                   IO3dArray -> (Int, Int, Int) -> Double -> IO ()
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) (VectorD -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
                   Int -> IO ()
loop (Int
iInt -> Int -> Int
forall 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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 = () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
                 | Bool
otherwise = do
                   IO3dArray -> (Int, Int, Int) -> Double -> IO ()
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) (VectorD -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
                   IO3dArray -> (Int, Int, Int) -> Double -> IO ()
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) (VectorD -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
                   Int -> IO ()
loop (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
      Int -> IO ()
loop Int
1

-- | Number of evaluations for each subregion
smpchc :: Int -> Int -> Int
smpchc :: Int -> Int -> Int
smpchc Int
nd Int
key
  | Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3
    = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div ((Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)) Int
6 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
  | Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
    = Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nd Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3
  | Int
key Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
    = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div ((Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)) Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
  | Bool
otherwise
    = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div ((Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)) Int
24 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
5 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int -> Int -> Int
forall a. Integral a => a -> a -> a
div ((Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Int
2

-- | Checks validity of parameters
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 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
er Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0
    = Int
5
  | Int
nf Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1
    = Int
4
  | Int
nd Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2
    = Int
3
  | Int
key Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Int
key Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
4
    = Int
2
  | Int
mxfs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
sbs Int -> Int -> Int
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)) <- IO3dArray -> IO ((Int, Int, Int), (Int, Int, Int))
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 Int -> Int -> Int
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 -> [Char] -> IO (VectorD, VectorD, Int, Bool)
forall a. HasCallStack => [Char] -> a
error [Char]
"integration rule must be between 1 and 4"
    Int
3 -> [Char] -> IO (VectorD, VectorD, Int, Bool)
forall a. HasCallStack => [Char] -> a
error [Char]
"dimension must be at least 2"
    Int
4 -> [Char] -> IO (VectorD, VectorD, Int, Bool)
forall a. HasCallStack => [Char] -> a
error [Char]
"number of components must be at least 1"
    Int
5 -> [Char] -> IO (VectorD, VectorD, Int, Bool)
forall a. HasCallStack => [Char] -> a
error [Char]
"requested errors must be positive"
    Int
_ -> [Char] -> IO (VectorD, VectorD, Int, Bool)
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 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ndInt -> Int -> Int
forall 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
  [[[Double]]]
simplices <- (Int -> IO [[Double]]) -> [Int] -> IO [[[Double]]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO3dArray -> Int -> Int -> IO [[Double]]
toSimplex IO3dArray
vrts (Int
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) [Int
1..Int
sbs]
  let vol :: Seq Double
vol = [Double] -> Seq Double
forall a. [a] -> Seq a
S.fromList ([Double] -> Seq Double) -> [Double] -> Seq Double
forall a b. (a -> b) -> a -> b
$ ([[Double]] -> Double) -> [[[Double]]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [[Double]] -> Double
simplexVolume [[[Double]]]
simplices
      nv :: Int
nv = Int
sbsInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
rcls
  Seq IOMatrix
matrices <- (Int -> IO IOMatrix) -> Seq Int -> IO (Seq IOMatrix)
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM
              (\Int
k -> ((Int, Int), (Int, Int))
-> ((Int, Int) -> (Int, Int, Int)) -> IO3dArray -> IO IOMatrix
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
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) (\(Int
i, Int
j) -> (Int
i, Int
j, Int
k)) IO3dArray
vrts)
              ([Int] -> Seq Int
forall a. [a] -> Seq a
S.fromList [Int
1..Int
sbs])
  Seq (IO1dArray, IO1dArray)
br <- ((IOMatrix, Double) -> IO (IO1dArray, IO1dArray))
-> Seq (IOMatrix, Double) -> IO (Seq (IO1dArray, IO1dArray))
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) (Seq IOMatrix -> Seq Double -> Seq (IOMatrix, Double)
forall a b. Seq a -> Seq b -> Seq (a, b)
S.zip Seq IOMatrix
matrices Seq Double
vol)
  Seq VectorD
aes <- ((IO1dArray, IO1dArray) -> IO VectorD)
-> Seq (IO1dArray, IO1dArray) -> IO (Seq VectorD)
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorD(IO1dArray -> IO VectorD)
-> ((IO1dArray, IO1dArray) -> IO1dArray)
-> (IO1dArray, IO1dArray)
-> IO VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(IO1dArray, IO1dArray) -> IO1dArray
forall a b. (a, b) -> b
snd) Seq (IO1dArray, IO1dArray)
br
  Seq VectorD
vls <- ((IO1dArray, IO1dArray) -> IO VectorD)
-> Seq (IO1dArray, IO1dArray) -> IO (Seq VectorD)
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorD(IO1dArray -> IO VectorD)
-> ((IO1dArray, IO1dArray) -> IO1dArray)
-> (IO1dArray, IO1dArray)
-> IO VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(IO1dArray, IO1dArray) -> IO1dArray
forall a b. (a, b) -> a
fst) Seq (IO1dArray, IO1dArray)
br
  let vl :: VectorD
vl = (VectorD -> VectorD -> VectorD) -> Seq VectorD -> VectorD
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) Seq VectorD
vls
      ae :: VectorD
ae = (VectorD -> VectorD -> VectorD) -> Seq VectorD -> VectorD
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
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'Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
dfcostInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
4Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
rcls Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
mxfs =
                     (VectorD, VectorD, Int, Bool) -> IO (VectorD, VectorD, Int, Bool)
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 = (VectorD -> Double) -> Seq VectorD -> Seq Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap VectorD -> Double
forall a. (Unbox a, Ord a) => Vector a -> a
UV.maximum Seq VectorD
aes'
                          imax :: Int
imax = Maybe Int -> Int
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Int -> Int) -> Maybe Int -> Int
forall a b. (a -> b) -> a -> b
$ (Double -> Bool) -> Seq Double -> Maybe Int
forall a. (a -> Bool) -> Seq a -> Maybe Int
S.findIndexL (Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Seq Double -> Double
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum Seq Double
maxs) Seq Double
maxs
                          vl0 :: VectorD
vl0 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
vl' (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
vls' Int
imax)
                          ae0 :: VectorD
ae0 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
ae' (Seq VectorD -> Int -> VectorD
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
imaxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
sbs' IO3dArray
vrts'
                      let vi :: Double
vi = Seq Double -> Int -> Double
forall a. Seq a -> Int -> a
index Seq Double
vol' Int
imax Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl Int
nregions
                          nv2 :: Int
nv2 = Int
nv' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nregionsInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
rcls Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
dfcost -- nregions*rcls ?
                          sbs2 :: Int
sbs2 = Int
sbs' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
nregionsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1
                      Seq IOMatrix
matrices2 <- (Int -> IO IOMatrix) -> Seq Int -> IO (Seq IOMatrix)
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM
                                   (\Int
k -> ((Int, Int), (Int, Int))
-> ((Int, Int) -> (Int, Int, Int)) -> IO3dArray -> IO IOMatrix
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
ndInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
                                          (\(Int
i, Int
j) -> (Int
i, Int
j, Int
k)) IO3dArray
vrts2)
                                   ([Int] -> Seq Int
forall a. [a] -> Seq a
S.fromList ((Int
imaxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[(Int
sbs'Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)..Int
sbs2]))
                      Seq (IO1dArray, IO1dArray)
br2 <- (IOMatrix -> IO (IO1dArray, IO1dArray))
-> Seq IOMatrix -> IO (Seq (IO1dArray, IO1dArray))
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 <- ((IO1dArray, IO1dArray) -> IO VectorD)
-> Seq (IO1dArray, IO1dArray) -> IO (Seq VectorD)
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorD(IO1dArray -> IO VectorD)
-> ((IO1dArray, IO1dArray) -> IO1dArray)
-> (IO1dArray, IO1dArray)
-> IO VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(IO1dArray, IO1dArray) -> IO1dArray
forall a b. (a, b) -> b
snd) Seq (IO1dArray, IO1dArray)
br2
                      Seq VectorD
basvals <- ((IO1dArray, IO1dArray) -> IO VectorD)
-> Seq (IO1dArray, IO1dArray) -> IO (Seq VectorD)
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorD(IO1dArray -> IO VectorD)
-> ((IO1dArray, IO1dArray) -> IO1dArray)
-> (IO1dArray, IO1dArray)
-> IO VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(IO1dArray, IO1dArray) -> IO1dArray
forall a b. (a, b) -> a
fst) Seq (IO1dArray, IO1dArray)
br2
                      let vl2 :: VectorD
vl2 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) VectorD
vl0
                                ((VectorD -> VectorD -> VectorD) -> Seq VectorD -> VectorD
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) Seq VectorD
basvals)
                          ae2 :: VectorD
ae2 = (Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) VectorD
ae0
                                ((VectorD -> VectorD -> VectorD) -> Seq VectorD -> VectorD
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 ((Double -> Double -> Double) -> VectorD -> VectorD -> VectorD
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) Seq VectorD
rgnerrs)
                          aes2 :: Seq VectorD
aes2 = Int -> VectorD -> Seq VectorD -> Seq VectorD
forall a. Int -> a -> Seq a -> Seq a
update Int
imax (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
rgnerrs Int
0) Seq VectorD
aes'
                                 Seq VectorD -> Seq VectorD -> Seq VectorD
forall a. Seq a -> Seq a -> Seq a
>< Int -> Seq VectorD -> Seq VectorD
forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
rgnerrs
                          vls2 :: Seq VectorD
vls2 = Int -> VectorD -> Seq VectorD -> Seq VectorD
forall a. Int -> a -> Seq a -> Seq a
update Int
imax (Seq VectorD -> Int -> VectorD
forall a. Seq a -> Int -> a
index Seq VectorD
basvals Int
0) Seq VectorD
vls'
                                 Seq VectorD -> Seq VectorD -> Seq VectorD
forall a. Seq a -> Seq a -> Seq a
>< Int -> Seq VectorD -> Seq VectorD
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 = Int -> Double -> Seq Double -> Seq Double
forall a. Int -> a -> Seq a -> Seq a
update Int
imax Double
vi Seq Double
vol'
                                 Seq Double -> Seq Double -> Seq Double
forall a. Seq a -> Seq a -> Seq a
>< Int -> Double -> Seq Double
forall a. Int -> a -> Seq a
S.replicate (Int
nregionsInt -> Int -> Int
forall 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
  -- dans le code R il fait rowSums mais ça me semble inutile
  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 = (Double -> Bool) -> VectorD -> Bool
forall a. Unbox a => (a -> Bool) -> Vector a -> Bool
UV.any (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
ea (VectorD -> Double
forall a. (Unbox a, Ord a) => Vector a -> a
UV.maximum ((Double -> Double) -> VectorD -> VectorD
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map ((Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
er)(Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Double
forall a. Num a => a -> a
abs) VectorD
v))) VectorD
a


toSimplex :: IO3dArray -> Int -> Int -> IO Simplex
toSimplex :: IO3dArray -> Int -> Int -> IO [[Double]]
toSimplex IO3dArray
m Int
n Int
k = do
  ((Int, Int, Int)
_, (Int
nrow, Int
_, Int
_)) <- IO3dArray -> IO ((Int, Int, Int), (Int, 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 = (Int, Int) -> (Int -> (Int, Int, Int)) -> IO3dArray -> IO IO1dArray
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 <- (Int -> IO IO1dArray) -> [Int] -> IO [IO1dArray]
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]
  (IO1dArray -> IO [Double]) -> [IO1dArray] -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM IO1dArray -> IO [Double]
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems [IO1dArray]
columns