kmn-programming (empty) → 0.9.1
raw patch · 6 files changed
+586/−0 lines, 6 filesdep +basedep +optparse-applicativedep +random
Dependencies added: base, optparse-applicative, random, random-shuffle, time, x86-64bit
Files
- ChangeLog.md +12/−0
- KMNProgramming.hs +224/−0
- LICENSE +30/−0
- Main.hs +283/−0
- SplitProblem.hs +9/−0
- kmn-programming.cabal +28/−0
+ ChangeLog.md view
@@ -0,0 +1,12 @@+# Revision history for kmn-programming++## 0.9.0.1 -- 2020-03-06++* Add Wmat44.txt+* Add SplitProblem.hs+* Add documentation in README.md++## 0.9.0.0 -- 2016-12-12++* First version.+
+ KMNProgramming.hs view
@@ -0,0 +1,224 @@+{-# language ViewPatterns #-}+{-# language RecursiveDo #-}+{-# language FlexibleInstances #-}+{-# language TypeSynonymInstances #-}+module KMNProgramming where++import Control.Monad hiding (unless)+import Control.Exception (evaluate)+import Foreign++import CodeGen.X86++{- register usage+eax: maximum value+rbx: level pointer+rcx: (aux) loop counter+rdx: (aux) loop pointer+rsp: (unused) stack pointer+rbp: (aux) loop pointer+rsi: (static) -16+rdi: (static) vector pointer & start level+r8: (static) trace level+r9: (static) stop level+r10: (static) debug flag+r11: (static) level length+r12: (static) loop length+r13: result bits (high)+r14: result bits (low)+r15: (aux) vector norm+-}++foreign import ccall "static stdio.h &fflush" fflush :: FunPtr a++type T = Ptr Int32 -> Int32+foreign import ccall "dynamic" callPBII :: FunPtr T -> T+instance Callable T where dynCCall = callPBII++foreign import ccall "dynamic" callIO :: FunPtr (IO ()) -> IO ()+instance Callable (IO ()) where dynCCall = callIO++umes_ :: Int -> Bool -> Int -> Int -> Int -> Int -> Int -> [Int32] -> [[Int32]] -> IO (Int32, [Int32])+umes_ guess trr ali ali' tr uroll lev levs m = umes guess trr ali ali' tr uroll lev levs (length m) (length $ head m) m++{-# NOINLINE umes #-}+umes :: Int -> Bool -> Int -> Int -> Int -> Int -> Int -> [Int32] -> Int -> Int -> [[Int32]] -> IO (Int32, [Int32])+umes guess trr ali ali' tr uroll lev levs rows_ d_ = fun+ where+ fun (align' (4*uroll) -> a: as)+ = withArrayAligned (concat $ a: zipWith (++) ((: [0,0,0]) <$> (replicate (rows-length levs) 0 ++ levs)) as) (2^8) $ \r -> do+ x <- evaluate $ code r+ xs <- forM [d+d1*lev, d+d1*(lev+1)..d+d1*(rows-1)] $ peekElemOff r+ return (x, xs)++ rows = rows_ - 1+ d = d_ `modd` (4*uroll)+ d1 = d+4++ modd a b = b * ((a - 1) `div` b + 1)++ code :: Ptr Int32 -> Int32+ code = compile $ saveNonVolatile $ mdo+ mov rdx $ fromIntegral $ 4*d1*(rows-length levs)+ xor_ rax rax+ mloop <- label+ cmp rdx $ fromIntegral $ 4*d1*lev+ unless Z $ do+ save [rdi, rdx] $ do+ xor_ r10 r10+ call $ ipRelValue main+ mov (addr $ rdi + rdx - 16) eax+ sub rdx $ fromIntegral $ 4*d1+ jmp mloop++ xor_ rdx rdx+ mov r10 $ fromIntegral $ fromEnum trr+ mov rax $ fromIntegral guess+ call $ ipRelValue main+ jmp end++ main <- label+ lea r9 $ addr8 $ rdi + fromIntegral (4 * (rows + 1) * d1)+ lea r8 $ addr8 $ rdi + fromIntegral (4 * tr * d1)+ mov r11 $ fromIntegral $ 4*d1+ mov r12 $ fromIntegral $ d `div` (4*uroll)+ mov rsi (-16)+ xor_ r13 r13+ xor_ r14 r14+ lea rdi $ addr8 $ rdi + rdx+ mov rbx rdi+ jmp chk++ trace <- label+ save [rax, rcx, rdx, rbx, rbp, rsi, rdi, r8, r9, r10, r11, r12, r13, r14, r15] $ do+ leaData arg1 $ CString "%ld.\n"+ mov arg2 r14+ printfflush+ ret++ trace2 <- label+ save [rax, rcx, rdx, rbx, rbp, rsi, rdi, r8, r9, r10, r11, r12, r13, r14, r15] $ do+ leaData arg1 $ CString "%016lx%016lx %d\n"+ mov arg2 r13+ mov arg3 r14+ mov (resizeOperand arg4) eax+ printfflush+ ret++ align ali+ chk <- label+ add rbx r11+ cmp rbx r9+ unless E $ mdo+ mov r15d $ addr $ rbx + rsi+ test r15d r15d+ if_ NZ (mdo+ mes+ cmp r15d eax+ j LE up) $ mdo+ cmp rbx r8+ unless NE $ call $ ipRelValue trace+ stc+ rcl r14 1+ when (rows > 64) $ rcl r13 1+ copyAdd+ jmp chk+ xor_ r15d r15d+ mes+ cmp r15d eax+ unless NA $ do+ mov eax r15d+ test r10 r10+ unless Z $ call $ ipRelValue trace2++ align ali+ up <- label+ sub rbx r11+ cmp rbx rdi+ unless E $ mdo+ test r14 1+ j Z zero+ xor_ r14 1+ copySub+ jmp chk+ zero <- label+ if rows > 64 then do+ shr r13 1+ rcr r14 1+ else do+ shr r14 1+ copyAdd+ jmp up+ ret++ end <- label+ return ()++ mes = mdo+ mov rcx r12+ mov rdx rdi+ pxor xmm2 xmm2+ align ali'+ doWhile NZ $ mdo+ replicateM_ uroll $ mdo+ movdqa xmm0 $ addr rdx+ sub rdx rsi+ movdqa xmm1 xmm0+ psrad xmm0 31+ pxor xmm1 xmm0+ psubd xmm2 xmm0+ paddd xmm2 xmm1+ sub rcx 1+ movdqa xmm1 xmm2+ psrldq xmm1 8+ paddd xmm2 xmm1+ movd ecx xmm2+ add r15d ecx+ psrldq xmm2 4+ movd ecx xmm2+ add r15d ecx++ copyAdd = copy $ mdo+ movdqa xmm0 $ addr rbp+ sub rbp rsi+ paddd xmm0 $ addr rdx+ movdqa (addr rdx) xmm0+ sub rdx rsi++ copySub = copy $ mdo+ movdqa xmm0 $ addr rbp+ sub rbp rsi+ movdqa xmm1 $ addr rdx+ paddd xmm0 xmm0+ psubd xmm1 xmm0+ movdqa (addr rdx) xmm1+ sub rdx rsi++ copy body = mdo+ mov rcx r12+ mov rbp rbx+ mov rdx rdi+ align ali'+ doWhile NZ $ mdo+ replicateM_ uroll body+ sub rcx 1++ printfflush = do+ xor_ rax rax+ callFun rbp printf+ xor_ arg1 arg1+ callFun rbp fflush++align' n xs = (++ pad) <$> xs+ where+ pad = replicate ((-length (head xs)) `mod` n) 0++withArrayAligned bs al cont = allocaBytesAligned (sizeOf (head bs) * length bs) al $ \r -> do+ forM_ (zip [0..] bs) $ uncurry $ pokeElemOff r+ cont r++save rs c = do+ sequence_ [push r | r <- rs]+ c+ sequence_ [pop r | r <- reverse rs]+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2016-2020, Péter Diviánszky++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Péter Diviánszky nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Main.hs view
@@ -0,0 +1,283 @@+{-# language PatternGuards #-}+{-# language ViewPatterns #-}+{-# language NoMonomorphismRestriction #-}+{-# language BlockArguments #-}+import Data.Function+import Data.List+import Data.Monoid+import Data.Maybe+import Data.Time.Clock+import Control.Arrow+import Control.Monad+import Control.Concurrent+import Control.Concurrent.Chan+import Options.Applicative hiding (asum)+import System.IO+import System.Random+import System.Random.Shuffle++import KMNProgramming++data Method = Masum | Mrandom (Maybe Int) | Mb | Mnosort+ deriving (Read, Show, Eq, Ord)++readMethod = eitherReader $ \s -> case s of+ "nosort" -> Right Mnosort+ "asum" -> Right Masum+ "b" -> Right Mb+ (splitAt 6 -> ("random", s))+ | "" <- s -> Right $ Mrandom Nothing+ | [(i, "")] <- reads s -> Right $ Mrandom $ Just i+ _ -> Left $ "available sort methods: asum random randomSEED b nosort"++main :: IO ()+main = do+ hSetBuffering stdout LineBuffering+ join $ execParser opts+ where+ opts = info (helper <*> options)+ ( fullDesc+ <> progDesc "Maximalize sum of the input matrix multiplied by tensor products of two vectors of +-1 elements"+ <> header "kmn-programming - specialized quadratic binary optimization")++ options :: Parser (IO ())+ options = + hsubparser+ ( command "sample" sampleCommand+ <> command "test" testCommand+ <> command "timerandom" randomCommand+ )+ <|> computeOptions++ sampleCommand = flip info (progDesc "print sample matrix") $ printTestMat+ <$> optional (strOption $ short 'o' <> long "output" <> metavar "FILE" <> help "output file" <> action "filenames")+ <*> argument auto (metavar "NAT")++ testCommand = flip info (progDesc "basic self-test") $ testTestMat+ <$> optional (argument auto $ metavar "NAT")++ randomCommand = flip info (progDesc "measure optimization time on random matrices") $ timeRandom+-- <$> optional (strOption $ short 'o' <> long "output" <> metavar "FILE" <> help "output file" <> action "filenames")+ <$> argument auto (metavar "NAT" <> help "biggest integer in random matrices")+ <*> argument auto (metavar "NAT" <> help "smallest width")+ <*> argument auto (metavar "NAT" <> help "biggest width")+ <*> argument auto (metavar "NAT" <> help "width step")+ <*> argument auto (metavar "NAT" <> help "smallest height")+ <*> argument auto (metavar "NAT" <> help "biggest height")+ <*> argument auto (metavar "NAT" <> help "height step")+ <*> argument auto (metavar "NAT" <> help "smallest level")+ <*> argument auto (metavar "NAT" <> help "biggest level")+ <*> argument auto (metavar "NAT" <> help "repeat computation")++ computeOptions = compute+ <$> switch (short 'd' <> long "delete" <> help "delete 0 rows and columns")+ <*> switch (short 't' <> long "transpose" <> help "transpose matrix if it has more rows than columns")+ <*> optional (option readMethod $ short 's' <> long "sort" <> metavar "SORTMETHOD" <> help "sort method - default is nosort" <> completeWith ["asum","random","b","nosort"])+ <*> switch (short 'm' <> long "multiply" <> help "multiply rows by +-1 to improve the first guess")+ <*> switch (short 'p' <> long "print" <> help "print matrix before optimization")+ <*> switch (short 's' <> long "silent" <> help "print just the result")+ <*> (fromMaybe 0 <$> optional (option auto $ short 'g' <> long "guess" <> metavar "NAT" <> help "guessed result - default is 0" <> completeWith ["0"]))+ <*> optional (strOption $ long "levelin" <> metavar "FILE" <> help "precomputed levels input file" <> action "filenames")+ <*> optional (strOption $ long "levelout" <> metavar "FILE" <> help "levels output file" <> action "filenames")+ <*> optional (option auto $ short 'l' <> long "level" <> metavar "NAT" <> help "level - default is number of rows / 4")+ <*> (fromMaybe 0 <$> optional (option auto $ long "trace" <> metavar "NAT" <> help "trace level - default is 0" <> completeWith ["0"]))+ <*> (fromMaybe 4 <$> optional (option auto $ short 'u' <> long "unroll" <> metavar "NAT" <> help "unroll cycles - default is 4" <> completeWith ["4"]))+ <*> (fromMaybe 8 <$> optional (option auto $ short 'a' <> long "align" <> metavar "NAT" <> help "code alignment - default is 8" <> completeWith ["8"]))+ <*> optional (strOption $ short 'o' <> long "output" <> metavar "FILE" <> help "output file" <> action "filenames")+ <*> ( {-Right . mkMatrix . filter (not . null) . concatMap splitSC <$ argument (readChars "matrix") (metavar "MATRIX") <*> many (argument (maybeReader readMat) (metavar "INT"))+ <|> -}Left <$> argument str (metavar "FILE" <> action "filenames")+ )+ <*> optional (option auto $ long "timeout" <> help "timeout in seconds")+ <*> (fromMaybe (Partial (Just 1) 0) <$> optional (option auto $ long "partial" <> help "do partial computation"))++readChars c = maybeReader f where+ f c' | c == c' = Just ()+ f _ = Nothing++readMat s@(_:_) = Just s+readMat _ = Nothing++mkMatrix s = case span (/= ",") s of+ (i, ",": s) -> i: mkMatrix s+ (i, []) -> [i]++splitSC s = case span (/= ',') s of+ (i, ',': s) -> i: ",": splitSC s+ (i, []) -> [i]++checkPred err p | p = return ()+checkPred err _ = error err++data Partial+ = Partial (Maybe Int) Int+ deriving (Eq, Ord, Show)++instance Read Partial where+ readsPrec _ ('*':'/':'2':'^': s)+ | [(j,s)] <- reads s+ = [(Partial Nothing j, s)]+ readsPrec _ s+ | [(i,'/':'2':'^':s)] <- reads s+ , [(j,s)] <- reads s+ = [(Partial (Just i) j, s)]+ readsPrec _ _ = []++compute+ :: Bool -> Bool -> Maybe Method -> Bool -> Bool -> Bool -> Int -> Maybe FilePath -> Maybe FilePath -> Maybe Int -> Int -> Int -> Int+ -> Maybe FilePath -> Either FilePath [[String]] -> Maybe Int -> Partial -> IO ()+compute del transp met mult printmat silent guess levelin levelout level_ tra uroll ali out fname timeout (Partial tasks splitPower) = do+ gen <- newStdGen+ s <- either (fmap (filter (not . null) . map words . lines) . readFile) return fname++ levs <- case levelin of+ Nothing -> return []+ Just f -> map read . lines <$> readFile f++ let vs = map (map read) s :: [[Integer]]++ mat_ = (if transp then transpose' else id)+ . (if del then transpose . filter (any (/=0)) . transpose . filter (any (/=0)) else id)+ $ map fromIntegral <$> vs++ transpose' xs = if length xs <= length t then xs else t+ where+ t = transpose xs++ ns@(~[width]) = nub $ length <$> mat_+ output = maybe putStrLn writeFile out++ method = case met of+ Nothing -> id+ Just Mnosort -> id+ Just Masum -> sortByASum+ Just Mb -> sortByASum'+ Just (Mrandom mi) -> \xs -> shuffle' xs (length xs) $ maybe gen mkStdGen mi++ mat = (if mult then multguess else id)+ . method+ $ mat_++ nx = length mat - splitPower+ level = fromMaybe (min (nx - length levs - 1) $ round $ fromIntegral nx / 4) level_++ checkPred "matix row lengths differ" $ length ns == 1+ checkPred "number of rows should be less than 129" $ nx <= 128+ checkPred "alignment should be 1,2,4,8,16,..." $ ali `elem` map (2^) [0..10]+ checkPred "unroll should be greater than 1" $ uroll >= 1+-- checkPred ("unroll " ++ show uroll ++ " less or equal then the width " ++ show width $ uroll <= width+ checkPred "level should not be negative" $ level >= 0+ checkPred "level should be less than the number of rows" $ level < nx+ checkPred "trace level should not be negative" $ tra >= 0+ checkPred "trace level should be less or equal than level" $ tra <= level++ when (not silent) $ do+ when (asum (concat vs) >= 2^31) $ output "!!! warning: overflow may happen !!!"+ output $ show nx ++ " x " ++ show width+ output $ "levels: " ++ show level++ let tasks' = maybe [1..2^splitPower] pure tasks++ resChan <- newChan++ tids <- forM tasks' \i -> forkIO $ do++ let mat' = mkPartial (2^splitPower-i) splitPower mat++ when (not silent && printmat && i == head tasks') $ output $ showMat mat'++ (res, levs') <- umes_ guess (not silent) ali ali tra uroll level levs mat'++ case levelout of+ Just f | i == head tasks' -> writeFile f $ unlines $ map show levs'+ _ -> return ()++ writeChan resChan res++ case timeout of+ Nothing -> return ()+ Just t -> void $ forkIO $ do+ threadDelay $ 1000000 * t+ forM_ tids killThread++ ress <- forM tasks' \_ -> readChan resChan+ output $ show $ maximum ress++mkPartial i 0 mat = mat+mkPartial i j (x: y: mat) = mkPartial (i`div`2) (j-1) (zipWith (if odd i then (+) else (-)) x y: mat)++timeRandom :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> IO ()+timeRandom r y y' sy x x' sx (fromIntegral -> z) (fromIntegral -> z') rep = forM_ (f $ zip [x, x+sx..x'] [y, y+sy..y']) $ \((i, j), z) -> do+ putStr $ show (i,j,z) ++ " " ++ show (log $ fromIntegral i) ++ " "+ mats <- replicateM rep $ replicateM i $ replicateM j $ fromIntegral <$> randomRIO (-r, r)+ t1 <- getCurrentTime+ let opt = umes 0 False 8 8 0 1 z [] i j+ forM_ mats $ \mat -> fst <$> opt mat+ t2 <- getCurrentTime+ putStrLn $ show $ log $ (realToFrac (diffUTCTime t2 t1) :: Double) / fromIntegral rep+ where+ f xs = zip xs $ map round [z, z + (z' - z)/(fromIntegral (length xs) - 1) ..]++printTestMat :: Maybe FilePath -> Int -> IO ()+printTestMat out n = maybe putStrLn writeFile out $ showMat $ testMat n++testTestMat :: Maybe Int -> IO ()+testTestMat n_ = do+ res <- fst <$> umes 0 False 8 8 0 1 (length m `div` 4) [] (length m) (length $ head m) m+ case () of+ _ | f n == res -> putStrLn "OK"+ | otherwise -> error "fatal error"+ where+ n = fromMaybe 6 n_+ m = testMat n+ f 2 = 4+ f 3 = 36+ f 4 = 120+ f 5 = 280+ f 6 = 540+ f 7 = 924++testMat n = [[(if i == j then 1 else 3) * inner i j | i<-fs] | j<-fs]+ where+ fs = [ins b 1 $ ins a 1 $ replicate (n-2) 0 | a <- [0..n-2], b <- [a+1..n-1]]+ ++ [ins b (-1) $ ins a 1 $ replicate (n-2) 0 | a <- [0..n-2], b <- [a+1..n-1]]++ ins i a (splitAt i -> (as, bs)) = as ++ a: bs++showMat = unlines . map unwords . transpose . map (pad . map (g . show)) . transpose+ where+ pad xs = map f xs+ where+ f = take n . (++ repeat ' ')+ n = maximum $ 0: map length xs+ g s@('-': _) = s+ g s = ' ': s++sortByASum = map snd . sortBy (flip compare `on` fst) . map (asum &&& id)++asum = sum . map abs++inner a b = sum $ zipWith (*) a b++sortByASum' [] = []+sortByASum' xs = reverse $ a: f a as+ where+ (a, as) = snd $ maximumBy (flip compare `on` fst) $ map (asum . fst &&& id) $ getOut xs++ f _ [] = []+ f x xs = a: f a as+ where+ (a, as) = snd $ maximumBy (flip compare `on` fst) $ map (abs . inner x . fst &&& id) $ getOut xs++ getOut :: [a] -> [(a, [a])]+ getOut xs = [(x, as ++ bs) | (as, x:bs) <- zip (inits xs) (tails xs)]++multguess (x: xs) = x: f x xs+ where+ f _ [] = []+ f v (x: xs)+ | asum a >= asum b = x: f a xs+ | otherwise = map negate x: f b xs+ where+ a = zipWith (+) v x+ b = zipWith (-) v x+
+ SplitProblem.hs view
@@ -0,0 +1,9 @@+import Control.Monad+import System.Environment++main = do+ (s: ss) <- getArgs+ let i = read s+ forM_ [1..2^i] $ \j -> do+ putStrLn $ unwords $ "kmn-programming": ss ++ ["--partial", show j ++ "/2^" ++ show i]+
+ kmn-programming.cabal view
@@ -0,0 +1,28 @@+name: kmn-programming+version: 0.9.1+synopsis: K_M,N quadratic programming+-- description: +license: BSD3+license-file: LICENSE+author: Péter Diviánszky+maintainer: divipp@gmail.com+-- copyright: +category: Math+build-type: Simple+extra-source-files: ChangeLog.md+ SplitProblem.hs+cabal-version: >=1.10++executable kmn-programming+ main-is: Main.hs+ ghc-options: -threaded -rtsopts=all -O2+ other-modules: KMNProgramming+ other-extensions:+ build-depends: base >=4.9 && <4.17+ , x86-64bit >=0.4 && <0.5+ , optparse-applicative >=0.13 && <0.17+ , random+ , random-shuffle+ , time+ -- hs-source-dirs: + default-language: Haskell2010