loopy-0.0: Stab.hs

{-# 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