Search algoritm for use with chord spaces
Donya Quick and Paul Hudak
Last modified: 13-Jan-2016

Implementation of a search algorithm for traversing chord spaces using
let-in constraints as well as progression level predicates.


> module Kulitta.Search where
> import Data.List
> import Kulitta.ChordSpaces
> import Control.DeepSeq
> import Control.Parallel.Strategies
> import System.Random
> import Kulitta.PTGG hiding (choose)
> import Kulitta.Grammars.MusicGrammars
> import Kulitta.PostProc
> import Kulitta.Constraints

> type Constraints = [[(Int, Int)]]
> type Index = Int
> type Bound = Int

> printIt = True -- flag for whether to show evidence of search progress
> printIter = 5000 -- number of solutions to check before reporting back

========= OBTAINING SOLUTIONS =========

The allSolns function finds all R-equivalent sequences given a quotient
space, qs, and a relation, r.

This is the obvious, but slightly less efficient way to obtain solutions.
The lazy evaluation in Haskell allows it to do aggressive pruning similarly
to how pairProg handles it.

> allSolns :: (Eq a, Show a) => QSpace a -> EqRel a -> [a] -> [[a]]
> allSolns qs r [] = [[]]
> allSolns qs r (x:xs) = 
>     [(y:ys) | y <- eqClass qs r x, ys <- allSolns qs r xs]

Unfortunately, allSolns won't work with really big problems. Large chord 
spaces and long progressions (even with small chord spaces) cause an 
exponential blowup in the search space. Two alternative options below
prune the solution space more aggressively.

The pairProg function applies a predicate as it generates the 
solutions.

> pairProg :: (Eq a, Show a) => QSpace a -> EqRel a -> Predicate (a,a) -> [a] -> [[a]]
> pairProg qs r c [] = []
> pairProg qs r c [x] = map (\a -> [a]) $ eqClass qs r x
> pairProg qs r c (x:xs) = 
>     let endSolns = pairProg qs r c xs -- map xs to a new progression
>         f y soln = c (y, head soln) -- filter using the classifier
>         newSolns = [(y:ys) | y<-eqClass qs r x, ys<-endSolns, f y ys]
>     in  if null newSolns then error "No solutions that satisfy the consraints!"
>         else newSolns

Now we define the greedy algorithm, greedyProg.

We define choose, a function to stochastically select an element
from a list.

> choose :: StdGen -> [a] -> (StdGen, a)
> choose g xs = 
>     let (r, g') = next g
>     in  (g', xs !! (r `mod` length xs))

And finally the recursive, greedy function, greedyProg, and its
"helper," greedyChord.


> type Fallback a = EqClass a -> StdGen -> a -> (StdGen, a)
		 	 
> greedyProg :: (Eq a, Show a) => QSpace a -> EqRel a -> 
>     Predicate (a,a) -> Fallback a -> StdGen -> [a] -> [a]
> greedyProg qs r c f g [] = []
> greedyProg qs r c f g (x:xs) =
>     let e = eqClass qs r x
>         (g', y0) = choose g e -- randomly choose the first AbsChord
>     in  greedyRec qs r c f g y0 xs where 
>         greedyRec qs r c f g y pts = 
>             let (g', yi) = greedyChord (eqClass qs r (head pts)) y c f g
>             in  if null pts then [y] 
>                 else y : greedyRec qs r c f g' yi (tail pts)

> greedyChord :: (Eq a, Show a) => EqClass a -> a -> Predicate (a,a) -> 
>     Fallback a -> StdGen -> (StdGen, a)
> greedyChord e yprev hpair f g = 
>     let (rand, g') = next g
>         yxs = zip (repeat yprev) e
>         ys = map snd $ filter hpair yxs
>     in  if null ys then f e g' yprev -- fallback case
>         else (g', ys !! (rand `mod` length ys)) -- ok case


Now we need some fall-back functions. We can construct them from
predicates.

> fallFun :: (Show a) => Predicate (a,a) -> EqClass a -> StdGen -> a -> (StdGen, a)
> fallFun c e g x = 
>     let ys = map snd $ filter c $ zip (repeat x) e
>         (i, g') = next g
>     in  if null ys then error ("Stuck at symbol "++show x++". No viable options left!")
>         else (g', ys !! (i `mod` length ys))

This is the default fallback function we used in our experiments:

> defFall = fallFun (maxClass 7)


And also a nearest neighbor fallback function that would always 
succeed if the equivalence class has at least one element (which 
should always be the case).

> nearFall ::  EqClass AbsChord -> StdGen -> AbsChord -> (StdGen, AbsChord)
> nearFall e g x = 
>     let ds = map (simpleDist x) e :: [Double]
>         y = e !! (head $ findIndices (==minimum ds) ds)
>     in  (g, y)


This version of greedyProg operates over a list of equivalence classes.

> greedyProg' :: (Eq a, Show a) => 
>     Predicate (a,a) -> Fallback a -> StdGen -> [EqClass a] -> [a]
> greedyProg' c f g [] = []
> greedyProg' c f g (e:es) =
>     let (g', y0) = choose g e -- randomly choose the first AbsChord
>     in  greedyRec' c f g y0 es where 
>         greedyRec' c f g y pts = 
>             let (g', yi) = greedyChord (head pts) y c f g
>             in  if null pts then [y] 
>                 else y : greedyRec' c f g' yi (tail pts)

> pairProg' :: (Eq a, Show a) => Predicate (a,a) -> [EqClass a] -> [[a]]
> pairProg' c [] = [[]]
> pairProg' c (e:es) = 
>     let newSolns = [(y:ys) | y<-e, ys<-pairProg' c es, 
>                     if null ys then True else c (y, head ys)]
>     in  if not $ null newSolns then newSolns
>         else error "No solutions that satisfy the consraints!"


======================

LET CONSTRAINTS

Let-in expressions impose constraints on the interpretation of
the results as music. These can be viewed as progression-level
constraints, or predicates.

mkCons fundtion looks through the let-in structure of a Term
and produces a progression-level predicate for use with the
ChordSpaces module.

> mkCons :: (Eq a) => [Term b c] -> Predicate [a]
> mkCons t xs = toCons (findInds [] t) xs where
>     toCons :: (Eq a) => [[(Int, Int)]] -> [a] -> Bool
>     toCons [] xs = True
>     toCons (c:cs) xs = 
>         let f (i,j) = take (j+1-i) $ drop i xs 
>         in  (and $ map (f (head c) == ) $ 
>             map f $ tail c) && toCons cs xs


The findInds function looks through a Term for let-in expressions.
When it finds Let x a exp, it calls findIndsSub on x and exp to 
determine which indices in the sequence are occupied by instances
of x. 

> findInds :: [(String, [Term a b])] -> [Term a b] -> [[(Int, Int)]]
> findInds e [] = []
> findInds e (t:ts) = let rest = findInds e ts in case t of 
>     Var x -> undefined
>     NT x -> map (map (pAdd 1)) $ findInds e ts
>     Let x a exp -> 
>         let a' = expand e a
>             exp' = expand ((x,a'):e) exp
>         in  findInds e a ++ findIndsSub x (length a') (expand' e exp) : 
>             map (map (pAdd $ length exp')) rest

> expand' :: [(String, Sentence a b)] -> Sentence a b -> Sentence a b
> expand' e [] = []
> expand' e (t:ts) = case t of 
>     Let x a exp -> expand' ((x, expand' e a) : e) exp ++ expand' e ts
>     Var x       -> case lookup x e of 
>                        Nothing -> Var x : expand' e ts
>                        Just a -> a ++ expand' e ts
>     x           -> x : expand' e ts


The findIndsSub function looks for instances of a variable and
determines what indices they occupy.

> findIndsSub :: String -> Int -> [Term a b] -> [(Int, Int)]
> findIndsSub x xLen [] = []
> findIndsSub x xLen (t:ts) = let rest = findIndsSub x xLen ts in case t of
>     Var y -> if x==y then (0, xLen-1) : map (pAdd xLen) rest else rest
>     NT y  -> map (pAdd 1) $ findIndsSub x xLen ts
>     Let y a e -> error "(find Instances) This point should be unreachable."

> pAdd amt (a,b) = (a+amt, b+amt)

========================

SEARCH IMPLEMENTATION


All of these functions assume that constraints are SORTED. This is used to  
"jump ahead" to the next solution that would meet constraints early-on in the 
depth-first traversal. We give preference to changing indices further right 
in the sequence. It is assumed that the constraints are fully sorted (inner 
and outer) by doing the following for k :: Constraints:

	sort (map sort k)
	
If this is not done, then the constraints cannot be guaranteed to be satisfied.
The findSoln function looks for the first solution satisfying two types of 
constraints: those from let-in statements and more generic, progression-level
constraints supplied as a Predicate function.

> findSoln :: (Eq a, Show a) => 
>     Constraints -> Predicate [a] -> [[a]] -> (Int, [a])
> findSoln k f ecs = 
>     let n = length ecs
>         bs = map (`elem` freeInds n k) [0..length ecs-1]
>         lens = map length ecs 
>         g i = zip3 bs i lens
>         frec j iCurr = 
>             let iNext = findNext2 k $ g iCurr 
>                 soln = zipWith (!!) ecs iCurr
>             in  if last iCurr < 0 then error "No more solutions." else
>                 if f soln then (j, soln) else frec (j+1) iNext
>     in  frec 0 $ take (length ecs) $ repeat 0 

> findSoln2 :: (NFData a, Eq a, Show a) => 
>     Constraints -> Predicate [a] -> [[a]] -> IO (Int, [a])
> findSoln2 k f ecs = 
>     let initVal = take (length ecs) $ repeat 0
>         bs = map (`elem` freeInds (length ecs) k) [0..length ecs-1]
>         lens = map length ecs 
>         g i = zip3 bs i lens
>         frec j iCurr = do
>             let iNext = findNext2 k (g iCurr)
>                 soln = makeSoln 0 ecs iCurr
>                 jNext = j+1
>             if last iCurr < 0 then error "No more solutions." else
>                 if f (force soln) then return (j, soln) else -- SEE NOTE BELOW
>                     if j `mod` printIter == 0 && printIt then 
>                             putStrLn ("Solutions examined: "++show j) >> 
>                             frec jNext iNext
>                     else frec jNext iNext
>     in  frec 0 initVal 

Note on use of force in findSoln2: if force is not used and used exactly where
it has been placed above, the function experiences a thunk leak. Printing iCurr, 
iNext, or soln can help to knock back memory usage, but it still leaks until 
printed unless force is used. Using seq can slow the leak to a trickle, but it 
still occurs and will cause long searches to crash. Using force appears to be
the best solution for this leak.


> findSolnPar :: (NFData a, Eq a, Show a) => 
>     Constraints -> Predicate [a] -> [[a]] -> Int -> IO (Int, [a])
> findSolnPar k f ecs parSize = 
>     let initVal = take (length ecs) $ repeat 0
>         bs = map (`elem` freeInds (length ecs) k) [0..length ecs-1]
>         lens = map length ecs 
>         g i = zip3 bs i lens
>         frec j iCurr = do
>             let iCurrs = filter ((>=0).last) (iCurr : findNextI k (g iCurr) (parSize - 1))
>                 iNext = findNext2 k $ g $ last iCurrs -- for next iteration
>                 solns = filter f $ seq (force iCurrs) 
>                         (map (makeSoln 0 ecs) iCurrs `using` parList rdeepseq)
>                 jNext = j+length iCurrs
>             if last iCurr < 0 then error "No more solutions." else
>                 if not $ null solns then return (j, head solns) else
>                     if printIt then 
>                         putStrLn ("Solutions examined: "++show j) >> frec jNext iNext
>                     else frec jNext iNext
>     in  frec 0 initVal


> makeSoln :: Int -> [[a]] -> [Int] -> [a]
> makeSoln j [] [] = []
> makeSoln j [] is = error "(makeSoln) Not enough equivalence relations!"
> makeSoln j ecs [] = error "(makeSoln) Not enough indices!"
> makeSoln j (e:ecs) (i:is) = 
>     if i < length e then (e !! i) : makeSoln (j+1) ecs is
>     else error ("Bad index at position "++show j++": "++show i++", classes="++show (length e))

> fetch xs i = if i >= length xs then error "(fetch) Index is too large!" else xs !! i


The findNext function performs the actual traversal of the space. It bypasses
a lot of irrelevant points by only looking at points that at least satisfy 
the let-in constraints. Those points may or may not satisfy the Predicate, f, 
and so the function may have to explore many solutions.


> findNext :: Constraints -> [Index] -> [Bound] -> [Index]
> findNext k is lens = 
>     let bs = map (`elem` freeInds (length is) k) [0..length is-1]
>         xs = zip3 bs (is) (lens)
>     in  foldl applyCons (incr xs) k 


This version is used for one of the recursive implementations.

> findNext2 :: Constraints -> [(Bool, Index, Bound)] -> [Int]
> findNext2 k xs = foldl' applyCons (incr2 xs) k

> findNextI :: Constraints -> [(Bool, Index, Bound)] -> Int -> [[Index]]
> findNextI k xs i = 
>     let xs' = findNext2 k xs
>         fixWith = zipWith (\(a,b,c) d -> (a,d,c))
>     in  xs' : take (i-1) (iterate (findNext2 k . fixWith xs) xs') where 

The applyCons function applys let-in constraints to an index list. Indices
on the left are given preference when satisfying constraints.

> applyCons :: [Int] -> [(Int, Int)] -> [Int]
> applyCons inds [] = inds
> applyCons inds ((i,j):ijs) = 
>     foldl (f val) inds (map fst ijs) where -- might need seq here
>     val = (take (j-i+1) $ drop i inds)
>     f val src i = take i src ++ val ++ 
>                   drop (i + length val) src
	
Solutions are explored by keeping track of a list of indices into equivalence
classes. This list is incremented by considering which indices can actually 
be changed (those unconstrained by let-in expressions). If the end of the 
search space is reached, an error is thrown since there is no solution.	



The following modification avoids a memory error if used for traversal
with something like iterate. It returns a "flag" at the end (-1) if 
no solutions exist. This must be picked up by the calling functions for
termination to occur.

> incr :: [(Bool, Int, Int)] -> [Int]
> incr ((b,i,l):xs) = let is = map (\(x,y,z) -> y) xs in
>     if b then if i >= l-1 then 0 : incr xs else (i+1) : is
>     else i : incr xs
> incr [] = error "No more solutions!"


The following version uses a flag instead of an error message:

> incr2 :: [(Bool, Int, Int)] -> [Int]
> incr2 ((b,i,l):xs) = let is = map (\(x,y,z) -> y) xs in
>     if b then if i >= l-1 then 0 : incr2 xs else (i+1) : is
>     else i : incr2 xs
> incr2 [] = [-1] -- no more solutions exist


The freeInds function turns let-in constraints into a list of indices that
can be updated or incremented.

> freeInds :: Int -> Constraints -> [Int]
> freeInds n k = 
>     let k' = map (map (\(i,j) -> [i..j])) k 
>         t = nub $ concat $ concatMap tail k'
>     in  filter (not . (`elem` t)) [0..n-1] 


==================================

THIS VERSION OF GREEDY PROG DOES NOT WORK - LETS ARE NOT DONE RIGHT

A version of greedyProg to support Let statements. Constraints are satisfied from 
left to right. Breaks are likely to occurr with repeats. For example, with the
progression

	let x = ... in x x

a break in voice-leading behavior is likely to occur between the last chord of the
first instance of x and the first chord of the second instance of x. There is 
currently no way around this using the greedy approach. The more exact searches
further up in this file are alternatives in such cases.

> greedyLet :: (Eq a, Show a) => Predicate (a,a) -> Fallback a -> Constraints -> 
>              [EqClass a] -> StdGen -> [a]
> greedyLet p f k es g = 
>     let n = length es
>         cs = greedyProg' p f g es
>         consPat = foldl applyCons [0..n-1] (sort k) 
>     in  map (cs !!) consPat

> greedyLetT :: QSpace AbsChord -> EqRel AbsChord -> Predicate (AbsChord,AbsChord) -> 
>               Fallback AbsChord -> Constraints -> [TChord] -> StdGen -> [TChord]
> greedyLetT q r p f k cs g = 
>     let justCs = map tnP cs
>         es = map (eqClass q r) justCs
>         justCs' = greedyLet p f k es g
>     in  zipWith newP cs justCs'