-- | Drawing of covariance model (http://www.tbi.univie.ac.at/software/cmcompare/) guide trees and highlighting comparison results
-- Drawing is done with the diagrams package
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies     #-}
{-# LANGUAGE RankNTypes #-}

module Bio.CMDraw
    (
     drawSingleCMComparisons,
     drawSingleCMs,
     drawCM,
     text',
     svgsize,
     diagramName,
     printCM,
     NodeIndices,
     buildRowIndexStructure,
     buildTreeIndexStructure,
     mergedSecondaryStructureVisualisation,
     perModelSecondaryStructureVisualisation
    ) where

import Diagrams.Prelude
import Bio.CMCompareResult
import Data.List
import qualified Data.Text as T
import qualified Data.Vector as V
import Bio.StockholmData
import Bio.StockholmDraw
import Diagrams.Backend.Cairo
import qualified Data.Vector.Unboxed as VU
import qualified Data.PrimitiveArray.Index.PhantomInt as PI
import Biobase.Types.Bitscore (Bitscore(..), score2Prob)
import Text.Printf
--import qualified Data.Colour.SRGB.Linear as R
import qualified Data.Colour.SRGB as R
import Data.Maybe
import Control.Monad.State
import qualified Data.Char as C
import Graphics.SVGFonts
import Bio.StockholmFont
import qualified Biobase.SElab.CM.Types as CM
import qualified Biobase.SElab.CM.ModelStructure as CM
import Data.Either.Unwrap
import qualified Data.Map as M
import Data.Function    

-- | Draw one or more CM
drawSingleCMComparisons :: String -> String -> Int -> Double -> String -> String -> Double -> Double -> [CM.CM] -> [Maybe StockholmAlignment] -> [CmcompareResult] -> [(QDiagram Cairo V2 Double Any,Maybe (QDiagram Cairo V2 Double Any))]
drawSingleCMComparisons layoutDirection modelDetail entryNumberCutoff transitionCutoff modelLayout emissiontype maxWidth scalef cms alns comparisons = diagrams
  where diagrams = map (drawCM layoutDirection modelDetail entryNumberCutoff transitionCutoff modelLayout emissiontype maxWidth scalef nameColorVector) zippedInput
        zippedInput = zip4 cms alns comparisonNodeLabels (V.toList colorVector)
        modelNumber = length cms
        comparisonNodeLabels = map (getComparisonNodeLabels comparisons nameColorVector) cms
        colorVector = makeColorVector modelNumber
        modelNames =  V.fromList (map (T.unpack . CM._name) cms)
        nameColorVector = V.zipWith (\a b -> (a,b)) modelNames colorVector

-- | Draw one or more CM
drawSingleCMs :: String -> String -> Int -> Double -> String -> String -> Double -> Double -> [CM.CM] -> [Maybe StockholmAlignment] -> [(QDiagram Cairo V2 Double Any,Maybe (QDiagram Cairo V2 Double Any))]
drawSingleCMs layoutDirection modelDetail entryNumberCutoff transitionCutoff modelLayout emissiontype maxWidth scalef cms alns = map (drawCM layoutDirection modelDetail entryNumberCutoff transitionCutoff modelLayout emissiontype maxWidth scalef emptyColorNameVector) zippedInput
    where zippedInput = zip4 cms alns comparisonNodeLabels colorList
          comparisonNodeLabels = map getBlankComparisonNodeLabels cms
          colorList = replicate (length cms) white
          emptyColorNameVector = V.empty

-- | Draw the guide Tree of a single CM
drawCM :: String -> String -> Int -> Double -> String -> String -> Double -> Double -> V.Vector (String,Colour Double) -> (CM.CM,Maybe StockholmAlignment,V.Vector (Int,V.Vector (Colour Double)), Colour Double) -> (QDiagram Cairo V2 Double Any,Maybe (QDiagram Cairo V2 Double Any))
drawCM layoutDirection modelDetail entryNumberCutoff transitionCutoff modelLayout emissiontype maxWidth scalef nameColorVector (inputCM,aln,comparisonNodeLabels,modelColor)
   | modelLayout == "tree" = ((applyAll ([bg white]) modelTreeLayout),alignmentDiagram)
   | modelLayout == "flat" = ((applyAll ([bg white]) modelFlatLayout),alignmentDiagram)
   | otherwise = ((applyAll ([bg white])modelTreeLayout),alignmentDiagram)
   where cm = fromLeft (CM._cm inputCM) -- select Flexible Model
         nodes = V.fromList (M.elems (CM._fmNodes cm))
         nodeNumber = V.length nodes
         allStates = CM._fmStates cm
         boxlength = fromIntegral (length alphabetSymbols) + 2
         alphabetSymbols = ['A','U','C','G']
         indices = V.toList (V.iterateN (nodeNumber-1) (1+) 0)
         (indexStructure,_)= runState (buildTreeIndexStructure 1 nodes indices) startState
         modelName = CM._name inputCM
         currentCat = if layoutDirection == "vertical" then vcat' with {_sep=0.01} else hcat' with {_sep=0.01}
         modelFlatLayout = alignTL (vcat' with {_sep=0.1} [modelHeader,nodeTransitions]) # scale scalef
         modelTreeLayout = alignTL (vcat' with {_sep=0.1} [modelHeader,nodeTreeTransitions]) #scale scalef
         nodeTreeTransitions = applyAll (arrowList ++ labelList) nodesTree
         nodeTransitions = applyAll (arrowList ++ labelList) nodesFlat
         firstInterval = fromJust (find (\(_,p,_,_,_) -> p == 0) (fst indexStructure))
         nodesTree = drawCMNodeTree layoutDirection modelDetail alphabetSymbols emissiontype boxlength allStates comparisonNodeLabels nodes (fst indexStructure) firstInterval
         modelHeader = makeModelHeader (T.unpack modelName) modelColor nameColorVector
         nodeIndices = V.iterateN nodeNumber (1+) 0
         nodesFlat = currentCat (V.toList (V.map (drawCMNode layoutDirection modelDetail alphabetSymbols emissiontype boxlength (0 :: Int) nodeNumber nodeNumber allStates comparisonNodeLabels nodes) nodeIndices))
         allConnectedStates = makeAllConnectedStates allStates
         highConnectedStates = V.filter (\(_,_,w) -> w >= transitionCutoff) allConnectedStates
         connectedStates = V.filter (\(stateId,targetStateIndex,_) -> stateId /= targetStateIndex) highConnectedStates
         selfConnectedStates = V.filter (\(stateId,targetStateIndex,_) -> stateId == targetStateIndex) highConnectedStates
         arrowList = case modelDetail of
                          "detailed" -> V.toList (V.map makeArrow connectedStates V.++ V.map makeSelfArrow selfConnectedStates)
                          "interval"-> map (makeArrow . indexStructureToConnections) (filter (\(acc,emit,_,_,_)-> acc /= emit)(fst indexStructure))
                          "minimal"-> V.toList (V.map makeArrow connectedStates)
                          "simple"-> V.toList (V.map makeArrow connectedStates)
                          _ -> []
         labelList = case modelDetail of
                          "detailed" -> V.toList (V.map makeLabel connectedStates V.++ V.map makeSelfLabel selfConnectedStates)
                          _ -> [] 
         alignmentDiagram = drawStockholmLinesComparisonLabel entryNumberCutoff maxWidth comparisonNodeLabels nodes aln

drawStockholmLinesComparisonLabel :: Int -> Double -> V.Vector (Int,V.Vector (Colour Double)) -> V.Vector CM.Node -> Maybe StockholmAlignment -> Maybe (QDiagram Cairo V2 Double Any)
drawStockholmLinesComparisonLabel entryNumberCutoff maxWidth comparisonNodeLabels nodes maybeAln
   | isJust maybeAln = Just alignmentVis
   | otherwise = Nothing
     where aln = fromJust maybeAln
           -- comparison labels do not include end nodes, but no root node.
           -- end nodes have no alignment columns associated
           columnComparisonLabels = getComparisonPerColumnLabels comparisonNodeLabels nodes
           alignmentVis = drawStockholmLines entryNumberCutoff maxWidth columnComparisonLabels aln
                                      
makeAllConnectedStates :: M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (String,String,Double)
makeAllConnectedStates allStates = allConnectedStates
  where indexStateTuples = M.assocs allStates
        allConnectedStates = V.fromList (concatMap makeStateConnections indexStateTuples)
  
makeStateConnections :: (PI.PInt () CM.StateIndex,CM.State) -> [(String,String,Double)]
makeStateConnections (pInt,currentState) = conns
  where stateId = show (PI.getPInt pInt)
        targetBitscoreVector = VU.toList (CM._stateTransitions currentState)
        conns = map (\(target,bitscore) ->  (stateId,show(PI.getPInt target),score2Prob 1 bitscore)) targetBitscoreVector


-- | Extracts consensus secondary structure from alignment and annotates cmcompare nodes for each model-model combination seperatly
perModelSecondaryStructureVisualisation :: String -> Double -> String -> [CM.CM] -> [Maybe StockholmAlignment] -> [CmcompareResult] -> [(String,String)]
perModelSecondaryStructureVisualisation selectedTool _ structureFilePath cms alns comparisons
  | selectedTool == "forna" = fornaVis
  | selectedTool == "r2r" = r2rVis
  | selectedTool == "fornaLink" = fornaLink
  | selectedTool == "r2rfornaLink" = fornaLink ++ r2rVis
  | selectedTool == "all" = fornaLink ++ r2rVis ++ fornaVis
  | otherwise = []
  where fornaVis = concatMap (buildFornaperModelInput structureFilePath) structureComparisonInfo
        fornaLink = concatMap (buildFornaLinksInput structureFilePath) structureComparisonInfo
        r2rVis = concatMap (buildR2RperModelInput structureFilePath) structureComparisonInfo
        modelNumber = length cms
        comparisonNodeLabels = map (getComparisonPerModelNodeLabels comparisons nameColorVector) cms
        colorVector = makeColorVector modelNumber
        modelNames = V.fromList (map (T.unpack . CM._name) cms)
        nameColorVector = V.zipWith (\a b -> (a,b)) modelNames colorVector
        structureComparisonInfo = zip3 cms alns comparisonNodeLabels

getComparisonPerModelNodeLabels :: [CmcompareResult] -> V.Vector (String, Colour Double) -> CM.CM -> V.Vector (String,Colour Double, V.Vector (Int,V.Vector (Colour Double)))
getComparisonPerModelNodeLabels comparsionResults colorVector model = modelComparisonLabels
   where modelName = T.unpack (CM._name model)
         relevantComparisons1 = filter ((modelName==) . model1Name) comparsionResults
         modelNodeInterval1 = map (\a -> (model2Name a,nub (model1matchednodes a)))  relevantComparisons1 
         relevantComparisons2 = filter ((modelName==) . model2Name) comparsionResults
         modelNodeInterval2 = map (\a -> (model1Name a,nub (model2matchednodes a)))  relevantComparisons2
         modelNodeIntervals =  V.fromList (modelNodeInterval1 ++ modelNodeInterval2)
         nodeNumber = (CM._nodesInModel model)
         modelComparisonLabels = V.map (getModelComparisonLabels modelName nodeNumber colorVector) modelNodeIntervals

getModelComparisonLabels :: String -> Int -> V.Vector (String, Colour Double) -> (String,[Int])-> (String,Colour Double,V.Vector (Int,V.Vector (Colour Double)))
getModelComparisonLabels _ nodeNumber colorVector (compModel,matchedNodes) = (compModel,modelColor,comparisonNodeLabels)
  where (modelColor,modelInterval) = modelToColor colorVector (compModel,matchedNodes)
        -- cm starts at node index 0 for root node and ends with end node
        -- cmcompare does not include root node, but end node 
        comparisonNodeLabels = V.generate (nodeNumber) (makeModelComparisonNodeLabel (modelColor,modelInterval))

makeModelComparisonNodeLabel :: (Colour Double,[Int]) -> Int -> (Int,V.Vector (Colour Double))
makeModelComparisonNodeLabel (modelColor, nodeInterval) nodeNumber 
  | elem nodeNumber nodeInterval = (nodeNumber,V.singleton modelColor)
  | otherwise = (nodeNumber,V.singleton white)

getComparisonPerColumnLabels :: V.Vector (Int,V.Vector (Colour Double)) -> V.Vector CM.Node -> V.Vector (Int, V.Vector (Colour Double))
getComparisonPerColumnLabels comparisonNodeLabels nodes = columnComparisonLabels   where 
         unsortedColumnComparisonLabel = concatMap (nodeToColumnComparisonLabel nodes) (V.toList comparisonNodeLabels)
         columnComparisonLabels = V.fromList (sortBy (compare `on` fst) unsortedColumnComparisonLabel)

nodeToColumnComparisonLabel:: V.Vector CM.Node -> (Int, V.Vector (Colour Double)) -> [(Int,V.Vector (Colour Double))]
nodeToColumnComparisonLabel nodes (nodeIndex,colors) = colLabels
  where currentNode = (V.!) nodes (nodeIndex)
        colIndices = nub [CM._nodeColL currentNode,CM._nodeColR currentNode]
        colLabels = map (\a->(a,colors)) colIndices
                                   
--
buildR2RperModelInput :: String -> (CM.CM, Maybe StockholmAlignment,V.Vector (String,Colour Double,V.Vector (Int,V.Vector (Colour Double)))) -> [(String,String)]
buildR2RperModelInput structureFilePath (inputCM,maybeAln,comparisonNodeLabels)
  | isNothing maybeAln = []
  | otherwise = if V.null comparisonNodeLabels then singler2rInput else V.toList r2rComparisonInputs
  where cm = fromLeft (CM._cm inputCM) -- select Flexible Model
        modelName = T.unpack (CM._name inputCM)
        nodes = V.fromList (M.elems (CM._fmNodes cm))
        aln = fromJust maybeAln
        r2rInputPrefix = sHeader ++ sConsensusStructure ++ sConsensusSequence ++ sConsensusSequenceColor ++ sCovarianceAnnotation 
        allColumnAnnotations = columnAnnotations aln
        consensusSequenceList = map annotation (filter (\annotEntry -> tag annotEntry == T.pack "RF") allColumnAnnotations)
        consensusSequence = if null consensusSequenceList then "" else T.unpack (head consensusSequenceList)
        gapFreeConsensusSequence = map C.toUpper (filter (not . isGap) consensusSequence)
        consensusStructureList = map (convertWUSStoDotBracket . annotation) (filter (\annotEntry -> tag annotEntry == T.pack "SS_cons") allColumnAnnotations)
        consensusStructure = if null consensusStructureList then "" else T.unpack (head consensusStructureList)
        indexedGapFreeConsensusStructure = extractGapfreeIndexedStructure consensusSequence consensusStructure
        consensusStructureColIndices = map ((+1) . fst) indexedGapFreeConsensusStructure
        gapFreeConsensusStructure = map snd indexedGapFreeConsensusStructure
        --maxEntryLength = length consensusStructure
        --convert node to column labels as needed for consensus secondary structure
        columnComparisonLabels = V.map (\(mname,mcolor,comparisonNodePerModelLabels) -> (mname,mcolor,getComparisonPerColumnLabels comparisonNodePerModelLabels nodes)) comparisonNodeLabels
        --filter for labels that are part of consensus secondary structure by index
        consensusStructureColumnComparisonLabels = V.map (\(mname,mcolor,colLabels) -> (mname,mcolor,V.filter (\(i,_) -> elem i consensusStructureColIndices) colLabels)) columnComparisonLabels
        --consensusStructureColumnComparisonLabels = V.filter (\(i,_) -> elem i consensusStructureColIndices) columnComparisonLabels
        sHeader =  "# STOCKHOLM 1.0\n"
        sConsensusStructure =     "#=GC SS_cons          " ++ gapFreeConsensusStructure ++ "\n"
        sConsensusSequence =      "#=GC cons             " ++ gapFreeConsensusSequence ++ "\n" -- ++ show consensusStructureColIndices ++ "\n" ++ show comparisonNodeLabels ++ "\n"
        sConsensusSequenceColor = "#=GC conss            " ++ replicate (length gapFreeConsensusSequence) '2' ++ "\n"
        sCovarianceAnnotation =   "#=GC cov_SS_cons      " ++ replicate (length gapFreeConsensusSequence) '.' ++ "\n"
        singleFilePath = structureFilePath ++ modelName ++ ".r2r"                        
        singler2rInput = [(singleFilePath,r2rInputPrefix)]
        -- for multiple comparisons we need to return different filenames and labels
        r2rComparisonInputs = V.map (buildR2RperModelComparisonInput modelName structureFilePath r2rInputPrefix) consensusStructureColumnComparisonLabels

buildR2RperModelComparisonInput :: String -> String -> String -> (String, Colour Double, V.Vector (Int, V.Vector (Colour Double))) -> (String,String)
buildR2RperModelComparisonInput modelName structureFilePath r2rInputPrefix (compModelName,modelColor,columnComparisonLabels) = (schemeFilePath,r2rInput)
  where schemeFilePath = structureFilePath ++ modelName ++ "." ++ compModelName ++ ".r2r"
        r2rLabels = map comparisonColLabelsToR2RLabel (V.toList columnComparisonLabels)
        sComparisonHighlight =    "#=GC R2R_LABEL        " ++ r2rLabels ++ "\n"
        backBoneColor = setBackboneColor modelColor
        sBackboneColorLabel =     "#=GF R2R shade_along_backbone s rgb:" ++ backBoneColor ++ "\n"
        r2rInput = r2rInputPrefix ++ sComparisonHighlight ++ sBackboneColorLabel

setBackboneColor :: Colour Double -> String
setBackboneColor modelColor = show ((R.channelRed currentColor)* 255) ++ "," ++ show ((R.channelGreen currentColor) * 255) ++ "," ++ show ((R.channelBlue currentColor) * 255)
  where currentColor = R.toSRGB modelColor

buildFornaperModelInput :: String -> (CM.CM,Maybe StockholmAlignment,V.Vector (String,Colour Double,V.Vector (Int,V.Vector (Colour Double)))) -> [(String, String)]
buildFornaperModelInput structureFilePath (inputCM,maybeAln,comparisonNodeLabelsPerModels)
  | isNothing maybeAln = []
  | otherwise = fornaInput:colorSchemes
  where cm = fromLeft (CM._cm inputCM) -- select Flexible Model
        nodes = V.fromList (M.elems (CM._fmNodes cm))
        aln = fromJust maybeAln
        fornaString = ">" ++ modelName ++ "\n" ++ gapfreeConsensusSequence ++ "\n" ++ gapFreeConsensusStructure
        fornaFilePath = structureFilePath ++ modelName ++ ".forna"
        fornaInput = (fornaFilePath,fornaString)
        allColumnAnnotations = columnAnnotations aln
        consensusSequenceList = map annotation (filter (\annotEntry -> tag annotEntry == T.pack "RF") allColumnAnnotations)
        consensusSequence = if null consensusSequenceList then "" else T.unpack (head consensusSequenceList)
        gapfreeConsensusSequence = map C.toUpper (filter (not . isGap) consensusSequence)
        consensusStructureList = map (convertWUSStoDotBracket . annotation) (filter (\annotEntry -> tag annotEntry == T.pack "SS_cons") allColumnAnnotations)
        consensusStructure = if null consensusStructureList then "" else T.unpack (head consensusStructureList)
        indexedGapFreeConsensusStructure = extractGapfreeIndexedStructure consensusSequence consensusStructure
        consensusStructureColIndices = map ((+1) . fst) indexedGapFreeConsensusStructure
        gapFreeConsensusStructure = map snd indexedGapFreeConsensusStructure
        modelName = T.unpack (CM._name inputCM)           
        columnComparisonLabels = V.map (\(mname,mcolor,comparisonNodePerModelLabels) -> (mname,mcolor,getComparisonPerColumnLabels comparisonNodePerModelLabels nodes)) comparisonNodeLabelsPerModels
        --filter for labels that are part of consensus secondary structure by index
        consensusStructureColumnComparisonLabels = V.map (\(mname,mcolor,colLabels) -> (mname,mcolor,V.filter (\(i,_) -> elem i consensusStructureColIndices) colLabels)) columnComparisonLabels
        colorSchemes = V.toList (V.map (makeColorScheme modelName structureFilePath) consensusStructureColumnComparisonLabels)
    
                       
buildFornaLinksInput :: String -> (CM.CM,Maybe StockholmAlignment,V.Vector (String,Colour Double,V.Vector (Int,V.Vector (Colour Double)))) -> [(String, String)]
buildFornaLinksInput structureFilePath (inputCM,maybeAln,comparisonNodeLabelsPerModels)
  | isNothing maybeAln = []
  | otherwise = if V.null comparisonNodeLabelsPerModels then singleFornaLink else fornaComparisons
  where cm = fromLeft (CM._cm inputCM) -- select Flexible Model
        nodes = V.fromList (M.elems (CM._fmNodes cm))
        aln = fromJust maybeAln
        --http://nibiru.tbi.univie.ac.at/forna/forna.html?id=fasta&file=%3Eheader\nAACGUUAGUU\n(((....)))&colors=%3Eheader\n0\n0.1\n0.2\n0.3\n0.4\n0.5\n0.6\n0.7\n0.8\n0.9\n1
        fornaURLPrefix = "http://rna.tbi.univie.ac.at/forna/forna.html?id=fasta&file=%3Eheader\\n" ++ gapfreeConsensusSequence ++ "\\n" ++ gapFreeConsensusStructure
        singleFornaLink = [(fornaFilePath,fornaURLPrefix)]
        fornaFilePath = structureFilePath ++ modelName ++ ".fornalink"
        allColumnAnnotations = columnAnnotations aln
        consensusSequenceList = map annotation (filter (\annotEntry -> tag annotEntry == T.pack "RF") allColumnAnnotations)
        consensusSequence = if null consensusSequenceList then "" else T.unpack (head consensusSequenceList)
        gapfreeConsensusSequence = map C.toUpper (filter (not . isGap) consensusSequence)
        modelName = T.unpack (CM._name inputCM)
        consensusStructureList = map (convertWUSStoDotBracket . annotation) (filter (\annotEntry -> tag annotEntry == T.pack "SS_cons") allColumnAnnotations)
        consensusStructure = if null consensusStructureList then "" else T.unpack (head consensusStructureList)
        indexedGapFreeConsensusStructure = extractGapfreeIndexedStructure consensusSequence consensusStructure
        consensusStructureColIndices = map ((+1) . fst) indexedGapFreeConsensusStructure
        gapFreeConsensusStructure = map snd indexedGapFreeConsensusStructure    
        columnComparisonLabels = V.map (\(mname,mcolor,comparisonNodePerModelLabels) -> (mname,mcolor,getComparisonPerColumnLabels comparisonNodePerModelLabels nodes)) comparisonNodeLabelsPerModels
        --filter for labels that are part of consensus secondary structure by index
        consensusStructureColumnComparisonLabels = V.map (\(mname,mcolor,colLabels) -> (mname,mcolor,V.filter (\(i,_) -> elem i consensusStructureColIndices) colLabels)) columnComparisonLabels
        fornaComparisons = V.toList (V.map (makeFornaComparisonLink modelName structureFilePath fornaURLPrefix) consensusStructureColumnComparisonLabels)

makeFornaComparisonLink ::  String -> String -> String -> (String,Colour Double,V.Vector (Int,V.Vector (Colour Double))) -> (String,String)
makeFornaComparisonLink modelName structureFilePath fornaURLPrefix (compModelName,_,comparisonColLabelsPerModel) = (comparisonPath,comparisonLink)
  where comparisonPath = structureFilePath ++ modelName ++ "." ++ compModelName ++ ".fornalink"
        comparisonLink = fornaURLPrefix ++ labelPrefix ++ singleColorLabels 
        labelPrefix = "&colors=%3Eheader\\n"            
        --forna only supports a single color per node, which has to be supplied as additional color scheme
        singleColorLabels = concatMap comparisonColLabelsToFornaLinkLabel (V.toList comparisonColLabelsPerModel)
        
comparisonColLabelsToFornaLinkLabel :: (Int, V.Vector (Colour Double)) -> String
comparisonColLabelsToFornaLinkLabel (_,colorVector)
  | V.null colorVector = ""
  | V.head colorVector /= white =  "1\\n"
  | otherwise = "0\\n"
    
makeColorScheme ::  String -> String -> (String,Colour Double,V.Vector (Int,V.Vector (Colour Double))) -> (String,String)
makeColorScheme modelName structureFilePath (compModelName,_,comparisonColLabelsPerModel) = (schemeFilePath,singleColorLabels)
  where schemeFilePath = structureFilePath ++ modelName ++ "." ++ compModelName ++ ".fornacolor"
        indexedComparisonColLabelsPerModel = V.indexed comparisonColLabelsPerModel
        --column indexes have to be mapped to gap free consensus sequence
        structureIndexedLabels = V.map (\(a,(_,c)) -> (a+1,c)) indexedComparisonColLabelsPerModel
        singleColorLabels = concatMap comparisonColLabelsToFornaLabel (V.toList structureIndexedLabels)
        
-- | Extracts consensus secondary structure from alignment and annotates cmcompare nodes for all comparisons in one merged output
mergedSecondaryStructureVisualisation :: String -> Double -> [CM.CM] -> [Maybe StockholmAlignment] -> [CmcompareResult] -> [(String,String)]
mergedSecondaryStructureVisualisation selectedTool _ cms alns comparisons
  | selectedTool == "forna" = fornaVis
  | selectedTool == "r2r" = r2rVis
  | otherwise = []
  where fornaVis = map buildMergedFornaInput structureComparisonInfo
        r2rVis = map buildMergedR2RInput structureComparisonInfo
        modelNumber = length cms
        comparisonNodeLabels = map (getComparisonNodeLabels comparisons nameColorVector) cms
        colorVector = makeColorVector modelNumber
        modelNames = V.fromList (map (T.unpack . CM._name) cms)
        nameColorVector = V.zipWith (\a b -> (a,b)) modelNames colorVector
        structureComparisonInfo = zip3 cms alns comparisonNodeLabels
        
buildMergedFornaInput :: (CM.CM,Maybe StockholmAlignment,V.Vector (Int, V.Vector (Colour Double))) -> (String, String)
buildMergedFornaInput (inputCM,maybeAln,comparisonNodeLabels)
  | isNothing maybeAln = ([],[])
  | otherwise = (fornaInput, colorScheme)
  where cm = fromLeft (CM._cm inputCM) -- select Flexible Model
        nodes = V.fromList (M.elems (CM._fmNodes cm))
        aln = fromJust maybeAln
        fornaInput = ">" ++ modelName ++ "\n" ++ gapfreeConsensusSequence ++ "\n" ++ gapFreeConsensusStructure
        allColumnAnnotations = columnAnnotations aln
        consensusSequenceList = map annotation (filter (\annotEntry -> tag annotEntry == T.pack "RF") allColumnAnnotations)
        consensusSequence = if null consensusSequenceList then "" else T.unpack (head consensusSequenceList)
        gapfreeConsensusSequence = map C.toUpper (filter (not . isGap) consensusSequence)
        modelName = T.unpack (CM._name inputCM)
        consensusStructureList = map (convertWUSStoDotBracket . annotation) (filter (\annotEntry -> tag annotEntry == T.pack "SS_cons") allColumnAnnotations)
        consensusStructure = if null consensusStructureList then "" else T.unpack (head consensusStructureList)
        indexedGapFreeConsensusStructure = extractGapfreeIndexedStructure consensusSequence consensusStructure
        --consensusStructureColIndices = map ((+1) . fst) indexedGapFreeConsensusStructure
        gapFreeConsensusStructure = map snd indexedGapFreeConsensusStructure
        columnComparisonLabels = getComparisonPerColumnLabels comparisonNodeLabels nodes
        --forna only supports a single color per node, which has to be supplied as additional color scheme
        singleColorLabels = concatMap comparisonColLabelsToFornaLabel (V.toList columnComparisonLabels)
        colorScheme = singleColorLabels

comparisonColLabelsToFornaLabel :: (Int, V.Vector (Colour Double)) -> String
comparisonColLabelsToFornaLabel (nodeNr,colorVector)
  | V.null colorVector = ""
  | V.head colorVector /= white =  " " ++ show nodeNr ++ ":blue "
  | otherwise = ""
        
buildMergedR2RInput :: (CM.CM, Maybe StockholmAlignment,V.Vector (Int,V.Vector (Colour Double))) -> (String,String)
buildMergedR2RInput (inputCM,maybeAln,comparisonNodeLabels)
   | isNothing maybeAln = ([],[])
   | otherwise = (r2rInput,[])
  where cm = fromLeft (CM._cm inputCM) -- select Flexible Model
        nodes = V.fromList (M.elems (CM._fmNodes cm))
        aln = fromJust maybeAln
        r2rInput = sHeader ++ sConsensusStructure ++ sConsensusSequence ++ sConsensusSequenceColor ++ sCovarianceAnnotation ++ sComparisonHighlight ++ sBackboneColorLabel
        allColumnAnnotations = columnAnnotations aln
        consensusSequenceList = map annotation (filter (\annotEntry -> tag annotEntry == T.pack "RF") allColumnAnnotations)
        consensusSequence = if null consensusSequenceList then "" else T.unpack (head consensusSequenceList)
        gapFreeConsensusSequence = map C.toUpper (filter (not . isGap) consensusSequence)
        consensusStructureList = map (convertWUSStoDotBracket . annotation) (filter (\annotEntry -> tag annotEntry == T.pack "SS_cons") allColumnAnnotations)
        consensusStructure = if null consensusStructureList then "" else T.unpack (head consensusStructureList)
        indexedGapFreeConsensusStructure = extractGapfreeIndexedStructure consensusSequence consensusStructure
        --consensusStructureColIndices = map ((+1) . fst) indexedGapFreeConsensusStructure
        gapFreeConsensusStructure = map snd indexedGapFreeConsensusStructure
        columnComparisonLabels = getComparisonPerColumnLabels comparisonNodeLabels nodes
        r2rLabels = map comparisonColLabelsToR2RLabel (V.toList columnComparisonLabels)
        sHeader =  "# STOCKHOLM 1.0\n"
        sConsensusStructure =     "#=GC SS_cons          " ++ gapFreeConsensusStructure ++ "\n"
        sConsensusSequence =      "#=GC cons             " ++ gapFreeConsensusSequence ++ "\n"
        sConsensusSequenceColor = "#=GC conss            " ++ replicate (length gapFreeConsensusStructure) '2' ++ "\n"
        sCovarianceAnnotation =   "#=GC cov_SS_cons      " ++ replicate (length gapFreeConsensusStructure) '.' ++ "\n"
        sComparisonHighlight =    "#=GC R2R_LABEL        " ++ r2rLabels ++ "\n"
        sBackboneColorLabel =     "#=GF R2R shade_along_backbone s rgb:200,0,0\n"

comparisonColLabelsToR2RLabel :: (Int, V.Vector (Colour Double)) -> Char
comparisonColLabelsToR2RLabel (_,colorVector)
  | V.null colorVector = '.'
  | V.head colorVector /= white = 's'
  | otherwise = '.'

-- nodeToColIndices :: (Int,(Int,V.Vector (Colour Double))) -> (Int,V.Vector (Colour Double))
-- nodeToColIndices (colIndex,(_,colors)) = (colIndex,colors)

-- fillComparisonColLabels :: Int -> V.Vector (Int, V.Vector (Colour Double)) ->  V.Vector (Int, V.Vector (Colour Double))
-- fillComparisonColLabels maxEntryLength sparseComparisonColLabels = fullComparisonColLabels
--    where fullComparisonColLabels = V.generate maxEntryLength (makeFullComparisonColLabel sparseComparisonColLabels)

-- makeFullComparisonColLabel :: V.Vector (Int, V.Vector (Colour Double)) -> Int -> (Int, V.Vector (Colour Double))
-- makeFullComparisonColLabel sparseComparisonColLabels colIndex = fullComparisonColLabel
--   where availableLabel = V.find (\(a,_)-> colIndex == a) sparseComparisonColLabels
--         fullComparisonColLabel = fromMaybe (colIndex,V.singleton white) availableLabel

indexStructureToConnections :: (Int, Int, String, Int, Int) -> (String,String,Double)
indexStructureToConnections (acc,emit,_,_,_) = (show emit,show acc,1)

data NodeIndices = S [Int] | L [Int] | R [Int]
  deriving (Show, Eq, Ord)

startState :: ([(Int,Int,String,Int,Int)],Int)
startState = ([],0::Int)

buildRowIndexStructure :: Int -> V.Vector CM.Node -> [Int] -> State ([(Int,Int,String,Int,Int)],Int) ([(Int,Int,String,Int,Int)],Int)
buildRowIndexStructure _ _ [] = do
  (a,b) <- get
  return (a,b)
buildRowIndexStructure row nodes (currentIndex:xs) = do
  (currentInterval,parentId) <- get
  let currentNode = nodes V.! currentIndex
  let currentEnd = getIndexEnd nodes (currentIndex:xs)
  let ntype = CM._nodeType currentNode
  case ntype of
    CM.Root -> put ((row,parentId,"S,",currentIndex,currentEnd):currentInterval,parentId) -- ROOT start tree             
    CM.BegL -> put ((row,parentId,"L,",currentIndex,currentEnd):currentInterval,parentId) -- BEGL set current label
    CM.BegR -> put ((row,parentId,"R,",currentIndex,currentEnd):currentInterval,parentId) -- BEGR set current label
    CM.Bif -> put (currentInterval,parentId+1)
    _ -> put (currentInterval,parentId)
  buildRowIndexStructure row nodes xs

buildTreeIndexStructure :: Int -> V.Vector CM.Node -> [Int] -> State ([(Int,Int,String,Int,Int)],Int) ([(Int,Int,String,Int,Int)],Int)
buildTreeIndexStructure intervalId nodes (currentIndex:xs) = do
  (currentInterval,parentId) <- get
  let currentNode = nodes V.! currentIndex
  let currentEnd = getIndexEnd nodes (currentIndex:xs)
  let ntype = CM._nodeType currentNode
  let maxId = if null currentInterval then 0 else maximum $ map (\(iid,_,_,_,_)-> iid) currentInterval
  let newId = maxId +1
  let nextId = setNextId ntype intervalId newId
  case ntype of
    CM.Root -> put ((intervalId,parentId,"S,",currentIndex,currentEnd):currentInterval,parentId)
    CM.BegL -> put ((newId,parentId,"L,",currentIndex,currentEnd):currentInterval,parentId)
    CM.BegR -> put ((newId,parentId,"R,",currentIndex,currentEnd):currentInterval,parentId)
    CM.Bif -> put (currentInterval,intervalId)
    _ -> put (currentInterval,parentId)
  buildTreeIndexStructure nextId nodes xs
buildTreeIndexStructure _ _ [] = do
  (a,b) <- get
  return (a,b)
  
setNextId :: CM.NodeType -> Int -> Int -> Int
setNextId ntype intervalId newId
  | ntype == CM.Root = newId
  | ntype == CM.BegL = newId
  | ntype == CM.BegR = newId
  | otherwise = intervalId

getIndexEnd :: V.Vector CM.Node -> [Int] -> Int
getIndexEnd nodes indices
  | null indices = length nodes -1
  | ntype == CM.End = currentIndex
  | ntype == CM.Bif = currentIndex
  | otherwise = getIndexEnd nodes remainingindices
   where currentIndex = head indices
         remainingindices = tail indices
         currentNode = nodes V.! currentIndex
         ntype = CM._nodeType currentNode

makeModelHeader :: String -> Colour Double -> V.Vector (String,Colour Double) -> QDiagram Cairo V2 Double Any
makeModelHeader mName modelColor nameColorVector = (strutX 2 ||| setModelName mName ||| strutX 1 ||| rect 12 12 # lw 0.1 # fc modelColor # translate (r2 (negate 0, 5))) === strutY 1 === (strutX 30 ||| modelLegend)
  where modelLegend = makeModelLegend otherModelsNameColorVector
        otherModelsNameColorVector = V.filter ((/=mName) . fst) nameColorVector

makeModelLegend :: V.Vector (String,Colour Double) -> QDiagram Cairo V2 Double Any
makeModelLegend nameColorVector
  | V.null nameColorVector = mempty
  | otherwise = (legendHead === legendBody) <> rect boxX boxY # lw 0.1 # translate (r2 ((boxX/2)-1, negate (boxY/2) + 6))
  where legendHead = setLegendLabel "Legend:"
        legendBody = vcat (V.toList (V.map makeLegendEntry nameColorVector))
        nameLengths = V.map (length . fst) nameColorVector
        maxNameLength = fromIntegral $ V.maximum nameLengths
        entryNumber = fromIntegral $ V.length nameColorVector
        boxX = maxNameLength * 6
        boxY = entryNumber * 15 + 10

makeLegendEntry :: (String,Colour Double) -> QDiagram Cairo V2 Double Any
makeLegendEntry (mName,mColor) = setLegendLabel mName ||| strutX 0.5 ||| rect 4 4 # lw 0.1 # fc mColor # translate (r2 (negate 0, 2))

setLabel :: String -> QDiagram Cairo V2 Double Any
setLabel t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 2 2) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 0.75, negate 0.25))
setTransition :: String -> QDiagram Cairo V2 Double Any
setTransition t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 2 2) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 0.75, negate 0.75))   # rotate (1/4 @@ turn)
setState :: String -> Double -> Double -> QDiagram Cairo V2 Double Any
setState t x y = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 3 3) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (x, y))
setNodeNumber :: String -> QDiagram Cairo V2 Double Any
setNodeNumber t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 5 5) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 2, 0))
setNodeLabel :: String -> QDiagram Cairo V2 Double Any
setNodeLabel t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 6 6) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 2, 0))
setLegendLabel :: String -> QDiagram Cairo V2 Double Any
setLegendLabel t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 10 10) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 0.75, negate 0.75))
setModelName :: String -> QDiagram Cairo V2 Double Any
setModelName t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 20 20) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 0.75, negate 0.75))

drawCMNodeTree :: String -> String -> String -> String -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double))-> V.Vector CM.Node -> [(Int, Int, String, Int, Int)] -> (Int,Int,String,Int,Int) -> QDiagram Cairo V2 Double Any
drawCMNodeTree layoutDirection modelDetail alphabetSymbols emissiontype boxlength allStates comparisonNodeLabels nodes indexStructure (intervalId,parentId,intervalType,currentIndex,currentEnd) = nodeTree
  where nodeTree = if layoutDirection == "vertical" then verticalNodeTree else horizontalNodeTree
        verticalNodeTree = currentIntervalDrawing === hcat' with {_sep = 20} (map (drawCMNodeTree layoutDirection modelDetail alphabetSymbols emissiontype boxlength allStates comparisonNodeLabels nodes indexStructure) nextIntervals)
        horizontalNodeTree = currentIntervalDrawing ||| vcat' with {_sep = 20} (map (drawCMNodeTree layoutDirection modelDetail alphabetSymbols emissiontype boxlength allStates comparisonNodeLabels nodes indexStructure) nextIntervals)
        nextIntervals = filter (\(_,p,_,_,_) -> intervalId == p) indexStructure
        currentIntervalDrawing = drawCMNodeInterval layoutDirection modelDetail alphabetSymbols emissiontype boxlength currentIndex currentEnd currentEnd allStates comparisonNodeLabels nodes (intervalId,parentId,intervalType,currentIndex,currentEnd)  --- ||| (text' (show intervalId ++ "I" ++ show indexStructure) <> rect 100 100)

--drawCMNodeRow :: String -> String -> String -> Int -> Int -> Int -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double))-> V.Vector CM.Node -> [(Int, Int, String, Int, Int)] -> QDiagram Cairo V2 Double Any
--drawCMNodeRow modelDetail alphabetSymbols emissiontype boxlength rowStart rowEnd lastIndex states comparisonNodeLabels nodes intervals = strutY 4 === hcat' with { _sep = 8 } (map (drawCMNodeInterval modelDetail alphabetSymbols emissiontype boxlength rowStart rowEnd lastIndex states comparisonNodeLabels nodes) intervals)

drawCMNodeInterval :: String -> String -> String -> String -> Int -> Int -> Int -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double))-> V.Vector CM.Node -> (Int, Int, String, Int, Int) -> QDiagram Cairo V2 Double Any
drawCMNodeInterval layoutDirection modelDetail alphabetSymbols emissiontype boxlength rowStart rowEnd lastIndex states comparisonNodeLabels nodes (intervalId,_,_,currentIndex,currentEnd)
  | modelDetail == "interval" = intervalVis
  | otherwise = nodeVis
  where intervalVis = if layoutDirection == "vertical" then verticalIntervalVis else horizontalIntervalVis
        verticalIntervalVis = rect 20 0 # named ("a" ++ intervalIdString)  # lw 0.0 === (rect 20 40 # lw 0.1 <> text' (show currentIndex ++ "-" ++ show currentEnd)) === rect 20 0 # named ("e" ++ intervalIdString)  # lw 0.0 === strutY 5.0
        horizontalIntervalVis = rect 20 0 # named ("a" ++ intervalIdString)  # lw 0.0 ||| (rect 20 40 # lw 0.1 <> text' (show currentIndex ++ "-" ++ show currentEnd)) ||| rect 20 0 # named ("e" ++ intervalIdString)  # lw 0.0 ||| strutX 5.0
        intervalIdString = show intervalId
        nodeVis = if layoutDirection == "vertical" then verticalNodeVis else horizontalNodeVis
        verticalNodeVis = strutY 4  === vcat' with { _sep = nodespacer } (map (drawCMNode layoutDirection modelDetail alphabetSymbols emissiontype boxlength rowStart rowEnd lastIndex states comparisonNodeLabels nodes) currentInterval)
        horizontalNodeVis = strutX 5 ||| hcat' with { _sep = nodespacer } (map (drawCMNode layoutDirection modelDetail alphabetSymbols emissiontype boxlength rowStart rowEnd lastIndex states comparisonNodeLabels nodes) currentInterval)
        currentInterval = [currentIndex..currentEnd]
        nodespacer = if modelDetail == "detailed" then (0 :: Double) else (0.5 :: Double)

getCMNodeType :: CM.Node -> String
getCMNodeType node
  | ntype == CM.Bif = "BIF"
  | ntype == CM.MatP = "MATP"
  | ntype == CM.MatL = "MATL"
  | ntype == CM.MatR = "MATR"
  | ntype == CM.BegL = "BEGL"
  | ntype == CM.BegR = "BEGR"
  | ntype == CM.Root = "ROOT"
  | ntype == CM.End = "END"
  | otherwise = "NA"
    where ntype = CM._nodeType node

text' :: String -> QDiagram Cairo V2 Double Any
text' t = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False 3 3) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate 0.75, negate 0.75))

--textWithSize' :: String -> Double -> QDiagram Cairo V2 Double Any
--textWithSize' t si = textSVG_ (TextOpts linLibertineFont INSIDE_H KERN False si si) t # fc black # fillRule EvenOdd # lw 0.0 # translate (r2 (negate siOffset, negate siOffset))
--  where siOffset = si/2

-- | Transform covariance model node labels to colors
--labelToColor :: String -> Colour Double
--labelToColor label
--   | label == "MATP" = sRGB24 211 211 211 -- P
--   | label == "MATL" = sRGB24 211 211 211 -- L
--   | label == "MATR" = sRGB24 211 211 211 -- R
--   | label == "BIF"  = sRGB24 255 069 064 -- B
--   | label == "ROOT" = sRGB24 245 245 245 -- S
--   | label == "BEGL" = sRGB24 211 211 211 -- S
--   | label == "BEGR" = sRGB24 211 211 211 -- S 
--   | label == "END"  = sRGB24 245 245 245 -- E
--labelToColor _ = sRGB24 245 245 245

drawCMNode :: String -> String -> String -> String -> Int -> Int -> Int -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double)) -> V.Vector CM.Node -> Int -> QDiagram Cairo V2 Double Any
drawCMNode layoutDirection modelDetail alphabetSymbols emissiontype boxlength _ _ _ states comparisonNodeLabels nodes nodeIndex
  | modelDetail == "minimal" = drawCMMinimalNodeBox layoutDirection alphabetSymbols emissiontype boxlength states comparisonNodeLabels node nodeIndex
  | modelDetail == "simple" = drawCMSimpleNodeBox layoutDirection alphabetSymbols emissiontype boxlength states comparisonNodeLabels node nodeIndex
  | otherwise = detailedNodeBox
  where node = nodes V.! nodeIndex
        --idNumber = nodeIndex
        --nId = show idNumber
        detailedNodeBox = drawCMDetailedNodeBox layoutDirection alphabetSymbols emissiontype boxlength states comparisonNodeLabels node nodeIndex
        --nodeType = getCMNodeType node
        --nodeLabels = V.toList (snd (comparisonNodeLabels V.! idNumber))

colorBox :: Double -> Colour Double -> QDiagram Cairo V2 Double Any
colorBox singleBoxYLength colColour = rect 5 singleBoxYLength # fc colColour # lw 0.1

drawCMMinimalNodeBox :: String -> String -> String -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double)) -> CM.Node -> Int -> QDiagram Cairo V2 Double Any
drawCMMinimalNodeBox layoutDirection alphabetSymbols emissiontype boxlength currentStates comparisonNodeLabels node nodeIndex
  | ntype == CM.Bif = currentCat minimalNode  splitStatesBox -- bifNode 
  | ntype == CM.BegL = currentCat splitStatesBox minimalNode -- begLNode 
  | ntype == CM.BegR = currentCat splitStatesBox minimalNode -- begRNode 
  | otherwise = minimalNode
    where ntype = CM._nodeType node
          idNumber = nodeIndex
          nId = show idNumber
          currentCat = if layoutDirection == "vertical" then (===) else (|||)
          stateIndices = V.toList (CM._nodeStates node)
          minimalNode = rect 5 5 # lw 0.1 # lc black <> text' nId # fontSize 2  <> wheel 2 nodeLabels # lw 0.1 # lc black
          splitStatesBox = hcat (map (drawCMSimpleStateBox nId alphabetSymbols emissiontype boxlength currentStates) stateIndices)
          --nodeType = getCMNodeType node
          nodeLabels = V.toList (snd (comparisonNodeLabels V.! idNumber))
          --boxNumber = fromIntegral $ length nodeLabels
          --totalBoxYlength = 5 
          --singleBoxYLength = totalBoxYlength / boxNumber
          --colourBoxes = vcat (map (colorBox singleBoxYLength) nodeLabels)

drawCMSimpleNodeBox :: String -> String -> String -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double)) -> CM.Node -> Int -> QDiagram Cairo V2 Double Any
drawCMSimpleNodeBox layoutDirection alphabetSymbols emissiontype boxlength currentStates comparisonNodeLabels node nodeIndex
  | ntype == CM.Bif = currentCat simpleNode splitStatesBox -- bifNode 
  | ntype == CM.BegL = currentCat splitStatesBox simpleNode -- begLNode 
  | ntype == CM.BegR = currentCat splitStatesBox simpleNode -- begRNode 
  | otherwise = simpleNode
    where ntype = CM._nodeType node
          idNumber = nodeIndex 
          nId = show idNumber
          currentCat = if layoutDirection == "vertical" then (===) else (|||)
          stateIndices = V.toList (CM._nodeStates node)
          simpleNode = rect 10 5 # lw 0.1 # lc black <>  ((text' nId # translate (r2 (negate 7.5,0)) <> colourBoxes # translate (r2 (negate 7.5, boxYoffset))) ||| text' nodeType # translate (r2 (14,0)))
          splitStatesBox = hcat (map (drawCMSimpleStateBox nId alphabetSymbols emissiontype boxlength currentStates) stateIndices)
          nodeType = getCMNodeType node
          nodeLabels = V.toList (snd (comparisonNodeLabels V.! idNumber))
          boxNumber = fromIntegral $ length nodeLabels
          totalBoxYlength = 5 
          singleBoxYLength = totalBoxYlength / boxNumber
          -- concatenated colorboxes are placed atop the simplenode box with the first colorbox
          boxYoffset = totalBoxYlength/2 - singleBoxYLength/2
          colourBoxes = vcat (map (colorBox singleBoxYLength) nodeLabels)

drawCMSimpleStateBox :: String -> String -> String -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> PI.PInt () CM.StateIndex -> QDiagram Cairo V2 Double Any
drawCMSimpleStateBox _ _ _ _ currentStates sIndex
  | stype == CM.S = sState 
  | stype == CM.B = bState 
  | otherwise = mempty
    where currentState = currentStates M.! sIndex
          stype = CM._stateType currentState
          stateIndx = show (PI.getPInt sIndex)
          sState = rect 1 0.001 # lw 0 # named ("a" ++ stateIndx)
          bState = rect 1 0.001 # lw 0 # named ("e" ++ stateIndx)

drawCMDetailedNodeBox :: String -> String -> String -> Int ->  M.Map (PI.PInt () CM.StateIndex) CM.State -> V.Vector (Int, V.Vector (Colour Double)) -> CM.Node -> Int -> QDiagram Cairo V2 Double Any
drawCMDetailedNodeBox layoutDirection alphabetSymbols emissiontype boxlength currentStates comparisonNodeLabels node nodeIndex
  | ntype == CM.Bif = bifNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.MatP = matPNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.MatL = matLNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.MatR = matRNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.BegL = begLNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.BegR = begRNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.Root = rootNode # translate (r2 (negate 25,25)) <> nodeBox
  | ntype == CM.End = endNode # translate (r2 (negate 25,25)) <> nodeBox
  | otherwise = endNode <> nodeBox
    where ntype = CM._nodeType node
          idNumber = nodeIndex
          nId = show idNumber
          nodeLabels = V.toList (snd (comparisonNodeLabels V.! idNumber))
          stateIndices = V.toList (CM._nodeStates node)
          splitStatesBox = hcat' with { _sep = 0.01 } (map (drawCMSplitStateBox nId alphabetSymbols emissiontype boxlength currentStates) stateIndices)
          insertStatesBox = hcat (map (drawCMInsertStateBox nId alphabetSymbols emissiontype boxlength currentStates) stateIndices)
          -- bif b
          bifNode = (idBox nId "BIF" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5 ||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- matP mp ml mr d il ir
          matPNode = (idBox nId "MATP" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- matL ml d il
          matLNode = (idBox nId "MATL" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5 ||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- matR mr d ir
          matRNode = (idBox nId "MATR" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- begL s
          begLNode = (idBox nId "BEGL" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5 ||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- begR s il
          begRNode = (idBox nId "BEGR" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- root s il ir
          rootNode = (idBox nId "ROOT" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5||| splitStatesBox) === strutY 5.0 === insertStatesBox
          -- end e
          endNode = (idBox nId "END" nodeLabels # rotate (1/4 @@ turn) # translate (r2 (1, negate 17)) ||| strutX 0.5 ||| splitStatesBox) === strutY 5.0 === insertStatesBox

idBox :: String -> String -> [Colour Double] -> QDiagram Cairo V2 Double Any
idBox nId nType nodeLabels = (setNodeNumber nId # translate (r2 (negate ((fromIntegral (length nId))/2), 0)) <>  wheel 4 nodeLabels # lw 0.1 # translate (r2 (0, 1)) <> rect 3 3 # lw 0) ||| strutX 1.0 ||| setNodeLabel nType
nodeBox :: QDiagram Cairo V2 Double Any
nodeBox = rect 60 60 # lw 0.1

wheel :: Double -> [Colour Double] -> QDiagram Cairo V2 Double Any
wheel wsize colors = wheel' # rotate r
   where
     wheel' = mconcat $ zipWith fc colors (iterateN n (rotate a) w)
     n = length colors
     a = 1 / fromIntegral n @@ turn
     w = wedge wsize xDir a # lwG 0
     r = (1/4 @@ turn)  ^-^  (1/(2*fromIntegral n) @@ turn)

drawCMSplitStateBox :: String -> String -> String -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> PI.PInt () CM.StateIndex -> QDiagram Cairo V2 Double Any
drawCMSplitStateBox _ _ emissiontype _ currentStates sIndex
  | stype == CM.D = dState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | stype == CM.MP = mpState # translate (r2 (negate 7,negate 1)) <> statebox 16.0 20.0 stateIndx
  | stype == CM.ML = mlState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | stype == CM.MR = mrState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | stype == CM.S = sState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | stype == CM.E = eState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | stype == CM.B = bState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | stype == CM.EL = elState # translate (r2 (negate 3,negate 1)) <> statebox 8.0 20.0 stateIndx
  | otherwise = mempty
    where currentState = currentStates M.! sIndex
          stype = CM._stateType currentState
          stateIndx = show (PI.getPInt sIndex)
          singleEmissionBitscores = CM._stateEmissions currentState
          singleEmissionEntries = setEmissions emissiontype 4 singleEmissionBitscores
          singleSymbolsAndEmissions = zip ["A","U","G","C"] (VU.toList singleEmissionEntries)
          pairEmissionBitscores = CM._stateEmissions currentState
          pairEmissionEntries = setEmissions emissiontype 16 pairEmissionBitscores
          pairSymbolsAndEmissions = zip ["AA","AU","AG","AC","UU","UA","UG","UC","GG","GA","GU","GC","CC","CA","CU","CG"] (VU.toList pairEmissionEntries)
          pairSymbolsAndEmissions1 = take 8 pairSymbolsAndEmissions
          pairSymbolsAndEmissions2 = drop 8 pairSymbolsAndEmissions
          dState = setState ("D" ++ stateIndx) (negate 0.5) (negate 1)  === strutY 1 
          mpState = setState ("MP" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1 === (vcat (map (emissionEntry emissiontype) pairSymbolsAndEmissions1) ||| strutX 0.5 ||| vcat (map (emissionEntry emissiontype) pairSymbolsAndEmissions2))
          mlState = setState ("ML" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1 === vcat (map (emissionEntry emissiontype) singleSymbolsAndEmissions)
          mrState = setState ("MR" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1 === vcat (map (emissionEntry emissiontype) singleSymbolsAndEmissions)
          sState = setState ("S" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1
          eState = setState ("E" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1
          bState = setState ("B" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1
          elState = setState ("EL" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1

drawCMInsertStateBox :: String -> String -> String -> Int -> M.Map (PI.PInt () CM.StateIndex) CM.State -> PI.PInt () CM.StateIndex -> QDiagram Cairo V2 Double Any
drawCMInsertStateBox _ _ emissiontype _ currentStates sIndex
  | stype == CM.IL = ((ilState # translate (r2 (negate 3,negate 1))) <> statebox 8.0 20.0 stateIndx) ||| strutX 38
  | stype == CM.IR = (irState # translate (r2 (negate 3,negate 1))) <> inverseStatebox 8.0 20.0 stateIndx
  | otherwise = mempty
    where currentState = currentStates M.! sIndex
          stype = CM._stateType currentState
          stateIndx = show (PI.getPInt sIndex)
          singleEmissionBitscores = CM._stateEmissions currentState
          singleEmissionEntries = setEmissions emissiontype 4 singleEmissionBitscores
          singleSymbolsAndEmissions = zip ["A","U","G","C"] (VU.toList singleEmissionEntries)
          ilState = setState ("IL" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1 === vcat (map (emissionEntry emissiontype) singleSymbolsAndEmissions)
          irState = setState ("IR" ++ stateIndx) (negate 0.5) (negate 1) === strutY 1 === vcat (map (emissionEntry emissiontype) singleSymbolsAndEmissions)

setEmissions :: String -> Double -> VU.Vector Bitscore -> VU.Vector Double
setEmissions emissiontype normalizationFactor emissions
  | emissiontype == "score" = scoreentries
  | emissiontype == "probability" = propentries
  | emissiontype == "bar" = barentries
  | otherwise = barentries
    where scoreentries = VU.map bitScore2Double emissions
          propentries = VU.map ((/normalizationFactor) . score2Prob 1) emissions
          barentries = VU.map ((/normalizationFactor) . score2Prob 1) emissions

--wrap :: a -> [a]
--wrap x = [x]

emissionEntry :: String -> (String,Double) -> QDiagram Cairo V2 Double Any
emissionEntry emissiontype (symbol,emission)
  | emissiontype == "probability" = textentry
  | emissiontype == "score" = textentry
  | emissiontype == "bar" = barentry
  | otherwise = barentry
    where textentry = setLabel (symbol ++ " " ++ printf "%.3f" emission)
          barentry = setLabel symbol  ||| strutX 0.2 ||| bar emission

bar :: Double -> QDiagram Cairo V2 Double Any
bar emission = rect (4 * emission) 1 # lw 0 # fc black # translate (r2 (negate (2 - (4 * emission/2)),0)) <> rect 4 1 # lw 0.03

--makeSingleEmissionIndices index = V.fromList [PAI.Z  PAI.:. index PAI.:. A,PAI.Z  PAI.:. index PAI.:. U,PAI.Z  PAI.:. index PAI.:. G,PAI.Z  PAI.:. index PAI.:. C]

--makePairEmissionIndices cindex = V.fromList [PAI.Z  PAI.:. cindex PAI.:. A PAI.:. A,PAI.Z  PAI.:. cindex PAI.:. A PAI.:. U,PAI.Z  PAI.:. cindex PAI.:. A PAI.:. G,PAI.Z  PAI.:. cindex PAI.:. A PAI.:. C,PAI.Z  PAI.:. cindex PAI.:. U PAI.:. U,PAI.Z  PAI.:. cindex PAI.:. U PAI.:. A,PAI.Z  PAI.:. cindex PAI.:. U PAI.:. G,PAI.Z  PAI.:. cindex PAI.:. U PAI.:. C,PAI.Z  PAI.:. cindex PAI.:. G PAI.:. G,PAI.Z  PAI.:. cindex PAI.:. G PAI.:. A,PAI.Z  PAI.:. cindex PAI.:. G PAI.:. U,PAI.Z  PAI.:. cindex PAI.:. G PAI.:. C,PAI.Z  PAI.:. cindex PAI.:. C PAI.:. C,PAI.Z  PAI.:. cindex PAI.:. C PAI.:. A,PAI.Z  PAI.:. cindex PAI.:. C PAI.:. U,PAI.Z  PAI.:. cindex PAI.:. C PAI.:.G]

statebox :: Double -> Double -> String -> QDiagram Cairo V2 Double Any
statebox x y si = (rect 0.05 0.1 # lw 0 # named ("s" ++ si) ||| rect 1 0.1 # lw 0 # named ("a" ++ si) ||| rect 2 0.1 # lw 0 ||| rect 0.05 0.1 # lw 0 # named ("z" ++ si)) === rect x y  # lw 0.1 === rect 1 0 # lw 0 # named ("e" ++ si)

inverseStatebox :: Double -> Double -> String -> QDiagram Cairo V2 Double Any
inverseStatebox x y si = (rect 0.05 0.1 # lw 0 # named ("s" ++ si) ||| rect 1 0.1 # lw 0 # named ("a" ++ si) ||| rect 2 0.1 # lw 0  ||| rect 0.05 0.1 # lw 0 # named ("z" ++ si)) === rect x y  # lw 0.1 === rect 1 0 # lw 0 # named ("e" ++ si)

--scaling
-- | Specifies the size of the diagram. Absolute adapts to overall size according to subdiagrams
svgsize :: SizeSpec V2 Double
svgsize = mkSizeSpec2D Nothing Nothing

-- | Check for available cairo output formats
diagramName :: String -> String -> Either String String
diagramName filename fileformat
  | fileformat == "pdf" = Right (filename ++ "." ++ fileformat )
  | fileformat == "svg" = Right (filename ++ "." ++ fileformat )
  | fileformat == "png" = Right (filename ++ "." ++ fileformat )
  | fileformat == "ps" = Right (filename ++ "." ++ fileformat )
  | otherwise = Left "Unsupported output format requested (use svg, pdf, ps, png)"
                
printCM :: FilePath -> SizeSpec V2 Double -> QDiagram Cairo V2 Double Any -> IO ()
printCM outputName = renderCairo outputName

getBlankComparisonNodeLabels :: CM.CM -> V.Vector (Int, V.Vector (Colour Double))
getBlankComparisonNodeLabels model = comparisonNodeLabels
   where comparisonNodeLabels = V.generate (nodeNumber ) makeBlankComparisonNodeLabel
         nodeNumber = (CM._nodesInModel model) 

makeBlankComparisonNodeLabel :: Int ->  (Int,V.Vector (Colour Double))
makeBlankComparisonNodeLabel nodeNumber = (nodeNumber,V.singleton white)

getComparisonNodeLabels :: [CmcompareResult] -> V.Vector (String, Colour Double) -> CM.CM -> V.Vector (Int, V.Vector (Colour Double))
getComparisonNodeLabels comparsionResults colorVector model = comparisonNodeLabels
   where modelName = T.unpack (CM._name model)
         relevantComparisons1 = filter ((modelName==) . model1Name) comparsionResults
         modelNodeInterval1 = map (\a -> (model2Name a,model1matchednodes a))  relevantComparisons1 
         relevantComparisons2 = filter ((modelName==) . model2Name) comparsionResults
         modelNodeInterval2 = map (\a -> (model1Name a,model2matchednodes a))  relevantComparisons2
         modelNodeIntervals =  V.fromList (modelNodeInterval1 ++ modelNodeInterval2)
         colorNodeIntervals = V.map (modelToColor colorVector) modelNodeIntervals
         nodeNumber = (CM._nodesInModel model)
         comparisonNodeLabels = V.generate (nodeNumber ) (makeComparisonNodeLabel colorNodeIntervals)

modelToColor :: V.Vector (String,Colour Double) ->  (String,[Int]) -> (Colour Double,[Int])
modelToColor colorVector (mName,nInterval) = nColorInterval
  where nColorInterval = (snd (fromJust entry),nInterval)
        --nColorInterval = maybe Nothing (\a -> Just (snd a,nInterval)) entry
        entry = V.find (\(a,_)-> mName == a) colorVector

makeComparisonNodeLabel :: V.Vector (Colour Double,[Int]) -> Int -> (Int,V.Vector (Colour Double))
makeComparisonNodeLabel colorNodeIntervals nodeNumber = comparisonNodeLabel
  where relevantColorNodeIntervals = V.filter (\(_,b) -> elem nodeNumber b) colorNodeIntervals
        modelColors = V.map fst relevantColorNodeIntervals
        comparisonNodeLabel = if null modelColors then (nodeNumber,V.singleton white) else (nodeNumber,modelColors)

-- First colors are picked from http://colorbrewer2.org,  scheme. Comparisons up to
-- 4 models are colorblind safe
makeColorVector :: Int -> V.Vector (Colour Double)
makeColorVector modelNumber = V.take modelNumber colorVector
   where colorVector = V.fromList [R.sRGB24read "#a6cee3", R.sRGB24read "#1f78b4", R.sRGB24read "#b2df8a", R.sRGB24read "#33a02c", R.sRGB24read "#ffff99", R.sRGB24read "#cab2d6", R.sRGB24read "#6a3d9a", R.sRGB24read "#b15928", R.sRGB24read "#fb9a99", R.sRGB24read "#e31a1c", R.sRGB24read "#fdbf6f", R.sRGB24read "#ff7f00", moccasin, lime, seagreen, aqua ,darkorange ,blue, blueviolet ,brown ,burlywood ,cadetblue ,chartreuse ,chocolate ,coral ,cornflowerblue ,cornsilk ,cyan ,darkblue ,darkcyan ,darkgoldenrod ,darkgray ,darkgreen ,darkgrey ,darkkhaki ,darkmagenta ,darkolivegreen ,darkorchid ,darkred ,darksalmon ,darkseagreen ,darkslateblue ,darkslategray ,darkslategrey ,darkturquoise ,darkviolet ,deeppink ,deepskyblue ,dimgray ,dimgrey ,dodgerblue ,firebrick ,forestgreen ,fuchsia ,gainsboro ,gold ,goldenrod ,gray ,grey ,green ,greenyellow ,honeydew ,hotpink ,indianred,indigo ,ivory ,khaki ,lavender ,lavenderblush ,lawngreen ,lemonchiffon ,lime ,limegreen ,linen ,magenta ,maroon ,mediumaquamarine ,mediumblue ,mediumorchid ,mediumpurple ,mediumseagreen ,mediumslateblue ,mediumspringgreen ,mediumturquoise ,mediumvioletred ,midnightblue ,mintcream ,mistyrose ,navy ,oldlace ,olive ,olivedrab ,orange ,orangered ,orchid ,papayawhip ,peachpuff ,peru ,pink ,plum ,powderblue ,purple ,red ,rosybrown ,royalblue ,saddlebrown ,salmon ,sandybrown ,seagreen]

makeArrow :: ([Char], [Char], Double) -> QDiagram Cairo V2 Double Any -> QDiagram Cairo V2 Double Any
makeArrow (lab1,lab2,weight) = connectOutside' arrowStyle1 ("e" ++ lab1) ("a" ++ lab2)
  where arrowStyle1 = with & arrowHead .~ spike & shaftStyle %~ lw (local 0.1) & headLength .~ local 0.01 & shaftStyle %~ dashingG [weight, 0.1] 0 & headStyle %~ fc black . opacity (weight * 2)

makeSelfArrow :: ([Char], [Char], Double) -> QDiagram Cairo V2 Double Any -> QDiagram Cairo V2 Double Any  
makeSelfArrow (lab1,_,weight) = connectPerim' arrowStyle ("s" ++ lab1) ("z" ++ lab1) (5/12 @@ turn) (8/12 @@ turn)
  where arrowStyle = with  & arrowHead .~ spike & arrowShaft .~ shaft' & arrowTail .~ lineTail & tailTexture .~ solid black  & shaftStyle %~ lw (local 0.1) & headLength .~ local 0.01  & tailLength .~ 0 & shaftStyle %~ dashingG [weight, 0.3] 0 & headStyle %~ fc black . opacity (weight * 2)
        shaft' = wedge 3 xDir (2/4 @@ turn)

makeLabel :: (String, String, Double) -> QDiagram Cairo V2 Double Any -> QDiagram Cairo V2 Double Any
makeLabel (n1,n2,weight) =
  withName ("a" ++ n1) $ \b1 ->
  withName ("s" ++ n2) $ \b2 ->
    let v = location b2 .-. location b1
        midpoint = location b1 .+^ (v ^/ 2)
        (xOffset,yOffset) = setLabelOffset (location b1 ^. _x) (location b2 ^. _x) (location b1 ^. _y) (location b2 ^. _y)
      in
        Diagrams.Prelude.atop (position [(midpoint # translateX xOffset # translateY yOffset, setTransition ((show (roundPos 3 weight))))])
        --Diagrams.Prelude.atop (position [(midpoint # translateX xOffset # translateY yOffset, setTransition (lclass ++"," ++ (show (roundPos 3 weight))))])         
        --Diagrams.Prelude.atop (position [(midpoint # translateX xOffset # translateY yOffset, setTransition (n1 ++"," ++ n2 ++"," ++lclass ++"," ++ (show (roundPos 3 weight))))]) 

setLabelOffset :: Double -> Double -> Double -> Double -> (Double,Double)
setLabelOffset x1 x2 y1 y2
  -- 
  | ydiff < 2 = (0,0)
  | xdiff < 2 = (negate 1,negate 1)
    -- 
  | x1 > x2 && (ydiff > 30) = (negate 1,negate 10)
  -- 
  | x1 < x2 && (ydiff > 30) = (1,negate 10)             
  -- between split and insert state of same node - left upper(A)
  | x1 > x2 && (ydiff < 30) = (negate 1,negate 12)
  -- between split and insert state of same node - right upper (B)
  | x1 < x2 && (ydiff < 30) = (0,negate 12)
  -- between same split states of different nodes (C)
  | otherwise = (0,0)
    where ydiff = abs (abs y1 - abs y2)
          xdiff = abs (abs x1 - abs x2)
                
makeSelfLabel :: (String, String, Double) -> QDiagram Cairo V2 Double Any -> QDiagram Cairo V2 Double Any
makeSelfLabel (n1,_,weight)
  | weight == 0 = mempty
  | otherwise = withName ("e" ++ n1) $ \b1 ->
                  let midpoint = location b1
                  in
                    Diagrams.Prelude.atop (position [(midpoint # translateX (negate 0.25)  # translateY 22, setTransition (show (roundPos 3 weight)))])

roundPos :: (RealFrac a) => Int -> a -> a
roundPos positions number  = fromInteger (round $ number * (10^positions)) / (10.0^^positions)

bitScore2Double :: Bitscore -> Double
bitScore2Double (Bitscore x) = x