{-# LANGUAGE BangPatterns #-}
module Numeric.Integration.SimplexCubature.Internal
  where
import           Data.Array.IO                       (IOUArray, getBounds, getElems
                                                     , mapIndices, newArray, newArray_
                                                     , newListArray, readArray
                                                     , writeArray)
import           Data.List                           (foldl', foldl1', sort)
import           Control.Monad                       ((<$!>), when)
import qualified Data.Foldable                       as DF
import           Data.Maybe                          (fromJust)
import           Data.Sequence                       (Seq, index, update, (><))
import qualified Data.Sequence                       as S
import           Data.Vector.Unboxed                 (Vector, fromList, toList
                                                     , unsafeFreeze)
import qualified Data.Vector.Unboxed                 as UV
import           Data.Vector.Unboxed.Mutable         (IOVector, new, unsafeRead
                                                     , unsafeWrite)
import qualified Data.Vector.Unboxed.Mutable         as UMV
import           Numeric.Integration.Simplex.Simplex (Simplex, simplexVolume)

type IOMatrix = IOUArray (Int, Int) Double
type IO1dArray = IOUArray Int Double
type IO3dArray = IOUArray (Int, Int, Int) Double
type IOVectorD = IOVector Double
type IOVectorI = IOVector Int
type VectorD = Vector Double

toDbl :: Int -> Double
toDbl :: Int -> Double
toDbl = forall a b. (Integral a, Num b) => a -> b
fromIntegral

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

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

square :: Double -> Double
square :: Double -> Double
square Double
x = Double
xforall a. Num a => a -> a -> a
*Double
x

smprms :: Int -> Int -> IO (IOMatrix, Seq VectorD, [Int])
smprms :: Int -> Int -> IO (IOMatrix, Seq VectorD, [Int])
smprms Int
n Int
key = do
  let (Int
rls, Int
gms, Int
wts) | Int
key forall a. Eq a => a -> a -> Bool
== Int
1 = (Int
3, Int
2, Int
3) :: (Int, Int, Int)
                      | Int
key forall a. Eq a => a -> a -> Bool
== Int
2 = (Int
5, Int
4, Int
6) :: (Int, Int, Int)
                      | Int
key forall a. Eq a => a -> a -> Bool
== Int
3 = (Int
7, Int
7, Int
11) :: (Int, Int, Int)
                      | Int
key forall a. Eq a => a -> a -> Bool
== Int
4 = (if Int
n forall a. Eq a => a -> a -> Bool
== Int
2 then (Int
7, Int
11, Int
20) else (Int
7, Int
12, Int
21))
                                    :: (Int, Int, Int)
                      | Bool
otherwise = forall a. HasCallStack => [Char] -> a
error [Char]
"this should not happen"
  IOMatrix
w <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
wts, Int
rls)) Double
0 :: IO IOMatrix
  IOVectorI
pts <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
UMV.replicate Int
wts Int
0 :: IO IOVectorI
  IOMatrix
g <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
nforall a. Num a => a -> a -> a
+Int
1, Int
wts)) Double
0 :: IO IOMatrix
  let np :: Int
np = Int
nforall a. Num a => a -> a -> a
+Int
1
      n2 :: Int
n2 = Int
np forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
2)
      n4 :: Int
n4 = Int
n2 forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
3) forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
4)
      n6 :: Int
n6 = Int
n4 forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
5) forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
6)
      n8 :: Int
n8 = Int
n6 forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
7) forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
+Int
8)
      o :: Int
o = forall a. Integral a => a -> a -> a
div (Int
nforall a. Num a => a -> a -> a
*Int
np) Int
2
      ndbl :: Double
ndbl = Int -> Double
toDbl Int
n
      sqrt15 :: Double
sqrt15 = Double
3.872983346207416885179265399782399611
      r1 :: Double
r1 = (Double
ndbl forall a. Num a => a -> a -> a
+ Double
4 forall a. Num a => a -> a -> a
- Double
sqrt15) forall a. Fractional a => a -> a -> a
/ (Double
ndblforall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
+ Double
8forall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
+ Double
1)
      s1 :: Double
s1 = Double
1 forall a. Num a => a -> a -> a
- Double
ndblforall a. Num a => a -> a -> a
*Double
r1
      l1 :: Double
l1 = Double
s1 forall a. Num a => a -> a -> a
- Double
r1
  forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
1) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
1))) [Int
1 .. Int
np]
  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
0 Int
1
  forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
1) Double
s1
  forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
1) Double
r1) [Int
2 .. Int
np]
  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
gms Int
np
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Ord a => a -> a -> Bool
< Int
4) forall a b. (a -> b) -> a -> b
$ do
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
1, Int
rls) Double
1
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
rlsforall a. Num a => a -> a -> a
-Int
1) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
1))
  let iw :: Int
iw = if Int
key forall a. Ord a => a -> a -> Bool
< Int
4 then Int
rlsforall a. Num a => a -> a -> a
-Int
2 else Int
rls
  forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
2) (Double
3forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3))
  forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
2) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3))) [Int
2 .. Int
np]
  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
1 Int
np
  let n2double :: Double
n2double = Int -> Double
toDbl Int
n2
  forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iw) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
2 forall a. Fractional a => a -> a -> a
/ (Double
4forall a. Num a => a -> a -> a
*Double
n2double))
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Ord a => a -> a -> Bool
> Int
1) forall a b. (a -> b) -> a -> b
$ do
    if Int
n forall a. Eq a => a -> a -> Bool
== Int
2
      then do
        let l2 :: Double
l2 = Double
0.62054648267200632589046034361711
            r1' :: Double
r1' = (Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt(Double
0.5 forall a. Num a => a -> a -> a
- Double
l2forall a. Num a => a -> a -> a
*Double
l2))forall a. Fractional a => a -> a -> a
/Double
3
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
1) (Double
1 forall a. Num a => a -> a -> a
- Double
2forall a. Num a => a -> a -> a
*Double
r1')
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
1) Double
r1') [Int
2 .. Int
np]
        forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
gms Int
3
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
1) (Double
1forall a. Fractional a => a -> a -> a
/Double
6)
        let r2 :: Double
r2 = (Double
1forall a. Num a => a -> a -> a
-Double
l2)forall a. Fractional a => a -> a -> a
/Double
3
            s2 :: Double
s2 = Double
1 forall a. Num a => a -> a -> a
- Double
2forall a. Num a => a -> a -> a
*Double
r2
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
s2
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
r2) [Int
2 .. Int
np]
        forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
1) Int
3
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
1) (Double
1forall a. Fractional a => a -> a -> a
/Double
6)
      else do
        let r2 :: Double
r2 = (Double
ndblforall a. Num a => a -> a -> a
+Double
4forall a. Num a => a -> a -> a
+Double
sqrt15) forall a. Fractional a => a -> a -> a
/ (Double
ndblforall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
8forall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
1)
            s2 :: Double
s2 = Double
1 forall a. Num a => a -> a -> a
- Double
ndblforall a. Num a => a -> a -> a
*Double
r2
            l2 :: Double
l2 = Double
s2 forall a. Num a => a -> a -> a
- Double
r2
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
s2
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
2) Double
r2) [Int
2 .. Int
np]
        forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
1) Int
np
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
1) ((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3)forall a. Num a => a -> a -> a
-Double
l1)forall a. Fractional a => a -> a -> a
/(Double
n2doubleforall a. Num a => a -> a -> a
*(Double
l2forall a. Num a => a -> a -> a
-Double
l1)forall a. Num a => a -> a -> a
*Double
l2forall a. Num a => a -> a -> a
*Double
l2))
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
1) ((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
3)forall a. Num a => a -> a -> a
-Double
l2)forall a. Fractional a => a -> a -> a
/(Double
n2doubleforall a. Num a => a -> a -> a
*(Double
l1forall a. Num a => a -> a -> a
-Double
l2)forall a. Num a => a -> a -> a
*Double
l1forall a. Num a => a -> a -> a
*Double
l1))
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
3) (Double
5forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
3) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))) [Int
2 .. Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
2 Int
np
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
4) (Double
3forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
4) (Double
3forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
4) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5))) [Int
3 .. Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
3 Int
o
    let tmp :: Double
tmp = Int -> Double
toDbl (Int
16forall a. Num a => a -> a -> a
*Int
n4)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
2) (- Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
5 forall a. Fractional a => a -> a -> a
/ Double
tmp)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
2) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
4 forall a. Fractional a => a -> a -> a
/ Double
tmp)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
2) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
4 forall a. Fractional a => a -> a -> a
/ Double
tmp)
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Ord a => a -> a -> Bool
> Int
2) forall a b. (a -> b) -> a -> b
$ do
    let tmp' :: Double
tmp' = Double
ndblforall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
+ Double
14forall a. Num a => a -> a -> a
*Double
ndbl forall a. Num a => a -> a -> a
- Double
11
        u1 :: Double
u1 = (Double
ndblforall a. Num a => a -> a -> a
+Double
7forall a. Num a => a -> a -> a
+Double
2forall a. Num a => a -> a -> a
*Double
sqrt15) forall a. Fractional a => a -> a -> a
/ Double
tmp'
        v1 :: Double
v1 = Double
0.5forall a. Num a => a -> a -> a
*(Double
1forall a. Num a => a -> a -> a
-(Double
ndblforall a. Num a => a -> a -> a
-Double
1)forall a. Num a => a -> a -> a
*Double
u1)
        d1 :: Double
d1 = Double
v1 forall a. Num a => a -> a -> a
- Double
u1
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
3) Double
v1
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
3) Double
v1
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
3) Double
u1) [Int
3 .. Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
2) Int
o
    let u2 :: Double
u2 = (Double
ndbl forall a. Num a => a -> a -> a
+ Double
7 forall a. Num a => a -> a -> a
- Double
2forall a. Num a => a -> a -> a
*Double
sqrt15) forall a. Fractional a => a -> a -> a
/ Double
tmp'
        v2 :: Double
v2 = (Double
1forall a. Num a => a -> a -> a
-(Double
ndblforall a. Num a => a -> a -> a
-Double
1)forall a. Num a => a -> a -> a
*Double
u2)forall a. Fractional a => a -> a -> a
/Double
2
        d2 :: Double
d2 = Double
v2 forall a. Num a => a -> a -> a
- Double
u2
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
4) Double
v2
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
4) Double
v2
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
j, Int
gmsforall a. Num a => a -> a -> a
+Int
4) Double
u2) [Int
3 .. Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
3) Int
o
    if Int
n forall a. Eq a => a -> a -> Bool
== Int
2
      then do
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
155forall a. Num a => a -> a -> a
-Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
1200)
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
155forall a. Num a => a -> a -> a
+Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
1200)
        forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
3) (Double
9forall a. Fractional a => a -> a -> a
/Double
40)
      else
        if Int
n forall a. Eq a => a -> a -> Bool
== Int
3
          then do
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
2665forall a. Num a => a -> a -> a
+Double
14forall a. Num a => a -> a -> a
*Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
37800)
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
3) ((Double
2665forall a. Num a => a -> a -> a
-Double
14forall a. Num a => a -> a -> a
*Double
sqrt15)forall a. Fractional a => a -> a -> a
/Double
37800)
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
3) (Double
10forall a. Fractional a => a -> a -> a
/Double
189)
            forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
3) Int
0
          else do
            let r2 :: Double
r2 = (Double
ndblforall a. Num a => a -> a -> a
+Double
4forall a. Num a => a -> a -> a
+Double
sqrt15) forall a. Fractional a => a -> a -> a
/ (Double
ndblforall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
8forall a. Num a => a -> a -> a
*Double
ndblforall a. Num a => a -> a -> a
+Double
1)
                l2 :: Double
l2 = Double
1 forall a. Num a => a -> a -> a
- (Double
ndblforall a. Num a => a -> a -> a
+Double
1)forall a. Num a => a -> a -> a
*Double
r2
                n4dbl :: Double
n4dbl = Int -> Double
toDbl Int
n4
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
1, Int
iwforall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2forall a. Num a => a -> a -> a
*(Double
27forall a. Num a => a -> a -> a
-Double
ndbl)forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
l2forall a. Num a => a -> a -> a
*(Double
13forall a. Num a => a -> a -> a
-Double
ndbl)) forall a. Fractional a => a -> a -> a
/
                         (Double -> Int -> Double
pow Double
l1 Int
4 forall a. Num a => a -> a -> a
* (Double
l1forall a. Num a => a -> a -> a
-Double
l2) forall a. Num a => a -> a -> a
* Double
n4dbl))
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2forall a. Num a => a -> a -> a
*(Double
27forall a. Num a => a -> a -> a
-Double
ndbl)forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
l1forall a. Num a => a -> a -> a
*(Double
13forall a. Num a => a -> a -> a
-Double
ndbl)) forall a. Fractional a => a -> a -> a
/
                         (Double -> Int -> Double
pow Double
l2 Int
4 forall a. Num a => a -> a -> a
*(Double
l2forall a. Num a => a -> a -> a
-Double
l1)forall a. Num a => a -> a -> a
*Double
n4dbl))
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
3, Int
iwforall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
d2) forall a. Fractional a => a -> a -> a
/ (Double
n4dbl forall a. Num a => a -> a -> a
* (Double
d1forall a. Num a => a -> a -> a
-Double
d2) forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow Double
d1 Int
4))
            forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
3)
                         ((Double
2forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
5)forall a. Num a => a -> a -> a
-Double
d1) forall a. Fractional a => a -> a -> a
/ (Double
n4dbl forall a. Num a => a -> a -> a
* (Double
d2forall a. Num a => a -> a -> a
-Double
d1) forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow Double
d2 Int
4))
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
5) (Double
7forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7))
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
5) (Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7))) [Int
2 .. Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
4 Int
np
    let invnp7 :: Double
invnp7 = Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
6) (Double
5forall a. Num a => a -> a -> a
*Double
invnp7)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
6) (Double
3forall a. Num a => a -> a -> a
*Double
invnp7)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
6) Double
invnp7) [Int
3 .. Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
5 (Int
npforall a. Num a => a -> a -> a
*Int
n)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
7) (Double
3forall a. Num a => a -> a -> a
*Double
invnp7)) [Int
1, Int
2, Int
3]
    forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
np forall a. Ord a => a -> a -> Bool
> Int
3) forall a b. (a -> b) -> a -> b
$
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
7) Double
invnp7) [Int
4..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
6 (forall a. Integral a => a -> a -> a
div ((Int
nforall a. Num a => a -> a -> a
-Int
1)forall a. Num a => a -> a -> a
*Int
nforall a. Num a => a -> a -> a
*Int
np) Int
6)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
4) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
7 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
128forall a. Num a => a -> a -> a
*Int
n4forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
+Int
5)))
    let tmp'' :: Double
tmp'' = Int -> Double
toDbl (Int
64forall a. Num a => a -> a -> a
*Int
n6)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
4) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
7 forall a. Fractional a => a -> a -> a
/ Double
tmp'')) [Int
3, Int
4]
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
4) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Double
tmp'')) [Int
5, Int
6, Int
7]
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
key forall a. Eq a => a -> a -> Bool
== Int
4) forall a b. (a -> b) -> a -> b
$ do
    let sg :: Double
sg = Double
1 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
23328forall a. Num a => a -> a -> a
*Int
n6)
        u5 :: Double
u5 = -Double
216 forall a. Num a => a -> a -> a
* Double
sg forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
52212 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*(Int
6353 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
1934forall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
*Int
27)))
        u6 :: Double
u6 = Double
1296 forall a. Num a => a -> a -> a
* Double
sg forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
7884 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*(Int
1541 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
9))
        u7 :: Double
u7 = -Double
7776 forall a. Num a => a -> a -> a
* Double
sg forall a. Num a => a -> a -> a
* Int -> Double
toDbl (Int
8292 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*(Int
1139 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
3))forall a. Fractional a => a -> a -> a
/(Double
ndbl forall a. Num a => a -> a -> a
+ Double
7)
        p0 :: Double
p0 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -Int
144 forall a. Num a => a -> a -> a
* (Int
142528 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
23073 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
115))
        p1 :: Double
p1 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -Int
12 forall a. Num a => a -> a -> a
* (Int
6690556 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
2641189 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
245378 forall a. Num a => a -> a -> a
- Int
nforall a. Num a => a -> a -> a
*Int
1495)))
        p2 :: Double
p2 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -Int
16 forall a. Num a => a -> a -> a
* (Int
6503401 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
4020794forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
787281forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
47323forall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
*Int
385))))
        p3 :: Double
p3 = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ -(Int
6386660 forall a. Num a => a -> a -> a
+ Int
nforall a. Num a => a -> a -> a
*(Int
4411997forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
951821forall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*(Int
61659forall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
*Int
665))))forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
+Int
7)
        a :: Double
a = Double
p2forall a. Fractional a => a -> a -> a
/(Double
3forall a. Num a => a -> a -> a
*Double
p3)
        p :: Double
p = Double
aforall a. Num a => a -> a -> a
*(Double
p1forall a. Fractional a => a -> a -> a
/Double
p2 forall a. Num a => a -> a -> a
- Double
a)
        q :: Double
q = Double
aforall a. Num a => a -> a -> a
*(Double
2forall a. Num a => a -> a -> a
*Double
aforall a. Num a => a -> a -> a
*Double
a forall a. Num a => a -> a -> a
- Double
p1forall a. Fractional a => a -> a -> a
/Double
p3) forall a. Num a => a -> a -> a
+ Double
p0forall a. Fractional a => a -> a -> a
/Double
p3
        th :: Double
th = forall a. Floating a => a -> a
acos(-Double
qforall a. Fractional a => a -> a -> a
/(Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt(-Double
pforall a. Num a => a -> a -> a
*Double
pforall a. Num a => a -> a -> a
*Double
p))) forall a. Fractional a => a -> a -> a
/ Double
3
        r :: Double
r = Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
sqrt(-Double
p)
        tp :: Double
tp = Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
3
        a1 :: Double
a1 = -Double
a forall a. Num a => a -> a -> a
+ Double
rforall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
th
        a3 :: Double
a3 = -Double
a forall a. Num a => a -> a -> a
+ Double
rforall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(Double
thforall a. Num a => a -> a -> a
+Double
tp)
        a2 :: Double
a2 = -Double
3forall a. Num a => a -> a -> a
*Double
a forall a. Num a => a -> a -> a
- Double
a1 forall a. Num a => a -> a -> a
- Double
a3 -- a3 = -a + r*(cos(th+tp))
        npdbl :: Double
npdbl = Int -> Double
toDbl Int
np
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
5) ((Double
1forall a. Num a => a -> a -> a
-Double
ndblforall a. Num a => a -> a -> a
*Double
a1)forall a. Fractional a => a -> a -> a
/Double
npdbl)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
5) ((Double
1forall a. Num a => a -> a -> a
+Double
a1)forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
4) Int
np
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
6) ((Double
1forall a. Num a => a -> a -> a
-Double
ndblforall a. Num a => a -> a -> a
*Double
a2)forall a. Fractional a => a -> a -> a
/Double
npdbl)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
6) ((Double
1forall a. Num a => a -> a -> a
+Double
a2)forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
5) Int
np
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
7) ((Double
1forall a. Num a => a -> a -> a
-Double
ndblforall a. Num a => a -> a -> a
*Double
a3)forall a. Fractional a => a -> a -> a
/Double
npdbl)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
7) ((Double
1forall a. Num a => a -> a -> a
+Double
a3)forall a. Fractional a => a -> a -> a
/Double
npdbl)) [Int
2..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
6) Int
np
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
5, Int
iwforall a. Num a => a -> a -> a
-Int
5)
               ((Double
u7forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
u6forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
*Double
u5)forall a. Fractional a => a -> a -> a
/(Double
a1forall a. Num a => a -> a -> a
*Double
a1forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
a1forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a3) forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a1 Int
5)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
6, Int
iwforall a. Num a => a -> a -> a
-Int
5)
               ((Double
u7forall a. Num a => a -> a -> a
-(Double
a1forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
u6forall a. Num a => a -> a -> a
+Double
a1forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
*Double
u5)forall a. Fractional a => a -> a -> a
/(Double
a2forall a. Num a => a -> a -> a
*Double
a2forall a. Num a => a -> a -> a
-(Double
a1forall a. Num a => a -> a -> a
+Double
a3)forall a. Num a => a -> a -> a
*Double
a2forall a. Num a => a -> a -> a
+Double
a1forall a. Num a => a -> a -> a
*Double
a3) forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a2 Int
5)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
7, Int
iwforall a. Num a => a -> a -> a
-Int
5)
               ((Double
u7forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a1)forall a. Num a => a -> a -> a
*Double
u6forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a1forall a. Num a => a -> a -> a
*Double
u5)forall a. Fractional a => a -> a -> a
/(Double
a3forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
-(Double
a2forall a. Num a => a -> a -> a
+Double
a1)forall a. Num a => a -> a -> a
*Double
a3forall a. Num a => a -> a -> a
+Double
a2forall a. Num a => a -> a -> a
*Double
a1) forall a. Fractional a => a -> a -> a
/ Double -> Int -> Double
pow Double
a3 Int
5)
    let invnp7' :: Double
invnp7' = Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
7)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
8) (Double
4forall a. Num a => a -> a -> a
*Double
invnp7')
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
8) (Double
4forall a. Num a => a -> a -> a
*Double
invnp7')
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
8) Double
invnp7') [Int
3..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
7) Int
o
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
8, Int
iwforall a. Num a => a -> a -> a
-Int
5) (Double
10 forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
729forall a. Num a => a -> a -> a
*Int
n6))
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
gmsforall a. Num a => a -> a -> a
+Int
9) (Double
5.5forall a. Num a => a -> a -> a
*Double
invnp7')
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
gmsforall a. Num a => a -> a -> a
+Int
9) (Double
2.5forall a. Num a => a -> a -> a
*Double
invnp7')
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
gmsforall a. Num a => a -> a -> a
+Int
9) Double
invnp7') [Int
3..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts (Int
gmsforall a. Num a => a -> a -> a
+Int
8) (Int
npforall a. Num a => a -> a -> a
*Int
n)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
gmsforall a. Num a => a -> a -> a
+Int
9, Int
iwforall a. Num a => a -> a -> a
-Int
5) (Double
64 forall a. Num a => a -> a -> a
* Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
6561forall a. Num a => a -> a -> a
*Int
n6))
    let tmp''' :: Double
tmp''' = Int -> Double
toDbl (Int
64forall a. Num a => a -> a -> a
*Int
n6)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
4, Int
iwforall a. Num a => a -> a -> a
-Int
5) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
7 forall a. Fractional a => a -> a -> a
/ Double
tmp''')
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
7, Int
iwforall a. Num a => a -> a -> a
-Int
5) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
6 forall a. Fractional a => a -> a -> a
/ Double
tmp''')
    let invnp9 :: Double
invnp9 = Double
1forall a. Fractional a => a -> a -> a
/(Double
ndblforall a. Num a => a -> a -> a
+Double
9)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
8) (Double
9forall a. Num a => a -> a -> a
*Double
invnp9)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
8) Double
invnp9) [Int
2..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
7 Int
np
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
9) (Double
7forall a. Num a => a -> a -> a
*Double
invnp9)
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
2, Int
9) (Double
3forall a. Num a => a -> a -> a
*Double
invnp9)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
9) Double
invnp9) [Int
3..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
8 (Int
npforall a. Num a => a -> a -> a
*Int
n)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
10) (Double
5forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
1, Int
2]
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
10) Double
invnp9) [Int
3..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
9 Int
o
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
1, Int
11) (Double
5forall a. Num a => a -> a -> a
*Double
invnp9)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
11) (Double
3forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
2, Int
3]
    forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
np forall a. Ord a => a -> a -> Bool
> Int
3) forall a b. (a -> b) -> a -> b
$
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
11) Double
invnp9) [Int
4..Int
np]
    forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
10 (Int
oforall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
-Int
1))
    forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
2, Int
iwforall a. Num a => a -> a -> a
-Int
6) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
3) Int
9 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
1536forall a. Num a => a -> a -> a
*Int
n6))
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
6)
                            (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
5) Int
9 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
512forall a. Num a => a -> a -> a
*Int
n6forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
+Int
7)))) [Int
3, Int
4]
    let tmp'''' :: Double
tmp'''' = Int -> Double
toDbl forall a b. (a -> b) -> a -> b
$ Int
256forall a. Num a => a -> a -> a
*Int
n8
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
6) (-Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
7) Int
9 forall a. Fractional a => a -> a -> a
/ Double
tmp'''')) [Int
5, Int
6, Int
7]
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
i, Int
iwforall a. Num a => a -> a -> a
-Int
6) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
9) Int
8 forall a. Fractional a => a -> a -> a
/ Double
tmp'''')) [Int
8..Int
11]
    forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
n forall a. Ord a => a -> a -> Bool
> Int
2) forall a b. (a -> b) -> a -> b
$ do
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
12) (Double
3forall a. Num a => a -> a -> a
*Double
invnp9)) [Int
1..Int
4]
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
i -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
g (Int
i, Int
12) Double
invnp9) [Int
5..Int
np]
      forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
pts Int
11 (forall a. Integral a => a -> a -> a
div (Int
npforall a. Num a => a -> a -> a
*Int
nforall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
-Int
1)forall a. Num a => a -> a -> a
*(Int
nforall a. Num a => a -> a -> a
-Int
2)) Int
24)
      forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
w (Int
12, Int
iwforall a. Num a => a -> a -> a
-Int
6) (Double -> Int -> Double
pow (Double
ndblforall a. Num a => a -> a -> a
+Double
9) Int
8 forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl (Int
256forall a. Num a => a -> a -> a
*Int
n8))
  Seq IO1dArray
rowsIO <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IOMatrix -> Int -> IO IO1dArray
extractRow IOMatrix
w) (forall a. [a] -> Seq a
S.fromList [Int
2..Int
wts])
  Seq VectorD
rows <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM IO1dArray -> IO VectorD
array1dToVectorD Seq IO1dArray
rowsIO
  let cols :: Seq VectorD
cols = Seq VectorD -> Seq VectorD
transpose Seq VectorD
rows
  Vector Int
pts_out <- forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
UV.unsafeFreeze IOVectorI
pts
  let ptsU :: VectorD
ptsU = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Int -> Double
toDbl Vector Int
pts_out
      row1 :: Seq Double
row1 = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\VectorD
col -> Double
1 forall a. Num a => a -> a -> a
- forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0
                           (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) (forall a. Unbox a => Vector a -> Vector a
UV.tail VectorD
ptsU) VectorD
col))
                  (forall a. Int -> Seq a -> Seq a
S.take Int
rls Seq VectorD
cols)
      wcols :: Seq VectorD
wcols = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Int
j -> forall a. Unbox a => a -> Vector a -> Vector a
UV.cons (forall a. Seq a -> Int -> a
index Seq Double
row1 Int
j) (forall a. Seq a -> Int -> a
index Seq VectorD
cols Int
j))
                   (forall a. [a] -> Seq a
S.fromList [Int
0..(Int
rlsforall a. Num a => a -> a -> a
-Int
1)])
      col1 :: VectorD
col1 = forall a. Seq a -> Int -> a
index Seq VectorD
wcols Int
0
      wcols2 :: Seq VectorD
wcols2 = forall a. a -> Seq a -> Seq a
(S.<|) VectorD
col1
                      (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\VectorD
col -> forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
col VectorD
col1) (forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
wcols))
      col2 :: VectorD
col2 = forall a. Seq a -> Int -> a
index Seq VectorD
wcols2 Int
1
      nb :: Double
nb = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
col1))
      ratio :: Double
ratio = Double
nb forall a. Fractional a => a -> a -> a
/ forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
col2))
      wcol2 :: VectorD
wcol2 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
ratio) VectorD
col2
      wcols3 :: Seq VectorD
wcols3 = forall a. a -> Seq a -> Seq a
(S.<|) VectorD
col1 (forall a. a -> Seq a -> Seq a
(S.<|) VectorD
wcol2 (forall a. Int -> Seq a -> Seq a
S.drop Int
2 Seq VectorD
wcols2))
  let updateW :: Seq VectorD -> Int -> Seq VectorD
      updateW :: Seq VectorD -> Int -> Seq VectorD
updateW Seq VectorD
cols' Int
k = forall a. Int -> a -> Seq a -> Seq a
update (Int
kforall a. Num a => a -> a -> a
-Int
1) VectorD
wknew Seq VectorD
cols'
         where
          ptsW :: VectorD
ptsW = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Fractional a => a -> a -> a
/Double
nb) (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a. Seq a -> Int -> a
index Seq VectorD
cols' (Int
kforall a. Num a => a -> a -> a
-Int
1)))
          slice :: Seq VectorD
slice = forall a. Int -> Seq a -> Seq a
S.drop Int
1 (forall a. Int -> Seq a -> Seq a
S.take (Int
kforall a. Num a => a -> a -> a
-Int
1) Seq VectorD
cols')
          prod1 :: VectorD
prod1 = (forall a. Unbox a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) forall a b. (a -> b) -> a -> b
$ -- ou alors pas de vecteurs
                  forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsW) Seq VectorD
slice
          rows' :: Seq VectorD
rows' = Seq VectorD -> Seq VectorD
transpose Seq VectorD
slice
          prod2 :: VectorD
prod2 = (forall a. Unbox a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) forall a b. (a -> b) -> a -> b
$
                  forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
prod1) Seq VectorD
rows'
          wk :: VectorD
wk = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) (forall a. Seq a -> Int -> a
index Seq VectorD
cols' (Int
kforall a. Num a => a -> a -> a
-Int
1)) VectorD
prod2
          ratio' :: Double
ratio' = Double
nb forall a. Fractional a => a -> a -> a
/
                   forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(*) VectorD
ptsU (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map Double -> Double
square VectorD
wk))
          wknew :: VectorD
wknew = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
ratio') VectorD
wk
      wcolsnew :: Seq VectorD
wcolsnew = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Seq VectorD -> Int -> Seq VectorD
updateW Seq VectorD
wcols3 [Int
3..Int
rls]
  forall (m :: * -> *) a. Monad m => a -> m a
return (IOMatrix
g, Seq VectorD -> Seq VectorD
transpose Seq VectorD
wcolsnew, forall a. Unbox a => Vector a -> [a]
toList forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => (a -> Bool) -> Vector a -> Vector Int
UV.findIndices (forall a. Eq a => a -> a -> Bool
/= Int
0) Vector Int
pts_out)

transpose :: Seq VectorD -> Seq VectorD
transpose :: Seq VectorD -> Seq VectorD
transpose Seq VectorD
cols =
  forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Int
i -> (forall a. Unbox a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> [a]
DF.toList) forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Unbox a => Vector a -> Int -> a
UV.!Int
i) Seq VectorD
cols)
       (forall a. [a] -> Seq a
S.fromList [Int
0..(forall a. Unbox a => Vector a -> Int
UV.length (forall a. Seq a -> Int -> a
index Seq VectorD
cols Int
0) forall a. Num a => a -> a -> a
- Int
1)])

matprod :: IOMatrix -> VectorD -> IO VectorD
matprod :: IOMatrix -> VectorD -> IO VectorD
matprod IOMatrix
mat VectorD
x = do
  ((Int, Int)
_, (Int
m, Int
n)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
mat
  IOVectorD
out <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
UMV.new Int
m :: IO IOVectorD
  let step :: Int -> IO VectorD
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
mforall a. Num a => a -> a -> a
+Int
1 = forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze IOVectorD
out
             | Bool
otherwise = do
              !Double
coef <- Int -> Double -> IO Double
innerstep Int
1 Double
0
              forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
out (Int
iforall a. Num a => a -> a -> a
-Int
1) Double
coef
              Int -> IO VectorD
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
              where
                innerstep :: Int -> Double -> IO Double
                innerstep :: Int -> Double -> IO Double
innerstep Int
j !Double
s | Int
j forall a. Eq a => a -> a -> Bool
== Int
nforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return Double
s
                               | Bool
otherwise = do
                                 Double
mat_ij <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
mat (Int
i, Int
j)
                                 Int -> Double -> IO Double
innerstep (Int
jforall a. Num a => a -> a -> a
+Int
1) (Double
s forall a. Num a => a -> a -> a
+ Double
mat_ij forall a. Num a => a -> a -> a
* (VectorD
x forall a. Unbox a => Vector a -> Int -> a
UV.! (Int
jforall a. Num a => a -> a -> a
-Int
1)))
  Int -> IO VectorD
step Int
1

smpsms :: IOMatrix -> Int -> (VectorD -> VectorD) -> IO1dArray
       -> Double -> IO IO1dArray
smpsms :: IOMatrix
-> Int
-> (VectorD -> VectorD)
-> IO1dArray
-> Double
-> IO IO1dArray
smpsms IOMatrix
vertex Int
nf VectorD -> VectorD
f IO1dArray
g Double
scalar = do
  [Double]
gAsList <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems IO1dArray
g
  [VectorD]
f_gPermuts <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap VectorD -> VectorD
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. IOMatrix -> VectorD -> IO VectorD
matprod IOMatrix
vertex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Unbox a => [a] -> Vector a
fromList)
                     (forall a. (Eq a, Ord a) => [a] -> [[a]]
permuteMultiset [Double]
gAsList)
  forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> [e] -> m (a i e)
newListArray (Int
1, Int
nf)
               (forall a. Unbox a => Vector a -> [a]
toList (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*Double
scalar) (forall a. (a -> a -> a) -> [a] -> a
foldl1' (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) [VectorD]
f_gPermuts)))

extractColumn :: IOMatrix -> Int -> IO IO1dArray
extractColumn :: IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
m Int
j = do
  ((Int, Int)
_, (Int
nrow, Int
_)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
  forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
nrow) (\Int
i -> (Int
i, Int
j)) IOMatrix
m

extractRow :: IOMatrix -> Int -> IO IO1dArray
extractRow :: IOMatrix -> Int -> IO IO1dArray
extractRow IOMatrix
m Int
i = do
  ((Int, Int)
_, (Int
_, Int
ncol)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
  forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
ncol) (\Int
j -> (Int
i, Int
j)) IOMatrix
m

outerProduct :: IO1dArray -> VectorD -> IO IOMatrix
outerProduct :: IO1dArray -> VectorD -> IO IOMatrix
outerProduct IO1dArray
x1 VectorD
x2 = do
  (Int
_, Int
n1) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO1dArray
x1
  let n2 :: Int
n2 = forall a. Unbox a => Vector a -> Int
UV.length VectorD
x2
  IOMatrix
out <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ ((Int
1, Int
1), (Int
n1, Int
n2)) :: IO IOMatrix
  let step :: Int -> IO IOMatrix
      step :: Int -> IO IOMatrix
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
n1forall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return IOMatrix
out
             | Bool
otherwise = do
                Double
x1_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
x1 Int
i
                Double -> Int -> IO IOMatrix
inner Double
x1_i Int
0
              where
                inner :: Double -> Int -> IO IOMatrix
inner !Double
x Int
j | Int
j forall a. Eq a => a -> a -> Bool
== Int
n2 = Int -> IO IOMatrix
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
                           | Bool
otherwise = do
                              forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
out (Int
i, Int
jforall a. Num a => a -> a -> a
+Int
1) (Double
x forall a. Num a => a -> a -> a
* (VectorD
x2 forall a. Unbox a => Vector a -> Int -> a
UV.! Int
j))
                              Double -> Int -> IO IOMatrix
inner Double
x (Int
jforall a. Num a => a -> a -> a
+Int
1)
  Int -> IO IOMatrix
step Int
1

sumMatrices :: [IOMatrix] -> IO IOMatrix
sumMatrices :: [IOMatrix] -> IO IOMatrix
sumMatrices [IOMatrix]
matrices = do
  ((Int, Int)
_, (Int
n1, Int
n2)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds (forall a. [a] -> a
head [IOMatrix]
matrices)
  IOMatrix
out <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ ((Int
1, Int
1), (Int
n1, Int
n2)) :: IO IOMatrix
  let step :: Int -> IO IOMatrix
      step :: Int -> IO IOMatrix
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
n1forall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return IOMatrix
out
             | Bool
otherwise = Int -> IO IOMatrix
inner Int
1
               where
                 inner :: Int -> IO IOMatrix
                 inner :: Int -> IO IOMatrix
inner Int
j | Int
j forall a. Eq a => a -> a -> Bool
== Int
n2forall a. Num a => a -> a -> a
+Int
1 = Int -> IO IOMatrix
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
                         | Bool
otherwise = do
                           [Double]
coefs <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\IOMatrix
m -> forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
m (Int
i, Int
j)) [IOMatrix]
matrices
                           forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
out (Int
i, Int
j) (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
coefs)
                           Int -> IO IOMatrix
inner (Int
jforall a. Num a => a -> a -> a
+Int
1)
  Int -> IO IOMatrix
step Int
1

smprul :: IOMatrix -> Int -> (VectorD -> VectorD) -> Double -> IOMatrix
       -> Seq VectorD -> [Int] -> IO (IO1dArray, IO1dArray)
smprul :: IOMatrix
-> Int
-> (VectorD -> VectorD)
-> Double
-> IOMatrix
-> Seq VectorD
-> [Int]
-> IO (IO1dArray, IO1dArray)
smprul IOMatrix
vrts Int
nf VectorD -> VectorD
f Double
vol IOMatrix
g Seq VectorD
w [Int]
pospts = do
  let rtmn :: Double
rtmn = Double
0.1
      small :: Double
small = Double
1e-12
      errcof :: Double
errcof = Double
8
      rls :: Int
rls = forall a. Unbox a => Vector a -> Int
UV.length (forall a. Seq a -> Int -> a
index Seq VectorD
w Int
0)
  [IOMatrix]
toSum <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\Int
k -> do
                         IO1dArray
g_colk <- IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
g (Int
kforall a. Num a => a -> a -> a
+Int
1)
                         IO1dArray
sms <- IOMatrix
-> Int
-> (VectorD -> VectorD)
-> IO1dArray
-> Double
-> IO IO1dArray
smpsms IOMatrix
vrts Int
nf VectorD -> VectorD
f IO1dArray
g_colk Double
vol
                         IO1dArray -> VectorD -> IO IOMatrix
outerProduct IO1dArray
sms (forall a. Seq a -> Int -> a
index Seq VectorD
w Int
k))
                [Int]
pospts
  IOMatrix
rule <- [IOMatrix] -> IO IOMatrix
sumMatrices [IOMatrix]
toSum
  IO1dArray
basval <- IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
rule Int
1
  IO1dArray
rgnerr <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray (Int
1, Int
nf) Double
0 :: IO IO1dArray
  let step :: Int -> IO ()
      step :: Int -> IO ()
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
nfforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
             | Bool
otherwise = do
                Double
basval_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
basval Int
i
                let nmbs :: Double
nmbs = forall a. Num a => a -> a
abs Double
basval_i
                (Double
rt, Double
nmcp) <- Int -> Double -> Double -> Double -> IO (Double, Double)
inner Int
rls Double
rtmn Double
0 Double
nmbs
                forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
rt forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Int
rls forall a. Ord a => a -> a -> Bool
> Int
3) forall a b. (a -> b) -> a -> b
$
                  forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (Double
rtforall a. Num a => a -> a -> a
*Double
nmcp)
                Double
rgnerr_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
rgnerr Int
i
                forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (forall a. Ord a => a -> a -> a
max (Double
errcofforall a. Num a => a -> a -> a
*Double
rgnerr_i) (Double
smallforall a. Num a => a -> a -> a
*Double
nmbs))
                Int -> IO ()
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
              where
               inner :: Int -> Double -> Double -> Double -> IO (Double, Double)
               inner :: Int -> Double -> Double -> Double -> IO (Double, Double)
inner Int
k !Double
x !Double
y !Double
z | Int
k forall a. Eq a => a -> a -> Bool
== Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return (Double
x, Double
y)
                                | Bool
otherwise = do
                                  Double
rule_ik <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
rule (Int
i, Int
k)
                                  Double
rule_ikm1 <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
rule (Int
i, Int
kforall a. Num a => a -> a -> a
-Int
1)
                                  let nmrl :: Double
nmrl = forall a. Ord a => a -> a -> a
max (forall a. Num a => a -> a
abs Double
rule_ik) (forall a. Num a => a -> a
abs Double
rule_ikm1)
                                  Double
rgnerr_i <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IO1dArray
rgnerr Int
i
                                  forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO1dArray
rgnerr Int
i (forall a. Ord a => a -> a -> a
max Double
nmrl Double
rgnerr_i)
                                  if Double
nmrl forall a. Ord a => a -> a -> Bool
> Double
smallforall a. Num a => a -> a -> a
*Double
z Bool -> Bool -> Bool
&& Int
k forall a. Ord a => a -> a -> Bool
< Int
rls
                                    then Int -> Double -> Double -> Double -> IO (Double, Double)
inner (Int
kforall a. Num a => a -> a -> a
-Int
2) (forall a. Ord a => a -> a -> a
max (Double
nmrlforall a. Fractional a => a -> a -> a
/Double
y) Double
x) Double
nmrl Double
z
                                    else Int -> Double -> Double -> Double -> IO (Double, Double)
inner (Int
kforall a. Num a => a -> a -> a
-Int
2) Double
x Double
nmrl Double
z
  Int -> IO ()
step Int
1
  forall (m :: * -> *) a. Monad m => a -> m a
return (IO1dArray
basval, IO1dArray
rgnerr)

rowMeans :: IOMatrix -> IO VectorD
rowMeans :: IOMatrix -> IO VectorD
rowMeans IOMatrix
m = do
  ((Int, Int)
_, (Int
nrow, Int
ncol)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IOMatrix
m
  IOVectorD
outIO <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
new Int
nrow :: IO IOVectorD
  let step :: Int -> IO ()
      step :: Int -> IO ()
step Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
nrowforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
             | Bool
otherwise = do
                !Double
sum_i <- Int -> Double -> IO Double
inner Int
1 Double
0
                forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
outIO (Int
iforall a. Num a => a -> a -> a
-Int
1) Double
sum_i
                Int -> IO ()
step (Int
iforall a. Num a => a -> a -> a
+Int
1)
              where
                inner :: Int -> Double -> IO Double
                inner :: Int -> Double -> IO Double
inner Int
j !Double
x | Int
j forall a. Eq a => a -> a -> Bool
== Int
ncolforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return (Double
x forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl Int
ncol)
                           | Bool
otherwise = do
                             Double
coef <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
m (Int
i, Int
j)
                             Int -> Double -> IO Double
inner (Int
jforall a. Num a => a -> a -> a
+Int
1) (Double
x forall a. Num a => a -> a -> a
+ Double
coef)
  Int -> IO ()
step Int
1
  forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
unsafeFreeze IOVectorD
outIO

array1dToVectorD :: IO1dArray -> IO VectorD
array1dToVectorD :: IO1dArray -> IO VectorD
array1dToVectorD IO1dArray
array = forall (m :: * -> *) a b. Monad m => (a -> b) -> m a -> m b
(<$!>) forall a. Unbox a => [a] -> Vector a
fromList (forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems IO1dArray
array)

getVectors :: Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD) -- 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 <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
n) (\Int
i -> (Int
i, Int
j1, Int
k)) IO3dArray
m)
  VectorD
out2U <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
n) (\Int
i -> (Int
i, Int
j2, Int
k)) IO3dArray
m)
  -- 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
  forall (m :: * -> *) a. Monad m => a -> m a
return (VectorD
out1U, VectorD
out2U)


smpdfs :: Int -> (VectorD -> VectorD) -> Int -> Int
       -> IO3dArray -> IO (Int, IO3dArray)
smpdfs :: Int
-> (VectorD -> VectorD)
-> Int
-> Int
-> IO3dArray
-> IO (Int, IO3dArray)
smpdfs Int
nd VectorD -> VectorD
f Int
top Int
sbs IO3dArray
vrts = do
  let cuttf :: Double
cuttf = Double
2.0
      cuttb :: Double
cuttb = Double
8.0
  IOMatrix
v <- forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1)) (\(Int
i, Int
j) -> (Int
i, Int
j, Int
top)) IO3dArray
vrts
  VectorD
cn <- IOMatrix -> IO VectorD
rowMeans IOMatrix
v
  let fc :: VectorD
fc = VectorD -> VectorD
f VectorD
cn
      dfmd :: Double
dfmd = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map forall a. Num a => a -> a
abs VectorD
fc)
  IOMatrix
frthdf <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1)) Double
0 :: IO IOMatrix
  IOVectorI
iejeitjtisjsls <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
new Int
7 :: IO IOVectorI
  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
1
  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
2
  IOVectorD
dfmxdfnx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
UMV.replicate Int
2 Double
0 :: IO IOVectorD
  let step :: Int -> Double -> IO ()
      step :: Int -> Double -> IO ()
step Int
i Double
x | Int
i forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
               | Bool
otherwise = do
                  !Double
emx <- Int -> Double -> IO Double
inner (Int
iforall a. Num a => a -> a -> a
+Int
1) Double
x
                  Int -> Double -> IO ()
step (Int
iforall a. Num a => a -> a -> a
+Int
1) Double
emx
                where
                  inner :: Int -> Double -> IO Double
                  inner :: Int -> Double -> IO Double
inner Int
j !Double
y | Int
j forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
2 = forall (m :: * -> *) a. Monad m => a -> m a
return Double
y
                             | Bool
otherwise = do
                              VectorD
vi <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
i)
                              VectorD
vj <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
j)
                              let h :: VectorD
h = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*(Double
2forall a. Fractional a => a -> a -> a
/(Double
5forall a. Num a => a -> a -> a
*(Int -> Double
toDbl Int
nd forall a. Num a => a -> a -> a
+Double
1))))
                                             (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
vi VectorD
vj)
                                  ewd :: Double
ewd = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map forall a. Num a => a -> a
abs VectorD
h)
                                  twoh :: VectorD
twoh = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*Double
2) VectorD
h
                                  t1 :: VectorD
t1 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
cn VectorD
twoh)
                                  t2 :: VectorD
t2 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
cn VectorD
twoh)
                                  t3 :: VectorD
t3 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (forall a. Num a => a -> a -> a
*Double
6) VectorD
fc
                                  t4 :: VectorD
t4 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
cn VectorD
h)
                                  t5 :: VectorD
t5 = VectorD -> VectorD
f (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
cn VectorD
h)
                                  t6 :: VectorD
t6 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*(-Double
4))forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
t4 VectorD
t5
                                  tsum :: VectorD
tsum = forall a. (a -> a -> a) -> [a] -> a
foldl1' (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) [VectorD
t1, VectorD
t2, VectorD
t3, VectorD
t6]
                                  dfr1 :: Double
dfr1 = forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
UV.foldr forall a. Num a => a -> a -> a
(+) Double
0 (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map forall a. Num a => a -> a
abs VectorD
tsum)
                                  dfr2 :: Double
dfr2 = if Double
dfmdforall a. Num a => a -> a -> a
+Double
dfr1forall a. Fractional a => a -> a -> a
/Double
8 forall a. Eq a => a -> a -> Bool
== Double
dfmd then Double
0 else Double
dfr1
                                  dfr3 :: Double
dfr3 = Double
dfr2forall a. Num a => a -> a -> a
*Double
ewd
                              Double
dfmx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
0
                              if Double
dfr3 forall a. Ord a => a -> a -> Bool
>= Double
dfmx
                                then do
                                 Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
                                 Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
                                 forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
is
                                 forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
3 Int
js
                                 forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
i
                                 forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
j
                                 forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
1 Double
dfmx
                                 forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
0 Double
dfr3
                                else do
                                 Double
dfnx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
1
                                 forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
dfr3 forall a. Ord a => a -> a -> Bool
>= Double
dfnx) forall a b. (a -> b) -> a -> b
$ do
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
i
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
3 Int
j
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
1 Double
dfr3
                              forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IOMatrix
frthdf (Int
i, Int
j) Double
dfr3
                              if Double
ewd forall a. Ord a => a -> a -> Bool
>= Double
y
                                then do
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
0 Int
i
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
1 Int
j
                                  Int -> Double -> IO Double
inner (Int
jforall a. Num a => a -> a -> a
+Int
1) Double
ewd
                                else Int -> Double -> IO Double
inner (Int
jforall a. Num a => a -> a -> a
+Int
1) Double
y
  Int -> Double -> IO ()
step Int
1 Double
0
  Double
dfmx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
0
  Double
dfnx <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorD
dfmxdfnx Int
1
  let nregions :: Int
nregions = if Double
dfnx forall a. Ord a => a -> a -> Bool
> Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttf then Int
4 else Int
3
  if Double
dfnx forall a. Ord a => a -> a -> Bool
> Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttf
    then forall (m :: * -> *) a. Monad m => a -> m a
return ()
    else
      if Double
dfmx forall a. Eq a => a -> a -> Bool
== Double
0
        then do
          Int
ie <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
0
          Int
je <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
1
          forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
ie
          forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
je
        else do
          let loop :: Int -> Double -> Int -> IO Int
              loop :: Int -> Double -> Int -> IO Int
loop Int
l !Double
x !Int
ls | Int
l forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
2 = forall (m :: * -> *) a. Monad m => a -> m a
return Int
ls
                            | Bool
otherwise = do
                              Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
                              Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
                              if (Int
l forall a. Eq a => a -> a -> Bool
/= Int
is) Bool -> Bool -> Bool
&& (Int
l forall a. Eq a => a -> a -> Bool
/= Int
js)
                                then do
                                  let it :: Int
it = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Int
l, Int
is, Int
js]
                                      jt :: Int
jt = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Int
l, Int
is, Int
js]
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
it
                                  forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
3 Int
jt
                                  let lt :: Int
lt = Int
isforall a. Num a => a -> a -> a
+Int
jsforall a. Num a => a -> a -> a
+Int
lforall a. Num a => a -> a -> a
-Int
itforall a. Num a => a -> a -> a
-Int
jt
                                  Double
dfr1 <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (Int
it, Int
lt)
                                  Double
dfr2 <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (Int
lt, Int
jt)
                                  let dfr :: Double
dfr = Double
dfr1 forall a. Num a => a -> a -> a
+ Double
dfr2
                                  if Double
dfr forall a. Ord a => a -> a -> Bool
>= Double
x
                                    then Int -> Double -> Int -> IO Int
loop (Int
lforall a. Num a => a -> a -> a
+Int
1) Double
dfr Int
l
                                    else Int -> Double -> Int -> IO Int
loop (Int
lforall a. Num a => a -> a -> a
+Int
1) Double
x Int
ls
                                else Int -> Double -> Int -> IO Int
loop (Int
lforall a. Num a => a -> a -> a
+Int
1) Double
x Int
ls
          !Int
ls <- Int -> Double -> Int -> IO Int
loop Int
1 Double
0 Int
0
          forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
6 Int
ls
          Int
is <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
          Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
          Double
difil <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (forall a. Ord a => a -> a -> a
min Int
is Int
ls, forall a. Ord a => a -> a -> a
max Int
is Int
ls)
          Double
diflj <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray IOMatrix
frthdf (forall a. Ord a => a -> a -> a
min Int
js Int
ls, forall a. Ord a => a -> a -> a
max Int
js Int
ls)
          let dfnx' :: Double
dfnx' = forall a. Ord a => a -> a -> a
max Double
difil Double
diflj
          forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorD
dfmxdfnx Int
1 Double
dfnx'
          forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttb forall a. Ord a => a -> a -> Bool
< Double
dfnx' Bool -> Bool -> Bool
&& Double
difil forall a. Ord a => a -> a -> Bool
> Double
diflj) forall a b. (a -> b) -> a -> b
$ do
            Int
is' <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
            Int
js' <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
            Int
it <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
2
            forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
2 Int
is'
            forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
4 Int
js'
            forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
unsafeWrite IOVectorI
iejeitjtisjsls Int
5 Int
it
  IO3dArray
vrts2 <- forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1, Int
sbsforall a. Num a => a -> a -> a
+Int
nregionsforall a. Num a => a -> a -> a
-Int
1))
                      (\(Int
i, Int
j, Int
k) -> (Int
i, Int
j, if Int
k forall a. Ord a => a -> a -> Bool
<= Int
sbs then Int
k else Int
top)) IO3dArray
vrts
  -- 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 <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
4
  Int
js <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
5
  VectorD
vti <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
is)
  VectorD
vtj <- forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) IO1dArray -> IO VectorD
array1dToVectorD (IOMatrix -> Int -> IO IO1dArray
extractColumn IOMatrix
v Int
js)
  if Int
nregions forall a. Eq a => a -> a -> Bool
== Int
4
    then do
      let vt :: VectorD
vt = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vti VectorD
vtj
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
top) (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vt
      IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vti
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
2) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vt
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vtj
      Int
it <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
2
      Int
jt <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
3
      (VectorD
vtit, VectorD
vtjt) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 Int
top Int
it Int
jt
      -- 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 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vtit VectorD
vtjt
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
jt, Int
top) (Int
it, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vtt
      IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
jt, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vtjt
      (VectorD
vti2, VectorD
vtj2) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 (Int
sbsforall a. Num a => a -> a -> a
+Int
2) Int
it Int
jt
      -- 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 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vti2 VectorD
vtj2
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
jt, Int
sbsforall a. Num a => a -> a -> a
+Int
2) (Int
it, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vt2
      IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
jt, Int
sbsforall a. Num a => a -> a -> a
+Int
3) VectorD
vtj2
    else do
      let vt :: VectorD
vt = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (\Double
x Double
y -> (Double
2forall a. Num a => a -> a -> a
*Double
xforall a. Num a => a -> a -> a
+Double
y)forall a. Fractional a => a -> a -> a
/Double
3) VectorD
vti VectorD
vtj
      IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
top) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
1) VectorD
vt
      if Double
dfmxforall a. Fractional a => a -> a -> a
/Double
cuttf forall a. Ord a => a -> a -> Bool
< Double
dfnx
        then do
          IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
1) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
          IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
          Int
ls <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
unsafeRead IOVectorI
iejeitjtisjsls Int
6
          (VectorD
vtj1, VectorD
vtl1) <- Int -> IO3dArray -> Int -> Int -> Int -> IO (VectorD, VectorD)
getVectors Int
nd IO3dArray
vrts2 (Int
sbsforall a. Num a => a -> a -> a
+Int
1) Int
js Int
ls
          -- 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 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (((forall a. Num a => a -> a -> a
*Double
0.5)forall b c a. (b -> c) -> (a -> b) -> a -> c
.)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a -> a
(+)) VectorD
vtj1 VectorD
vtl1
          IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
ls, Int
sbsforall a. Num a => a -> a -> a
+Int
1) (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vt1
          IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
ls, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtl1
        else do
          let vv :: VectorD
vv = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (\Double
x Double
y -> (Double
xforall a. Num a => a -> a -> a
+Double
2forall a. Num a => a -> a -> a
*Double
y)forall a. Fractional a => a -> a -> a
/Double
3) VectorD
vti VectorD
vtj
          IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
1) (Int
is, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vv
          IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
vrts2 (Int
js, Int
sbsforall a. Num a => a -> a -> a
+Int
2) VectorD
vtj
  forall (m :: * -> *) a. Monad m => a -> m a
return (Int
nregions, IO3dArray
vrts2)
  where
    replaceDimension :: IO3dArray -> (Int, Int) -> VectorD -> IO ()
    replaceDimension :: IO3dArray -> (Int, Int) -> VectorD -> IO ()
replaceDimension IO3dArray
m (Int
j, Int
k) VectorD
v = do
      let loop :: Int -> IO ()
          loop :: Int -> IO ()
loop Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
                 | Bool
otherwise = do
                   forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO3dArray
m (Int
i, Int
j, Int
k) (forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iforall a. Num a => a -> a -> a
-Int
1))
                   Int -> IO ()
loop (Int
iforall a. Num a => a -> a -> a
+Int
1)
      Int -> IO ()
loop Int
1
    replaceDimensions :: IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
    replaceDimensions :: IO3dArray -> (Int, Int) -> (Int, Int) -> VectorD -> IO ()
replaceDimensions IO3dArray
m (Int
j1, Int
k1) (Int
j2, Int
k2) VectorD
v = do
      let loop :: Int -> IO ()
          loop :: Int -> IO ()
loop Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
ndforall a. Num a => a -> a -> a
+Int
1 = forall (m :: * -> *) a. Monad m => a -> m a
return ()
                 | Bool
otherwise = do
                   forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO3dArray
m (Int
i, Int
j1, Int
k1) (forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iforall a. Num a => a -> a -> a
-Int
1))
                   forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray IO3dArray
m (Int
i, Int
j2, Int
k2) (forall a. Unbox a => Vector a -> Int -> a
(UV.!) VectorD
v (Int
iforall a. Num a => a -> a -> a
-Int
1))
                   Int -> IO ()
loop (Int
iforall a. Num a => a -> a -> a
+Int
1)
      Int -> IO ()
loop Int
1

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

-- | 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 forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
er forall a. Ord a => a -> a -> Bool
< Double
0
    = Int
5
  | Int
nf forall a. Ord a => a -> a -> Bool
< Int
1
    = Int
4
  | Int
nd forall a. Ord a => a -> a -> Bool
< Int
2
    = Int
3
  | Int
key forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Int
key forall a. Ord a => a -> a -> Bool
> Int
4
    = Int
2
  | Int
mxfs forall a. Ord a => a -> a -> Bool
< Int
sbs forall a. Num a => a -> a -> a
* Int -> Int -> Int
smpchc Int
nd Int
key
    = Int
1
  | Bool
otherwise
    = Int
0

adsimp :: Int -> Int -> Int -> (VectorD -> VectorD) -> Double -> Double
       -> Int -> IO3dArray -> IO (VectorD, VectorD, Int, Bool)
adsimp :: Int
-> Int
-> Int
-> (VectorD -> VectorD)
-> Double
-> Double
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
adsimp Int
nd Int
nf Int
mxfs VectorD -> VectorD
f Double
ea Double
er Int
key IO3dArray
vrts = do
  ((Int, Int, Int)
_, (Int
_, Int
_, Int
sbs)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO3dArray
vrts
  case Int -> Int -> Int -> Double -> Double -> Int -> Int -> Int
check Int
nd Int
nf Int
mxfs Double
ea Double
er Int
sbs Int
key of
    Int
0 -> Int
-> Int
-> (VectorD -> VectorD)
-> Int
-> Double
-> Double
-> Int
-> Int
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
smpsad Int
nd Int
nf VectorD -> VectorD
f Int
mxfs Double
ea Double
er Int
key (Int -> Int -> Int
smpchc Int
nd Int
key) Int
sbs IO3dArray
vrts
    Int
1 -> Int
-> Int
-> (VectorD -> VectorD)
-> Int
-> Double
-> Double
-> Int
-> Int
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
smpsad Int
nd Int
nf VectorD -> VectorD
f (Int
sbs forall a. Num a => a -> a -> a
* Int -> Int -> Int
smpchc Int
nd Int
key) Double
ea Double
er Int
key (Int -> Int -> Int
smpchc Int
nd Int
key) Int
sbs IO3dArray
vrts
    Int
2 -> forall a. HasCallStack => [Char] -> a
error [Char]
"integration rule must be between 1 and 4"
    Int
3 -> forall a. HasCallStack => [Char] -> a
error [Char]
"dimension must be at least 2"
    Int
4 -> forall a. HasCallStack => [Char] -> a
error [Char]
"number of components must be at least 1"
    Int
5 -> forall a. HasCallStack => [Char] -> a
error [Char]
"requested errors must be positive"
    Int
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"this should not happen"

type Params = (Bool, Int, Int, Seq VectorD, Seq VectorD, 
               VectorD, VectorD, IO3dArray, Seq Double)

smpsad :: Int -> Int -> (VectorD -> VectorD) -> Int -> Double -> Double -> Int
       -> Int -> Int -> IO3dArray -> IO (VectorD, VectorD, Int, Bool)
smpsad :: Int
-> Int
-> (VectorD -> VectorD)
-> Int
-> Double
-> Double
-> Int
-> Int
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
smpsad Int
nd Int
nf VectorD -> VectorD
f Int
mxfs Double
ea Double
er Int
key Int
rcls Int
sbs IO3dArray
vrts = do
  let dfcost :: Int
dfcost = Int
1 forall a. Num a => a -> a -> a
+ Int
2forall a. Num a => a -> a -> a
*Int
ndforall a. Num a => a -> a -> a
*(Int
ndforall a. Num a => a -> a -> a
+Int
1)
  (IOMatrix
g, Seq VectorD
w, [Int]
pospts) <- Int -> Int -> IO (IOMatrix, Seq VectorD, [Int])
smprms Int
nd Int
key
  [Simplex]
simplices <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO3dArray -> Int -> Int -> IO Simplex
toSimplex IO3dArray
vrts (Int
ndforall a. Num a => a -> a -> a
+Int
1)) [Int
1..Int
sbs]
  let vol :: Seq Double
vol = forall a. [a] -> Seq a
S.fromList forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map Simplex -> Double
simplexVolume [Simplex]
simplices
      nv :: Int
nv = Int
sbsforall a. Num a => a -> a -> a
*Int
rcls
  Seq IOMatrix
matrices <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM
              (\Int
k -> forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1)) (\(Int
i, Int
j) -> (Int
i, Int
j, Int
k)) IO3dArray
vrts)
              (forall a. [a] -> Seq a
S.fromList [Int
1..Int
sbs])
  Seq (IO1dArray, IO1dArray)
br <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\(IOMatrix
m, Double
v) -> IOMatrix
-> Int
-> (VectorD -> VectorD)
-> Double
-> IOMatrix
-> Seq VectorD
-> [Int]
-> IO (IO1dArray, IO1dArray)
smprul IOMatrix
m Int
nf VectorD -> VectorD
f Double
v IOMatrix
g Seq VectorD
w [Int]
pospts) (forall a b. Seq a -> Seq b -> Seq (a, b)
S.zip Seq IOMatrix
matrices Seq Double
vol)
  Seq VectorD
aes <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> b
snd) Seq (IO1dArray, IO1dArray)
br
  Seq VectorD
vls <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> a
fst) Seq (IO1dArray, IO1dArray)
br
  let vl :: VectorD
vl = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
vls
      ae :: VectorD
ae = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
aes
      fl :: Bool
fl = VectorD -> VectorD -> Bool
getFL VectorD
ae VectorD
vl
  let loop :: Params -> IO (VectorD, VectorD, Int, Bool)
      loop :: Params -> IO (VectorD, VectorD, Int, Bool)
loop !Params
params | Bool -> Bool
not Bool
fl' Bool -> Bool -> Bool
|| Int
nv'forall a. Num a => a -> a -> a
+Int
dfcostforall a. Num a => a -> a -> a
+Int
4forall a. Num a => a -> a -> a
*Int
rcls forall a. Ord a => a -> a -> Bool
> Int
mxfs =
                     forall (m :: * -> *) a. Monad m => a -> m a
return (VectorD
vl', VectorD
ae', Int
nv', Bool
fl')
                   | Bool
otherwise = do
                      let maxs :: Seq Double
maxs = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. (Unbox a, Ord a) => Vector a -> a
UV.maximum Seq VectorD
aes'
                          imax :: Int
imax = forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> Seq a -> Maybe Int
S.findIndexL (forall a. Eq a => a -> a -> Bool
== forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum Seq Double
maxs) Seq Double
maxs
                          vl0 :: VectorD
vl0 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
vl' (forall a. Seq a -> Int -> a
index Seq VectorD
vls' Int
imax)
                          ae0 :: VectorD
ae0 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith (-) VectorD
ae' (forall a. Seq a -> Int -> a
index Seq VectorD
aes' Int
imax)
                      (Int
nregions, IO3dArray
vrts2) <- Int
-> (VectorD -> VectorD)
-> Int
-> Int
-> IO3dArray
-> IO (Int, IO3dArray)
smpdfs Int
nd VectorD -> VectorD
f (Int
imaxforall a. Num a => a -> a -> a
+Int
1) Int
sbs' IO3dArray
vrts'
                      let vi :: Double
vi = forall a. Seq a -> Int -> a
index Seq Double
vol' Int
imax forall a. Fractional a => a -> a -> a
/ Int -> Double
toDbl Int
nregions
                          nv2 :: Int
nv2 = Int
nv' forall a. Num a => a -> a -> a
+ Int
nregionsforall a. Num a => a -> a -> a
*Int
rcls forall a. Num a => a -> a -> a
+ Int
dfcost -- nregions*rcls ?
                          sbs2 :: Int
sbs2 = Int
sbs' forall a. Num a => a -> a -> a
+ Int
nregionsforall a. Num a => a -> a -> a
-Int
1
                      Seq IOMatrix
matrices2 <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM
                                   (\Int
k -> forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices ((Int
1, Int
1), (Int
nd, Int
ndforall a. Num a => a -> a -> a
+Int
1))
                                          (\(Int
i, Int
j) -> (Int
i, Int
j, Int
k)) IO3dArray
vrts2)
                                   (forall a. [a] -> Seq a
S.fromList ((Int
imaxforall a. Num a => a -> a -> a
+Int
1)forall a. a -> [a] -> [a]
:[(Int
sbs'forall a. Num a => a -> a -> a
+Int
1)..Int
sbs2]))
                      Seq (IO1dArray, IO1dArray)
br2 <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\IOMatrix
m -> IOMatrix
-> Int
-> (VectorD -> VectorD)
-> Double
-> IOMatrix
-> Seq VectorD
-> [Int]
-> IO (IO1dArray, IO1dArray)
smprul IOMatrix
m Int
nf VectorD -> VectorD
f Double
vi IOMatrix
g Seq VectorD
w [Int]
pospts) Seq IOMatrix
matrices2
                      Seq VectorD
rgnerrs <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> b
snd) Seq (IO1dArray, IO1dArray)
br2
                      Seq VectorD
basvals <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO1dArray -> IO VectorD
array1dToVectorDforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a b. (a, b) -> a
fst) Seq (IO1dArray, IO1dArray)
br2
                      let vl2 :: VectorD
vl2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
vl0
                                (forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
basvals)
                          ae2 :: VectorD
ae2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+) VectorD
ae0
                                (forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
UV.zipWith forall a. Num a => a -> a -> a
(+)) Seq VectorD
rgnerrs)
                          aes2 :: Seq VectorD
aes2 = forall a. Int -> a -> Seq a -> Seq a
update Int
imax (forall a. Seq a -> Int -> a
index Seq VectorD
rgnerrs Int
0) Seq VectorD
aes'
                                 forall a. Seq a -> Seq a -> Seq a
>< forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
rgnerrs
                          vls2 :: Seq VectorD
vls2 = forall a. Int -> a -> Seq a -> Seq a
update Int
imax (forall a. Seq a -> Int -> a
index Seq VectorD
basvals Int
0) Seq VectorD
vls'
                                 forall a. Seq a -> Seq a -> Seq a
>< forall a. Int -> Seq a -> Seq a
S.drop Int
1 Seq VectorD
basvals
                          fl2 :: Bool
fl2 = VectorD -> VectorD -> Bool
getFL VectorD
ae2 VectorD
vl2
                          vol2 :: Seq Double
vol2 = forall a. Int -> a -> Seq a -> Seq a
update Int
imax Double
vi Seq Double
vol'
                                 forall a. Seq a -> Seq a -> Seq a
>< forall a. Int -> a -> Seq a
S.replicate (Int
nregionsforall a. Num a => a -> a -> a
-Int
1) Double
vi
                      Params -> IO (VectorD, VectorD, Int, Bool)
loop (Bool
fl2, Int
nv2, Int
sbs2, Seq VectorD
aes2, Seq VectorD
vls2, VectorD
ae2, VectorD
vl2, IO3dArray
vrts2, Seq Double
vol2)
                    where
                      (Bool
fl', Int
nv', Int
sbs', Seq VectorD
aes', Seq VectorD
vls', VectorD
ae', VectorD
vl', IO3dArray
vrts', Seq Double
vol') = Params
params
  -- 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 = forall a. Unbox a => (a -> Bool) -> Vector a -> Bool
UV.any (forall a. Ord a => a -> a -> Bool
> forall a. Ord a => a -> a -> a
max Double
ea (forall a. (Unbox a, Ord a) => Vector a -> a
UV.maximum (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map ((forall a. Num a => a -> a -> a
*Double
er)forall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a
abs) VectorD
v))) VectorD
a


toSimplex :: IO3dArray -> Int -> Int -> IO Simplex
toSimplex :: IO3dArray -> Int -> Int -> IO Simplex
toSimplex IO3dArray
m Int
n Int
k = do
  ((Int, Int, Int)
_, (Int
nrow, Int
_, Int
_)) <- forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds IO3dArray
m
  let getColumn :: Int -> IO IO1dArray
      getColumn :: Int -> IO IO1dArray
getColumn Int
col = forall (a :: * -> * -> *) e (m :: * -> *) i j.
(MArray a e m, Ix i, Ix j) =>
(i, i) -> (i -> j) -> a j e -> m (a i e)
mapIndices (Int
1, Int
nrow) (\Int
i -> (Int
i, Int
col, Int
k)) IO3dArray
m
  [IO1dArray]
columns <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM Int -> IO IO1dArray
getColumn [Int
1..Int
n]
  forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m [e]
getElems [IO1dArray]
columns