-- | Clarence Barlow. \"Two Essays on Theory\".
-- /Computer Music Journal/, 11(1):44-60, 1987.
-- Translated by Henning Lohner.
module Music.Theory.Meter.Barlow_1987 where

import Data.List {- base -}
--import Debug.Trace

import qualified Data.Numbers.Primes as P {- primes -}

import qualified Music.Theory.Math as T {- hmt-base -}

traceShow :: a -> b -> b
traceShow :: forall a b. a -> b -> b
traceShow a
_ b
x = b
x

-- | One indexed variant of 'genericIndex'.
--
-- > map (at1 [11..13]) [1..3] == [11,12,13]
at1 :: Integral n => [a] -> n -> a
at1 :: forall n a. Integral n => [a] -> n -> a
at1 [a]
x n
i = [a]
x forall n a. Integral n => [a] -> n -> a
`genericIndex` (n
i forall a. Num a => a -> a -> a
- n
1)

-- | Variant of 'at1' with boundary rules and specified error message.
--
-- > map (at1_bnd_err 'x' [11..13]) [0..4] == [1,11,12,13,1]
-- > at1_bnd_err 'x' [0] 3 == undefined
at1_bnd_err :: (Num a,Show a,Integral n,Show n,Show m) => m -> [a] -> n -> a
at1_bnd_err :: forall a n m.
(Num a, Show a, Integral n, Show n, Show m) =>
m -> [a] -> n -> a
at1_bnd_err m
m [a]
x n
i =
    let n :: n
n = forall i a. Num i => [a] -> i
genericLength [a]
x
    in if n
i forall a. Eq a => a -> a -> Bool
== n
0 Bool -> Bool -> Bool
|| n
i forall a. Eq a => a -> a -> Bool
== n
n forall a. Num a => a -> a -> a
+ n
1
       then a
1 -- error (show ("at':==",m,x,i))
       else if n
i forall a. Ord a => a -> a -> Bool
< n
0 Bool -> Bool -> Bool
|| n
i forall a. Ord a => a -> a -> Bool
> n
n forall a. Num a => a -> a -> a
+ n
1
            then forall a. HasCallStack => String -> a
error (forall a. Show a => a -> String
show (String
"at1_bnd_err",m
m,[a]
x,n
i))
            else [a]
x forall n a. Integral n => [a] -> n -> a
`genericIndex` (n
i forall a. Num a => a -> a -> a
- n
1)

-- | Variant of 'mod' with input constraints.
--
-- > mod_pos_err (-1) 2 == 1
-- > mod_pos_err 1 (-2) == undefined
mod_pos_err :: (Integral a,Show a) => a -> a -> a
mod_pos_err :: forall a. (Integral a, Show a) => a -> a -> a
mod_pos_err a
a a
b =
    let r :: a
r = forall a. Integral a => a -> a -> a
mod a
a a
b
    in if a
r forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
|| a
r forall a. Ord a => a -> a -> Bool
>= a
b
       then forall a. HasCallStack => String -> a
error (forall a. Show a => a -> String
show (String
"mod_pos_err",a
a,a
b,a
r))
       else a
r

-- | Type-specialised variant of 'fromIntegral'.
to_r :: Integral n => n -> Double
to_r :: forall n. Integral n => n -> Double
to_r = forall a b. (Integral a, Num b) => a -> b
fromIntegral

-- | Variant on 'div' with input constraints.
div_pos_err :: (Integral a,Show a) => String -> a -> a -> a
div_pos_err :: forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
m a
i a
j =
    if a
i forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
|| a
j forall a. Ord a => a -> a -> Bool
< a
0
    then forall a. HasCallStack => String -> a
error (forall a. Show a => a -> String
show (String
"div_pos_err",String
m,a
i,a
j))
    else forall a b. (RealFrac a, Integral b) => a -> b
truncate (forall n. Integral n => n -> Double
to_r a
i forall a. Fractional a => a -> a -> a
/ forall n. Integral n => n -> Double
to_r a
j)

-- | A stratification is a tree of integral subdivisions.
type Stratification t = [t]

-- | Indispensibilities from stratification.
--
-- > indispensibilities [3,2,2] == [11,0,6,3,9,1,7,4,10,2,8,5]
-- > indispensibilities [2,3,2] == [11,0,6,2,8,4,10,1,7,3,9,5]
-- > indispensibilities [2,2,3] == [11,0,4,8,2,6,10,1,5,9,3,7]
-- > indispensibilities [3,5] == [14,0,9,3,6,12,1,10,4,7,13,2,11,5,8]
indispensibilities :: (Integral n,Show n) => Stratification n -> [n]
indispensibilities :: forall n.
(Integral n, Show n) =>
Stratification n -> Stratification n
indispensibilities Stratification n
x = forall a b. (a -> b) -> [a] -> [b]
map (forall a. (Integral a, Show a) => Stratification a -> a -> a -> a
lower_psi Stratification n
x (forall i a. Num i => [a] -> i
genericLength Stratification n
x)) [n
1 .. forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product Stratification n
x]

-- | The indispensibility measure (ψ).
--
-- > map (lower_psi [2] 1) [1..2] == [1,0]
-- > map (lower_psi [3] 1) [1..3] == [2,0,1]
-- > map (lower_psi [2,2] 2) [1..4] == [3,0,2,1]
-- > map (lower_psi [5] 1) [1..5] == [4,0,3,1,2]
-- > map (lower_psi [3,2] 2) [1..6] == [5,0,3,1,4,2]
-- > map (lower_psi [2,3] 2) [1..6] == [5,0,2,4,1,3]
lower_psi :: (Integral a,Show a) => Stratification a -> a -> a -> a
lower_psi :: forall a. (Integral a, Show a) => Stratification a -> a -> a -> a
lower_psi Stratification a
q a
z a
n =
    let s8 :: a -> a
s8 a
r =
            let s1 :: a
s1 = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product Stratification a
q
                s2 :: a
s2 = (a
n forall a. Num a => a -> a -> a
- a
2) forall a. (Integral a, Show a) => a -> a -> a
`mod_pos_err` a
s1
                s3 :: a
s3 = let f :: a -> a
f a
k = forall a n m.
(Num a, Show a, Integral n, Show n, Show m) =>
m -> [a] -> n -> a
at1_bnd_err String
"s3" Stratification a
q (a
z forall a. Num a => a -> a -> a
+ a
1 forall a. Num a => a -> a -> a
- a
k)
                     in forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map a -> a
f [a
0 .. a
r])
                s4 :: a
s4 = a
1 forall a. Num a => a -> a -> a
+ forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
"s4" a
s2 a
s3
                c :: a
c = forall a n m.
(Num a, Show a, Integral n, Show n, Show m) =>
m -> [a] -> n -> a
at1_bnd_err String
"c" Stratification a
q (a
z forall a. Num a => a -> a -> a
- a
r)
                s5 :: a
s5 = a
s4 forall a. (Integral a, Show a) => a -> a -> a
`mod_pos_err` a
c
                s6 :: a
s6 = forall a. (Integral a, Show a) => a -> a -> a
upper_psi a
c (a
1 forall a. Num a => a -> a -> a
+ a
s5)
                s7 :: a
s7 = let f :: a -> a
f = forall a n m.
(Num a, Show a, Integral n, Show n, Show m) =>
m -> [a] -> n -> a
at1_bnd_err String
"s7" Stratification a
q
                     in forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map a -> a
f [a
0 .. a
z forall a. Num a => a -> a -> a
- a
r forall a. Num a => a -> a -> a
- a
1])
            in forall a b. a -> b -> b
traceShow (String
"lower_psi:s",a
s1,a
s2,a
s3,a
s4,a
s5,a
s6,a
s7) (a
s7 forall a. Num a => a -> a -> a
* a
s6)
    in forall a b. a -> b -> b
traceShow (String
"lower_psi",Stratification a
q,a
z,a
n) (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map a -> a
s8 [a
0 .. a
z forall a. Num a => a -> a -> a
- a
1]))

-- | The first /n/ primes, reversed.
--
-- > reverse_primes 14 == [43,41,37,31,29,23,19,17,13,11,7,5,3,2]
-- > length (reverse_primes 14) == 14
reverse_primes :: Integral n => n -> [n]
reverse_primes :: forall n. Integral n => n -> [n]
reverse_primes n
n = forall a. [a] -> [a]
reverse (forall i a. Integral i => i -> [a] -> [a]
genericTake n
n forall int. Integral int => [int]
P.primes)

-- | Generate prime stratification for /n/.
--
-- > map prime_stratification [2,3,5,7,11] == [[2],[3],[5],[7],[11]]
-- > map prime_stratification [6,8,9,12] == [[3,2],[2,2,2],[3,3],[3,2,2]]
-- > map prime_stratification [22,10,4,1] == [[11,2],[5,2],[2,2],[]]
-- > map prime_stratification [18,16,12] == [[3,3,2],[2,2,2,2],[3,2,2]]
prime_stratification :: (Integral n,Show n) => n -> Stratification n
prime_stratification :: forall n. (Integral n, Show n) => n -> Stratification n
prime_stratification =
    let go :: [t] -> t -> [t]
go [t]
x t
k =
            case [t]
x of
              t
p:[t]
x' -> if t
k forall a. Integral a => a -> a -> a
`rem` t
p forall a. Eq a => a -> a -> Bool
== t
0
                      then t
p forall a. a -> [a] -> [a]
: [t] -> t -> [t]
go [t]
x (forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
"ps" t
k t
p)
                      else [t] -> t -> [t]
go [t]
x' t
k
              [] -> []
    in forall {t}. (Integral t, Show t) => [t] -> t -> [t]
go (forall n. Integral n => n -> [n]
reverse_primes n
14)

-- | Fundamental indispensibilities for prime numbers (Ψ).
--
-- > map (upper_psi 2) [1..2] == [1,0]
-- > map (upper_psi 3) [1..3] == [2,0,1]
-- > map (upper_psi 5) [1..5] == [4,0,3,1,2]
-- > map (upper_psi 7) [1..7] == [6,0,4,2,5,1,3]
-- > map (upper_psi 11) [1..11] == [10,0,6,4,9,1,7,3,8,2,5]
-- > map (upper_psi 13) [1..13] == [12,0,7,4,10,1,8,5,11,2,9,3,6]
upper_psi :: (Integral a,Show a) => a -> a -> a
upper_psi :: forall a. (Integral a, Show a) => a -> a -> a
upper_psi a
p a
n =
    if a
p forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` forall n. Integral n => n -> [n]
reverse_primes a
14
    then forall a. HasCallStack => String -> a
error (forall a. Show a => a -> String
show (String
"upper_psi",String
"not prime",a
p,a
n))
    else if a
p forall a. Eq a => a -> a -> Bool
== a
2
         then a
p forall a. Num a => a -> a -> a
- a
n
         else if a
n forall a. Eq a => a -> a -> Bool
== a
p forall a. Num a => a -> a -> a
- a
1
              then forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
"upper_psi" a
p a
4
              else let n' :: a
n' = a
n forall a. Num a => a -> a -> a
- forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
"n'" a
n a
p
                       s :: Stratification a
s = forall n. (Integral n, Show n) => n -> Stratification n
prime_stratification (a
p forall a. Num a => a -> a -> a
- a
1)
                       q :: a
q = forall a. (Integral a, Show a) => Stratification a -> a -> a -> a
lower_psi Stratification a
s (forall i a. Num i => [a] -> i
genericLength Stratification a
s) a
n'
                       q' :: Double
q' = forall n. Integral n => n -> Double
to_r a
q
                       p' :: Double
p' = forall n. Integral n => n -> Double
to_r a
p
                   in forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double
q' forall a. Num a => a -> a -> a
+ Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt ((Double
q' forall a. Num a => a -> a -> a
+ Double
1) forall a. Fractional a => a -> a -> a
/ Double
p'))

-- | Table such that each subsequent row deletes the least
-- indispensibile pulse.
--
-- > thinning_table [3,2] == [[True,True,True,True,True,True]
-- >                         ,[True,False,True,True,True,True]
-- >                         ,[True,False,True,False,True,True]
-- >                         ,[True,False,True,False,True,False]
-- >                         ,[True,False,False,False,True,False]
-- >                         ,[True,False,False,False,False,False]]
thinning_table :: (Integral n,Show n) => Stratification n -> [[Bool]]
thinning_table :: forall n. (Integral n, Show n) => Stratification n -> [[Bool]]
thinning_table Stratification n
s =
    let x :: Stratification n
x = forall n.
(Integral n, Show n) =>
Stratification n -> Stratification n
indispensibilities Stratification n
s
        n :: n
n = forall i a. Num i => [a] -> i
genericLength Stratification n
x
        true :: i -> [Bool]
true i
i = forall i a. Integral i => i -> a -> [a]
genericReplicate i
i Bool
True
        false :: i -> [Bool]
false i
i = forall i a. Integral i => i -> a -> [a]
genericReplicate i
i Bool
False
        f :: n -> [Bool]
f n
i = forall {i}. Integral i => i -> [Bool]
true (n
i forall a. Num a => a -> a -> a
+ n
1)  forall a. [a] -> [a] -> [a]
++ forall {i}. Integral i => i -> [Bool]
false (n
n forall a. Num a => a -> a -> a
- n
i forall a. Num a => a -> a -> a
- n
1)
    in forall a. [[a]] -> [[a]]
transpose (forall a b. (a -> b) -> [a] -> [b]
map n -> [Bool]
f Stratification n
x)

-- | Trivial pretty printer for 'thinning_table'.
--
-- > putStrLn (thinning_table_pp [3,2])
-- > putStrLn (thinning_table_pp [2,3])
--
-- > ******   ******
-- > *.****   *.****
-- > *.*.**   *.**.*
-- > *.*.*.   *..*.*
-- > *...*.   *..*..
-- > *.....   *.....
thinning_table_pp :: (Integral n,Show n) => Stratification n -> String
thinning_table_pp :: forall n. (Integral n, Show n) => Stratification n -> String
thinning_table_pp Stratification n
s =
    let f :: Bool -> Char
f Bool
x = if Bool
x then Char
'*' else Char
'.'
    in [String] -> String
unlines (forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map Bool -> Char
f) (forall n. (Integral n, Show n) => Stratification n -> [[Bool]]
thinning_table Stratification n
s))

-- | Scale values against length of list minus one.
--
-- > relative_to_length [0..5] == [0.0,0.2,0.4,0.6,0.8,1.0]
relative_to_length :: (Real a, Fractional b) => [a] -> [b]
relative_to_length :: forall a b. (Real a, Fractional b) => [a] -> [b]
relative_to_length [a]
x =
    let n :: Int
n = forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
x forall a. Num a => a -> a -> a
- Int
1
    in forall a b. (a -> b) -> [a] -> [b]
map ((forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Real a, Fractional b) => a -> b
realToFrac) [a]
x

-- | Variant of 'indispensibilities' that scales value to lie in
-- @(0,1)@.
--
-- relative_indispensibilities [3,2] == [1,0,0.6,0.2,0.8,0.4]
relative_indispensibilities :: (Integral n,Show n) => Stratification n -> [Double]
relative_indispensibilities :: forall n. (Integral n, Show n) => Stratification n -> [Double]
relative_indispensibilities = forall a b. (Real a, Fractional b) => [a] -> [b]
relative_to_length forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall n.
(Integral n, Show n) =>
Stratification n -> Stratification n
indispensibilities

-- | Align two meters (given as stratifications) to least common
-- multiple of their degrees.  The 'indispensibilities' function is
-- given as an argument so that it may be relative if required.  This
-- generates Table 7 (p.58).
--
-- > let r = [(5,5),(0,0),(2,3),(4,1),(1,4),(3,2)]
-- > in align_meters indispensibilities [2,3] [3,2] == r
--
-- > let r = [(1,1),(0,0),(0.4,0.6),(0.8,0.2),(0.2,0.8),(0.6,0.4)]
-- > in align_meters relative_indispensibilities [2,3] [3,2] == r
--
-- > align_meters indispensibilities [2,2,3] [3,5]
-- > align_meters relative_indispensibilities [2,2,3] [3,5]
align_meters :: (t -> [b]) -> t -> t -> [(b,b)]
align_meters :: forall t b. (t -> [b]) -> t -> t -> [(b, b)]
align_meters t -> [b]
f t
s1 t
s2 =
    let i1 :: [b]
i1 = t -> [b]
f t
s1
        i2 :: [b]
i2 = t -> [b]
f t
s2
        n1 :: Int
n1 = forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
i1
        n2 :: Int
n2 = forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
i2
        n :: Int
n = forall a. Integral a => a -> a -> a
lcm Int
n1 Int
n2
        i1' :: [b]
i1' = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (forall a. Int -> a -> [a]
replicate (Int
n forall a. Integral a => a -> a -> a
`div` Int
n1) [b]
i1)
        i2' :: [b]
i2' = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (forall a. Int -> a -> [a]
replicate (Int
n forall a. Integral a => a -> a -> a
`div` Int
n2) [b]
i2)
    in forall a b. [a] -> [b] -> [(a, b)]
zip [b]
i1' [b]
i2'

-- | Type pairing a stratification and a tempo.
type S_MM t = ([t],t)

-- | Variant of 'div' that requires 'mod_pos_err be @0@.
whole_div :: Integral a => a -> a -> a
whole_div :: forall a. Integral a => a -> a -> a
whole_div a
i a
j =
    case a
i forall a. Integral a => a -> a -> (a, a)
`divMod` a
j of
      (a
k,a
0) -> a
k
      (a, a)
_ -> forall a. HasCallStack => String -> a
error String
"whole_div"

-- | Variant of 'quot' that requires 'rem' be @0@.
whole_quot :: Integral a => a -> a -> a
whole_quot :: forall a. Integral a => a -> a -> a
whole_quot a
i a
j =
    case a
i forall a. Integral a => a -> a -> (a, a)
`quotRem` a
j of
      (a
k,a
0) -> a
k
      (a, a)
_ -> forall a. HasCallStack => String -> a
error String
"whole_quot"

-- | Rule to prolong stratification of two 'S_MM' values such that
-- pulse at the deeper level are aligned.  (Paragraph 2, p.58)
--
-- > let x = ([2,2,2],1)
-- > in prolong_stratifications x x == (fst x,fst x)
--
-- > let r = ([2,5,3,3,2],[3,2,5,5])
-- > in prolong_stratifications ([2,5],50) ([3,2],60) == r
--
-- > prolong_stratifications ([2,2,3],5) ([3,5],4) == ([2,2,3],[3,5])
prolong_stratifications :: (Integral n,Show n) => S_MM n -> S_MM n -> ([n],[n])
prolong_stratifications :: forall n. (Integral n, Show n) => S_MM n -> S_MM n -> ([n], [n])
prolong_stratifications ([n]
s1,n
v1) ([n]
s2,n
v2) =
    let t1 :: n
t1 = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [n]
s1 forall a. Num a => a -> a -> a
* n
v1
        t2 :: n
t2 = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [n]
s2 forall a. Num a => a -> a -> a
* n
v2
        t :: n
t = forall a. Integral a => a -> a -> a
lcm n
t1 n
t2
        s1' :: [n]
s1' = [n]
s1 forall a. [a] -> [a] -> [a]
++ forall n. (Integral n, Show n) => n -> Stratification n
prime_stratification (n
t forall a. Integral a => a -> a -> a
`whole_div` n
t1)
        s2' :: [n]
s2' = [n]
s2 forall a. [a] -> [a] -> [a]
++ forall n. (Integral n, Show n) => n -> Stratification n
prime_stratification (n
t forall a. Integral a => a -> a -> a
`whole_div` n
t2)
    in ([n]
s1',[n]
s2')

-- | Composition of 'prolong_stratifications' and 'align_meters'.
--
-- > align_s_mm indispensibilities ([2,2,3],5) ([3,5],4)
align_s_mm :: (Integral n,Show n) => ([n] -> [t]) -> S_MM n -> S_MM n -> [(t,t)]
align_s_mm :: forall n t.
(Integral n, Show n) =>
([n] -> [t]) -> S_MM n -> S_MM n -> [(t, t)]
align_s_mm [n] -> [t]
f ([n]
s1,n
v1) ([n]
s2,n
v2) =
    let ([n]
s1',[n]
s2') = forall n. (Integral n, Show n) => S_MM n -> S_MM n -> ([n], [n])
prolong_stratifications ([n]
s1,n
v1) ([n]
s2,n
v2)
    in forall t b. (t -> [b]) -> t -> t -> [(b, b)]
align_meters [n] -> [t]
f [n]
s1' [n]
s2'

-- | An attempt at Equation 5 of the /CMJ/ paper.  When /n/ is /h-1/
-- the output is incorrect (it is the product of the correct values
-- for /n/ at /h-1/ and /h/).
--
-- > map (upper_psi' 5) [1..5] /= [4,0,3,1,2]
-- > map (upper_psi' 7) [1..7] /= [6,0,4,2,5,1,3]
-- > map (upper_psi' 11) [1..11] /= [10,0,6,4,9,1,7,3,8,2,5]
-- > map (upper_psi' 13) [1..13] /= [12,0,7,4,10,1,8,5,11,2,9,3,6]
upper_psi' :: (Integral a,Show a) => a -> a -> a
upper_psi' :: forall a. (Integral a, Show a) => a -> a -> a
upper_psi' a
h a
n =
    if a
h forall a. Ord a => a -> a -> Bool
> a
3
    then let omega :: a -> a
omega a
x = if a
x forall a. Eq a => a -> a -> Bool
== a
0 then a
0 else a
1
             h4 :: a
h4 = forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
"h4" a
h a
4
             n' :: a
n' = a
n forall a. Num a => a -> a -> a
- a
1 forall a. Num a => a -> a -> a
+ forall {a} {a}. (Eq a, Num a, Num a) => a -> a
omega (a
h forall a. Num a => a -> a -> a
- a
n)
             p :: Stratification a
p = forall n. (Integral n, Show n) => n -> Stratification n
prime_stratification (a
h forall a. Num a => a -> a -> a
- a
1)
             x0 :: a
x0 = forall a. (Integral a, Show a) => Stratification a -> a -> a -> a
lower_psi Stratification a
p (forall i a. Num i => [a] -> i
genericLength Stratification a
p) a
n'
             x1 :: a
x1 = a
x0 forall a. Num a => a -> a -> a
+ forall {a} {a}. (Eq a, Num a, Num a) => a -> a
omega (forall a. (Integral a, Show a) => String -> a -> a -> a
div_pos_err String
"z" a
x0 a
h4)
             x2 :: a
x2 = forall {a} {a}. (Eq a, Num a, Num a) => a -> a
omega (a
h forall a. Num a => a -> a -> a
- a
n forall a. Num a => a -> a -> a
- a
1)
             x3 :: a
x3 = a
x2 forall a. Num a => a -> a -> a
+ a
h4 forall a. Num a => a -> a -> a
* (a
1 forall a. Num a => a -> a -> a
- a
x2)
         in forall a b. a -> b -> b
traceShow (String
"upper_psi'",a
h,a
n,a
n',a
x0,a
x1,a
x2,a
x3) (a
x1 forall a. Num a => a -> a -> a
* a
x3)
    else (a
h forall a. Num a => a -> a -> a
+ a
n forall a. Num a => a -> a -> a
- a
2) forall a. (Integral a, Show a) => a -> a -> a
`mod_pos_err` a
h

-- | The /MPS/ limit equation given on p.58.
--
-- > mps_limit 3 == 21 + 7/9
mps_limit :: Floating a => a -> a
mps_limit :: forall a. Floating a => a -> a
mps_limit a
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a
n forall a. Floating a => a -> a -> a
** a
4 forall a. Fractional a => a -> a -> a
/ a
9
                  ,a
n forall a. Floating a => a -> a -> a
** a
3 forall a. Fractional a => a -> a -> a
/ a
3
                  ,a
13 forall a. Num a => a -> a -> a
* (a
n forall a. Floating a => a -> a -> a
** a
2 ) forall a. Fractional a => a -> a -> a
/ a
36
                  ,a
n forall a. Fractional a => a -> a -> a
/ a
6
                  ,a
1 forall a. Fractional a => a -> a -> a
/ a
36]

-- | The square of the product of the input sequence is summed, then
-- divided by the square of the sequence length.
--
-- > mean_square_product [(0,0),(1,1),(2,2),(3,3)] == 6.125
-- > mean_square_product [(2,3),(4,5)] == (6^2 + 20^2) / 2^2
mean_square_product :: Fractional n => [(n,n)] -> n
mean_square_product :: forall n. Fractional n => [(n, n)] -> n
mean_square_product [(n, n)]
x =
    let f :: (n, n) -> n
f = forall a. Num a => a -> a
T.square forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a. Num a => a -> a -> a
(*)
        n :: n
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length [(n, n)]
x)
    in forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (n, n) -> n
f [(n, n)]
x) forall a. Fractional a => a -> a -> a
/ forall a. Num a => a -> a
T.square n
n

-- | An incorrect attempt at the description in paragraph two of p.58
-- of the /CMJ/ paper.
--
-- > let p ~= q = abs (p - q) < 1e-4
-- > metrical_affinity [2,3] 1 [3,2] 1 ~= 0.0324
-- > metrical_affinity [2,2,3] 20 [3,5] 16 ~= 0.0028
metrical_affinity :: (Integral n,Show n) => [n] -> n -> [n] -> n -> Double
metrical_affinity :: forall n. (Integral n, Show n) => [n] -> n -> [n] -> n -> Double
metrical_affinity [n]
s1 n
v1 [n]
s2 n
v2 =
    let ([n]
s1',[n]
s2') = forall n. (Integral n, Show n) => S_MM n -> S_MM n -> ([n], [n])
prolong_stratifications ([n]
s1,n
v1) ([n]
s2,n
v2)
        i1 :: [Double]
i1 = forall n. (Integral n, Show n) => Stratification n -> [Double]
relative_indispensibilities [n]
s1'
        i2 :: [Double]
i2 = forall n. (Integral n, Show n) => Stratification n -> [Double]
relative_indispensibilities [n]
s2'
        v :: n
v = forall a. Integral a => a -> a -> a
lcm n
v1 n
v2
        i1' :: [Double]
i1' = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (forall i a. Integral i => i -> a -> [a]
genericReplicate (n
v forall a. Integral a => a -> a -> a
`div` n
v1) [Double]
i1)
        i2' :: [Double]
i2' = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (forall i a. Integral i => i -> a -> [a]
genericReplicate (n
v forall a. Integral a => a -> a -> a
`div` n
v2) [Double]
i2)
    in forall n. Fractional n => [(n, n)] -> n
mean_square_product (forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
i1' [Double]
i2')

-- | An incorrect attempt at Equation 6 of the /CMJ/ paper, see
-- omega_z.
--
-- > let p ~= q = abs (p - q) < 1e-4
-- > metrical_affinity' [2,2,2] 1 [2,2,2] 1 ~= 1.06735
-- > metrical_affinity' [2,2,2] 1 [2,2,3] 1 ~= 0.57185
-- > metrical_affinity' [2,2,2] 1 [2,3,2] 1 ~= 0.48575
-- > metrical_affinity' [2,2,2] 1 [3,2,2] 1 ~= 0.45872
--
-- > metrical_affinity' [3,2,2] 3 [2,2,3] 2 ~= 0.10282
metrical_affinity' :: (Integral t,Show t) => [t] -> t -> [t] -> t -> Double
metrical_affinity' :: forall n. (Integral n, Show n) => [n] -> n -> [n] -> n -> Double
metrical_affinity' [t]
s1 t
v1 [t]
s2 t
v2 =
    let ([t]
s1',[t]
s2') = forall n. (Integral n, Show n) => S_MM n -> S_MM n -> ([n], [n])
prolong_stratifications ([t]
s1,t
v1) ([t]
s2,t
v2)
        ix :: (Integer -> x) -> Integer -> x
        ix :: forall x. (Integer -> x) -> Integer -> x
ix Integer -> x
f Integer
i = case Integer
i of
                   Integer
1 -> Integer -> x
f Integer
1
                   Integer
2 -> Integer -> x
f Integer
2
                   Integer
_ -> forall a. HasCallStack => String -> a
error (forall a. Show a => a -> String
show (String
"ix",Integer
i))
        s :: Integer -> [t]
s = forall x. (Integer -> x) -> Integer -> x
ix (forall n a. Integral n => [a] -> n -> a
at1 [[t]
s1,[t]
s2])
        v :: Integer -> t
v = forall x. (Integer -> x) -> Integer -> x
ix (forall n a. Integral n => [a] -> n -> a
at1 [t
v1,t
v2])
        u :: Integer -> Int
u = forall x. (Integer -> x) -> Integer -> x
ix (forall i a. Num i => [a] -> i
genericLength forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> [t]
s)
        s' :: Integer -> [t]
s' = forall x. (Integer -> x) -> Integer -> x
ix (forall n a. Integral n => [a] -> n -> a
at1 [[t]
s1',[t]
s2'])
        z :: Integer -> t
z = forall x. (Integer -> x) -> Integer -> x
ix (forall i a. Num i => [a] -> i
genericLength forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> [t]
s')
        q :: Integer -> n -> t
q Integer
i n
j = Integer -> [t]
s Integer
i forall n a. Integral n => [a] -> n -> a
`at1` n
j
        omega_u :: Integer -> t
omega_u Integer
i = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map (forall {n}. Integral n => Integer -> n -> t
q Integer
i) [Int
1::Int .. Integer -> Int
u Integer
i])
        omega_z :: p -> t
omega_z p
_ = forall a. Integral a => a -> a -> a
lcm (Integer -> t
v Integer
1 forall a. Num a => a -> a -> a
* Integer -> t
omega_u Integer
1) (Integer -> t
v Integer
2 forall a. Num a => a -> a -> a
* Integer -> t
omega_u Integer
2)
        omega_0 :: t
omega_0 = forall a. Integral a => a -> a -> a
lcm (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (Integer -> [t]
s' Integer
1)) (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (Integer -> [t]
s' Integer
2))
        x0 :: t -> Integer -> t
x0 t
n Integer
i = forall a. (Integral a, Show a) => Stratification a -> a -> a -> a
lower_psi (Integer -> [t]
s' Integer
i) (Integer -> t
z Integer
i) (t
1 forall a. Num a => a -> a -> a
+ ((t
n forall a. Num a => a -> a -> a
- t
1) forall a. (Integral a, Show a) => a -> a -> a
`mod_pos_err` forall {p}. p -> t
omega_z Integer
i))
        x1 :: t -> t
x1 t
n = forall a. Num a => a -> a
T.square (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map (t -> Integer -> t
x0 t
n) [Integer
1,Integer
2]))
        x2 :: t
x2 = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map t -> t
x1 [t
1 .. t
omega_0])
        x3 :: t
x3 = t
18 forall a. Num a => a -> a -> a
* t
x2 forall a. Num a => a -> a -> a
- t
2
        x4 :: p -> t
x4 p
i = forall a. Num a => a -> a
T.square (forall {p}. p -> t
omega_z p
i forall a. Num a => a -> a -> a
- t
1)
        x5 :: t
x5 = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map forall {p}. p -> t
x4 [Integer
1::Integer,Integer
2])
        x6 :: t
x6 = t
7 forall a. Num a => a -> a -> a
* t
omega_0 forall a. Num a => a -> a -> a
* t
x5
        x7 :: Double
x7 = forall n. Integral n => n -> Double
to_r t
x3 forall a. Fractional a => a -> a -> a
/ forall n. Integral n => n -> Double
to_r t
x6
        x8 :: Double
x8 = Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
x7
        x9 :: Double
x9 = forall a. Num a => a -> a
negate (forall a. Fractional a => a -> a
recip Double
x8)
    in forall a b. a -> b -> b
traceShow (forall {p}. p -> t
omega_z,t
omega_0,t
x2,t
x3,t
x5,t
x6,Double
x7,Double
x8,Double
x9) Double
x9