# Standard interface for PNA and MNA algorithms
using NetalignMeasures, NetalignUtils
using StatsBase

# rename nodes in aln bassed on nts and nodesets
function maln_mapped(aln::Vector,nts,nodesets)
    nodemap = Dict{String,String}()
    for i = 1:length(nts), j = 1:length(nodesets[i])
        nodemap[nodesets[i][j]] = nts[i].nodes[j]
    end
    alnout = Vector{Vector{String}}(length(aln))
    for i = 1:length(aln)
        vals = filter(x -> x!="-", aln[i])
        alnout[i] = String[nodemap[x] for x in vals] 
    end 
    alnout
end
function readmaln_mapped(alnfile::String,nts,nodesets)
    aln = open(alnfile,"r") do fd
        map(line->split(strip(line),[' ','\t','\n']),
            filter(line->!isempty(strip(line)),readlines(fd)))
    end
    maln_mapped(aln,nts,nodesets)
end

function readpaln_mapped(alnfile,nt1,nt2,nodes1,nodes2)
    aln = open(alnfile,"r") do fd
        map(line->split(strip(line),[' ','\t','\n']),
            filter(line->!isempty(strip(line)),readlines(fd)))
    end
    nodemap1 = Dict{String,String}()
    nodemap2 = Dict{String,String}()
    for j = 1:length(nodes1) nodemap1[nodes1[j]] = nt1.nodes[j] end
    for j = 1:length(nodes2) nodemap2[nodes2[j]] = nt2.nodes[j] end
    
    alnout = Vector{Vector{String}}(length(aln))
    for i = 1:length(aln)
        alnout[i] = String[nodemap1[aln[i][1]], nodemap2[aln[i][2]]]
    end 
    alnout
end

function write_wavecpp_network(fd::IO,G,nodes=0:size(G,1)-1)
    m, n = size(G)
    println(fd, n, " ", div(nnz(G),2))
    writeedgelist(fd,Network(G,Vector{String}(nodes)))
end
function write_wavecpp_network(file::String,G,nodes=0:size(G,1)-1)
    open(file,"w") do fd write_wavecpp_network(fd,G,nodes) end
end

function runwave(nt1::Network,nt2::Network,S::AbstractMatrix{Float64})
    G1 = nt1.G
    G2 = nt2.G
    nodes1 = map(string,randperm(size(G1,1))-1)
    nodes2 = map(string,randperm(size(G2,1))-1)
    file1 = "a.txt"
    file2 = "b.txt"
    simfile = "a-b.sim"
    alnfile = "aln.out"
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do 
        write_wavecpp_network(file1,G1,nodes1)
        write_wavecpp_network(file2,G2,nodes2)
        writelistmat(simfile,S,nodes1,nodes2)
        waveapp = joinpath(curdirname(),"namethods/WAVE/WAVE")
        cmd = `$waveapp $file1 $file2 $simfile $alnfile`
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readpaln_mapped(alnfile,nt1,nt2,nodes1,nodes2)
        aln,runtime
    end
    #writescorefile("wave.sim",S,nodes1,nodes2,:evals)        
    rm(workdir,recursive=true)
    aln,runtime
end

function runwavetop(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix)
    G1 = nt1.G
    G2 = nt2.G
    T = NodeSimMeasure(:hgraalgdvs, G1,G2,gdv1,gdv2).S
    runwave(nt1,nt2,T)
end
function runwaveseq(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix,
                    S::AbstractMatrix{Float64})
    G1 = nt1.G
    G2 = nt2.G
    T = NodeSimMeasure(:hgraalgdvs, G1,G2,gdv1,gdv2).S
    runwave(nt1, nt2, 0.5*T+0.5*S)
    #runwave(nt1, nt2, S)
end

function runghost(nt1::Network,nt2::Network,S::AbstractMatrix{Float64}=zeros(Float64,0,0))
    dosim = !isempty(S)
    nprocs = 8
    G1 = nt1.G
    G2 = nt2.G
    nodes1 = map(string,randperm(size(G1,1))-1)
    nodes2 = map(string,randperm(size(G2,1))-1)
    file1 = "a.txt"
    file2 = "b.txt"
    simfile = "a-b.sim"
    gfile1 = "a.gexf"
    gfile2 = "b.gexf"
    cfgfile = "a-b.cfg"
    if dosim simfile = tempname() end
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do 
        writeedgelist(file1, Network(G1,nodes1))
        writeedgelist(file2, Network(G2,nodes2))
        a2g = joinpath(curdirname(),"namethods/GHOST/adj_to_gexf.py")
        run(`python $a2g $file1 $gfile1`)
        run(`python $a2g $file2 $gfile2`)
        if dosim writelistmat(simfile,S,nodes1,nodes2) end # have to be BLAST e-values   
        ghostapp = joinpath(curdirname(),"namethods/GHOST/GHOST")
        open(cfgfile,"w") do fd
            println(fd,"[main]")
            println(fd,"network1: $gfile1")
            println(fd,"network2: $gfile2")
            if dosim
                println(fd,"alpha: 0.5")
                println(fd,"sequencescores: $simfile")
                println(fd,"beta: 1e-10")
            end
            println(fd,"nneighbors: 4")
            println(fd,"processors: $nprocs")
        end
        cmd = `$ghostapp -c $cfgfile`
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readpaln_mapped("a_vs_b.af",nt1,nt2,nodes1,nodes2)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end
function runghosttop(nt1::Network,nt2::Network)
    runghost(nt1,nt2)
end
function runghostseq(nt1::Network,nt2::Network,S::AbstractMatrix{Float64})
    runghost(nt1,nt2,S)
end

function runmagna(nt1::Network,nt2::Network,S::AbstractMatrix{Float64},alpha::Float64)
    G1 = nt1.G
    G2 = nt2.G
    nprocs = 8
    nodes1 = map(string,randperm(size(G1,1))-1)
    nodes2 = map(string,randperm(size(G2,1))-1)
    file1 = "a.gw"
    file2 = "b.gw"
    simfile = "a-b.evals"
    alnfile = "aln.txt"
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do 
        writegw(file1, Network(G1,nodes1))
        writegw(file2, Network(G2,nodes2))
        writelistmat(simfile,S,nodes1,nodes2)
        magnaapp = joinpath(curdirname(),"namethods/MAGNA++/magna")
        cmd = `$magnaapp -G $file1 -H $file2 -d $simfile -o $alnfile -a $alpha -p 15000 -n 10000 -f 100 -t 8 -m S3`
        #cmd = `$magnaapp -G $file1 -H $file2 -d $simfile -o $alnfile -a $alpha -p 15 -n 100 -f 100 -t 8 -m S3`
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readpaln_mapped(string(alnfile,"_final_alignment.txt"),nt1,nt2,nodes1,nodes2)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end
function runmagnatop(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix)
    G1 = nt1.G
    G2 = nt2.G
    T = NodeSimMeasure(:hgraalgdvs, G1,G2,gdv1,gdv2).S
    runmagna(nt1,nt2,T,0.5)
end
function runmagnaseq(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix,S::AbstractMatrix{Float64})
    G1 = nt1.G
    G2 = nt2.G
    T = NodeSimMeasure(:hgraalgdvs, G1,G2,gdv1,gdv2).S
    runmagna(nt1, nt2, 0.5*S+0.25*T, 0.25)
end

function runlgraal(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix,S::AbstractMatrix{Float64}=zeros(Float64,0,0))
    dosim = !isempty(S)
    G1 = nt1.G
    G2 = nt2.G
    nodes1 = map(string,randperm(size(G1,1))-1)
    nodes2 = map(string,randperm(size(G2,1))-1)
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do 
        file1 = joinpath(workdir,"a.gw")
        file2 = joinpath(workdir,"b.gw")
        gfile1 = joinpath(workdir,"a.gdv")
        gfile2 = joinpath(workdir,"b.gdv")
        simfile = joinpath(workdir,"a-b.sim")
        alnfile = joinpath(workdir,"aln.txt")
        writegw(file1, Network(G1,nodes1))
        writegw(file2, Network(G2,nodes2))
        writegdv(gfile1,gdv1,nodes1)
        writegdv(gfile2,gdv2,nodes2)
        if dosim
            writelistmat(simfile,S,nodes1,nodes2) # (1-evalues, or -log(evalues), or bitscores)
        else
            open(print,simfile,"w")
        end
        lgraalapp = joinpath(curdirname(),"namethods/L-GRAAL/L-GRAAL")
        if dosim
            cmd = `$lgraalapp -Q $file1 -T $file2 -q $gfile1 -t $gfile2 -B $simfile -o $alnfile -a 0.5`
        else
            cmd = `$lgraalapp -Q $file1 -T $file2 -q $gfile1 -t $gfile2 -B $simfile -o $alnfile -a 0.0`
        end
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readpaln_mapped(alnfile,nt1,nt2,nodes1,nodes2)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end
function runlgraaltop(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix)
    runlgraal(nt1,nt2,gdv1,gdv2)
end
function runlgraalseq(nt1::Network,nt2::Network,gdv1::Matrix,gdv2::Matrix,S::AbstractMatrix{Float64})
    runlgraal(nt1,nt2,gdv1,gdv2,S)
end

function runisorankn(nts::Vector{Network},Ss::Matrix,alpha::Float64)
    Rs = Matrix{AbstractMatrix}(size(Ss))
    maxveclen = 5000000
    scorefile = "scores.evals"
    for i = 1:length(nts), j = i+1:length(nts)
        S = copy(Ss[i,j])
        vcutoff = select(vec(S), min(maxveclen,length(S)), rev=true)
        S[S .< vcutoff] = 0
        println("vcutoff $i $j: $vcutoff")
        Rs[i,j] = sparse(S)
    end
    println("Starting IsoRankN setup")
    isoranknapp = joinpath(curdirname(),"namethods/IsoRankN/isorank-n-v3-64")
    nodesets = Vector(length(nts))
    files = Vector(length(nts))
    c = 0
    for i = 1:length(nts)
        nodesets[i] = map(string, c + randperm(size(nts[i].G,1))-1)
        c += size(nts[i].G,1)
        files[i] = string(('a':'z')[i], ".tab") 
    end
    simfiles = Matrix(size(Rs)...)
    for i = 1:length(nts), j = i+1:length(nts)
        simfiles[i,j] = string(('a':'z')[i], '-', ('a':'z')[j], ".evals")
    end
    alnfile = "aln.out"
    cfgfile = "aln.cfg"
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do
        println("Writing to temporary files")
        for i = 1:length(nts)
            open(files[i],"w") do fd
                println(fd, "INTERACTOR_A INTERACTOR_B")
                writeedgelist(fd, Network(nts[i].G,nodesets[i]))
            end
        end
        for i = 1:length(nts), j = i+1:length(nts)
            writelistmat(simfiles[i,j], Rs[i,j], nodesets[i], nodesets[j])
        end        
        open(cfgfile,"w") do fd
            println(fd, ".\n-\n$(length(nts))")
            for i = 1:length(nts) println(fd, ('a':'z')[i]) end
        end
        if isapprox(alpha, 0.0, atol=1e-5)
            open(scorefile, "w") do fd
                for i = 1:length(nts), j = i+1:length(nts)
                    writelistmat(fd, Rs[i,j], nodesets[i], nodesets[j])
                end
            end
            cmd = `$isoranknapp --scorefile $scorefile --K 3 --thresh 1e-4 --alpha $alpha --maxveclen $maxveclen --o $alnfile $cfgfile`
        else
            cmd = `$isoranknapp --K 20 --thresh 1e-4 --alpha $alpha --maxveclen $maxveclen --o $alnfile $cfgfile`
        end
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readmaln_mapped(alnfile,nts,nodesets)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end

function runisorankntop(nts::Vector{Network},gdvs::Vector{Matrix{Float64}})
    Cs = Matrix{AbstractMatrix{Float64}}(length(nts),length(nts))
    for i = 1:length(nts), j = i+1:length(nts)
        #Cs[i,j] = hgraalgdvs(nts[i].G, nts[j].G, gdvs[i], gdvs[j])
        Cs[i,j] = ones(Float64,size(nts[i].G,1),size(nts[j].G,1))
    end
    runisorankn(nts,Cs,1.0)
end    
function runisoranknseq(nts::Vector{Network},Ss::Matrix)
    runisorankn(nts,Ss,0.5)
end 

function migraalconf(Xs::AbstractMatrix...)
    C = zeros(Float64,size(Xs[1]))
    for X in Xs
        for i = 1:size(X,1)
            C[i,:] += competerank(X[i,:]) ./ size(X,2)
        end
    end
    C ./= length(Xs)
    C
end

function runmiisotop(nts::Vector{Network}, gdvs::Vector)
    Cs = Matrix{AbstractMatrix{Float64}}(length(nts),length(nts))
    for i = 1:length(nts), j = i+1:length(nts)
        Cs[i,j] = migraalconf(NodeSimMeasure(:hgraalgdvs, nts[i].G, nts[j].G, gdvs[i], gdvs[j]).S)
        #Cs[i,j] = NodeSimMeasure(:hgraalgdvs, nts[i].G, nts[j].G, gdvs[i], gdvs[j]).S
    end
    runisorankn(nts,Cs,0.0)
end
function runmiisoseq(nts::Vector{Network}, gdvs::Vector, Ss::Matrix)
    Cs = Matrix{AbstractMatrix{Float64}}(size(Ss)...)
    for i = 1:length(nts), j = i+1:length(nts)
        Cs[i,j] = migraalconf(Ss[i,j], NodeSimMeasure(:hgraalgdvs, nts[i].G, nts[j].G, gdvs[i], gdvs[j]).S)
        #Cs[i,j] = 0.5 .* Ss[i,j] .+ 0.5 .* NodeSimMeasure(:hgraalgdvs, nts[i].G, nts[j].G, gdvs[i], gdvs[j]).S
    end
    runisorankn(nts,Cs,0.0)
end

function runmultimagna(nts::Vector{Network},Ss=Matrix(), alpha::Float64=1.0)
    dosim = !isempty(Ss)
    magnaapp = joinpath(curdirname(),"namethods/multiMAGNA++/magna")
    nodesets = Vector(length(nts))
    files = Vector(length(nts))
    c = 0
    for i = 1:length(nts)
        nodesets[i] = map(string, c + randperm(size(nts[i].G,1))-1)
        c += size(nts[i].G,1)
        files[i] = string(('a':'z')[i], ".txt")        
    end
    if dosim
        simfiles = Matrix(size(Ss)...)
        for i = 1:length(nts), j = i+1:length(nts)
            simfiles[i,j] = string(('a':'z')[i], '-', ('a':'z')[j], ".evals")
        end
    end
    alnfile = "aln.out"
    cfgfile = "aln.cfg"
    cfgsimfile = "alnsim.cfg"
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do
        for i = 1:length(nts)
            open(files[i],"w") do fd
                writeedgelist(fd, Network(nts[i].G,nodesets[i]))
            end
        end
        if dosim
            for i = 1:length(nts), j = i+1:length(nts)
                writelistmat(simfiles[i,j], Ss[i,j], nodesets[i], nodesets[j]) # have to be BLAST bitscores
            end
        end
        open(cfgfile,"w") do fd
            for i = 1:length(nts) println(fd, files[i]) end
        end
        if dosim
            open(cfgsimfile,"w") do fd
                for i = 1:length(nts), j = i+1:length(nts) println(fd, simfiles[i,j]) end
            end
        end
        npop = 15000; ngen = 100000
        #npop = 15; ngen = 100
        if dosim
            cmd = `$magnaapp -G $cfgfile -d $cfgsimfile -o $alnfile -a $alpha -m CIQ -p $npop -n $ngen -f 100 -t 8`
        else
            cmd = `$magnaapp -G $cfgfile -o $alnfile -m CIQ -p $npop -n $ngen -f 100 -t 8`
        end
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readmaln_mapped("$(alnfile)_final_alignment.txt",nts,nodesets)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end
function runmultimagnatop(nts::Vector{Network},gdvs::Vector{Matrix{Float64}})
    Cs = Matrix{AbstractMatrix{Float64}}(length(nts),length(nts))
    for i = 1:length(nts), j = i+1:length(nts)
        Cs[i,j] = NodeSimMeasure(:hgraalgdvs, nts[i].G, nts[j].G, gdvs[i], gdvs[j]).S
    end
    runmultimagna(nts,Cs,0.5)
end
function runmultimagnaseq(nts::Vector{Network},gdvs::Vector{Matrix{Float64}},
                          Ss::Matrix)
    Cs = Matrix{AbstractMatrix{Float64}}(size(Ss)...)
    for i = 1:length(nts), j = i+1:length(nts)
        T = NodeSimMeasure(:hgraalgdvs, nts[i].G, nts[j].G, gdvs[i], gdvs[j]).S
        Cs[i,j] = (2/3) .* Ss[i,j] .+ (1/3) .* T
    end
    runmultimagna(nts,Cs,0.25)
end

function rungedevom(nts::Vector{Network}, gdvs=Vector{Matrix}(), Ss=Matrix())
    dogdvs = !isempty(gdvs)
    alnapp = joinpath(curdirname(),"namethods/GEDEVO-M/gedevo-m")
    nodesets = Vector(length(nts))
    files = Vector(length(nts))
    gfiles = Vector(length(nts))
    c = 0
    for i = 1:length(nts)
        nodesets[i] = map(string, c + randperm(size(nts[i].G,1))-1)
        c += size(nts[i].G,1)
        files[i] = string(('a':'z')[i], ".txt")
        gfiles[i] = string(('a':'z')[i], ".gdv")        
    end
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do
        for i = 1:length(nts)
            writeedgelist(files[i], Network(nts[i].G,nodesets[i]))
            if dogdvs writegdv(gfiles[i], gdvs[i], nodesets[i]) end
        end
        grps = ('a':'z')[1:length(nts)]
        strfiles = vcat(map(i -> ["--edgelist", string(grps[i]), files[i]], 1:length(nts))...)
        strgfiles = dogdvs ? vcat(map(i -> ["--grsig", string(grps[i]), gfiles[i]], 1:length(nts))...) : []
        maxiter = 100000
        cmd = `$alnapp --grnum $(length(nts)) --groups $grps $strfiles $strgfiles --pop 500 --maxsame 3000 --maxiter $maxiter --no-prematch --no-workfiles --undirected ---threads 8 --no-save`
        println(cmd)
        tic(); run(cmd); runtime = toq()
        alnfile = filter(x -> splitext(x)[2]==".mmapping", readdir())[1]
        aln = open(alnfile,"r") do fd
            lines = filter(line -> !(line[1]=='#'), map(strip,readlines(fd)))
            aln = map(lines) do line
                filter(val -> val!="-", split(line, [' ','\t','\n'])[1:length(nts)])
            end
            filter(x->length(x)>1, aln)            
        end
        aln = maln_mapped(aln,nts,nodesets)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end

function runfuse(nts::Vector{Network}, Ss::Matrix)
    fusedir = joinpath(curdirname(),"namethods/FUSE")
    fuseapp = joinpath(fusedir,"FUSE.exe")
    nodesets = Vector(length(nts))
    gwfiles = Vector(length(nts))
    files = Vector(length(nts))
    c = 0
    for i = 1:length(nts)
        nodesets[i] = map(string, c + randperm(size(nts[i].G,1))-1)
        c += size(nts[i].G,1)
        files[i] = string(('a':'z')[i], ".fuselist")        
        gwfiles[i] = string(('a':'z')[i], ".gw")        
    end
    alnfile = "aln.out"
    simfile = "blast.fuselst"
    nmtffile = "nmtf.fuselist"
    cfgfile = "aln.cfg"
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do
        for i = 1:length(nts)
            nt = Network(nts[i].G,nodesets[i])
            writeedgelist(files[i], nt, suffix=" 1.0")
            writegw(gwfiles[i], nt)
        end
        open(simfile,"w") do fd
            println(fd, "Gene1 Gene2 Association_strength")
            for i = 1:length(nts), j = i+1:length(nts)
                writelistmat(fd, Ss[i,j], nodesets[i], nodesets[j]) # have to be BLAST bitscores
            end
        end
        niter = 1000
        nthreads = 8
        matlabcmd = string("maxNumCompThreads($nthreads); addpath('$(fusedir)'); run_nmtf({",
                           join(map(x->"'$(x)'",files),", "),
                           "}, '$(simfile)', ",
                           "[",join(fill(80,length(nts)),", "),"], $niter, 0.01, ",
                           "'$(nmtffile)'); exit")
        open(cfgfile,"w") do fd
            for i = 1:length(nts) println(fd, gwfiles[i]) end
        end
        cmd = `$fuseapp -a 0.5 -n $cfgfile -s $simfile -t $nmtffile -o $alnfile`
        println(matlabcmd)
        println(cmd)
        tic(); run(`matlab -nodisplay -nojvm -r $matlabcmd`); run(cmd); runtime = toq()
        aln = open(alnfile,"r") do fd
            map(line->split(strip(line),[' ','\t','\n'])[3:end],
                filter(line->!isempty(strip(line)),readlines(fd)))
        end
        aln = maln_mapped(aln,nts,nodesets)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end

function runbeams(nts::Vector{Network},Ss::Matrix)
    beamsapp = joinpath(curdirname(),"namethods/BEAMS/beams")
    nodesets = Vector(length(nts))
    files = Vector(length(nts))
    c = 0
    for i = 1:length(nts)
        nodesets[i] = map(string, c + randperm(size(nts[i].G,1))-1)
        c += size(nts[i].G,1)
        files[i] = string(('a':'z')[i], ".tab")        
    end
    simfiles = Matrix(size(Ss)...)
    for i = 1:length(nts), j = i+1:length(nts)
        simfiles[i,j] = string(('a':'z')[i], '-', ('a':'z')[j], ".evals")
    end
    alnfile = "aln.out"
    cfgfile = "aln.cfg"
    workdir = tempname(); println("Work dir: $workdir"); mkdir(workdir)
    aln,runtime = cd(workdir) do
        for i = 1:length(nts)
            open(files[i],"w") do fd
                println(fd, "INTERACTOR_A INTERACTOR_B")
                writeedgelist(fd, Network(nts[i].G,nodesets[i]))
            end
        end
        for i = 1:length(nts), j = i+1:length(nts)
            writelistmat(simfiles[i,j], Ss[i,j], nodesets[i], nodesets[j]) # have to be BLAST bitscores
        end
        open("empty.evals","w") do fd print(fd,"") end
        open(cfgfile,"w") do fd
            println(fd, "$(length(nts))")
            for i = 1:length(nts) println(fd, files[i]) end
            for i = 1:length(nts) println(fd, "empty.evals") end
            for i = 1:length(nts), j = i+1:length(nts) println(fd, simfiles[i,j]) end
        end
        cmd = `$beamsapp $cfgfile 0.5 0.4 aln.backbones $alnfile`
        println(cmd)
        tic(); run(cmd); runtime = toq()
        aln = readmaln_mapped(alnfile,nts,nodesets)
        aln,runtime
    end
    rm(workdir,recursive=true)
    aln,runtime
end

function align(algo::String,
                      nts::Vector{Network}, gdvfiles::Vector, evalfiles::Matrix)
    if startswith(algo,"PW_") algo = algo[4:end] end

    if algo in ["LGRAAL_top","LGRAAL_seq"]
        p = sortperm(nts, by=x -> nnz(x.G))
    else
        p = sortperm(nts, by=x -> size(x.G,1))
    end
    nts = nts[p]
    gdvfiles = gdvfiles[p]
    evalfiles = evalfiles[p,p]
    
    if algo in ["WAVE_top","MAGNA_top","LGRAAL_top",
                "WAVE_seq","MAGNA_seq","LGRAAL_seq",
                "IsoRankN","MI-Iso_top","MI-IsoN","MI-Iso_seq","MI-IsoN_seq",
                "MAGNA-M_top","MAGNA-M_plusgdvs","MAGNA-M_seq","MAGNA-M_graalnlogevals",
                "GEDEVO-M_plusgdvs"]
        println("Reading GDVs")
        gdvs = readgdvfiles(nts,gdvfiles)
    end
    if algo in ["WAVE_seq","GHOST_seq","MAGNA_seq","LGRAAL_seq",
                "IsoRankN_a0.5","BEAMS","FUSE","MI-Iso_seq","MI-IsoN_seq",
                "MAGNA-M_seq","MAGNA-M_graalnlogevals"]
        println("Reading E-values")
        simmats = readevaltosimfiles(nts,evalfiles)
    end
    println("Starting aligner")
    aln,runtime =
        if algo == "WAVE_top"
            runwavetop(nts[1], nts[2], gdvs[1], gdvs[2])
        elseif algo == "WAVE_seq"
            runwaveseq(nts[1], nts[2], gdvs[1], gdvs[2], simmats[1,2])
        elseif algo == "GHOST_top"
            runghost(nts[1], nts[2])
        elseif algo == "GHOST_seq"
            runghost(nts[1], nts[2],simmats[1,2])
        elseif algo == "MAGNA_top"
            runmagnatop(nts[1], nts[2], gdvs[1], gdvs[2])
        elseif algo == "MAGNA_seq"
            runmagnaseq(nts[1], nts[2], gdvs[1], gdvs[2], simmats[1,2])
        elseif algo == "LGRAAL_top"
            runlgraaltop(nts[1], nts[2], gdvs[1], gdvs[2])
        elseif algo == "LGRAAL_seq"
            runlgraalseq(nts[1], nts[2], gdvs[1], gdvs[2], simmats[1,2])
        elseif algo == "IsoRankN"
            runisorankntop(nts,gdvs)
        elseif algo == "IsoRankN_a0.5"
            runisoranknseq(nts,simmats)
        elseif algo in ["MI-Iso_top","MI-IsoN"]
            runmiisotop(nts,gdvs)
        elseif algo in ["MI-Iso_seq","MI-IsoN_seq"]
            runmiisoseq(nts,gdvs,simmats)
        elseif algo == "BEAMS"
            runbeams(nts,simmats)
        elseif algo in ["MAGNA-M_top","MAGNA-M_plusgdvs"]
            runmultimagnatop(nts,gdvs)
        elseif algo in ["MAGNA-M_seq","MAGNA-M_graalnlogevals"]
            runmultimagnaseq(nts,gdvs,simmats)
        elseif algo == "GEDEVO-M"
            rungedevom(nts)
        elseif algo == "GEDEVO-M_plusgdvs"
            rungedevom(nts,gdvs)
        elseif algo == "FUSE"
            runfuse(nts,simmats)
        end
    aln,runtime
end

function readevaltosimfiles(nts,simfiles,normalize=false)
    k = length(nts)
    Es = Matrix{AbstractMatrix{Float64}}(k,k)
    for i = 1:k, j = i+1:k
        simfilex,simfiley,colnum,issparse,nodenameskip = simfiles[i,j]
        gnodes = map(x->x[nodenameskip+1:end], nts[i].nodes)
        hnodes = map(x->x[nodenameskip+1:end], nts[j].nodes)
        if isfile(simfilex)
            Es[i,j] = NodeSimMeasure(:evalues, readevalues(simfilex, gnodes, hnodes)).S
        elseif isfile(simfiley)
            Es[i,j] = NodeSimMeasure(:evalues, readevalues(simfiley, hnodes, gnodes)').S
        else
            error("Files ", simfilex, " and ", simfiley, " do not exist")
        end
    end
    Es
end

function readgdvfiles(nts,gdvfiles)
    map(nts,gdvfiles) do nt,file
        readgdv(file, nt.nodes)
    end
end

function NetalignUtils.writelistmat(fd::IO, S::Matrix, nodes1::Vector{<:AbstractString},
                      nodes2::Vector{<:AbstractString})
    for i = 1:size(S,1), j = 1:size(S,2)
        println(fd, "$(nodes1[i]) $(nodes2[j]) $(S[i,j])")
    end
end

function NetalignUtils.writelistmat(fd::IO, S::Matrix, nodes1::Vector{<:AbstractString},
                      nodes2::Vector{<:AbstractString}, pred = _ -> true)
    for i = 1:size(S,1), j = 1:size(S,2)
        if pred(S[i,j])
            println(fd, "$(nodes1[i]) $(nodes2[j]) $(S[i,j])")
        end
    end
end

function nu_writelistmat(fd::IO, S::SparseMatrixCSC, nodes1::Vector{<:AbstractString},
                         nodes2::Vector{<:AbstractString}, pred = _ -> true)
    I,J,V = findnz(S)
    for (i,j,v) in zip(I,J,V)
        if pred(v)
            println(fd, "$(nodes1[i]) $(nodes2[j]) $v")
        end
    end
end
nu_writelistmat(file::AbstractString, args...;vargs...) =
    open(fd -> nu_writelistmat(fd, args...;vargs...), file, "w")

curdirname() = if isinteractive() "." else dirname(Base.source_path()) end

function writemaln(fd::IO,aln)
    for i = 1:length(aln)
        for j = 1:length(aln[i])
            print(fd, aln[i][j])
            if j < length(aln[i]) print(fd," ") else println(fd) end
        end
    end
end
writemaln(file::String,aln) = open(fd -> writemaln(fd,aln), file, "w")

function readgwlist(nwlistfile)
    nts = Network[]
    open(nwlistfile,"r") do fd
	for line = eachline(fd)
	    nwfile = joinpath(dirname(nwlistfile), strip(line))
            push!(nts, readgw(nwfile))
        end
    end
    return nts
end

function getevalfiles(df,netsgroup)
    simfilemap = Dict()
    for i = 1:size(df,1), j = 1:size(df,1)
        g = df[i,:network_file]
        h = df[j,:network_file]
        # simfile,format,issparse,nodenameskip
        simfilemap[(g,h)]  = (joinpath(curdirname(),"../data/evals",
                                       string(df[i,:evalue_file_prefix],"_",df[j,:evalue_file_prefix],"_",netsgroup,df[i,:evalue_file_suffix])),
                              joinpath(curdirname(),"../data/evals",
                                       string(df[j,:evalue_file_prefix],"_",df[i,:evalue_file_prefix],"_",netsgroup,df[i,:evalue_file_suffix])),
                              df[i,:evalue_colnum], true, df[i,:node_name_skip])
    end
    simfilemap
end

function getgdvfiles(df)
    gdvfilemap = Dict()
    k = size(df,1)
    for i = 1:k
        file = df[i,:network_file]
        gdvfilemap[file] = joinpath(curdirname(),"../data",df[i,:gdv_file])
    end
    gdvfilemap
end   
