packages feed

sigma-ij 0.2 → 0.2.0.1

raw patch · 3 files changed

+381/−386 lines, 3 filesdep +sigma-ij

Dependencies added: sigma-ij

Files

+ main/sigmaij.hs view
@@ -0,0 +1,377 @@++-- | Calculates Thom polynomial of Sigma^{ij} with localization +-- and the substituion trick+--+-- Some example usages:+--+-- > sigma-ij -h                           # help+-- > sigma-ij -i3 -j2 -n7                  # compute @Tp(Sigma^{3,2}(7))@ in the default ring+-- > sigma-ij -i3 -j2 -n7  -rZp            # compute @Tp(Sigma^{3,2}(7))@ in the hard-coded prime field+-- > sigma-ij -i3 -j2 -n10 -rZp -b3 -B10   # compute the 3rd part of @Tp(Sigma^{3,2}(10))@ divided into 10 pieces+--+-- The task can be parallezied using the @-B@ and @-b@ options+--++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, PackageImports, PatternGuards #-}+module Main where++--------------------------------------------------------------------------------++import Data.Char+import Data.List+import Data.Ratio+import Data.Monoid++import Control.Monad+import Control.Applicative+import Control.Concurrent+import Control.Concurrent.MVar++import System.Environment+import System.Mem+import System.IO+import System.Exit++import "time" Data.Time.Clock.POSIX++import Math.Combinat.Numbers.Primes++import Math.FreeModule.Symbol+import Math.FreeModule.SortedList+import Math.FreeModule.PrettyPrint+import Math.FreeModule.PP+-- import FreeModule.Parser++import Math.Algebra.ModP+import Math.Algebra.Schur++import Math.ThomPoly.Shared+import Math.ThomPoly.SigmaI   as SigmaI+import Math.ThomPoly.SigmaIJ  as SigmaIJ+import Math.ThomPoly.Formulae as Formulae++import Options.Applicative+    +--------------------------------------------------------------------------------++data Config = Config+  { _problem   :: !AnyProblem  +  , _tgtDim    :: !(Maybe Int)+  , _ring      :: !Ring +  , _outFile   :: !(Maybe FilePath)+  , _batch     :: !(Maybe Batch)+  , _printStat :: !Bool+  , _dry       :: !Bool+  , _timeout   :: !(Maybe Int)+  }+  deriving Show++--------------------------------------------------------------------------------+  +run :: Config -> IO ()+run config = do+  -- print config+  void $ mbTimeout (_timeout config) $ do++    let problem = _problem config+        batch   = maybe defaultBatch id (_batch config)+        ring    = selectRing (_ring config)++    when (_printStat config) $ do+      print $ case problem of+        PI  si  -> calcStats si +        PIJ sij -> calcStats sij++    let fname = case _outFile config of+          Just fname -> fname+          Nothing    -> case problem of+            PI  si  -> fullFName batch si +            PIJ sij -> fullFName batch sij        ++    let answer = case problem of+          PI  si  -> solveAny ring batch si +          PIJ sij -> solveAny ring batch sij++    let text   = pretty answer++    answer `seq` do+      unless (_dry config) $ writeFile fname text++--------------------------------------------------------------------------------+-- * configuration++data AnyProblem +  = PI  !SigmaI+  | PIJ !SigmaIJ+  | PI1 !SigmaI1            -- ^ we have an explicit formula for @Sigma^{i,1}@+  deriving Show++-- | Coefficient ring +data Ring+  = Integers+  | Rationals+  | HardCodedZp             -- ^ temporary+  | PrimeField !Integer+  | SpecPrime  !Int         -- ^ special primes just below @2^k@ for @k=7,15,31,63@+  deriving Show++selectRing :: Ring -> AnyRing+selectRing r = case r of+  Integers    -> ringZZ+  Rationals   -> ringQQ+  HardCodedZp -> ringZp++-- | Primes close to the bounds of (signed) machine words.+specPrimes :: [(Int,Integer)]+specPrimes = +  [ ( 7   , 2^7   -  1 )+  , ( 15  , 2^15  - 19 )+  , ( 31  , 2^31  -  1 )+  , ( 63  , 2^63  - 25 )+--  , ( 127 , 2^127 -  1 )+--  , ( 255 , 2^255 - 19 )+  ]++--------------------------------------------------------------------------------++maybeRead :: Read a => String -> Maybe a+maybeRead s = case reads s of +  [(x,"")] -> Just x+  _        -> Nothing++parseRing :: String -> Either String Ring+parseRing str0 +  | str `elem` ["zz","integer" ,"integers"  ]   =  Right Integers+  | str `elem` ["qq","rational","rationals" ]   =  Right Rationals+  | str `elem` ["zp","primefield" ]             =  Right HardCodedZp+--  | take 2+  where+    str = map toLower str0 ++--------------------------------------------------------------------------------++class Validate a where+  isValid :: a -> Maybe String+  +instance Validate Batch where+  isValid (Batch a b) +    | b < 1           =  Just "the number of batches B should be at least 1"+    | a < 1 || a > b  =  Just "batch index b should be between 1 and B"+    | otherwise       =  Nothing++instance Validate Ring where+  isValid r = case r of+    PrimeField p -> if isProbablyPrime p +      then Nothing +      else Just "order of the finite field should be a prime"+    SpecPrime  q -> case lookup q specPrimes of +      Nothing     -> Just "unimplemented special prime field (BITS should be one of 7, 15, 31 or 63)"+      _           -> Nothing+    _ -> Nothing       ++instance Validate AnyProblem where++  isValid problem = case problem of+    PI (SigmaI i n)+      | i < 1     -> Just "the index I should be at least 1"+      | n < 1     -> Just "the source dimension N should be at least 1"+      | otherwise -> Nothing+    PIJ (SigmaIJ i j n)+      | i < 1          -> Just "the index I should be at least 1"+      | j < 1 || j > i -> Just "the index J should be between 1 and I"+      | n < 1          -> Just "the source dimension N should be at least 1"+      | otherwise      -> Nothing+    PI1 (SigmaI1 i n)+      | i < 1     -> Just "the index I should be at least 1"+      | n < 1     -> Just "the source dimension N should be at least 1"+      | otherwise -> Nothing++--------------------------------------------------------------------------------+-- * option parsing++configOpt :: Parser Config+configOpt = Config +  <$> problemOpt +  <*> mOpt +  <*> ringNameOpt+  <*> outOpt +  <*> batchOpt +  <*> statFlag +  <*> dryFlag +  <*> timeoutOpt++problemOpt :: Parser AnyProblem+problemOpt = f <$> iOpt <*> jOpt <*> nOpt where+  f i mbj n = case mbj of+    Nothing -> PI $ SigmaI i n +    Just j  -> case j of +      0 -> PI  $ SigmaI  i   n+      _ -> PIJ $ SigmaIJ i j n++batchOpt :: Parser (Maybe Batch)+batchOpt = f <$> whichBatchOpt <*> nbatchOpt where+  f a b +    | a >= 1 && a <= b  =  if b > 1 then Just (Batch a b) else Nothing+    | otherwise         =  error "the batch index should be between 1 and B"++ringNameOpt :: Parser Ring+ringNameOpt = option (eitherReader parseRing)+  (  long    "ring"+  <> short   'r'+  <> metavar "R"+  <> value   Rationals+  <> help    "The coefficient ring (or field) R we compute in, for example a prime field" +  <> completeWith [ "Integers" , "Rationals" , "ZZ" , "QQ"+                  , "PrimeField" , "Zp"+--                  , "Zp7bit" , "Zp15bit" , "Zp31bit" , "Zp63bit"+                  ]+  <> showDefault+  )++timeoutOpt :: Parser (Maybe Int)+timeoutOpt = option (Just <$> auto)+  (  long    "timeout"+  <> short   't'+  <> metavar "TIMEOUT"+  <> value   Nothing+  <> help    "Timeout (specified in minutes)" +  )++primeOpt :: Parser Integer+primeOpt = option auto+  (  long    "prime"+  <> short   'p'+  <> metavar "P"+  <> value   1000000007+  <> help    "The order of the prime field" +  )++bitsOpt :: Parser Int+bitsOpt = option auto+  (  long    "bits"+  <> short   'q'+  <> metavar "BITS"+  <> value   63+  <> help    "Number of bits in the order of a special prime fields" +  )++statFlag :: Parser Bool+statFlag = switch+  (  long  "stats"+  <> short 's'+  <> help  "print \"statistics\" (codimension, algebraic multiplicity)" +  )++dryFlag :: Parser Bool+dryFlag = switch+  (  long  "dry"+  <> help  "do not write the result into a file" +  )++outOpt :: Parser (Maybe FilePath)+outOpt = option (Just <$> str)+  (  long    "output"+  <> short   'o'+  <> metavar "FILE"+  <> value   Nothing+  <> help    "Write output to FILE (use --dry to skip)" +  )++nbatchOpt :: Parser Int+nbatchOpt = option auto+  (  long    "nbatches"+  <> short   'B'+  <> metavar "B"+  <> value   1+  <> help    "number of batches" +  )++whichBatchOpt :: Parser Int+whichBatchOpt = option auto+  (  long "batch"+  <> short   'b'+  <> metavar "b"+  <> value   1+  <> help    "which batch to run (from 1 to B)" +  )++iOpt :: Parser Int+iOpt = option auto+  (  short   'i'+  <> metavar "I"+  <> help    "first Thom-Boardman index (I)"+  <> noArgError (ErrorMsg "specifying I is mandatory")+  )++jOpt :: Parser (Maybe Int)+jOpt = option (Just <$> auto)+  (  short   'j'+  <> metavar "J"+  <> value   Nothing+  <> help    "second Thom-Boardman index (J)"+  )++nOpt :: Parser Int+nOpt = option auto+  (  short   'n'+  <> metavar "N"+  <> help    "source dimension (N)"+  )++mOpt :: Parser (Maybe Int)+mOpt = option (Just <$> auto)+  (  short   'm'+  <> metavar "N"+  <> value   Nothing+  <> help    "target dimension (M, optional)"+  <> hidden+  )++--------------------------------------------------------------------------------++main :: IO ()+main = execParser opts >>= run where+  opts = info (helper <*> configOpt)+    (  fullDesc+    <> progDesc shortDesc+    <> header   longDesc+    )+  shortDesc = "Thom polynomials of second order Thom-Boardman singularities"+  longDesc  = "A program computing Thom polynomials of second order Thom-Boardman singularities"++--------------------------------------------------------------------------------+-- * timeout++-- | argument: number of minutes+mbTimeout :: Maybe Int -> IO a -> IO (Maybe a)+mbTimeout mb action = +  case mb of +    Nothing      -> Just <$> action+    Just minutes -> do+      mv <- newEmptyMVar +      t0 <- getPOSIXTime+      threadid <- forkIO $ do+        y <- action +        putMVar mv $! y+      wait mv t0 minutes threadid++  where      +    wait mv t0 minutes threadid = do+      let seconds = minutes * 60+      let loop = do+            threadDelay 1000000     -- wait 1 sec+            mb <- tryTakeMVar mv+            case mb of+              Just y  -> return $ Just y+              Nothing -> do+                t <- getPOSIXTime+                if t - t0 < fromIntegral seconds+                  then loop+                  else do+                    putStrLn $ "timeout after " ++ show minutes ++ " minutes"+                    killThread threadid+                    return Nothing+      loop+      +--------------------------------------------------------------------------------+
sigma-ij.cabal view
@@ -1,6 +1,6 @@  Name:                sigma-ij-Version:             0.2+Version:             0.2.0.1 Synopsis:            Thom polynomials of second order Thom-Boardman singularities Description:         A program to compute Thom polynomials of second order Thom-Boardman                       singularities, using the localization method described in the@@ -27,16 +27,11 @@  Executable sigma-ij -  hs-source-dirs:      src+  hs-source-dirs:      main   main-is:             sigmaij.hs -  Build-Depends:       base >= 4 && < 5, array >= 0.5, containers >= 0.5, random,-                       time, parsec2, optparse-applicative, -                       combinat >= 0.2.8--  -- cabal gets confused if the executable is in the same source tree...-  c-sources:           src/cbits/c_det.c    -  cc-options:          -std=c99 +  Build-Depends:       base >= 4 && < 5, time, optparse-applicative,+                       sigma-ij >= 0.2.0.1, combinat >= 0.2.8    Default-Language:    Haskell2010 
− src/sigmaij.hs
@@ -1,377 +0,0 @@---- | Calculates Thom polynomial of Sigma^{ij} with localization --- and the substituion trick------ Some example usages:------ > sigma-ij -h                           # help--- > sigma-ij -i3 -j2 -n7                  # compute @Tp(Sigma^{3,2}(7))@ in the default ring--- > sigma-ij -i3 -j2 -n7  -rZp            # compute @Tp(Sigma^{3,2}(7))@ in the hard-coded prime field--- > sigma-ij -i3 -j2 -n10 -rZp -b3 -B10   # compute the 3rd part of @Tp(Sigma^{3,2}(10))@ divided into 10 pieces------ The task can be parallezied using the @-B@ and @-b@ options-----{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, PackageImports, PatternGuards #-}-module Main where------------------------------------------------------------------------------------import Data.Char-import Data.List-import Data.Ratio-import Data.Monoid--import Control.Monad-import Control.Applicative-import Control.Concurrent-import Control.Concurrent.MVar--import System.Environment-import System.Mem-import System.IO-import System.Exit--import "time" Data.Time.Clock.POSIX--import Math.Combinat.Numbers.Primes--import Math.FreeModule.Symbol-import Math.FreeModule.SortedList-import Math.FreeModule.PrettyPrint-import Math.FreeModule.PP--- import FreeModule.Parser--import Math.Algebra.ModP-import Math.Algebra.Schur--import Math.ThomPoly.Shared-import Math.ThomPoly.SigmaI   as SigmaI-import Math.ThomPoly.SigmaIJ  as SigmaIJ-import Math.ThomPoly.Formulae as Formulae--import Options.Applicative-    -----------------------------------------------------------------------------------data Config = Config-  { _problem   :: !AnyProblem  -  , _tgtDim    :: !(Maybe Int)-  , _ring      :: !Ring -  , _outFile   :: !(Maybe FilePath)-  , _batch     :: !(Maybe Batch)-  , _printStat :: !Bool-  , _dry       :: !Bool-  , _timeout   :: !(Maybe Int)-  }-  deriving Show-----------------------------------------------------------------------------------  -run :: Config -> IO ()-run config = do-  -- print config-  void $ mbTimeout (_timeout config) $ do--    let problem = _problem config-        batch   = maybe defaultBatch id (_batch config)-        ring    = selectRing (_ring config)--    when (_printStat config) $ do-      print $ case problem of-        PI  si  -> calcStats si -        PIJ sij -> calcStats sij--    let fname = case _outFile config of-          Just fname -> fname-          Nothing    -> case problem of-            PI  si  -> fullFName batch si -            PIJ sij -> fullFName batch sij        --    let answer = case problem of-          PI  si  -> solveAny ring batch si -          PIJ sij -> solveAny ring batch sij--    let text   = pretty answer--    answer `seq` do-      unless (_dry config) $ writeFile fname text------------------------------------------------------------------------------------- * configuration--data AnyProblem -  = PI  !SigmaI-  | PIJ !SigmaIJ-  | PI1 !SigmaI1            -- ^ we have an explicit formula for @Sigma^{i,1}@-  deriving Show---- | Coefficient ring -data Ring-  = Integers-  | Rationals-  | HardCodedZp             -- ^ temporary-  | PrimeField !Integer-  | SpecPrime  !Int         -- ^ special primes just below @2^k@ for @k=7,15,31,63@-  deriving Show--selectRing :: Ring -> AnyRing-selectRing r = case r of-  Integers    -> ringZZ-  Rationals   -> ringQQ-  HardCodedZp -> ringZp---- | Primes close to the bounds of (signed) machine words.-specPrimes :: [(Int,Integer)]-specPrimes = -  [ ( 7   , 2^7   -  1 )-  , ( 15  , 2^15  - 19 )-  , ( 31  , 2^31  -  1 )-  , ( 63  , 2^63  - 25 )---  , ( 127 , 2^127 -  1 )---  , ( 255 , 2^255 - 19 )-  ]------------------------------------------------------------------------------------maybeRead :: Read a => String -> Maybe a-maybeRead s = case reads s of -  [(x,"")] -> Just x-  _        -> Nothing--parseRing :: String -> Either String Ring-parseRing str0 -  | str `elem` ["zz","integer" ,"integers"  ]   =  Right Integers-  | str `elem` ["qq","rational","rationals" ]   =  Right Rationals-  | str `elem` ["zp","primefield" ]             =  Right HardCodedZp---  | take 2-  where-    str = map toLower str0 ------------------------------------------------------------------------------------class Validate a where-  isValid :: a -> Maybe String-  -instance Validate Batch where-  isValid (Batch a b) -    | b < 1           =  Just "the number of batches B should be at least 1"-    | a < 1 || a > b  =  Just "batch index b should be between 1 and B"-    | otherwise       =  Nothing--instance Validate Ring where-  isValid r = case r of-    PrimeField p -> if isProbablyPrime p -      then Nothing -      else Just "order of the finite field should be a prime"-    SpecPrime  q -> case lookup q specPrimes of -      Nothing     -> Just "unimplemented special prime field (BITS should be one of 7, 15, 31 or 63)"-      _           -> Nothing-    _ -> Nothing       --instance Validate AnyProblem where--  isValid problem = case problem of-    PI (SigmaI i n)-      | i < 1     -> Just "the index I should be at least 1"-      | n < 1     -> Just "the source dimension N should be at least 1"-      | otherwise -> Nothing-    PIJ (SigmaIJ i j n)-      | i < 1          -> Just "the index I should be at least 1"-      | j < 1 || j > i -> Just "the index J should be between 1 and I"-      | n < 1          -> Just "the source dimension N should be at least 1"-      | otherwise      -> Nothing-    PI1 (SigmaI1 i n)-      | i < 1     -> Just "the index I should be at least 1"-      | n < 1     -> Just "the source dimension N should be at least 1"-      | otherwise -> Nothing------------------------------------------------------------------------------------- * option parsing--configOpt :: Parser Config-configOpt = Config -  <$> problemOpt -  <*> mOpt -  <*> ringNameOpt-  <*> outOpt -  <*> batchOpt -  <*> statFlag -  <*> dryFlag -  <*> timeoutOpt--problemOpt :: Parser AnyProblem-problemOpt = f <$> iOpt <*> jOpt <*> nOpt where-  f i mbj n = case mbj of-    Nothing -> PI $ SigmaI i n -    Just j  -> case j of -      0 -> PI  $ SigmaI  i   n-      _ -> PIJ $ SigmaIJ i j n--batchOpt :: Parser (Maybe Batch)-batchOpt = f <$> whichBatchOpt <*> nbatchOpt where-  f a b -    | a >= 1 && a <= b  =  if b > 1 then Just (Batch a b) else Nothing-    | otherwise         =  error "the batch index should be between 1 and B"--ringNameOpt :: Parser Ring-ringNameOpt = option (eitherReader parseRing)-  (  long    "ring"-  <> short   'r'-  <> metavar "R"-  <> value   Rationals-  <> help    "The coefficient ring (or field) R we compute in, for example a prime field" -  <> completeWith [ "Integers" , "Rationals" , "ZZ" , "QQ"-                  , "PrimeField" , "Zp"---                  , "Zp7bit" , "Zp15bit" , "Zp31bit" , "Zp63bit"-                  ]-  <> showDefault-  )--timeoutOpt :: Parser (Maybe Int)-timeoutOpt = option (Just <$> auto)-  (  long    "timeout"-  <> short   't'-  <> metavar "TIMEOUT"-  <> value   Nothing-  <> help    "Timeout (specified in minutes)" -  )--primeOpt :: Parser Integer-primeOpt = option auto-  (  long    "prime"-  <> short   'p'-  <> metavar "P"-  <> value   1000000007-  <> help    "The order of the prime field" -  )--bitsOpt :: Parser Int-bitsOpt = option auto-  (  long    "bits"-  <> short   'q'-  <> metavar "BITS"-  <> value   63-  <> help    "Number of bits in the order of a special prime fields" -  )--statFlag :: Parser Bool-statFlag = switch-  (  long  "stats"-  <> short 's'-  <> help  "print \"statistics\" (codimension, algebraic multiplicity)" -  )--dryFlag :: Parser Bool-dryFlag = switch-  (  long  "dry"-  <> help  "do not write the result into a file" -  )--outOpt :: Parser (Maybe FilePath)-outOpt = option (Just <$> str)-  (  long    "output"-  <> short   'o'-  <> metavar "FILE"-  <> value   Nothing-  <> help    "Write output to FILE (use --dry to skip)" -  )--nbatchOpt :: Parser Int-nbatchOpt = option auto-  (  long    "nbatches"-  <> short   'B'-  <> metavar "B"-  <> value   1-  <> help    "number of batches" -  )--whichBatchOpt :: Parser Int-whichBatchOpt = option auto-  (  long "batch"-  <> short   'b'-  <> metavar "b"-  <> value   1-  <> help    "which batch to run (from 1 to B)" -  )--iOpt :: Parser Int-iOpt = option auto-  (  short   'i'-  <> metavar "I"-  <> help    "first Thom-Boardman index (I)"-  <> noArgError (ErrorMsg "specifying I is mandatory")-  )--jOpt :: Parser (Maybe Int)-jOpt = option (Just <$> auto)-  (  short   'j'-  <> metavar "J"-  <> value   Nothing-  <> help    "second Thom-Boardman index (J)"-  )--nOpt :: Parser Int-nOpt = option auto-  (  short   'n'-  <> metavar "N"-  <> help    "source dimension (N)"-  )--mOpt :: Parser (Maybe Int)-mOpt = option (Just <$> auto)-  (  short   'm'-  <> metavar "N"-  <> value   Nothing-  <> help    "target dimension (M, optional)"-  <> hidden-  )------------------------------------------------------------------------------------main :: IO ()-main = execParser opts >>= run where-  opts = info (helper <*> configOpt)-    (  fullDesc-    <> progDesc shortDesc-    <> header   longDesc-    )-  shortDesc = "Thom polynomials of second order Thom-Boardman singularities"-  longDesc  = "A program computing Thom polynomials of second order Thom-Boardman singularities"------------------------------------------------------------------------------------- * timeout---- | argument: number of minutes-mbTimeout :: Maybe Int -> IO a -> IO (Maybe a)-mbTimeout mb action = -  case mb of -    Nothing      -> Just <$> action-    Just minutes -> do-      mv <- newEmptyMVar -      t0 <- getPOSIXTime-      threadid <- forkIO $ do-        y <- action -        putMVar mv $! y-      wait mv t0 minutes threadid--  where      -    wait mv t0 minutes threadid = do-      let seconds = minutes * 60-      let loop = do-            threadDelay 1000000     -- wait 1 sec-            mb <- tryTakeMVar mv-            case mb of-              Just y  -> return $ Just y-              Nothing -> do-                t <- getPOSIXTime-                if t - t0 < fromIntegral seconds-                  then loop-                  else do-                    putStrLn $ "timeout after " ++ show minutes ++ " minutes"-                    killThread threadid-                    return Nothing-      loop-      ----------------------------------------------------------------------------------