loopy (empty) → 0.0
raw patch · 12 files changed
+2306/−0 lines, 12 filesdep +GoogleChartdep +basedep +cmdargssetup-changed
Dependencies added: GoogleChart, base, cmdargs, containers, directory, filepath, hmatrix, process, random
Files
- Avalon.hs +279/−0
- DotFile.hs +172/−0
- FeedingRates.hs +100/−0
- Graph.hs +182/−0
- LICENSE +30/−0
- LoopProp.hs +162/−0
- Main.hs +510/−0
- Mlw.hs +321/−0
- Setup.hs +2/−0
- Stab.hs +336/−0
- Types.hs +181/−0
- loopy.cabal +31/−0
+ Avalon.hs view
@@ -0,0 +1,279 @@+{-# LANGUAGE ScopedTypeVariables #-} +module Avalon where + + +import System.Environment +import System.Random +import Control.Monad +import System.Directory +import System.Cmd +import Data.Char +import Data.Maybe +import Data.List +import Numeric +import System.FilePath +import Numeric.LinearAlgebra.LAPACK +import Data.Packed.Matrix hiding (Matrix) +import Data.Packed.Vector +import Data.Complex +import Types +import FeedingRates +import Stab + + +--not finished yet +changeChemoProp :: Foodweb -> Double -> Double -> Foodweb --first is frondose, second is rangeomorph prop osmo +changeChemoProp (Foodweb names info wij cm) cf cr = Foodweb names (changeInfo info) (changeWij wij) cm + where + changeInfo (fl:fm:fs:rl:rm:rs:rest) = frond fl: frond fm : frond fs : rangeo rl: rangeo rm: rangeo rs: rest + where + frond x = x{aeff =(cf*0.3+(1-cf)*0.8) , peff =(cf*0.8+(1-cf)*0.3) } + rangeo x = x{aeff =(cr*0.3+(1-cr)*0.8) , peff =(cr*0.8+(1-cr)*0.3) } + + changeWij wij = mapMatrixInd f wij + where + f i j x|i==11 && j >= 1 && j <=3 =cf*100 -- |i==12 && j >= 1 && j <=3 =cr*100 + |i==11 && j >= 4 && 6 <=3 =cf*100 -- |i==12 && j >= 1 && j <=3 =cr*100 + |otherwise = x + +createChemoRange :: FilePath -> Int-> IO () +createChemoRange fp n = do + fw <- readFoodweb fp + createDirectoryIfMissing True (dropExtension fp ++ "ChemoRange") + sequence_ [writeFoodweb file $ changeChemoProp fw x y + | x <- [0,1/(fromIntegral n)..1], y <- [0,1/(fromIntegral n)..1], let file = dropExtension fp ++ "ChemoRange"</> dropExtension fp ++ "_" ++ show3dp x ++ "_" ++ show3dp y ++ ".fw"] + + +createOsmoRange :: FilePath -> Int-> IO () +createOsmoRange fp n = do + fw <- readFoodweb fp + createDirectoryIfMissing True (dropExtension fp ++ "OsmoRange") + sequence_ [writeFoodweb file $ changeFeedingStrategy fw x y + | x <- [0,1/(fromIntegral n)..1], y <- [0,1/(fromIntegral n)..1], let file = dropExtension fp ++ "OsmoRange"</> dropExtension fp ++ "_" ++ show3dp x ++ "_" ++ show3dp y ++ ".fw"] + +createOsmoRange2 :: FilePath -> Int-> IO () +createOsmoRange2 fp n = do + fw <- readFoodweb fp + createDirectoryIfMissing True (dropExtension fp ++ "OsmoRange2") + sequence_ [writeFoodweb file $ changeFS2 fw x y + | x <- [0,1/(fromIntegral n)..1], y <- [0,1/(fromIntegral n)..1], let file = dropExtension fp ++ "OsmoRange2"</> dropExtension fp ++ "_" ++ show3dp x ++ "_" ++ show3dp y ++ ".fw"] + + +{- +createFeedingRange :: FilePath -> [Int]-> IO () +createFeedingRange fp = do + fw <- return$ fileToFoodweb fp + --createDirectoryIfMissing True (dropExtension fp ++ "SpeciesRange") + let file = [dropExtension fp ++ "SpeciesRange"</> dropExtension fp ++ ".fw"] + writeFoodweb file $ changeFeeding fw [2,2,2,3,3,3] + return() +-} +changeFeeding :: Foodweb -> [Int] -> Foodweb +changeFeeding (Foodweb names info wij cm) fts = Foodweb names (changeInfoSp info fts) (changeWijSp wij fts ((length names))) cm -- note that cm will be wrong + +generateFeedingTypes :: Int -> [[Int]] +generateFeedingTypes n = undefined --[k|k<-(replicate n i) && i<-[1..3] ] + +changePlanktonicBiomass :: Foodweb -> [Double] -> Foodweb +changePlanktonicBiomass (Foodweb names info wij cm) bm = Foodweb names (changeBiomass info) wij cm + where + changeBiomass (a:b:c:d:e:f:p:hb:ppp:rest) = a:b:c:d:e:f:pro p:het hb:prim ppp:rest + where + pro x = x{biomass =bm!!1} + het x =x{biomass =bm!!2} + prim x =x{biomass =bm!!3} +----update to include chemo stuff +changeInfoSp :: [Info] -> [Int] -> [Info] +changeInfoSp is [] = basicInfo is +changeInfoSp (i:is) (ft:fts) = [change1 i ft] ++ changeInfoSp is fts + where + change1 x t = if t==1 then x{aeff =0.3 , peff =0.8 } else + if t==2 then x{aeff = 0.8 , peff = 0.11 } else + x{aeff =1 , peff =0.3 } + +----update to include chemo stuff +changeWijSp :: Matrix Double -> [Int] -> Int -> Matrix Double +changeWijSp ws [] det = ws +changeWijSp wij fts det= mapMatrixInd f wij -- det is detrius row, found using getDetRow + where + n = 5 + length fts + stPl = length fts + enPl = length fts + 3 + chemo = [i| (i,j)<- zip [1..] fts, j ==1] + sf = [i| (i,j)<- zip [1..] fts, j ==2] + osmo = [i| (i,j)<- zip [1..] fts, j ==3] + f i j x|i >stPl && i <=enPl && j `elem` sf =33.33 + |i ==det && j `elem` osmo =100 + |otherwise = x + + +--returns just a list of the last 5 info +basicInfo :: [Info] -> [Info] +basicInfo (i:is) = if (length (i:is))> 5 then basicInfo is else (i:is) +{----- + +-} + + +outputAlpha :: FilePath -> IO () +outputAlpha fp = do + fw1<- readFoodweb fp + let fw = checkFoodweb fw1 + let alpha = getCommat fw + writeFile (dropExtension fp ++".alpha") (outputMatrix $ matrixDoubleToString alpha) + return() + -- print $ alphaijElem info wij res 14 0 + + +checkFoodweb :: Foodweb -> Foodweb +checkFoodweb fw = deleteSps fw b0 + where b0 = [i | (i,j) <- zip [0..] (info fw), biomass j == 0]--if biomass ==0 output col number + + +checkFoodweb2 :: Foodweb->Foodweb +checkFoodweb2 fw = fw{cm=getCommat fw} + + +getDetRow :: Foodweb -> [Int] +getDetRow fw = [i | (i,j) <- zip [0..] (info fw), deathRate j==0] +{-----} + + +changeInfo :: Foodweb -> Maybe Double ->Maybe Double -> Maybe [Double] -> Maybe [Double] -> Maybe Double-> Foodweb +changeInfo (Foodweb names info wij cm) marealCoverage mrFRatio msizeClassesF msizeClassesR mareaBiomassRatio = Foodweb names (change info) wij cm + where change = if mrFRatio==Just 0 then changeF else + if mrFRatio==Just 1 then changeR else + changeInfo + where + changeInfo (fl:fm:fs:rl:rm:rs:rest) = frond fl 0: frond fm 1: frond fs 2: rangeo rl 0: rangeo rm 1: rangeo rs 2: rest--numbers refer to sizeclasses + changeF (fl:fm:fs:rest) = frond fl 0: frond fm 1: frond fs 2: rest--numbers refer to sizeclasses + changeR (rl:rm:rs:rest) = rangeo rl 0: rangeo rm 1: rangeo rs 2: rest--numbers refer to sizeclasses + frond x y= x{biomass = (arealCoverage*(1-rFRatio)*(sizeClassesF!!y)*areaBiomassRatio)} + rangeo x y= x{biomass = (arealCoverage*rFRatio*(sizeClassesR!!y)*areaBiomassRatio)} + arealCoverage = fromMaybe 0.1 marealCoverage + rFRatio = fromMaybe 0.5 mrFRatio + sizeClassesR = fromMaybe [0.95,0.045,0.005] msizeClassesR -- need to extend to rangeo and fronds + sizeClassesF = fromMaybe [0.9,0.09,0.01] msizeClassesF -- need to extend to rangeo and fronds + areaBiomassRatio = fromMaybe 1000 mareaBiomassRatio +changeFS2 :: Foodweb -> Double -> Double -> Foodweb --first is frondose, second is rangeomorph prop osmo +changeFS2 (Foodweb names info wij cm) cf cr = Foodweb names (changeInfo info) (changeWij wij) cm + where + changeInfo (sp1:sp2:rest) = species1 sp1: species2 sp2 : rest + where + species1 x = x{aeff =(cf*1+(1-cf)*0.8) , peff =(cf*0.3+(1-cf)*0.11) } + species2 x = x{aeff =(cr*1+(1-cr)*0.8) , peff =(cr*0.3+(1-cr)*0.11) } + + changeWij wij = mapMatrixInd f wij + where + f i j x|i >=5 && i <=7 && j == 1 =(1-cf)*100 + |i >=5 && i <=7 && j == 2 =(1-cr)*100 + |i==9 && j ==1 =cf*100 -- |i==12 && j >= 1 && j <=3 =cr*100 + |i==9 && j ==2 =cr*100 -- |i==12 && j >= 1 && j <=3 =cr*100 + |otherwise = x + +--needs create a food web before going in +changeFeedingStrategy :: Foodweb -> Double -> Double -> Foodweb --first is frondose, second is rangeomorph prop osmo +changeFeedingStrategy (Foodweb names info wij cm) cf cr = Foodweb names (changeInfo info) (changeWij wij) cm + where + changeInfo (fl:fm:fs:rl:rm:rs:rest) = frond fl: frond fm : frond fs : rangeo rl: rangeo rm: rangeo rs: rest + where + frond x = x{aeff =(cf*1+(1-cf)*0.8) , peff =(cf*0.3+(1-cf)*0.11) } + rangeo x = x{aeff =(cr*1+(1-cr)*0.8) , peff =(cr*0.3+(1-cr)*0.11) } + + changeWij wij = mapMatrixInd f wij + where + f i j x|i >=7 && i <=9 && j >= 1 && j <=3 =(1-cf)*100 + |i >=7 && i <=9 && j >= 4 && j <=6 =(1-cr)*100 + |i==11 && j >= 1 && j <=3 =cf*100 -- |i==12 && j >= 1 && j <=3 =cr*100 + |i==11 && j >= 4 && j <=6 =cr*100 -- |i==12 && j >= 1 && j <=3 =cr*100 + |otherwise = x +---- +--- changing the biomass in some way +-- three things influence the biomasses: +--(1) total biota coverage of bedding plane - default 0.10 +--(2) prop of rangeo:frond - default 0.5 +--(3) prop of size classes in each tier - default [0.9,0.09,0.01] +--(4) the area to biomass ratio - default 1000 +--each size class is given by 1*2*3 + +changeBiomass :: Foodweb -> Maybe Double ->Maybe Double -> Maybe [Double] -> Maybe [Double] -> Maybe Double-> Foodweb +changeBiomass (Foodweb names info wij cm) marealCoverage mrFRatio msizeClassesF msizeClassesR mareaBiomassRatio = Foodweb names (changeInfo info) wij cm + where + changeInfo (fl:fm:fs:rl:rm:rs:rest) = frond fl 0: frond fm 1: frond fs 2: rangeo rl 0: rangeo rm 1: rangeo rs 2: rest--numbers refer to sizeclasses + frond x y= x{biomass = (arealCoverage*(1-rFRatio)*(sizeClassesF!!y)*areaBiomassRatio)} + rangeo x y= x{biomass = (arealCoverage*rFRatio*(sizeClassesR!!y)*areaBiomassRatio)} + arealCoverage = fromMaybe 0.1 marealCoverage + rFRatio = fromMaybe 0.5 mrFRatio + sizeClassesR = fromMaybe [0.95,0.045,0.005] msizeClassesR -- need to extend to rangeo and fronds + sizeClassesF = fromMaybe [0.9,0.09,0.01] msizeClassesF -- need to extend to rangeo and fronds + areaBiomassRatio = fromMaybe 1000 mareaBiomassRatio + +changeBiomass2 :: Foodweb -> Maybe Double ->Maybe Double -> Maybe [Double] -> Maybe [Double] -> Maybe Double-> Foodweb +changeBiomass2 (Foodweb names info wij cm) marealCoverage mrFRatio msizeClassesF msizeClassesR mareaBiomassRatio = changeInfo (changeFoodweb (Foodweb names info wij cm) mrFRatio) marealCoverage mrFRatio msizeClassesF msizeClassesR mareaBiomassRatio + +changeFoodweb :: Foodweb -> Maybe Double-> Foodweb +changeFoodweb fw r = if r == Just 0 --r is the ratio of rangeomorphs to fronds + then + deleteSp (deleteSp (deleteSp fw 3) 3) 3 + else + if r == Just 1 + then + deleteSp (deleteSp (deleteSp fw 0) 0) 0 + else + fw + + + +fileToFoodweb2 :: FilePath -> IO Foodweb +fileToFoodweb2 fp = do + y <- fileToFoodweb fp + return $ y{cm=getCommat y} -- FIXME: Reorder to remove this + +{- +-} +--takes away any zero biomass species NEED TO FINISH THIS OFF PROPERLY + + + + +--needs to be an ordered list +deleteSps :: Foodweb -> [Int] -> Foodweb +deleteSps fw [] = fw +deleteSps fw spNu = deleteSps (deleteSp fw (head spNu)) (reduceListNu (tail spNuO)) + where spNuO = spNu + +reduceListNu :: [Int] -> [Int] +reduceListNu [] = [] +reduceListNu (x:xs) = [x-1] ++ reduceListNu xs + +deleteSp :: Foodweb -> Int -> Foodweb +deleteSp (Foodweb names info wij cm) spNu = Foodweb (delListLine spNu names) (delListLine spNu info) (delListListLine spNu wij) cm + +delListLine :: Int->[a]->[a] +delListLine i xs =take i xs ++ drop (i+1) xs -- deleting line i from xs, with numbers starting at 0 + +delListListLine :: Int -> [[a]] -> [[a]] +delListListLine i xs = map (delListLine i) (delListLine i xs) + +createFrondRange :: FilePath -> Int -> Int -> IO () -- assumes fronds and rangeomorphs behave in the same way +createFrondRange fp n1 n2 = do + fw <- readFoodweb fp + createDirectoryIfMissing True (dropExtension fp ++ "FrondRange") + sequence_ [writeFoodweb file $ changeFeedingStrategy (changeBiomass2 fw Nothing (Just y) Nothing Nothing Nothing) x x + | x <- [0,1/(fromIntegral n1)..1], y <- [0,1/(fromIntegral n2)..1], let file = dropExtension fp ++ "FrondRange"</> dropExtension fp ++ "_" ++ show3dp y ++ "_" ++ show3dp x ++ ".fw"] + + + + + + + + + + + + + + + +
+ DotFile.hs view
@@ -0,0 +1,172 @@+module DotFile where + + +import System.Environment +import System.Random +import Control.Monad +import System.Directory +import System.Cmd +import Data.Char +import Data.Maybe +import Data.List +import Numeric +import System.FilePath +import Numeric.LinearAlgebra.LAPACK +import Data.Packed.Matrix hiding (Matrix) +import Data.Packed.Vector +import Data.Complex +import Types +import FeedingRates +import Stab + +-------------------------------------------------------------------------------- +--generating a dot file from a fw +funct :: [Int] -> Int -> [Int] +funct [] a = [] +funct m rw = if (head m) ==0 + then [0] ++ funct (tail m) rw + else [rw] ++ funct (tail m) rw + +wijSimple :: Matrix Double -> Matrix Int +wijSimple wij = map rowSimple wij + +rowSimple [] = [] +rowSimple rij = if (head rij) == 0 then [0] ++ rowSimple (tail rij) else [1] ++ rowSimple (tail rij) + +findPP :: Matrix Double -> [Int] +findPP [] = [] +findPP m = if (sum (head m)) == 0 then [1]++ findPP (tail m) else [0] ++ findPP (tail m) +--where m is transpose $ wijSimple $ wij fw + +--outputs the row numbers of the trophic level put in +ppRowNumbers :: [(Int,Int)] -> [Int] +ppRowNumbers [] = [] +ppRowNumbers zpp = if snd (head zpp) == 1 then [fst (head zpp)]++ ppRowNumbers (tail zpp) else ppRowNumbers (tail zpp) + --where zpp = zip [1::Int ..] (findPP wij) + +trophicElem :: Matrix Double -> Matrix Double -> Int -> Int -> Double +trophicElem z w i j | w!!i!!j ==0 = 0 + | otherwise = 1+sum0 [z!!m!!i | m<-[0..(length w -1)]] + + + +trophicMatrixZ :: Matrix Double -> Matrix Double +trophicMatrixZ w = z + where len = length w + z = [[trophicElem z w i j | j <- [0..len-1]]| i <-[0..len-1] ] + + + +drawFoodweb :: FilePath -> IO () +drawFoodweb fp = do + fw<-fileToFoodweb fp + let nms = createNodeNames fw + let ranks = createNodeTL fw + let edges = concat $ createEdges fw + let pp = createPPRank fw + let res = "digraph{\nrankdir=BT;\ngraph[label=" ++ show (dropExtension fp) ++ "];\n" ++ nms ++ "\n" ++ ranks ++ "\n" ++ pp ++" \n" ++ edges ++"}" -- error at ranks + writeFile (dropExtension fp ++ ".dot") res + + + +createEdges :: Foodweb -> [String] +createEdges fw = [if (wij fw)!!i!!j /= 0 then (show i) ++ "->" ++ (show j) ++ ";" else "" |i<-[0..(length (names fw) -1)], j<-[0..(length (names fw) -1)]] + + +createNodeTL :: Foodweb -> String +createNodeTL fw = concat $ map outputTL (nodeTL fw) + +createPPRank :: Foodweb -> String +createPPRank fw = concat $ ["{rank=source;"] ++[if tls!!i == 1 then show i ++ ";" else "" |i<-[0..(length (names fw) -1)]] ++ ["}"] + where tls=map (+1) $ map sum0 $transpose (trophicMatrixZ (wij fw)) + +outputTL :: (Integer, Double) -> String +outputTL tp = "{rank=" ++ show (snd tp) ++ ";" ++ show (fst tp) ++ ";}" + +nodeTL :: Foodweb -> [(Integer,Double)] +nodeTL fw = zip [0..] tls + where tls = map (+1) $ map sum0 $transpose (trophicMatrixZ (wij fw)) + +createNodeNames :: Foodweb -> String +createNodeNames fw = concat $ map nameLabel (zip [0..] (names fw)) + +nameLabel :: (Integer,String) -> String +nameLabel tp = show (fst tp) ++ "[label=" ++ (snd tp) ++ "];" + +writeTLMatrix :: FilePath -> IO () +writeTLMatrix fp = do + fw<-fileToFoodweb fp + let res = trophicMatrixZ (wij fw) + writeFile ( dropExtension fp ++ ".tl") (outputMatrix $matrixDoubleToString $ res) + +count0 :: [Double] -> Double +count0 [] = 0 +count0 (x:xs) = if x == 0 then count0 xs else 1+ count0 xs + +sum0 :: [Double] -> Double +sum0 x = if count0 x == 0 then 0 else sum(x)/count0 x + +-------------------------------------------------------------------------------- +--generating fig 2 from anje 1995 paper + +fig2 :: FilePath -> IO () +fig2 fp = do + fw1<-fileToFoodweb fp + let fw = putCM fw1 + f = fij (info fw) (wij fw) + a = cm fw + n = names fw + title = "Predator \t Fij \t aij \t aji \t Prey \n" + contents = concat$ [if (wij fw)!!i!!j /= 0 then n!!i ++ "\t" ++ show ( f!!i!!j) ++ "\t" ++ show (a!!i!!j) ++ "\t" ++ show( a!!j!!i) ++ "\t" ++ n!!j ++ "\n" else "" |i<-[0..(length (names fw) -1)], j<-[0..(length (names fw) -1)]] + res = title ++ contents + writeFile (dropExtension fp ++ ".fig2") res +-------------------------------------------------------------------------------- +--generating biomass ratios of mlw3 loop for a directory of files +{--} -- problem is generating them for a directory +getBiomassRatio :: FilePath -> Int -> IO () +getBiomassRatio fp ll1 = do + fw<-fileToFoodweb fp + --getLoopyInfo fp + x <-readLoopPropFile (dropExtension fp ++ ".log") + let ll = if ll1 >(1+(length x)) then 2 else ll1 + let loop = read (x!!(ll-2)!!8)::[Int] + tls = map snd $ nodeTL fw + i = info fw + b =map biomass i + zp =zip [tls!!(loop!!i) | i<-[0..(ll-1)]] [b!!(loop!!i) | i<-[0..(ll-1)]] + pmax = snd $ maximum zp + pmin = snd $ minimum zp + writeFile (dropExtension fp ++ ".br" ) (show (pmax/pmin)) + --putStr (show (pmax/pmin)) + return ( ) + + +writeBiomassRatios :: FilePath -> IO () +writeBiomassRatios fp= do + stabs<- readStabs fp + brs<-readBrs fp + let res = zip stabs brs + writeFile (dropExtension fp ++ ".stabBr") (tableTuple res) + +writeBrs :: FilePath -> Int -> IO () +writeBrs fp ll= do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".fw") allfiles + res<- sequence [getBiomassRatio (fp </> s) ll| s <- ss] + return() +{--} + +readBrs :: FilePath -> IO [Double] +readBrs fp = do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".br") allfiles + res<- sequence [readBr (fp </> s)| s <- ss] + stabs<-return $ res + return(stabs) + + +readBr :: FilePath -> IO Double -- checked ok +readBr x = do + s<-readFile $x--fp</>x + let y =read s :: Double + return(y)
+ FeedingRates.hs view
@@ -0,0 +1,100 @@+module FeedingRates where + +import System.Environment +import System.Random +import Control.Monad +import System.Directory +import System.Cmd +import Data.Char +import Data.Maybe +import Data.List +import Numeric +import Types +import System.FilePath + + +getCommat :: Foodweb -> Matrix Double +getCommat inputFW = functionToMatrix (alphaijElem inputFW res) (length (wij inputFW)) + + where + res = fij (info inputFW) (wij inputFW) + +{- writeFile (x ++ ".alpha.txt") $ outputMatrix $ map (map show) alphaij +-} +--splits a txt file with matrices in, into lines, then drops any extra spaces, then drops empty lists from start of lists +splitFile file = do + src <- readFile file + src <- return $ dropWhile null $ map (dropWhile isSpace) $ lines src + let (s1,_:rest) = break null src -- breaks the first matrix from the rest + (s2,_:s3) = break null $ dropWhile null rest -- separates out the last two matrices from each other + return $ map (unlines . map tabify) [s1,s2,s3] + + + +putCM :: Foodweb -> Foodweb +putCM fw = fw{cm=getCommat fw} + +tableTuple stabs1 = unlines $ [(show (map fst stabs1!!n)) ++ "\t" ++ (show $ map snd stabs1!!n)|n<-[0..(length stabs1)-1]] + + +------------change matrix type downwards +sumT :: [Info] -> Matrix Double -> Int -> Double +sumT info wij j= sum $ zipWith (*) (transpose wij !! j) (map biomass info) --checked + +fij :: [Info] -> Matrix Double -> Matrix Double +fij info wij = res + where + is = [0..length wij-1] + res = [[fijElem info wij res j i | i <- is] | j <- is] + +fijElem :: [Info] -> Matrix Double -> Matrix Double -> Int -> Int -> Double +fijElem info wij fij i j = + if sumT info wij j == 0 then 0 + else + ((wij!!i!!j) * biomass infoi * ((deathRate infoj * biomass infoj) + mi j)) / + (aeff infoj * peff infoj * sumT info wij j) + where + infoi = info!!i + infoj = info!!j + mi j = sum [fij !! j !! i | i <- [0..j-1]] + +--fij checked and is correct and working + + +alphaijElem :: Foodweb->Matrix Double-> Int -> Int -> Double +alphaijElem (Foodweb names info wij _) res i j + | i > j && i /= d = -resij/(biomassj) + | i < j && i /= d = aeffi*peffi*resji/biomassj --error $ show (aeffj,peffj,resij,biomassi) -- + | i == j && i /= d = -s*(deathRate infoi) -- will need to link up so we can insert real value of s in here + | i == j && i == d = -(sumTdd)/(biomassd) -- need to check correct + | i /= j && i == d = (deathRate infoj) + (1-aeffj)*(sumTdi1) + (sumTdi2)/(biomassj) - res!!d!!j/biomassj + + where + s = 1 + d = fromMaybe (-1) $ findIndex (== 0) (map deathRate info) + infoi = info!!i + infoj = info!!j + infod = info!!d + resij = res!!i!!j + resji = res!!j!!i + resdj = res!!d!!j + biomassd = biomass infod + biomassj = biomass infoj + biomassi = biomass infoi + peffi = peff infoi + peffj = peff infoj + aeffi = aeff infoi + aeffj = aeff infoj + sumTdd = sum [ (res!!d!!k)* ( (aeff (info!!k) ) * ( biomass (info!!k)))/(biomass (info!!k)) | k <-[0..length wij-1]] --need to check correct + sumTdi1 = sum [(res!!k!!j)/(biomass (info!!j))| k <-[0..length wij-1]] + sumTdi2 =sum [ (res!!j!!k)* ((1-(aeff (info!!k)))) | k <-[0..length wij-1] ] + +--perhaps create a sum over k function? + + + +{-wrap "table" . concatMap (wrap "tr" . concatMap (wrap "td" . unshow . showFFloat (Just 3))) + where wrap tag str = "<" ++ tag ++ ">" ++ str ++ "</" ++ tag ++ ">"-} + + +unshow f = f ""
+ Graph.hs view
@@ -0,0 +1,182 @@+--need to change matrix types +module Graph where + +import Types +import System.IO.Unsafe(unsafePerformIO) + +import Data.List +import Data.IntMap(IntMap) +import qualified Data.IntMap as IntMap +import Control.Monad +import Data.Bits + +type Graph = [(Int,Int)] --list of edges +type Graph2 = IntMap (IntMap Chunk, IntMap Chunk) +type Loop = [Int] +type Size = Int + +graphFromMatrix :: Int -> Matrix Double -> Graph +graphFromMatrix n m = [(i,j) | i <- [0..n-1], j <- [0..n-1], i /= j, let v = m!!i!!j, v /= 0] + + +-- Index is the lowest of src and dest link +-- min (src, dest) + +data Chunk = Nil + | Atom Int + | CrossL Chunk Int + | CrossR Int Chunk + | Cross Chunk Int Chunk + | Or Chunk Chunk + deriving Show + +cross Nil i Nil = Atom i +cross Nil i x = CrossR i x +cross x i Nil = CrossL x i +cross x i y = Cross x i y + + +loops :: Graph -> [Loop] +loops g = deleteNodes $ graph2 g + + +graph2 :: Graph -> Graph2 +graph2 xs = addLinks IntMap.empty [(a,Nil,b) | (a,b) <- xs] + + +addLinks :: Graph2 -> [(Int,Chunk,Int)] -> Graph2 +addLinks = foldl addLink + +addLink :: Graph2 -> (Int,Chunk,Int) -> Graph2 +addLink g (src,i,dest) | src < dest = IntMap.insertWith add src (IntMap.empty, IntMap.singleton dest i) g + where add _ (x,y) = (x, IntMap.insertWith Or dest i y) +addLink g (src,i,dest) | otherwise = IntMap.insertWith add dest (IntMap.singleton src i , IntMap.empty) g + where add _ (x,y) = (IntMap.insertWith Or src i x , y) + + +deleteNodes :: Graph2 -> [Loop] +deleteNodes g | IntMap.null g = [] + | otherwise = a ++ deleteNodes b + where (a,b) = deleteNode g + + +deleteNode :: Graph2 -> ([Loop], Graph2) +deleteNode g = (concatMap flatten [cross Nil i c | (i,c,_) <- lop], addLinks g2 new) + where + (lop,new) = partition (\(i,_,j) -> i == j) [(src,cross srci k desti,dest) | (src,srci) <- lhss, (dest,desti) <- rhss] + (lhss,rhss) = (IntMap.toList lhs, IntMap.toList rhs) + Just ((k,(lhs,rhs)),g2) = IntMap.minViewWithKey g + + + +flatten :: Chunk -> [Loop] +flatten x = map snd . f $ x + where + -- the bits in the Integer represent which indecies are present + f :: Chunk -> [(Integer, Loop)] + f Nil = [(0, [])] + f (Or a b) = f a ++ f b + f (Cross as b cs) = + [(aci, a ++ c) + | let ass = [(ai,a) | (ai,a) <- f as, not $ testBit ai b] + , (ciPre,cPre) <- f cs, let c = b : cPre + , not $ testBit ciPre b, let ci = setBit ciPre b + , (ai,a) <- ass + , ai .&. ci == 0, let aci = ai .|. ci] + + -- specialisations of Cross + f (Atom b) = [(bit b, [b])] + f (CrossR b cs) = + [(ci, c) + | (ciPre,cPre) <- f cs, let c = b : cPre + , not $ testBit ciPre b, let ci = setBit ciPre b] + f (CrossL as b) = + [(ai, a) + | (aiPre,aPre) <- f as, let a = aPre ++ [b] + , not $ testBit aiPre b, let ai = setBit aiPre b] +--removes any loops containing detritus +detLoops :: Int -> [Loop] -> [Loop] +detLoops detRow [] = [] +detLoops detRow ls = if (elem detRow (head ls)) then (detLoops detRow (tail ls)) else (ls ++ (detLoops detRow $ tail ls)) + +{- +detLoops2 :: Int -> [String] -> [String] +detLoops2 detRow [] = [] +detLoops2 detRow ls = if (elem detRow (l)) then (detLoops2 detRow (tail ls)) else (ls ++ (detLoops2 detRow $ tail ls)) + where l = read l2 :: Loop + l2 = (head ls)!!0 + + +flatten :: Chunk -> [Loop] +flatten = map snd . f + where + f :: Chunk -> [([Int], Loop)] + f Nil = [([], [])] + f (Or a b) = f a ++ f b + f (Cross as b cs) = + [(aci, a ++ c) + | let ass = filter (notElemInt b . fst) $ f as + , (ciPre,cPre) <- f cs, let c = b : cPre + , b `notElemInt` ciPre, let ci = insertInt b ciPre + , (ai,a) <- ass + , disjoint ai ci, let aci = merge ai ci] + + +notElemInt :: Int -> [Int] -> Bool +notElemInt x [] = True +notElemInt x (y:ys) = x /= y && notElemInt x ys + + +insertInt :: Int -> [Int] -> [Int] +insertInt x (y:ys) | x < y = x : y : ys + | otherwise = y : insertInt x ys +insertInt x [] = [x] + + +merge :: [Int] -> [Int] -> [Int] +merge (x:xs) (y:ys) + | x < y = x : merge xs (y:ys) + | otherwise = y : merge (x:xs) ys +merge xs [] = xs +merge [] ys = ys + + +disjoint :: [Int] -> [Int] -> Bool +disjoint xs [] = True +disjoint [] ys = True +disjoint (x:xs) (y:ys) = case compare x y of + EQ -> False + LT -> disjoint xs (y:ys) + GT -> disjoint (x:xs) ys +-} + + + + + +loopsOld :: Graph -> [Loop] +loopsOld g = f [(x,[],y) | (x,y) <- g] + where + f [] = out "done" [] + f xs = out (length res, length xs) $ res ++ f xs4 + where + res = [a:b | (a,b,c) <- loops] + kill = fst3 $ head xs + (from,xs2) = partition ((==) kill . fst3) xs + xs3 = concatMap add xs2 + (loops,xs4) = partition isLoop xs3 + + add (a,b,c) | c == kill = [(a,b ++ d:e,f) | (d,e,f) <- from, disjointed b e] + add x = [x] + + isLoop (a,b,c) = a == c + + +fst3 (a,b,c) = a + +disjointed x y = null $ x `intersect` y + + + +{-# NOINLINE out #-} +out x y = unsafePerformIO $ print x >> return y
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Emily Mitchell 2009-2011.+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 Neil Mitchell 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.
+ LoopProp.hs view
@@ -0,0 +1,162 @@+ +--need to change matrix types +-- trace (show (a,b,c)) +module LoopProp where + +import Types +import Graph + +import Data.Function +import Control.Arrow +import Data.List +import Data.Maybe +import Graphics.Google.Chart +import Data.Map(Map) +import qualified Data.Map as Map + +import Debug.Trace + +{- +Options: + +divide by d or not +showing mlw loop, option to choose any of the other loops to show, + +generate graph with r^2 value, or just output r^2 val +-} + + +-- A loop (propLoop) and properties about that loop +data Prop = Prop + {propLoop :: !Loop + ,propLength :: !Int + ,propWeight :: !Double + ,propPositive :: !Bool + ,propEven :: !Int + ,propNet :: !Double + } + deriving Show + +-- The properties on a set of loops +data PropSum = PropSum + {propsLength :: !Int + ,propsTotal :: !Int + ,propsAll :: !(Maybe Prop) -- need to check that it is just positive and negative, not All including net + ,propsPositive :: !(Maybe Prop) + ,propsNegative :: !(Maybe Prop) + ,propsEven :: !(Maybe Prop) + ,propsOmni:: !(Maybe Prop) + ,propsNet:: !(Maybe Prop) -- the addition of the forward and back loops + } + deriving Show + +---------------- + +propSum :: Prop -> PropSum +propSum p = PropSum (propLength p) 1 + (test True) (test $ propPositive p) (test $ not $ propPositive p) (test $ propEven p == 0) -- changed first bracket from just (test True) $ propEven p /= propLength p + (test $ propEven p == 2 - propLength p || propEven p ==propLength p - 2) (test True) + where test b = if b then Just p else Nothing + +propSumOutput :: Prop -> String +propSumOutput p = (show (propLoop p))++ "\t"++(show (propLength p)) ++ "\t"++ (show3dp (propWeight p) )++ "\t"++ (if (propPositive p) then "+" else "-" ) ++ "\t"++ (show (propEven p)) ++ "\t"++ (show3dp (propNet p)) ++ "\n" + +--adds two propSums together, taking the maximum values +propSumPlus :: PropSum -> PropSum -> PropSum +propSumPlus a b = PropSum (propsLength a) (propsTotal a + propsTotal b) + (both propsAll) (both propsPositive) (both propsNegative) (both propsEven) (both propsOmni) (bothr propsNet) + where + both f = join (f a) (f b) + join Nothing x = x + join x Nothing = x + join (Just a) (Just b) | propWeight a > propWeight b = Just a + | otherwise = Just b + bothr f = joinr (f a) (f b) + joinr Nothing x = x + joinr x Nothing = x + joinr (Just a) (Just b) | propNet a > propNet b = Just a + | otherwise = Just b + +--check that code written this week for signum is correct +-- | Finds the properties for a loop, m is the interaction matrix +classify :: Matrix Double -> Loop -> Prop +classify m loop = Prop loop n -- the loop and the length of the loop + (abs (pas / product ds) ** (1 / fromIntegral n)) -- weight of loop-- / product ds taken out dividing by diag for richard law things.... + (pas > 0) -- is the loop positive + (sum $ map (round . signum) as) -- the sum of the loop + (abs ( (abs (pas / product ds) ** (1 / fromIntegral n)) - (abs (pasr / product ds) ** (1 / fromIntegral n)) ))-- reverse weight of loop-- + where + pas = product as + as = map (\(x,y) -> m !! x !! y) $ links loop + ds = [m!!i!!i | i <- loop] + n = length loop + pasr = product asr + asr = map (\(x,y) -> m !! x !! y) $ links $ reverse loop + + + +-- sorted by the Int +separate :: [Prop] -> [(Int, [Prop])] +separate = map (fst . head &&& map snd) . groupBy ((==) `on` fst) . sortBy (compare `on` fst) . map (propLength &&& id) + + +summarize :: [Prop] -> [PropSum] +summarize = Map.elems . foldl' f Map.empty + where + f :: Map Int PropSum -> Prop -> Map Int PropSum + f mp x = Map.insertWith' propSumPlus (propLength x) (propSum x) mp + + +propSumRow :: PropSum -> [String] +propSumRow (PropSum a b c d e f h i) = [show a, show b, g c, g d, g e, g f, g h, g3 i, g2 c] + where g = maybe "NA" (show3dp . propWeight) + g2 (Just x) = show $ propLoop x + g3 = maybe "NA" (show3dp . propNet) + +---propSumRowOmni currently doesn't work 25th aug 10 +propSumRowOmni :: PropSum -> [String] +propSumRowOmni (PropSum a b c d e f h i) = [show a, show b, g c, g d, g e, g f, g h, g3 i, g2 h] + where g = maybe "NA" (show3dp . propWeight) + g2 = maybe "NA" (show . propLoop) + g3 = maybe "NA" (show3dp . propNet) + + +propSumRowChoice :: Int -> PropSum -> [String] +propSumRowChoice i1 (PropSum a b c d e f h i) = [show a, show b, g c, g d, g e, g f, g h, g3 i, g2 l1] + where g = maybe "NA" (show3dp . propWeight) + g3 = maybe "NA" (show3dp . propNet) + g2 = maybe "NA" (show . propLoop) + l1 = if i1 == 1 then c else + if i1 == 2 then d else + if i1 == 3 then e else + if i1 == 4 then f else + if i1 == 5 then h else i +{--} +-- [1,2,3] = [(1,2),(2,3),(3,1)] + +links :: [Int] -> [(Int,Int)] +links (x:xs) = f x xs + where f y [] = [(y,x)] + f y (z:zs) = (y,z) : f z zs + + +outputChart :: [PropSum] -> String +outputChart xs = chartURL $ + setSize 400 257 $ + setTitle "Maximum Loop Weights" $ + setData (encodeDataSimple $ map f [propsAll,propsPositive,propsNegative,propsEven]) $ + setDataColors ["000000","FF0000","0000FF","00000088"] $ + setAxisLabels [map (show . propsLength) xs, ["0",show3dp mx]] $ + setAxisTypes [AxisBottom, AxisLeft] $ + setLegend ["MLW","MLW+","MLW-","MLWe"] $ + newLineChart + where + mx = maximum $ map propWeight $ catMaybes $ concat [[a,b,c,d,e,f] | PropSum _ _ a b c d e f<- xs] + f g = [floor $ (i * 61) / mx | y <- map g xs, let i = maybe 0 propWeight y] +--type LoopPropFile = [PropSum] + +--convert :: [String] -> PropSum +--convert [a,b,c,d,e,f]=PropSum (read a) (read b) (read c) (read d) (read e) (read f) +{- + +-}
+ Main.hs view
@@ -0,0 +1,510 @@+{-# LANGUAGE DeriveDataTypeable, RecordWildCards #-} +--runhaskell Main.hs *.fw + + +module Main(main) where + +import System.Environment + +--import Matrix +import Graph +import LoopProp +import FeedingRates +import Types +import Stab +import DotFile +import Data.Char +import System.FilePath +import System.Directory +import Avalon +import Data.List +import Mlw +import Graphics.Plot +import Data.Maybe +import System.Console.CmdArgs + +--import LinearRegression +--import Modern +--import Numeric.LinearAlgebra + +data LoopyOpt = LoopyOpt + {filename :: String + ,loopOutput :: Bool + ,omni :: Bool + } deriving (Data,Typeable) + +loopyOpt = cmdArgsMode $ LoopyOpt + {filename = "check.fw" &= args &= typ "FILE/DIR" + ,loopOutput = False &= help "Output All loops found. Slows program" + ,omni = False &= help "Loops shown under MLW are Omnivorous Loops"} + &= summary "\nEmily's Loop Finder" + +--this is the command line directory +main:: IO () +main = do + o@LoopyOpt{..} <-cmdArgsRun loopyOpt + --[fp] <- getArgs + --x<-return fp + if (length (takeExtension filename)) == 0 then mainPath o else mainFile o -- + --mainPath fp# + -- if input alpha, don't stability +--if fw then output alpha + + + +mainPath :: LoopyOpt-> IO () +mainPath o@LoopyOpt{filename=fp,..} = do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".fw" ) allfiles + tt1 = filter (\xx -> takeExtension xx == ".alpha" ) allfiles + tt = myMerge tt1 ss + res<- sequence [mainFile o{filename=(fp </> s)} | s <- tt] --fp <\> creates a dir + --writeBrs fp 3 + --writeBiomassRatios fp + --putStr (show tt) + return() + +mainFile :: LoopyOpt -> IO () +mainFile LoopyOpt{filename=fp,..} = do + l<- returnComponents fp + if length l==1 then getLoopyInfoA loopOutput fp + else do + outputAlpha fp + if omni then do + getLoopyInfoOmni loopOutput fp + putStr "Only finding Omnivorous Loops\n" + else getLoopyInfo loopOutput fp + getStabInfo fp + --putStr ("stab and loopy done\n") + drawFoodweb fp + --putStr("next\n") + --fig2 fp + writeEV fp + writeEV0 fp + summarizeProps fp 3 + + return() + +mainFileOmni :: LoopyOpt -> IO () +mainFileOmni LoopyOpt{filename=fp,..} = do + l<- returnComponents fp + if length l==1 then getLoopyInfoA loopOutput fp + else do + outputAlpha fp + getLoopyInfoOmni loopOutput fp + getStabInfo fp + --putStr ("stab and loopy done") + drawFoodweb fp + fig2 fp + summarizeProps fp 3 + return() + +getLoopyInfo :: Bool ->FilePath -> IO () +getLoopyInfo loopOutput x = do + putStr ("Loopy running on: " ++ x ++ "\n") + z <- fileToFoodweb x + let y = checkFoodweb z + let t = getCommat y + let n = length t + m = matrixFromList t + g = graphFromMatrix n m + ps = summarize $ map (classify m) (loops g) --[propSum] + let titles = ["LL", "Num", "MLW","MLW+","MLW-", "MLWe","MLWo","MLWnet","MLW Loop"] -- , + tbl = titles : map propSumRow ps + res = outputMatrix tbl + ps2 = ("Loop\tLL\tLW\tSign\tCount\tNet\n") ++ (concat $ map propSumOutput $ [(classify m) ((loops g)!!i) | i <-[0..(length (loops g)-1)]]) + -- writeFile (dropExtension x ++ ".html") $ prefix ++ img (outputChart ps) ++ "<br/>" ++ htmlMatrix tbl ++ suffix + if loopOutput then writeFile (dropExtension x ++ ".loops") ps2 else putStr "Loops not outputted\n" + writeFile (dropExtension x ++ ".log") res + putStr (res)--from res + +--this currently only shows the omni loops in the log file, and that is the only difference +getLoopyInfoOmni :: Bool -> FilePath -> IO () +getLoopyInfoOmni loopOutput x = do + putStr ("Loopy running on: " ++ x ++ "\n") + z <- fileToFoodweb x + let y = checkFoodweb z + let t = getCommat y + let n = length t + m = matrixFromList t + g = graphFromMatrix n m + ps = summarize $ map (classify m) (loops g) --[propSum] + let titles = ["LL", "Num", "MLW","MLW+","MLW-", "MLWe","MLWo","MLWnet","MLWo Loop"] -- , + tbl = titles : map propSumRowOmni ps + res = outputMatrix tbl + ps2 = ("Loop\tLL\tLW\tSign\tCount\tNet\n") ++ (concat $ map propSumOutput $ [(classify m) ((loops g)!!i) | i <-[0..(length (loops g)-1)]]) + -- writeFile (dropExtension x ++ ".html") $ prefix ++ img (outputChart ps) ++ "<br/>" ++ htmlMatrix tbl ++ suffix + if loopOutput then writeFile (dropExtension x ++ ".loops") ps2 else putStr "Loops not outputted\n" + writeFile (dropExtension x ++ ".log") res + putStr (res)--from res + + + + +getLoopyInfoA :: Bool -> FilePath -> IO () +getLoopyInfoA loopOutput y = do + putStr ("Loopy running on: " ++ y ++ "\n") + input <-readAlpha y + let t = input + n = length t + m = matrixFromList t + g = graphFromMatrix n m + ps = summarize $ map (classify m) (loops g) --[propSum] + titles = ["LL", "Num", "MLW","MLW+","MLW-", "MLWe","MLWo","MLWnet","MLW Loop"] -- , + tbl = titles : map propSumRow ps + res = outputMatrix tbl + ps2 = ("Loop\tLL\tLW\tSign\tCount\tNet\n") ++ (concat $ map propSumOutput $ [(classify m) ((loops g)!!i) | i <-[0..(length (loops g)-1)]]) + + if loopOutput then writeFile (dropExtension y ++ ".loops") ps2 else putStr "Loops not outputted\n" + + writeFile (dropExtension y ++ ".log") res + putStr (res) + {- -} +{- +getLoopyInfoDetDir :: FilePath -> IO () +getLoopyInfoDetDir fp = do + allfiles <- getDirectoryContents fp + let allLoops = filter (\x -> takeExtension x == ".loops" ) allfiles + sequence [getLoopyInfoDet (fp </> s3) | s3 <- allLoops] --fp <\> creates a dir + return() +--loopsG = (detLoops (fromJust detRow) (loops g) +getLoopyInfoDet :: FilePath -> IO () +getLoopyInfoDet ll = do + loops <- readLoopList ll + fw<-fileToFoodweb ((dropExtension fp) ++ ".fw") + let ds = map deathRate $ info fw + let detRow = findIndex (==0) ds +-} + +writeEV :: FilePath -> IO [Double] +writeEV fp = do + fw<-fileToFoodweb fp + let cm = getCommat fw + res = getEVals cm + writeFile (dropExtension fp ++ ".evs") (unlines $ map show res) + return (res) + +writeEV0 :: FilePath -> IO [Double] +writeEV0 fp = do + fw<-fileToFoodweb2 fp + let res = getEVals $ diagS 0 fw + putStr(show res ++ "\n") + writeFile (dropExtension fp ++ ".evs0") (unlines $ map show res) + return (res) + +writeEV0a :: FilePath -> IO [Double] +writeEV0a fp = do + fw<-fileToFoodweb fp + let d = findD fw + cm = getCommat fw + cm0 = zerodiagd cm (fromJust d) + res = getEVals cm0 + writeFile (dropExtension fp ++ ".evs0") (unlines $ map show res) + return (res) + +zerodiagd :: Matrix Double -> Int -> Matrix Double +zerodiagd m d= mapMatrixInd f m + where + f i j x | i==j&&i/=(d) = 0 + | otherwise = x + +zerodiag :: Matrix Double -> Matrix Double +zerodiag m= mapMatrixInd f m + where + f i j x | i==j = 0 + | otherwise = x + +getLoopyInfoChoice :: FilePath -> Int-> IO () +getLoopyInfoChoice x ll = do + putStr ("Loopy running on: " ++ x ++ "\n") + z <- fileToFoodweb x + let y = checkFoodweb z + let t = getCommat y + let n = length t + m = matrixFromList t + g = graphFromMatrix n m + ps = summarize $ map (classify m) (loops g) --[propSum] + let titles = ["LL", "Num", "MLW","MLW+","MLW-", "MLWe","MLWo","MLWnet","MLWo Loop"] -- , + tbl = titles : map (propSumRowChoice ll ) ps + res = outputMatrix tbl + ps2 = ("Loop\tLL\tLW\tSign\tCount\tNet\n") ++ (concat $ map propSumOutput $ [(classify m) ((loops g)!!i) | i <-[0..(length (loops g)-1)]]) + -- writeFile (dropExtension x ++ ".html") $ prefix ++ img (outputChart ps) ++ "<br/>" ++ htmlMatrix tbl ++ suffix + --writeFile (dropExtension x ++ ".loops") ps2 + writeFile (dropExtension x ++ ".log") res + putStr (res)--from res + +getStabInfo :: FilePath -> IO () +getStabInfo x = do + putStr ("Stability for " ++ dropExtension x++": " ) + y <- fileToFoodweb x + let z = checkFoodweb y + z <- return $ z{cm=getCommat z} -- FIXME: Reorder to remove this + let s = findStab z + writeFile (dropExtension x ++ ".stab") (show s) + let evals = getEVals $ diagS s z + putStr (show evals ++ "\n") + putStr (show s ++ "\n") + +img x = "<img src=\"" ++ x ++ "\" /><br/>" + + +myMerge ::[FilePath] -> [FilePath] -> [FilePath] +myMerge [] ys = ys +myMerge (x:xs) ys = if (elem (dropExtension x) (map dropExtension ys)) then myMerge xs ys else [x]++ myMerge xs ys + +removeDet :: Foodweb -> Foodweb +removeDet f = deleteSps f (getDetRow f) + +writeFwND :: FilePath -> IO () +writeFwND fp = do + fw<-fileToFoodweb fp + let f = removeDet fw + writeFoodweb (dropExtension fp ++ "_ND.fw") f --dropExtension fp ++ "_ND" </> + return() + +--working, but not putting int hte new directory +writeFwNDDir :: FilePath -> IO () +writeFwNDDir fp = do + allfiles <- getDirectoryContents fp + createDirectoryIfMissing True $ dropExtension fp ++ "_ND" + let fs = filter (\x -> takeExtension x == ".fw") allfiles + putStr ("fs done\n") + sequence [writeFwND (fp </> s )| s<-fs] + return() + +outputGraphsDir :: FilePath -> IO () +outputGraphsDir fp = do + allfiles <- getDirectoryContents fp + let allLoops = filter (\x -> takeExtension x == ".loops" ) allfiles + sequence [writeTableLW (fp </> s1) | s1<- allLoops] --fp <\> creates a dir + sequence [writeTableLWo (fp </> s2) | s2 <- allLoops] --fp <\> creates a dir + sequence [writeTableLWd (fp </> s3) | s3 <- allLoops] --fp <\> creates a dir + --sequence [writeTableLWdNotONot (fp </> s3) | s3 <- allLoops] --fp <\> creates a dir + sequence [writeTableLWdNot (fp </> s3) | s3 <- allLoops] --fp <\> creates a dir + sequence [writeSummaryTableTypes (fp </> s3) | s3 <- allLoops] --fp <\> creates a dir + + + --putStr (show tt) + return() + +changeToFoodweb :: FilePath -> IO () +changeToFoodweb fp = do + fw<-fileToFoodweb fp + writeFoodweb (dropExtension fp ++ ".fw") fw + putStr (show fp) + return() + +changeToFoodwebDir :: FilePath -> IO () +changeToFoodwebDir fp = do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".DAT" ) allfiles + res<- sequence [changeToFoodweb (fp </> s) | s <- ss] --fp <\> creates a dir + return() + + +foodwebList :: [[Int]] -> Foodweb -> [Foodweb] +foodwebList [] fw = [] +foodwebList (x:xs) fw =[changeFeeding fw x] ++ foodwebList xs fw + +stabList :: [Foodweb] -> [Double] +stabList [] = [] +stabList (f:fs) =[findStab f] ++ stabList fs + + + +breakList :: Int -> [a] -> [[a]] +breakList n as =if length as <= n + then [as] + else fst y : breakList n ((snd y)) + where y = splitAt n as + + + + +{- +everything :: IO () +everything = do + main --on all the files + getStabVsProp fp rn cn = do +-} +--run main of a directory of files +{--} + +--create a foodweb +--vary biomasses if appropriate +--then run the below + +{- + +writeSpeciesChangeTable :: FilePath -> Int -> IO [Double] +writeSpeciesChangeTable fp n = do + createSFtoOsmoRange fp n + main (dropExtension fp ++ "SFtoOsmo") + stabs<-readStabs (dropExtension fp ++ "SFtoOsmo") + ss<- return $breakList (n+1) stabs + writeFile (dropExtension fp ++ "SFtoOsmo.txt") $ unlines $ map (intercalate "\t") $map2 show ss + return(stabs) +-} + +cf :: LoopyOpt -> IO () +cf o@LoopyOpt{filename=fp,..} = do + fw2<-fileToFoodweb fp + let g2 = checkFoodweb fw2 + let n = (length $ info g2)-5 + let g22=changeFeeding g2 (replicate n 2) + let g33=changeFeeding g2 (replicate n 3) + writeFoodweb (dropExtension fp ++ "33a.fw") g33 + writeFoodweb (dropExtension fp ++ "22a.fw") g22 + getLoopyInfo loopOutput (dropExtension fp ++ "33a.fw") + getStabInfo (dropExtension fp ++ "33a.fw") + getLoopyInfo loopOutput (dropExtension fp ++ "22a.fw") + getStabInfo (dropExtension fp ++ "22a.fw") + + +funct23 :: [[Int]] -> [[Int]] +funct23 [] = [] +funct23 x = [2:(head x), 3:(head x)]++ funct23 (tail x) + +generate23 n = iterate funct23 [[]]!!n + +funct13 :: [[Int]] -> [[Int]] +funct13 [] = [] +funct13 x = [1:(head x), 3:(head x)]++ funct13 (tail x) + +generate13 n = iterate funct13 [[]]!!n + +{-something going wrong here -} + +allFSDir :: FilePath -> IO () +allFSDir fp =do + allfiles <- getDirectoryContents fp + let tt = filter (\x -> takeExtension x == ".fw" ) allfiles + res<- sequence [allFS s | s <- tt] --fp <\> creates a dir + return() + +allFS :: FilePath -> IO ([([Int],Double)]) +allFS fp = do + fw2<-fileToFoodweb fp + let fw = checkFoodweb fw2 + let n = (length $ info fw)-5 + fws = foodwebList (generate23 n) fw + fwsA = map putCM fws + res<- sequence [writeFoodweb (dropExtension fp ++ show s ++ ".fw") (fwsA!!s) | s <- [0..(length (generate23 n)-1)]] + let stabs = stabList fwsA + res = zip (generate23 n) stabs + writeFile (dropExtension fp ++ ".stabLog") (tableTuple res) + return(res) + +---chemo fs stuff +{- +writeChemoRangeTable :: FilePath -> Int -> IO [Double] +writeChemoRangeTable fp n= do + createChemoRange fp n + mainPath (dropExtension fp ++ "ChemoRange") + sumAlphas<-readAlphaSum (dropExtension fp ++ "ChemoRange") + ss<- return $breakList (n+1) sumAlphas + writeFile (dropExtension fp ++ "ChemoRange.txt") $ unlines $ map (intercalate "\t") $map2 show ss + return(sumAlphas) +-} + + +allFSChemo :: FilePath -> IO ([([Int],Double)]) +allFSChemo fp = do + fw2<-fileToFoodweb fp + let fw = checkFoodweb fw2 + let n = (length $ info fw)-5 + fws = foodwebList (generate13 n) fw + fwsA = map putCM fws + res<- sequence [writeFoodweb (dropExtension fp ++ show s ++ ".fw") (fwsA!!s) | s <- [0..(length (generate23 n)-1)]] + let cms = map cm fwsA + let aphs = map sumAlpha cms -- change to alphas + res = zip (generate13 n) aphs + writeFile (dropExtension fp ++ ".sumAlphaLog") (tableTuple res) + return(res) + +--moved putCM tableTuple to FeedingRates.hs 20/08/2010 + +--------modern fish stuff +changeNoFish :: Foodweb -> Double -> Foodweb +changeNoFish (Foodweb names info wij cm) osProp= Foodweb names (info) (changeWij wij) cm + where + changeWij wij = mapMatrixInd f wij + where + f i j x|i ==6 && j ==3 =(1-osProp)*100 --sepenson + |i ==5 && j ==3 =osProp*100 --eating off detritus + |otherwise = x + +createNFRange :: FilePath -> Int-> IO () +createNFRange fp n = do + fw <- fileToFoodweb fp + createDirectoryIfMissing True (dropExtension fp ++ "PhyRange") + sequence_ [writeFoodweb file $ changeNoFish fw x + | x <- [0,1/(fromIntegral n)..1], let file = dropExtension fp ++ "PhyRange"</> dropExtension fp ++ "_" ++ show3dp x ++ ".fw"] + + + +------------osmo vs sf +writeOsmoRangeTable :: LoopyOpt -> Int -> IO [Double] +writeOsmoRangeTable o@LoopyOpt{filename=fp,..} n= do + createOsmoRange fp n + mainPath o{filename=(dropExtension fp ++ "OsmoRange")} + stabs<-readStabs (dropExtension fp ++ "OsmoRange")--sequence [readStabs (fp ++ "OsmoRange")| file<-files] + ss<- return $breakList (n+1) stabs + writeFile (dropExtension fp ++ "OsmoRange.txt") $ unlines $ map (intercalate "\t") $map2 show ss + return(stabs) + +writeOsmoRange2Table :: LoopyOpt -> Int -> IO [Double] +writeOsmoRange2Table o@LoopyOpt{filename=fp,..} n= do + createOsmoRange2 fp n + mainPath o{filename=(dropExtension fp ++ "OsmoRange2")} + stabs<-readStabs (dropExtension fp ++ "OsmoRange2")--sequence [readStabs (fp ++ "OsmoRange")| file<-files] + ss<- return $breakList (n+1) stabs + writeFile (dropExtension fp ++ "OsmoRange2.txt") $ unlines $ map (intercalate "\t") $map2 show ss + return(stabs) + + +writeChemoRangeTable :: LoopyOpt -> Int -> IO [Double] +writeChemoRangeTable o@LoopyOpt{filename=fp,..} n= do + createChemoRange fp n + mainPath o{filename=(dropExtension fp ++ "ChemoRange")} + sumAlphas<-readAlphaSum (dropExtension fp ++ "ChemoRange") + ss<- return $breakList (n+1) sumAlphas + writeFile (dropExtension fp ++ "ChemoRange.txt") $ unlines $ map (intercalate "\t") $map2 show ss + return(sumAlphas) + +--readAlphaSum is equivlanf stabs + +writeFrondRangeTable :: LoopyOpt -> Int -> Int -> IO [Double] +writeFrondRangeTable o@LoopyOpt{filename=fp,..} n1 n2 = do + createFrondRange fp n1 n2 + mainPath o{filename=(dropExtension fp ++ "FrondRange")} + stabs<-readStabs (dropExtension fp ++ "FrondRange")--sequence [readStabs (fp ++ "FrondRange")| file<-files] + ss<- return $breakList (n1+1) stabs + writeFile (dropExtension fp ++ "FrondRange.txt") $ unlines $ map (intercalate "\t") $map2 show ss + return(stabs) + +--------------extra stuff + +prefix = unlines + ["<html>" + ,"<head>" + ,"<style type='text/css'>" + ,"table {border: 1px solid black;}" + ,"body {font-family: sans-serif; font-size: small;}" + ,"td {padding-right: 10px;}" + ,"</style>" + ,"</head>" + ,"<body>" + ] + +suffix = "</body></html>" + +htmlMatrix :: Matrix String -> String +htmlMatrix = wrap "table" . concatMap (wrap "tr" . concatMap (wrap "td")) + where wrap tag str = "<" ++ tag ++ ">" ++ str ++ "</" ++ tag ++ ">" + + + + + +-- loops g = the list of all loops for graph g
+ Mlw.hs view
@@ -0,0 +1,321 @@+module Mlw where + + +import System.Environment +import System.Random +import Control.Monad +import System.Directory +import System.Cmd +import Data.Char +import Data.Maybe +import Data.List +import Numeric +import System.FilePath +import Numeric.LinearAlgebra.LAPACK +import Data.Packed.Matrix hiding (Matrix) +import Data.Packed.Vector +import Data.Complex +import Types +import FeedingRates +import Stab +import Graph + + +-------------------------------------------------------------------------------- + +readLoopList :: FilePath -> IO [[String]] +readLoopList fp = do + x<-readFile fp + let z = map words (lines x) + return(tail(z)) + + +-- output table of LW vs LL for all points + +writeTableLW :: FilePath -> IO () +writeTableLW fp = do + loops<-readLoopList fp + let x = map createCoord loops + writeFile (dropExtension fp ++ ".LWvsLL") (unlines x) + return() + + +createCoord :: [String] -> String +createCoord l = concat $ intersperse "\t" [l!!1,l!!2,l!!3] + + +-- output table of LWo vs LL for all points + +writeTableLWo :: FilePath -> IO () +writeTableLWo fp = do + loops<-readLoopList fp + let loopsO = filter (/=[]) (map filterOmni loops) + let x = map createCoordo loopsO + writeFile (dropExtension fp ++ ".LWovsLL") (unlines x) + return() + + +filterOmni :: [String] -> [String] +filterOmni [] = [] +filterOmni l = if ((readAsInt (l!!4)) == 2-(readAsInt (l!!1))) then l else [] + + +filterOmniNot :: [String] -> [String] +filterOmniNot [] = [] +filterOmniNot l = if ((readAsInt (l!!4)) == 2-(readAsInt (l!!1))) then [] else l + + +readAsInt :: String -> Int +readAsInt i = read i :: Int + +createCoordo :: [String] -> String +createCoordo l = concat $ intersperse "\t" [l!!1,l!!2,l!!3,l!!4] + {- read first column as-} +-- read in loop row with props and det row number +--convert first to list +--if list contains +--if elem detRow loop then output else Nothing + + +writeSummaryTableDet :: FilePath -> IO () +writeSummaryTableDet fp = do + let lw = ((dropExtension fp) ++ ".LWvsLL") + lwd = ((dropExtension fp) ++ ".LWdvsLL") + lwnd = ((dropExtension fp) ++ ".LWdNotvsLL") + a<-findMLW3 lw + b<-findMLWpos lw + c<-findMLWneg lw + d<-findMLW3o lw + a2<-findMLW3 lwd + b2<-findMLWpos lwd + c2<-findMLWneg lwd + d2<-findMLW3o lwd + a3<-findMLW3 lwnd + b3<-findMLWpos lwnd + c3<-findMLWneg lwnd + d3<-findMLW3o lwnd + + let alll = show a ++ "\t" ++ show b ++ "\t" ++ show c ++ "\t" ++ show d ++ "\n" --("\tLW\tLWo\tLWd\tLWnd\nMLW\t") ++ + let det = show a2 ++ "\t" ++ show b2 ++ "\t" ++ show c2 ++ "\t" ++ show d2 ++ "\n" --("MLW3\t") ++ + let tro = show a3 ++ "\t" ++ show b3 ++ "\t" ++ show c3 ++ "\t" ++ show d3 ++ "\n" --("MLW3\t") ++ + + writeFile ((dropExtension fp) ++ ".sumLWsDET")(alll ++ det ++ tro)-- + return() + + +writeSummaryTableTypes :: FilePath -> IO () +writeSummaryTableTypes fp = do + let lw = ((dropExtension fp) ++ ".LWvsLL") + lwo = ((dropExtension fp) ++ ".LWovsLL") + lwd = ((dropExtension fp) ++ ".LWdvsLL") + lwnd = ((dropExtension fp) ++ ".LWdNotvsLL") + a<- findMLWn2 lw + b<-findMLWn2 lwo + c<-findMLWn2 lwd + d<-findMLWn2 lwnd + a2<- findMLW3 lw + b2<-findMLW3 lwo + c2<-findMLW3 lwd + d2<-findMLW3 lwnd + let mlws = show a ++ "\t" ++ show b ++ "\t" ++ show c ++ "\t" ++ show d ++ "\n" --("\tLW\tLWo\tLWd\tLWnd\nMLW\t") ++ + let mlw3s = show a2 ++ "\t" ++ show b2 ++ "\t" ++ show c2 ++ "\t" ++ show d2 ++ "\n" --("MLW3\t") ++ + writeFile ((dropExtension fp) ++ ".sumLWs")(mlws ++ mlw3s)-- + return() +{--} +--LWd vs LL [(findMLWn2 lw),(findMLWn2 lwo),(findMLWn2 lwd),(findMLWn2 lwnd)] +writeTableLWd :: FilePath -> IO () +writeTableLWd fp = do + loops<-readLoopList ((dropExtension fp) ++ ".loops") + fw<-fileToFoodweb ((dropExtension fp) ++ ".fw") + --putStr("fw\n") + let ds = map deathRate $ info fw + --putStr(show ds ++ "\n") + let detRow = findIndex (==0) ds + --putStr(show detRow ++ "\n") + + let loopsD = filter (/= []) $ map (filterDetritus (fromJust detRow)) loops + loopsDV = filter (/=[]) $ map (filterDetrivoreNotStr fw) loopsD + putStr(show (length (loopsDV)) ++ " loopsDV\n"++ show (length (loops)) ++ " loops\n" ++ show (length (loopsD)) ++ " loopsD\n") + let x = map createCoord loopsDV + --putStr(show x++"\n") + writeFile ( (dropExtension fp) ++ ".LWdvsLL") (unlines x) + return() +filterDetrivoreNot :: Foodweb -> Loop -> Loop +filterDetrivoreNot fw l = if (groupElem de es) then [] else l -- if (length de) >=1 then () else l + where + e = getEdges fw + de = getDetrEdges n e + n = fromJust $ detRow fw + es = loopToEdges l + + +filterDetrivoreNotStr :: Foodweb -> [String] -> [String] +filterDetrivoreNotStr fw []= [] +filterDetrivoreNotStr fw l = if (groupElem de es) then [] else l -- [f (head (head (convertToList ll!!0))) detRow] ++ [filterDetritus f (tail ll) detRow] + where + de = getDetrEdges n (getEdges fw) + n = fromJust $ detRow fw + es = loopToEdges $ (convertToLoop (head l)) +--not working properly, i think it is still considering reverse end edges +writeTableLWdNot :: FilePath -> IO () +writeTableLWdNot fp = do + loops<-readLoopList ((dropExtension fp) ++ ".loops") + fw<-fileToFoodweb ((dropExtension fp) ++ ".fw") + --putStr("fw\n") + let ds = map deathRate $ info fw + --putStr(show ds ++ "\n") + let detRow = findIndex (==0) ds + --putStr(show detRow ++ "\n") + + let loopsD = filter (/= []) $ map (filterDetritusNot (fromJust detRow)) loops + loopsDV = filter (/=[]) $ map (filterDetrivoreStr fw) loops + loopsAll = loopsD -- ++ loopsDV + --putStr(show loopsD ++ "loopsD\n") + let x = map createCoord loopsAll + --putStr(show x++"\n") + writeFile ( (dropExtension fp) ++ ".LWdNotvsLL") (unlines x) + return() + + +filterDetrivore :: Foodweb -> Loop -> Loop +filterDetrivore fw l = if (length de) >=1 then (if (groupElem de es) then l else []) else [] + where + e = getEdges fw + de = getDetrEdges n e + n = fromJust $ detRow fw + es = loopToEdges l + +filterDetrivoreStr :: Foodweb -> [String] -> [String] +filterDetrivoreStr fw []= [] +filterDetrivoreStr fw l = if (groupElem de es) then l else [] -- [f (head (head (convertToList ll!!0))) detRow] ++ [filterDetritus f (tail ll) detRow] + where + de = getDetrEdges n (getEdges fw) + n = fromJust $ detRow fw + es = loopToEdges $ (convertToLoop (head l)) + +filterDV :: Foodweb -> [String] -> [String] +filterDV fw []= [] +filterDV fw l = undefined + + +filterDetritus :: Int -> [String] -> [String] +filterDetrius detRow []= [] +filterDetritus detRow l = if (elem detRow (convertToList (head l))) then l else [] -- [f (head (head (convertToList ll!!0))) detRow] ++ [filterDetritus f (tail ll) detRow] + +getEdges :: Foodweb -> Graph +getEdges fw = edges + where + edges = graphFromMatrix n (wij fw) + n = length (wij fw) --fromJust $ detRow fw + +getDetrEdges :: Int ->Graph-> Graph +getDetrEdges a [] = [] +getDetrEdges n e = if ((fst (head e)) == n) then [head e] ++ (getDetrEdges n (tail e)) else (getDetrEdges n (tail e)) + +getDetrEdgesFw :: Foodweb -> Graph +getDetrEdgesFw fw =getDetrEdges (fromJust (detRow fw)) (getEdges fw) + + +detRow :: Foodweb -> Maybe Int +detRow fw = findIndex (==0)$ map deathRate $ info fw + +loopToEdges2 :: Loop -> Graph +loopToEdges2 [] = [] +loopToEdges2 l = if ((length l) > 1) then [(head l ,head (tail l))] ++ loopToEdges2 ((tail l)) else [] + +loopToEdges :: Loop -> Graph +loopToEdges l = loopToEdges2 l ++ [(last l, head l)] + +justLoops :: [String] -> [String] +justLoops ll =undefined -- ( map head ll )-- :: [Loop] + + +--if any of l1 are contained in l2 returns l2 else returns empty list +groupElem :: (Eq a) => [a] -> [a] -> Bool +groupElem _ [] = False +groupElem [] l2 = False +groupElem l1 l2 = if elem (head l1) l2 then True else (groupElem (tail l1) l2) + +notGroupElem :: (Eq a) => [a] -> [a] -> Bool +notGroupElem [] _ = False +notGroupElem l1 l2 = if notElem (head l1) l2 then True else (notGroupElem (tail l1) l2) + +--NOTE the below has not been updated +writeTableLWdNotONot :: FilePath -> IO () +writeTableLWdNotONot fp = do + loops<-readLoopList ((dropExtension fp) ++ ".loops") + fw<-fileToFoodweb ((dropExtension fp) ++ ".fw") + --putStr("fw\n") + let ds = map deathRate $ info fw + --putStr(show ds ++ "\n") + let detRow = findIndex (==0) ds + --putStr(show detRow ++ "\n") + + let loopsD = filter (/= []) $ map (filterNoDetritusNorOmni (fromJust detRow)) loops + --putStr(show loopsD ++ "loopsD\n") + let x = map createCoord loopsD + --putStr(show x++"\n") + writeFile ( (dropExtension fp) ++ ".LWdNotOvsLL") (unlines x) + return() + +filterNoDetritusNorOmni :: Int -> [String] -> [String] +filterNoDetritusNorOmni detRow []= [] +filterNoDetritusNorOmni detRow l = filterOmniNot $ filterDetritusNot detRow l + +filterDetritusNot :: Int -> [String] -> [String] +filterDetritusNot detRow []= [] +filterDetritusNot detRow l = if (notElem detRow (convertToList (head l))) then l else [] -- [f (head (head (convertToList ll!!0))) detRow] ++ [filterDetritus f (tail ll) detRow] + + + +convertToList :: String -> [Int] +convertToList str = read str :: [Int] + +convertToLoop :: String -> Loop +convertToLoop str = read str :: Loop +--test $ propEven p == 2 - propLength p + +--output table of LWdet vs LL showing det vs non det + +--take in table showing LL, LW, sign output + +findMLW :: FilePath -> IO (Double) +findMLW fp =do + ls<-readLoopList fp + let ts = ((transpose ls)!!1) + let weights = map read ts:: [Double] + let res = maximum (weights) + return (res) + +findMLWn2 :: FilePath -> IO (Double) +findMLWn2 fp = do + ls <-readLoopList fp -- + let lsn2 = (map filterNot2s ls) + res = maximum (lsn2) + return(res) +findMLW3 :: FilePath -> IO (Double) +findMLW3 fp = do + ls <-readLoopList fp -- + let lsn2 = (map filter3s ls) + res = maximum (lsn2) + return(res) +findMLW3o :: FilePath -> IO (Double) +findMLW3o fp = undefined +findMLWpos :: FilePath -> IO (Double) +findMLWpos fp = undefined +findMLWneg :: FilePath -> IO (Double) +findMLWneg fp = undefined + + +filteros :: [String] -> Double +filteros ls = undefined +filterpos :: [String] -> Double +filterpos ls = undefined +filterneg :: [String] -> Double +filterneg ls = undefined + +filterNot2s :: [String] -> Double +filterNot2s ls = if ((ls!!0) /="2") then (read (ls!!1) :: Double) else 0 + +filter3s :: [String] -> Double +filter3s ls = if ((ls!!0) =="3") then (read (ls!!1) :: Double) else 0
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ Stab.hs view
@@ -0,0 +1,336 @@+{-# LANGUAGE ScopedTypeVariables #-} +module Stab where + + +import System.Environment +import System.Random +import Control.Monad +import System.Directory + +import System.Cmd +import Data.Char +import Data.Maybe +import Data.List +import Numeric +import System.FilePath +import Numeric.LinearAlgebra.LAPACK +import Data.Packed.Matrix hiding (Matrix) +import Data.Packed.Vector +import Data.Complex +import Types +import FeedingRates +import LoopProp +import Debug.Trace + +diagS :: Double -> Foodweb -> Matrix Double +diagS s fw = mapMatrixInd0 f $ cm fw + where + d = findD fw + f i j x|i == j && Just i /= d = s*x -- + |otherwise = x + +findD :: Foodweb -> Maybe Int +findD fw = findIndex (\x -> deathRate x == 0) $ info fw + +getEVals :: Matrix Double -> [Double] +getEVals m = map realPart $ toList $ fst $ eigR $ fromLists m + +filterPos :: [Double] -> Int +filterPos evals = length $ filter (>0) evals + +isHigher :: Foodweb -> Double -> Bool +isHigher fw smed = if zeroE ==0 && medE /= 0 || zeroE /=0 && medE == 0 then False else True + where + medE = filterPos $ getEVals $ diagS smed fw + zeroE = filterPos $ getEVals $ diagS 0 fw + +findF :: (Double -> Bool) -> (Double,Double) -> Double +findF f (lo,hi) | hi-lo < 0.00000001 = med + | f med = findF f (med, hi) + | otherwise = findF f (lo, med) + where + med = (lo + hi)/(fromIntegral 2) + +findStab :: Foodweb -> Double +findStab fw = if zeroE == 0 then 0 else findF (isHigher fw) (0,2) + where + zeroE = filterPos $ getEVals $ diagS 0 fw + +-- error $ show (v, stable m v, stable m (v+0.001), stable m (v-0.001)) +-- findStab fw = findF (isHigher fw) (0,2)-- error $ show (v, stable m v, stable m (v+0.001), stable m (v-0.001)) +{- +stable :: Matrix Double -> Double -> Bool +stable commat s = all (<= 0) $ getEVals $ diagS s commat +-} + +output :: [(Maybe Double, Double)] -> String +output xs = unlines [show a ++ "\t" ++ show b | (Just a,b) <- xs]--just shows things if they aren't NA + +showJust (Just x) = show x +showJust Nothing = "NA" + +nothing0 (Just x) = x +nothing0 Nothing = 0 + +sumAlpha :: Matrix Double -> Double +sumAlpha alpha = sum $ map sum alpha + +readAlpha :: FilePath -> IO (Matrix Double) +readAlpha fp = do + input <- readFile fp + let alpha = matrixStringToDouble $ toMatrix input + return(alpha) + +readAlphas :: FilePath -> IO [Matrix Double] +readAlphas fp = do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".alpha") allfiles + res<- sequence [readAlpha (fp </> s)| s <- ss] + alphas<-return $ res + return(alphas) + +readAlphaSum :: FilePath -> IO [Double] +readAlphaSum fp = do + alphas<-readAlphas fp + let sums = map sumAlpha alphas + return(sums) + +readStabs :: FilePath -> IO [Double] +readStabs fp = do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".stab") allfiles + res<- sequence [readStab (fp </> s)| s <- ss] + stabs<-return $ res + return(stabs) +{- +testStabs :: FilePath -> IO [Double] +testStabs fp = do + --res1<-return $readStabs + old<-return $ [0.0002,0.006,0.0005,0.0003,0.277,0.374,0.259,0.134,0.0092,0.0126,0.0107,0.0354,0.317,0.14,0.223,0.205,0.0001,0.0001,0,0.0001,0.359,0.251,0.309,0.29,0.235,0.178,0.124,0.24,0.204,0.42,0.459,0.285] :: [Double] + return(old) +-} + +readProps :: FilePath -> Int-> Int-> Int -> IO [Maybe Double] -- if want mlw not mlw3 etc use 0 or 1 as rn lm is min ll wanted, need to take intoaccoutn NAs +readProps fp ll cn lm= do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".log") allfiles + props<- sequence [readProp (fp </> s) cn lm | s <- ss] + -- + res<-if ll<=1 + then return $ [maximum proplist | proplist<-props] + else return $ [llProp ll lm proplist | proplist<-props] + return(res) + --not working due to integer/int confusion +--find max for each thing, total for num and not mlw loop +summarizeProps :: FilePath -> Int -> IO () +summarizeProps fw lmin = do + let fp=(dropExtension fw ++".log") + t2<-readProp fp (fromIntegral 2) (fromIntegral 2) + let lm = if length t2< lmin then 2 else lmin + t<-readProp fp (fromIntegral 2) (fromIntegral lm) + props<-sequence [readProp fp i (fromIntegral lm)|i<-[1..8]] + let num = sum $map fromJust t + maxtps = map (maxProp (fromIntegral lm)) props + maxll = concat $ intersperse "\t" $ map show $ map snd maxtps + maxs = concat $ intersperse "\t" $ map show $ map fst maxtps + title = "\n\nTable of Maximums for each loop catagory with the corresponding LL \n \t LL \t Num \t MLW \t MLW+ \t MLW- \t MLWe \t MLWo \t MLWnet \n" + res = "Total Number of Loops = " ++ show num ++ title ++ "LL\t" ++ maxll ++ "\n" ++"maxs\t" ++ maxs ++ "\n" + writeFile (dropExtension fp ++ ".sum") res + --putStr(res) + +findMaxLL :: FilePath -> IO () +findMaxLL fp = do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".sum") allfiles + mlws<-sequence [takeMLW (fp </> s) | s <- ss] + writeFile (dropExtension fp ++ ".mlws") (unlines mlws) + +takeMLW :: FilePath -> IO String +takeMLW sumfile = do + t<-readFile sumfile + let res = ((map words (lines t))!!4!!3)-- :: Double + return(res) +--get the proplist using readProp +maxProp :: Integer ->[Maybe Double] -> (Double,Integer) +maxProp lm proplist = maximum $ zip (map nothing0 proplist) [lm..] + + +readPropsT :: FilePath -> Int-> Int-> Int -> IO [Maybe Double] -- if want mlw not mlw3 etc use 0 or 1 as rn lm is min ll wanted, need to take intoaccoutn NAs +readPropsT fp ll cn lm= do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".log") allfiles + props<- sequence [readProp (fp </> s) cn lm | s <- ss] + -- + res<-if ll<=1 + then return $ [maximum proplist | proplist<-props] + else return $ [llProp ll lm proplist | proplist<-props] + return(res) + +--writeStabLoops + +writeStabProps :: FilePath -> Int -> Int -> Int-> IO () +writeStabProps fp ll cn lm= do + stabs<-readStabs fp + props<-readProps fp ll cn lm + res<-return $ zip props stabs + writeFile (fp ++ "Stabvs" ++ show ll ++ "_" ++ show cn ++ ".txt") $ output res + return() + +maybeRead :: Read a => String -> Maybe a +maybeRead = fmap fst . listToMaybe + . reads +readProp :: FilePath -> Int -> Int -> IO [Maybe Double]--need to adjust so returns [] if too big +readProp fp cn lm = do + x <-readLoopPropFile fp + n<- return $length x + prop<-return$ [x!!y!!(cn-1) | y<-[(lm-2)..(n-1)]]--check this is correct regarding the n-1 bit + prop2<-return $map maybeRead prop + return(prop2) + +--mapD :: String -> Maybe Double +--mapD x = read x :: Maybe Double + +llProp :: Int -> Int-> [Maybe Double] -> Maybe Double +llProp ll lm proplist = if (ll-lm+1)<= length proplist then proplist!!(ll-lm) else Nothing -- not working properly, not returning nothing + + +readStab :: FilePath -> IO Double -- checked ok +readStab x = do + s<-readFile $x--fp</>x + let y =read s :: Double + return(y) + +readLoopPropFile :: FilePath -> IO [[String]] +readLoopPropFile file = do + x <- readFile file + return $ tail $ map words (lines x) + + +readLoops :: FilePath -> Int -> IO [[Int]] +readLoops fp lm= do + x <-readLoopPropFile fp + n<- return $length x + prop<-return$ [x!!y!!8 | y<-[(lm-2)..(n-1)]]--check this is correct regarding the n-1 bit + --error $ show $ head prop + prop2<-return $map read prop + return(prop2) + +{--} +readLoop :: FilePath -> Int -> Int -> IO [Int] +readLoop fp lm ll = do + x <-readLoops fp lm + return(x!!(ll-lm)) + +{---assume that loopy has already been run +loopBiomassRatio :: Foodweb ->[Int] -> [Int] +loopBiomassRatio fw loop1= [ biomassRatio loop1!!x loop1!!y | x<-[0..(length loop1-1)],y<-(if x==(length loop1-1) then 0 else x+1) ] +-} + +biomassRatio :: Int -> Int -> Double +biomassRatio x y = (fromIntegral x)/(fromIntegral y) + + + +{- +------------------------------------------------------------------------------------------------------------------------------------------ +getSpProp :: FilePath -> Int -> Int ->IO String -- rn is the LL and cn is the prop col number +getSpProp file rn cn = do + x <-readLoopPropFile file + return $ x!!(rn-2)!!(cn-1) + +getStabVsProp :: FilePath -> Int -> Int -> IO () +getStabVsProp fp rn cn = do -- ll then prop + getAllStab fp + getAllProp fp rn cn + s<-readFile (fp ++ "Stabilities.txt") + s2 <- return $ read s :: IO [(String,Double)] + p <- readFile (fp ++ "Props.txt") -- if equal nothing need to remove the s val + p2 <- return $ read p :: IO [String] --need to read as double if number + stabs<- return $ map snd s2 + res<- return $ zip p2 stabs + writeFile (fp ++ "Stabvs" ++ show rn++ "_" ++ show cn ++ ".txt") $ show res + return () + + if head e1 == '[' then return $ Just $ Left e1 + else if e1 == "NA" then return Nothing + else return $ Just $ Right (read e1 :: Double) + +getAllStab :: FilePath -> IO () +getAllStab fp= do + --read list of files and on each one make list of name and stab - tuples + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".stab") allfiles + res<- sequence [getStabFromFile (fp </> s)| s <- ss] + writeFile (fp ++ "Stabilities.txt") $ show res + return() + + + + + +getAllProp:: FilePath-> Int -> Int -> IO () +getAllProp fp rn cn =do + allfiles <- getDirectoryContents fp + let ss = filter (\x -> takeExtension x == ".log") allfiles + res<- sequence [getSpProp (fp </> s) rn cn | s <- ss] --fp <\> creates a dir + writeFile (fp ++ "Props.txt") $ show res + return() + +-} + + +--data Either a b = Left a | Right b +--data Maybe a = Nothing | Just a + +{---get a certain result from loop matrix, eg mlw3 +getSpProp :: FilePath ->String -> Int -> IO () +getSpProp fp pr n = do + --read list of files and on each one make list of name and stab - tuples + allfiles <- getDirectoryContents fp + let loopInfo = filter (\x -> takeExtension x == ".log") allfiles + res<- sequence $ map getProp loopInfo pr + writeFile (fp ++ "Properties.txt") $ show res + return() +-} +--read file +--match string to column number +--get row from row number (-1) and match to column number +{- +read file, and choose a column + +with column, choose row number or max/min +y needs to read, then convert to matrix. match top column to pr, then n columns to loop length or max, min, with or without 2 loops + +-} +--plot s against different props loop weight + +--add on further loops prop - either direction of loop, omnivorous loops, divide by death rate or not +-- +-- +-- +-- +{- +NEED TO REDO THE BELOW STUFF, SO DON'T OUTPUT STAB AND PROP SUMMARY FILES, JUST THE OVERALL ONES +readStabDir :: FilePath -> IO [StabInfo] +readStabDir = + + +process path = do + file <- getdirectorycontents path + res <- sequence $ map readStabLogFile file + .... + writeFile "summary.txt" (output res2) + + + +processMe :: ([StabInfo] -> String) -> FilePath -> IO () + + +processMLEnet stabs = .... + + + + + +-} +--need to have extra function for looking at the actual loops
+ Types.hs view
@@ -0,0 +1,181 @@+{-# LANGUAGE ScopedTypeVariables #-} +module Types where + + +import System.Environment +import System.Random +import Control.Monad +import System.Directory +import System.Cmd +import Data.Char +import Data.Maybe +import Data.List +import Numeric +import System.FilePath +--maybe also contain basic utilities + +type Matrix a = [[a]] + +data Info = Info -- this says what is in the second input matrix + {deathRate :: Double + ,aeff :: Double + ,peff :: Double + ,biomass :: Double + } deriving Show + +data Foodweb = Foodweb + {names :: [String] + , info :: [Info] + , wij :: Matrix Double + ,cm :: Matrix Double + } deriving Show + + + +fileToFoodweb :: FilePath -> IO Foodweb +fileToFoodweb = readFoodweb -- FIXME: Inline this everywhere + + +readFoodweb :: FilePath -> IO Foodweb +readFoodweb file = do + input <-readFile file + let [names,info,wij] :: [Matrix String] = map toMatrix (splitString (lines (formatString input))) + return $ Foodweb (matrixStringToListString names) (map toInfo (matrixStringToDouble info)) (matrixStringToDouble wij) undefined + + +returnComponents :: FilePath -> IO [String] +returnComponents fp = do + input <-readFile fp + output <- return $ splitString (lines (formatString input)) + return(output) + +writeFoodweb :: FilePath -> Foodweb -> IO () +writeFoodweb file (Foodweb names info wij _) = writeFile file $ unlines $ + names ++ + [""] ++ + map (intercalate "\t" . map show . infoToDouble) info ++ + [""] ++ + [outputMatrix (map2 show wij)] + + +formatString :: String -> String +formatString input = unlines $ dropWhile null $ map (dropWhile isSpace) $ lines input +-- unlines takes list of strings to string +splitString :: [String] -> [String] +splitString input = + if null xs then + [] + else + (unlines y:splitString ys ) + where xs = dropWhile null input + (y,ys) = break null xs + +--readFile :: FilePath -> IO String +-- show stringyfys ints, doubles, lists of stuff, tuples , boolens, strings to strings +splitAndWriteFile :: FilePath -> IO () +splitAndWriteFile file = do + input <- readFile file + output :: [Matrix String]<- return (map toMatrix (splitString (lines (formatString input)))) + createDirectoryIfMissing True ("input") + names <- return [ ("input/" ++ dropExtension file ++ show i ++".txt")| i <- [1..] ]--return (dropExtension file ++ matrixNumber) + sequence $ map (uncurry writeFile) (zip names (map show output)) + return () + --uncurry takes a 2 input function and makes it accept a tuple instead + + +--replace any non empty sequences of space with tab +tabify (x:xs) | isSpace x = '\t' : tabify (dropWhile isSpace xs) + | otherwise = x : tabify xs +tabify [] = [] + +matrixFromList :: [[a]] -> Matrix a +matrixFromList x = x + + +toInfo :: [Double] -> Info +toInfo [_,a,b,c,d] = Info a b c d + +show3dp :: Double -> String +show3dp x = showFFloat (Just 3) x "" + +infoToMatrixDouble :: [Info] -> Matrix Double +infoToMatrixDouble input = map infoToDouble input + +infoToDouble :: Info -> [Double] +infoToDouble (Info a b c d) = [0,a,b,c,d] + +mean :: [Double] -> Double +mean a = sum(a)/(fromIntegral (length(a))) -- cant' divide a double by an int + +functionToMatrix :: (Int -> Int -> Double) -> Int -> Matrix Double +functionToMatrix f len = [[f i j | j <- [0..(len-1)]]| i <-[0..(len-1)]] + +outputMatrix :: Matrix String -> String +outputMatrix = unlines . map (concat . intersperse "\t") + +toMatrix :: String -> Matrix String +toMatrix input = map words (lines input) + +matrixStringToListString :: Matrix String -> [String] +matrixStringToListString = map head + +matrixStringToDouble :: Matrix String -> Matrix Double +matrixStringToDouble input = map (map read) input + +matrixDoubleToString :: Matrix Double -> Matrix String +matrixDoubleToString input = map (map show) input + +stringsToMatrixDouble :: [String] -> Matrix Double +stringsToMatrixDouble input = matrixStringToDouble $ map words input + +--applys a function to a matrix element, then returns the changed function +mapInd :: (Int->a->b) -> [a] -> [b] +mapInd f a = zipWith f [1..] a + +changeElemList :: Int -> (Double -> Double) -> [Double] -> [Double] +changeElemList i f a = mapInd (\inew x -> if inew ==i then f x else x) a + +mapMatrixInd ::(Int->Int->a->b) -> Matrix a -> Matrix b +mapMatrixInd f xs = mapInd (\i x -> mapInd (f i) x) xs + +changeElem :: Int -> Int -> (Double -> Double) ->Matrix Double -> Matrix Double +changeElem i j f m = mapMatrixInd (\inew jnew x -> if inew== i && jnew == j then f x else x ) m + +--applys a function to a matrix element, then returns the changed function +mapInd0 :: (Int->a->b) -> [a] -> [b] +mapInd0 f a = zipWith f [0..] a + +mapMatrixInd0 ::(Int->Int->a->b) -> Matrix a -> Matrix b +mapMatrixInd0 f xs = mapInd0 (\i x -> mapInd0 (f i) x) xs + + + +------------------------------------------------------------------- + +--matrix :: Size -> (Int -> Int -> a) -> Matrix a +--matrix _ m = m + +------------------------------- +--vectorToList :: Size -> Vector a -> [a] +--vectorToList n v = [v i | i <- [0..n-1]] + +--matrixToList :: Size -> Matrix a -> [[a]] +--matrixToList n m = [vectorToList n $ m i | i <- [0..n-1]] + +--operation + +transposeMatrix m i j = m j i + +rowMatrix m i = m i + +colMatrix m j = \i -> m i j + +-- General + +map2 :: (a -> b) -> [[a]] -> [[b]] +map2 f = map (map f) + + +table :: String -> [[String]] +table = map words . lines +
+ loopy.cabal view
@@ -0,0 +1,31 @@+cabal-version: >= 1.6+build-type: Simple+name: loopy+version: 0.0+stability: Alpha+license: BSD3+license-file: LICENSE+category: Unclassified+author: Emily Mitchell <emily.g.h.mitchell@gmail.com>+maintainer: Emily Mitchell <emily.g.h.mitchell@gmail.com>+copyright: Emily Mitchell 2009-2011+synopsis: Find all biological feedback loops within an ecosystem graph.+description: Find all biological feedback loops within an ecosystem graph.+homepage: http://www.esc.cam.ac.uk/people/research-students/emily-king++executable loopy+ main-is:+ Main.hs++ build-depends:+ base == 4.*, hmatrix, GoogleChart, cmdargs, containers, filepath, process, directory, random++ other-modules:+ Avalon+ DotFile+ FeedingRates+ Graph+ LoopProp+ Mlw+ Stab+ Types