text:
#for the sampling
include("Mantzaris.jl")
#create csv of the graphs
include("dot2csvNetTotalCollusion.jl")

#needed for the plots
using Plots
#use pyplot()
pyplot()

function demo(stYr=1975,endYr=1985,windowSize=5)
    southWest = ["Portugal","Spain","Malta","SanMarino","Andorra","Monaco","Morocco","Italy"]
    north = ["Iceland","Denmark","Norway","Sweden","Finland"]
    northWest = ["UnitedKingdom","Ireland","Belgium","France","Luxembourg"]
    central = ["Germany","Austria","TheNetherlands","Switzerland","Slovenia","CzechRepublic","Hungary"]
    southEast = ["Greece","Montenegro","Cyprus","Albania","Bulgaria","Croatia","BosniaHerzegovina","Turkey","FYRMacedonia","Romania","Serbia","Israel","Yugoslavia"]
    east = ["Russia","Ukraine","Moldova","Belarus","Poland","Georgia","Armenia","Azerbaijan","Estonia","Lithuania","Latvia"]
    yr=stYr
    collusionDict = Mantzaris(stYr,endYr,windowSize)
    print("\n");print(collusionDict);print("\n")
    ratioTwoWay = []
    totalTwoWay = []
    totalOneWay = []
    while( (yr+windowSize) <= endYr )
        oneWays = collusionDict["1way:$(yr)-$(yr + windowSize)"]
        twoWays = collusionDict["2way:$(yr)-$(yr + windowSize)"]
        network = "digraph{"
        yrs = string("years","$(yr)","to","$(yr + windowSize)")# to "# $(yr + windowSize)"
        network = string(network,string(" graph [label=","$(yrs)",", fontsize=34]; "))

        seenC = []
        for pair in oneWays
            pair = replace(pair,"&","")
            pair = replace(pair,".","")
            pairTmp = split(pair,'-')
            pairTmp[1] =  replace(pairTmp[1]," ","")
            pairTmp[2] =  replace(pairTmp[2]," ","")
            if(!(pairTmp[1] in seenC))
                if(pairTmp[1] in east)
                    network = string(network,string(pairTmp[1]," [style=filled,fillcolor=olivedrab3]; "))

                elseif(pairTmp[1] in southWest)
                    network = string(network,string(pairTmp[1]," [style=filled,fillcolor=indianred1]; "))
                elseif(pairTmp[1] in north)
                    network = string(network,string(pairTmp[1]," [style=filled,fillcolor=dodgerblue]; "))
                elseif(pairTmp[1] in northWest)
                    network = string(network,string(pairTmp[1]," [style=filled,fillcolor=darkslategray2]; "))
                elseif(pairTmp[1] in central)
                    network = string(network,string(pairTmp[1]," [style=filled,fillcolor=gray]; "))
                elseif(pairTmp[1] in southEast)
                    network = string(network,string(pairTmp[1]," [style=filled,fillcolor=darkgoldenrod1]; "))
                end
                seenC = [seenC,pairTmp[1]]
            end
            if(!(pairTmp[2] in seenC))
                if(pairTmp[2] in east)
                    network = string(network,string(pairTmp[2]," [style=filled,fillcolor=olivedrab3]; "))
                elseif(pairTmp[2] in southWest)
                    network = string(network,string(pairTmp[2]," [style=filled,fillcolor=indianred1]; "))
                elseif(pairTmp[2] in north)
                    network = string(network,string(pairTmp[2]," [style=filled,fillcolor=dodgerblue]; "))
                elseif(pairTmp[2] in northWest)
                    network = string(network,string(pairTmp[2]," [style=filled,fillcolor=darkslategray2]; "))
                elseif(pairTmp[2] in central)
                    network = string(network,string(pairTmp[2]," [style=filled,fillcolor=gray]; "))
                elseif(pairTmp[2] in southEast)
                    network = string(network,string(pairTmp[2]," [style=filled,fillcolor=darkgoldenrod1]; "))
                end
                seenC = [seenC,pairTmp[2]]
            end
        end

        collusion1 = []
        for (ind,pair) in enumerate(twoWays)
            print(pair);print("\n")
            pair = replace(pair,"&","")
            pair = replace(pair,".","")
            pairTmp = split(pair,'-')
            pairTmp[1] =  replace(pairTmp[1]," ","")
            pairTmp[2] =  replace(pairTmp[2]," ","")
            edge = string(pairTmp[1],"->",pairTmp[2]," [dir=both color=red penwidth=3];")
            network = string(network,edge)
            append!(collusion1,[string(pairTmp[1],pairTmp[2])])
            append!(collusion1,[string(pairTmp[2],pairTmp[1])])
        end
        print(collusion1)
        for pair in oneWays
            pair = replace(pair,"&","")
            pair = replace(pair,".","")
            pairTmp = split(pair,'-')
            pairTmp[1] =  replace(pairTmp[1]," ","")
            pairTmp[2] =  replace(pairTmp[2]," ","")
            edge = string(pairTmp[1],"->",pairTmp[2],";")
            if(!( string(pairTmp[1],pairTmp[2]) in collusion1))
                network = string(network,edge)
            end

        end
        network = string(network,"}")
        fileName = string("network",":$(yr)-$(yr + windowSize)",".dot")
        print(fileName)
        filePNG = string("network",":$(yr)-$(yr + windowSize)",".png")
        print(filePNG)
        writedlm(string("./",fileName), [network])

        run(`dot $fileName -Tpng -o $filePNG`)
        print(network)
        yr = yr + windowSize

        append!(ratioTwoWay,[length(twoWays) / (length(oneWays) - length(twoWays))])
        append!(totalOneWay,length(oneWays)-length(twoWays))
        append!(totalTwoWay,length(twoWays))
    end

yr = stYr
    repeatWinds = Dict()
    while( (yr+windowSize) <= endYr )
        oneWays = collusionDict["1way:$(yr)-$(yr + windowSize)"]
        twoWays = collusionDict["2way:$(yr)-$(yr + windowSize)"]
        for pair in oneWays
            pair = replace(pair,"&","")
            pair = replace(pair,".","")
            pairTmp = split(pair,'-')
            pairTmp[1] =  replace(pairTmp[1]," ","")
            pairTmp[2] =  replace(pairTmp[2]," ","")
            edgeName = string(pairTmp[1],"->",pairTmp[2])
            print(edgeName)
            if(haskey(repeatWinds, string(pairTmp[1],"->",pairTmp[2]) ))
                repeatWinds[string(pairTmp[1],"->",pairTmp[2])] = 1 + repeatWinds[string(pairTmp[1],"->",pairTmp[2])]
            else
                repeatWinds[string(pairTmp[1],"->",pairTmp[2])] = 1
            end
        end
        yr = yr + windowSize
    end
    print("\n the aggregrate of repeated edges \n")
print(repeatWinds)
#####
network = "digraph{"
yrs = string("years","$(stYr)","to","$(endYr)")

network = string(network,string(" graph [label=","$(yrs)",", fontsize=34]; "))
seenC = []
repeatWinds2 = repeatWinds
for entry in repeatWinds2
    keyTmp = entry[1]
    weight = entry[2]
    pairTmp = split(keyTmp,"->")

    if(!(pairTmp[1] in seenC))
        if(pairTmp[1] in east)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=olivedrab3]; "))
        elseif(pairTmp[1] in southWest)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=indianred1]; "))

        elseif(pairTmp[1] in north)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=dodgerblue]; "))
        elseif(pairTmp[1] in northWest)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=darkslategray2]; "))
        elseif(pairTmp[1] in central)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=gray]; "))
        elseif(pairTmp[1] in southEast)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=darkgoldenrod1]; "))
        end
        seenC = [seenC,pairTmp[1]]
    end

    if(!(pairTmp[2] in seenC))
       if(pairTmp[2] in east)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=olivedrab3]; "))

       elseif(pairTmp[2] in southWest)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=indianred1]; "))
       elseif(pairTmp[2] in north)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=dodgerblue]; "))
       elseif(pairTmp[2] in northWest)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=darkslategray2]; "))

       elseif(pairTmp[2] in central)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=gray]; "))
       elseif(pairTmp[2] in southEast)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=darkgoldenrod1]; "))
       end
       seenC = [seenC,pairTmp[2]]
   end

end

for entry in repeatWinds
    keyTmp = entry[1]
   # if(entry[2] >= ceil((0.33)*windowSize))
    weight = entry[2]*1.5
    edgeName = string(keyTmp," [penwidth=$(weight)];")
    network = string(network,edgeName)
    #end
end
    network = string(network,"}")
    fileName = string("networkTotal",":$(stYr)-$(endYr)","windowSize$(windowSize)",".dot")
    print(fileName)
    filePNG = string("networkTotal",":$(stYr)-$(endYr)","windowSize$(windowSize)",".png")
    print(filePNG)
    writedlm(string("./",fileName), [network])

    run(`dot $fileName -Tpng -o $filePNG`)
    print(network)

#######collusion now
yr = stYr
    repeatWinds = Dict()
    while( (yr+windowSize) <= endYr )

        twoWays = collusionDict["2way:$(yr)-$(yr + windowSize)"]
        for pair in twoWays
            pair = replace(pair,"&","")
            pair = replace(pair,".","")
            pairTmp = split(pair,'-')
            pairTmp[1] =  replace(pairTmp[1]," ","")
            pairTmp[2] =  replace(pairTmp[2]," ","")
            edgeName = string(pairTmp[1],"->",pairTmp[2])
            print(edgeName)
            if(haskey(repeatWinds, string(pairTmp[1],"->",pairTmp[2]) ))
                repeatWinds[string(pairTmp[1],"->",pairTmp[2])] = 1 + repeatWinds[string(pairTmp[1],"->",pairTmp[2])]
            else
                repeatWinds[string(pairTmp[1],"->",pairTmp[2])] = 1
            end
        end
        yr = yr + windowSize
    end
    print("\n the aggregrate of repeated collusion \n")
    print(repeatWinds)
#####
network = "digraph{"
yrs = string("years","$(stYr)","to","$(endYr)")

network = string(network,string(" graph [label=","$(yrs)",", fontsize=34]; "))
seenC = []
repeatWinds2 = repeatWinds
for entry in repeatWinds2
    keyTmp = entry[1]
    weight = entry[2]
    pairTmp = split(keyTmp,"->")

    if(!(pairTmp[1] in seenC))
        if(pairTmp[1] in east)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=olivedrab3]; "))
        elseif(pairTmp[1] in southWest)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=indianred1]; "))

        elseif(pairTmp[1] in north)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=dodgerblue]; "))
        elseif(pairTmp[1] in northWest)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=darkslategray2]; "))
        elseif(pairTmp[1] in central)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=gray]; "))
        elseif(pairTmp[1] in southEast)
            network = string(network,string(pairTmp[1]," [style=filled,fillcolor=darkgoldenrod1]; "))
        end
        seenC = [seenC,pairTmp[1]]
    end

    if(!(pairTmp[2] in seenC))
       if(pairTmp[2] in east)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=olivedrab3]; "))

       elseif(pairTmp[2] in southWest)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=indianred1]; "))
       elseif(pairTmp[2] in north)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=dodgerblue]; "))
       elseif(pairTmp[2] in northWest)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=darkslategray2]; "))

       elseif(pairTmp[2] in central)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=gray]; "))
       elseif(pairTmp[2] in southEast)
           network = string(network,string(pairTmp[2]," [style=filled,fillcolor=darkgoldenrod1]; "))
       end
       seenC = [seenC,pairTmp[2]]
   end

end

for entry in repeatWinds
    keyTmp = entry[1]
    weight = entry[2]*3
    edgeName = string(keyTmp," [dir=both color=red penwidth=3 penwidth=$(weight)];")
    network = string(network,edgeName)
end
    network = string(network,"}")
    fileName = string("networkTotalCollusion",":$(stYr)-$(endYr)","windowSize$(windowSize)",".dot")
    print(fileName)
    filePNG = string("networkTotalCollusion",":$(stYr)-$(endYr)","windowSize$(windowSize)",".png")
    print(filePNG)
    writedlm(string("./",fileName), [network])

    run(`dot $fileName -Tpng -o $filePNG`)
    print(network)

#print("\n two way $(yr)-$(yr + windowSize) ratio=")
print("\n the ratio of the two Way\n")
println(ratioTwoWay)
fileNameTmp = string("netRatios",":$(stYr)-$(endYr)","windowSize$(windowSize)",".txt")
writedlm(string("./",fileNameTmp), ratioTwoWay)
drawTrend(string("netRatiosTwoWayToOneWay",":$(stYr)-$(endYr)","windowSize$(windowSize)"),stYr,endYr,windowSize, ratioTwoWay)

print("\n the total of the two Way\n")
println(totalTwoWay)
fileNameTmp = string("netTotalTwoWays",":$(stYr)-$(endYr)","windowSize$(windowSize)",".txt")
writedlm(string("./",fileNameTmp), totalTwoWay)
drawTrend(string("netTotalTwoWays",":$(stYr)-$(endYr)","windowSize$(windowSize)"),stYr,endYr,windowSize, totalTwoWay)

print("\n the total of the one Way\n")
println(totalOneWay)
fileNameTmp = string("netTotalOneWays",":$(stYr)-$(endYr)","windowSize$(windowSize)",".txt")
writedlm(string("./",fileNameTmp), totalOneWay)
drawTrend(string("netTotalOneWays",":$(stYr)-$(endYr)","windowSize$(windowSize)"),stYr,endYr,windowSize, totalOneWay)

cntryNumYrs = countryYrNumber()

twoWayCntryRatio = []
oneWayCntryRatio = []
yr=stYr
ii = 1
while( (yr) < endYr )
    mNum = mean([cntryNumYrs[tmp1] for tmp1 in (yr:(yr+windowSize))])
    #println(totalTwoWay[ii] / mNum)
    println(mNum)
    append!(twoWayCntryRatio,totalTwoWay[ii] / mNum)
    append!(oneWayCntryRatio,totalOneWay[ii] / mNum)
    yr = yr + windowSize
    ii += 1
end
print("\n the total of the two Way country Ratio\n")
println(twoWayCntryRatio)
fileNameTmp = string("netTotalTwoWaysCntryRatio",":$(stYr)-$(endYr)","windowSize$(windowSize)",".txt")
writedlm(string("./",fileNameTmp), twoWayCntryRatio)
drawTrend(string("netTotalTwoWaysCntryRatio",":$(stYr)-$(endYr)","windowSize$(windowSize)"),stYr,endYr,windowSize, twoWayCntryRatio)

print("\n the total of the one Way country Ratio\n")
println(oneWayCntryRatio)
fileNameTmp = string("netTotalOneWaysCntryRatio",":$(stYr)-$(endYr)","windowSize$(windowSize)",".txt")
writedlm(string("./",fileNameTmp), oneWayCntryRatio)
drawTrend(string("netTotalOneWaysCntryRatio",":$(stYr)-$(endYr)","windowSize$(windowSize)"),stYr,endYr,windowSize, oneWayCntryRatio)

#convert to get the CSVs from the graphs
dot2csvConvert()

end

#run(`diplsay network.png`)

function countryYrNumber()
    countryYearsNum = Dict{Integer,Integer}()
    resultsFile = readdir("./dataTables/")
    yrMin = 100000
    yrMax = -1
    for rf in resultsFile
        fileTmp = open(string("./dataTables/",rf))
        linesTmp = readlines(fileTmp)         #readfile lines
        yrTmp = parse(Int,((split(rf,"."))[1]))
        countryNumTmp = length(split(linesTmp[1],",")) - 1#COLUMN NAME LIST (still the same regarding large set mapping to small set in cols) OF WHO CAN RECEIVE VOTES
        countryYearsNum[yrTmp] = countryNumTmp #NUM that receives votes
        close(fileTmp)
        if(yrTmp < yrMin)
            yrMin = yrTmp
        end
        if(yrTmp > yrMax)
            yrMax = yrTmp
        end
    end
    countryYearsNum
end

function drawTrend(titleStr,yrStr,yrEnd,windowSize, ypnts)
    x = [string(yr,"-",yr+windowSize)  for yr in (yrStr:windowSize:(yrEnd-windowSize))];
    y = ypnts
    p = bar(x,y,title=titleStr,leg=false,xrotation=60)
    #p = plot(x,y,title=titleStr,lw=3)
    xlabel!("years")
    savefig(p,"$(titleStr).png")
    p2 = plot(x,y,title=titleStr,lw=4,leg=false,xrotation=60)
    #p = plot(x,y,title=titleStr,lw=3)
    xlabel!("years")
    savefig(p2,"$(titleStr)line.png")
end

max_stars_repo_name: mantzaris/eurovision

max_issues_repo_name: mantzaris/eurovision

max_forks_repo_name: mantzaris/eurovision


text:
# =======================================================================
# 2UMFLP

using vOptSpecific

# ==========================================
# instance

id = load2UMFLP("../examples/F50-51.txt")

# ==========================================
# vOptSpecific : UMFLP_Delmee2017

Z1, Z2, Y, X, isEdge = vSolve(id)

for i = 1:length(z1)
    println("(" , Z1[i], " | ", Z2[i], ") ")
end

ext: jl

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
## Modules call
using Revise, FLPDevelopment, BrowseTables
## Plot settings
# Font size and type
xyfont = font(18, "Bookman Light")
legfont = font(14, "Bookman Light")

# Figure/Plot image size
pixel_factor = 600
ps_x = pixel_factor * 1
ps_y = pixel_factor * 1
fig_size = (ps_x, ps_y)

theme(:default)
pgfplotsx(size = fig_size,
    tick_orientation = :out,
    grid = false,
    markerstrokecolor = :black,
    markersize = 8,
    thickness_scaling = 1.5,
    titlefont = xyfont,
    guidefont = xyfont,
    legendfont = legfont)
## Load and adjust data
include("Young_to_run2.jl")
include("Caspase_to_run.jl")
for df in (Age_p, Age_b, Age_s, Cas_p, Cas_b, Cas_s)
    filter!(r -> r.Protocol == "90/90" &&
    r.MouseID != "CD09" && # biting, see B1_CD09_2020-07-13 minute30
    r.MouseID != "RJ58" && # blind
    r.MouseID != "RJ67" && # biting, see B3_RJ67_2020-09-28 minute 7:33
    !(r.MouseID in first_females_group) &&
    r.ProtocolSession == 1
    ,df)
end
for df in (Cas_p, Cas_b, Cas_s)
    filter!(r -> r.Gen == "Rbp4-cre", df) # exclude wild type animals
end
#adjustments for pokes DataFrames
for pokedf in [Age_p, Cas_p]
    # transform time in log 10 scale
    pokedf.LogOut = log10.(pokedf.Out)
    pokedf.LogInterpoke = log10.(pokedf.PostInterpoke)
    # zscore values for logistic regression
    transform!(pokedf, [:Streak, :Out, :LogOut] .=> zscore)
end
#adjustments for trials DataFrames
for streakdf in [Age_s, Cas_s]
    # transform time in log 10 scale
    streakdf.LogDuration = log10.(streakdf.Trial_Duration)
    streakdf.ElapsedTime = log10.(streakdf.AfterLast_Duration .+ 0.1)
    streakdf.LogTravel = log10.(streakdf.Travel_to)
    # bin trials
    # streakdf.BinnedStreak = bin_axis(streakdf.Streak; unit_step = 4)
    # streakdf.BlockBinnedStreak = bin_axis(streakdf.Streak; unit_step = 15)
    streakdf.LongBinnedStreak = bin_axis(streakdf.Streak; unit_step = 20)
    transform!(streakdf, :Streak .=> zscore)
end
open_html_table(FLPDevelopment.summarydf(Age_s,Age_p))
open_html_table(FLPDevelopment.summarydf(Cas_s,Cas_p))
########################### Example Session Plots ######################################
## Adult example session
lowertrial = 20
uppertrial = 45
tt = filter(r -> r.MouseID == "RJ23" && lowertrial <= r.Streak <= uppertrial, Age_p)
session_plot(tt)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","AdultSession.pdf"))
## Juvenile example session
tt = filter(r -> r.MouseID == "RJ02" && lowertrial <= r.Streak <= uppertrial, Age_p)
session_plot(tt)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","JuvenileSession.pdf"))
## Well trained example session
comp_df = CSV.read("/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Stimulations/DRN_Opto_again/pokesDRN_Opto_again.csv", DataFrame)
comp90 = filter(r-> r.Protocol == "90/90", comp_df)
sessions = union(comp90.Session)
tt = filter(r-> r.Session == sessions[99] && lowertrial <= r.Streak <= uppertrial, comp90)
transform!(groupby(tt, [:Session, :Streak]), [:PokeIn, :PokeOut] => ((i,o) -> (In = i .- i[1], Out = o .- i[1])) => AsTable)
session_plot(tt)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","TrainedSession.pdf"))
########################### Poking rate ######################################
# Age experiment
# calculate the total amount of time spent in poking and the trial duration
Age_df = combine(groupby(Age_p, [:MouseID,:Streak, :Age]), :PokeDur => sum => :PokingTime,
    [:PokeIn, :PokeOut] => ((i,o) -> o[end] - i[1]) => :TrialDuration)
# calculate the fraction of time spent poking per trial: Poking rate
Age_df.PokingRate = Age_df.PokingTime./Age_df.TrialDuration
# Test the difference in the poking rate for age group
Age_Poking = Difference(Age_df, :Age, :PokingRate; ind_summary = median, ylabel = "Poking rate", ylims = (0,1))
Age_Poking.plt
Age_Poking.groupdf
Age_Poking.test
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","DutyCycle_test.pdf"))
## Same for ablation experiment
Cas_df = combine(groupby(Cas_p, [:MouseID,:Streak, :Virus]), :PokeDur => sum => :PokingTime, [:PokeIn, :PokeOut] => ((i,o) -> o[end] - i[1]) => :TrialDuration)
Cas_df.PokingRate = Cas_df.PokingTime ./ Cas_df.TrialDuration
Cas_Poking = Difference(Cas_df, :Virus, :PokingRate; ind_summary = median, ylabel = "Poking rate", ylims = (0,1))
Cas_Poking.plt
Cas_Poking.groupdf
Cas_Poking.test
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","DutyCycle_test.pdf"))
###########################  Leving rate ######################################
Age_LevingAn, Age_LevingAn_df = function_analysis(Age_s,:LogDuration, cumulative_algorythm; grouping = :Age, calc = :bootstrapping)
xprop = ("Poke Time (log10 s)", xyfont,(log10.([0.1,1,10,100,1000]),["0.1","1","10","100","1000"]))
yprop = ("Probablity of leaving", xyfont)
plot!(Age_LevingAn, xaxis = xprop, yaxis = yprop, legend = false)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","LevingRate.pdf"))
##
Cas_SurvivalAn, Cas_SurvivalAn_df = function_analysis(Cas_s,:LogDuration, cumulative_algorythm; grouping = :Virus, calc = :bootstrapping)
plot!(Cas_SurvivalAn, xaxis = xprop, yaxis = yprop, legend = false)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","LevingRate.pdf"))
###########################  Model bootstrap ######################################
Age_Basic_verb = @formula(Leave ~ 1 + Streak_zscore + LogOut_zscore +  (1+Streak_zscore+LogOut_zscore|MouseID));
Age_Basic = fit(MixedModel,Age_Basic_verb, Age_p, Bernoulli())
Age_Full_verb = @formula(Leave ~ 1 + Streak_zscore * Age + LogOut_zscore * Age +  (1+Streak_zscore+LogOut_zscore|MouseID));
Age_Full = fit(MixedModel,Age_Full_verb, Age_p, Bernoulli())
MixedModels.likelihoodratiotest(Age_Basic,Age_Full)
Age_BootDf = CSV.read(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","1000AgeBootstrap.csv"),DataFrame)
# Age_BootDf = bootstrapdf(Age_p, Age_Full; n = 1000)
# CSV.write(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","1000AgeBootstrap.csv"),Age_BootDf)
Age_btdf = Age_BootDf[[1,4,2,3,5,6],:]
Age_btdf.err = [tuple(parse.(Float64,split(err[2:end-1],", "))...) for err in Age_btdf.err]
Age_btdf.variable
yprop = ("",font(10, "Bookman Light"),(collect(1:nrow(Age_btdf)),
    [L"Intercept",L"PokeTime",L"Trial",L"Juveniles",
    L"Trial \&",
    L"PokeTime \&"]))
    xprop = ("Coefficient estimate", xyfont, (-1.33,1.3))
    Titolo = L"+ \Leftarrow Behavioural Control \Rightarrow -"
    @df Age_btdf scatter(:coef ,1:nrow(Age_btdf), xerror = :err,legend = false,
    xaxis = xprop, yaxis = yprop, markercolor = :gray75, title = Titolo)
    vline!([0], linecolor = :red, legend = false, linestyle = :dash)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","1000AgeModelPlot.pdf"))
# leave_modelplt(Age_BootDf)
##
Cas_Basic_verb = @formula(Leave ~ 1 + Streak_zscore + LogOut_zscore +  (1+Streak_zscore+LogOut_zscore|MouseID));
Cas_Basic = fit(MixedModel,Cas_Basic_verb, Cas_p, Bernoulli())
Cas_Full_verb = @formula(Leave ~ 1 + Streak_zscore * Virus + LogOut_zscore * Virus +  (1+Streak_zscore+LogOut_zscore|MouseID));
Cas_Full = fit(MixedModel,Cas_Full_verb, Cas_p, Bernoulli())
MixedModels.likelihoodratiotest(Cas_Basic,Cas_Full)
# Cas_BootDf = bootstrapdf(Cas_p, Cas_Full, n = 1000)
# CSV.write(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","1000CasBootstrap.csv"),Cas_BootDf)
Cas_BootDf = CSV.read(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","1000CasBootstrap.csv"), DataFrame)
Cas_btdf = Cas_BootDf[[1,4,2,3,5,6],:]
Cas_btdf.err = [tuple(parse.(Float64,split(err[2:end-1],", "))...) for err in Cas_btdf.err]
Cas_btdf.variable
yprop = ("",font(10, "Bookman Light"),(collect(1:nrow(Cas_btdf)),
    [L"Intercept",L"PokeTime",L"Trial",L"Caspase",
    L"Trial \&",
    L"PokeTime\&"]))
    xprop = ("Coefficient estimate", xyfont, (-1.4,1.35))
    @df Cas_btdf scatter(:coef ,1:nrow(Cas_btdf), xerror = :err,legend = false,
    xaxis = xprop, yaxis = yprop, markercolor = :gray75, title = Titolo)
    vline!([0], linecolor = :red, legend = false, linestyle = :dash)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","1000ModelPlot.pdf"))
# leave_modelplt(Cas_BootDf)

########################## Incorrect scatter ######################################
Age_Inc = Difference(Age_s, :Age, :IncorrectLeave; ind_summary = incorrect_fraction, ylabel = "Fraction of premature leaving", ylims = (0,1))
Age_Inc.plt
Age_Inc.groupdf
Age_Inc.test
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","Incorrect_test.pdf"))
##
Cas_Inc = Difference(Cas_s, :Virus, :IncorrectLeave; ind_summary = incorrect_fraction, ylabel = "Fraction of premature leaving", ylims = (0,1))
Cas_Inc.plt
Cas_Inc.groupdf
Cas_Inc.test
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","Incorrect_test.pdf"))
########################### Age Sex Model bootstrap ######################################
###########################  Leving rate ######################################
Sex_LevingAn, Sex_LevingAn_df = function_analysis(Age_s,:LogDuration, cumulative_algorythm;
    grouping = :Sex, calc = :bootstrapping, color = [:violetred3 :sandybrown])
    xprop = ("Poke Time (log10 s)", xyfont,(log10.([0.1,1,10,100,1000]),["0.1","1","10","100","1000"]))
    yprop = ("Probablity of leaving", xyfont)
    plot!(Sex_LevingAn, xaxis = xprop, yaxis = yprop, legend = (0.8, 0.25))
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","SFig1","SexLevingRate.pdf"))

########################### Age Sex Model bootstrap #####################################
Age_Sex_verb1 = @formula(Leave ~ 1 + Streak_zscore + LogOut_zscore + Age +  (1+Streak_zscore+LogOut_zscore|MouseID));
Age_Sex_verb2 = @formula(Leave ~ 1 + Streak_zscore + LogOut_zscore + Age*Sex +  (1+Streak_zscore+LogOut_zscore|MouseID));
Age_Sex1 = fit(MixedModel,Age_Sex_verb1, Age_p, Bernoulli())
Age_Sex2 = fit(MixedModel,Age_Sex_verb2, Age_p, Bernoulli())
MixedModels.likelihoodratiotest(Age_Sex1,Age_Sex2)
# Sex_BootDf2 = bootstrapdf(Age_p, Age_Sex2; n = 1000)
# CSV.write(joinpath(replace(path,basename(path)=>""),"Development_Figures","SFig1","1000SexBootstrap.csv"),Sex_BootDf2)
Sex_BootDf2 = CSV.read(joinpath(replace(path,basename(path)=>""),"Development_Figures","SFig1","1000SexBootstrap.csv"),DataFrame)
Sex_btdf = Sex_BootDf2[[1,3,2,4,5,6],:]
Sex_btdf.err = [tuple(parse.(Float64,split(err[2:end-1],", "))...) for err in Sex_btdf.err]
yprop = ("",font(10, "Bookman Light"),(collect(1:nrow(Sex_btdf)),
    [L"Intercept",L"PokeTime",L"Trial",L"Juveniles",
    L"Male",
    L"Juveniles \&"]))
    xprop = ("Coefficient estimate", xyfont, (-1.33,1.3))
    @df Sex_btdf scatter(:coef ,1:nrow(Sex_btdf), xerror = :err,legend = false,
    xaxis = xprop, yaxis = yprop, markercolor = :gray75, title = Titolo)
    vline!([0], linecolor = :red, legend = false, linestyle = :dash)
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","SFig1","1000FullSexModelPlot.pdf"))
########################### Over trials median Survival rate ######################################
Age_TrialMedianSurvival=plot();
att = combine(groupby(Age_s,:Age)) do dd
    mediansurvival_analysis(dd,:LogDuration,:LongBinnedStreak; plt = Age_TrialMedianSurvival)
end
Age_TrialMedianSurvival
xprop = ("Trial", (10,90), xyfont)
yprop = ("Median laeving time (log10 s)", xyfont,(log(1),log10(30)),(log10.([1,5,10,15,20, 30]),string.([1,5,10,15,20,30])))
plot!(xaxis = xprop, yaxis = yprop, grid = true, size = (600,600))
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig1","PartialMedianRatePerTrial.pdf"))
##
Cas_TrialMedianSurvival = plt = plot()
combine(groupby(Cas_s,:Virus)) do dd
    mediansurvival_analysis(dd,:LogDuration,:LongBinnedStreak; plt = Cas_TrialMedianSurvival)
end
Cas_TrialMedianSurvival
plot!(xaxis = xprop, yaxis = yprop, grid = true, size = (600,600))
savefig(joinpath(replace(path,basename(path)=>""),"Development_Figures","Fig3","PartialMedianRatePerTrial.pdf"))


text:
"""
    SegmentedMax{V <: AbstractVector{<:Number}} <: AbstractAggregation

[`AbstractAggregation`](@ref) implementing segmented max aggregation:

``
f(\\{x_1, \\ldots, x_k\\}) = \\max_{i = 1, \\ldots, k} x_i
``

Stores a vector of parameters `ψ` that are filled into the resulting matrix in case an empty bag is encountered.

See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`max_aggregation`](@ref),
    [`SegmentedMean`](@ref), [`SegmentedSum`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
"""
struct SegmentedMax{V <: AbstractVector{<:Number}} <: AbstractAggregation
    ψ::V
end

Flux.@functor SegmentedMax

SegmentedMax(T::Type, d::Int) = SegmentedMax(zeros(T, d))
SegmentedMax(d::Int) = SegmentedMax(Float32, d)

Flux.@forward SegmentedMax.ψ Base.getindex, Base.length, Base.size, Base.firstindex, Base.lastindex,
        Base.first, Base.last, Base.iterate, Base.eltype

Base.vcat(as::SegmentedMax...) = reduce(vcat, as |> collect)
function Base.reduce(::typeof(vcat), as::Vector{<:SegmentedMax})
    SegmentedMax(reduce(vcat, [a.ψ for a in as]))
end

function (a::SegmentedMax)(x::Maybe{AbstractMatrix{T}}, bags::AbstractBags,
                              w::Optional{AbstractVecOrMat{T}}=nothing) where T
    _check_agg(a, x)
    segmented_max_forw(x, a.ψ, bags)
end

segmented_max_forw(::Missing, ψ::AbstractVector, bags::AbstractBags) = repeat(ψ, 1, length(bags))
function segmented_max_forw(x::AbstractMatrix, ψ::AbstractVector, bags::AbstractBags)
    T = promote_type(eltype(x), eltype(ψ))
    y = Matrix{T}(fill(_typemin(T), size(x, 1), length(bags)))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                y[i, bi] = ψ[i]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    y[i, bi] = max(y[i, bi], x[i, j])
                end
            end
        end
    end
    y
end

function segmented_max_back(Δ, y, x, ψ, bags)
    dx = zero(x)
    dψ = zero(ψ)
    v = similar(x, size(x, 1))
    idxs = zeros(Int, size(x, 1))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                dψ[i] += Δ[i, bi]
            end
        else
            fi = first(b)
            v .= x[:,fi]
            idxs .= fi
            for j in b
                for i in 1:size(x,1)
                    if v[i] < x[i, j]
                        idxs[i] = j
                        v[i] = x[i, j]
                    end
                end
            end
            for i in 1:size(x, 1)
                dx[i, idxs[i]] += Δ[i, bi]
            end
        end
    end
    dx, dψ, NoTangent()
end

function segmented_max_back(Δ, y, x::Missing, ψ, bags)
    dψ = zero(ψ)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(ψ)
            dψ[i] += Δ[i, bi]
        end
    end
    ZeroTangent(), dψ, NoTangent()
end

function ChainRulesCore.rrule(::typeof(segmented_max_forw), args...)
    y = segmented_max_forw(args...)
    grad = Δ -> (NoTangent(), segmented_max_back(Δ, y, args...)...)
    y, grad
end


text:
module Runoffs
"""
Module for the modelling of waterborne pollutants.
"""

using ArchGDAL
using Dates

include(".\\Utils\\Functions.jl")

const agd = ArchGDAL

function getFeatureByFid( features::Vector{agd.Feature}, fid )
    for f in features
        if agd.getfid(f) == fid
            return f
        end
    end
    return nothing
end

"""
    hash_adjacent( t::Tuple{Int64, Int64} )

Given a tuple with both values in `-1:1` (except `(0, 0)`) return a value in `1:8`, matching the cell displacement with its numeric direction reppresentation.
"""
function hash_adjacent( t::Tuple{Int64, Int64} )
    # 3x3 matrix linearization
     # Its necessary to sum 2 to both values of t as they reppresent index displacements from the center of the matrix
      # the linearization process is based on starting index 0, hence the need to sum 2 instead of 3
      # for the same reason 1 needs to be substracted from the first index, ending up with "t[1] + 1"
    res = 3(t[1] + 1) + (t[2] + 2)
    # The cell of index 5 is the center and it's not needed
    return res >= 5 ? res - 1 : res
end

"""
    hash_adjacent( r::Int64, c::Int64 )

Given a two values in `-1:1` (except `(0, 0)`) return a value in `1:8`, matching the cell displacement with its numeric direction reppresentation.
"""
function hash_adjacent( r::Int64, c::Int64 )
    # 3x3 matrix linearization
     # Its necessary to sum 2 to both values of t as they reppresent index displacements from the center of the matrix
      # the linearization process is based on starting index 0, hence the need to sum 2 instead of 3
      # for the same reason 1 needs to be substracted from the first index, ending up with "t[1] + 1"
    res = 3(r + 1) + (c + 2)
    # The cell of index 5 is the center and it's not needed
    return res >= 5 ? res - 1 : res
end

"""
    hash_adjacent( i::Int64 )

Given a value in `1:8` return the corresponding tuple with both values in `-1:1`, except `(0, 0)`, matching the direction reppresentaion to the cell displacement.
"""
function hash_adjacent( i::Int64 )
    # To account for the fact that "i" = 5 for the tuple "(0, 0)"
    if i > 4
        i += 1
    end
    r = ceil(Int64, i / 3)
    c = ((i - 1) % 3) + 1
    # Substract 2 to obtain the index displacement from the center
    return r - 2, c - 2
end

"""
    connectivity_batch!( mat, dem_band::AbstractArray{T}, noDataValue::Real ) where {T <: Number}

Given NxMx8 (N and M must be the same between `mat` and `dem_band`) matrix `mat`, for each group of 8 cells with indexes r and c,
store in them the differences in values between cell (r, c) of dem_band and its 8 adjacent cells, unless the value is equal to `noDataValue`.
"""
function connectivity_batch!( mat, dem_band::AbstractArray{T}, noDataValue::Real ) where {T <: Number}
    rows, cols = size(dem_band)
    # For each cell of the dem's band
    @inbounds for r in 1:rows, c in 1:cols
        if dem_band[r, c] != noDataValue
            # Indexes of adjacent cells
            for (i, j) in keys(indexes)
                if ( r+i >= 1 && r+i <= rows ) && ( c+j >= 1 && c+j <= cols ) && dem_band[r+i, c+j] != noDataValue
                    if mat[ r, c, hash_adjacent(i, j) ] == noDataValue
                        mat[ r, c, hash_adjacent(i, j) ] = dem_band[r, c] - dem_band[r+i, c+j]
                    end
                    if mat[ r+i, c+j, hash_adjacent(-i, -j) ] == noDataValue
                        mat[ r+i, c+j, hash_adjacent(-i, -j) ] = dem_band[r+i, c+j] - dem_band[r, c]
                    end
                end
            end
        end
    end
end

"""
    connectivity( dem_band::Matrix{T}, batch_size::Integer, noDataValue::Real ) where {T <: Number}

Generate a matrix with the same size as `dem_band` and 8 layers, containing for each cell (r, c) of dem_band the difference between the value of that cell and
its 8 adjacent cells, unless the have value equal to `noDataValue`.
The operation will be executed one `batch_size`x`batch_size` section of the matrix at a time.
"""
function connectivity( dem_band::Matrix{T}, batch_size::Integer, noDataValue::Real ) where {T <: Number}
    rows, cols = size(dem_band)
    mat = fill( convert(Float32, noDataValue) , rows, cols, 8 )
    if batch_size == 1
        connectivity_batch!(mat, dem_band, noDataValue)
        return mat
    end
    n, m = ceil.( Int64, [rows, cols] ./ (batch_size - 1) )
    @inbounds for i in 1:n, j in 1:m
        # Find the starting and ending indexes for the current slice of the matrix
         # if the batch is one of the ending ones its ending index will be the size of the matrix for that dimension
        rows_range::UnitRange{Int64} = ( (batch_size - 1) * (i - 1) + 1 ) : ( i != n ? (batch_size - 1) * i + 1 : rows )
        cols_range::UnitRange{Int64} = ( (batch_size - 1) * (j - 1) + 1 ) : ( j != m ? (batch_size - 1) * j + 1 : cols )
        # Use the indexes to run the function on a view of the matrix (passing also the corresponding view of the dem)
        connectivity_batch!( view(mat, rows_range, cols_range, :), view(dem_band, rows_range, cols_range), noDataValue )
    end
    return mat
end

"""
    createRasterizedConnectivity( file::AbstractString, dtm_path::AbstractString )

Create and save as `file` a raster with the same size as the one pointed at by `dtm_path` and 8 layers, containing for each cell (r, c) of  the latter the difference
between the value of that cell and its 8 adjacent cells, unless the have value equal to `noDataValue`.
The operation will be executed one `batch_size`x`batch_size` section of the matrix at a time.
"""
function createRasterizedConnectivity( file::AbstractString, dtm_path::AbstractString )
    mat = connectivity(band_mat, 2048, ndv)
    dtm = agd.read(dtm_path)
    rows, cols = size(mat)
    res = agd.create( file, driver=agd.getdriver("GTiff"), width=rows, height=cols, nbands=8, dtype=Float32 )
    agd.setgeotransform!( res, agd.getgeotransform(dtm) )
    agd.setproj!( res, agd.getproj(dtm) )
    valNoData = agd.getnodatavalue(dtm)
    bands = [ agd.getband( res, i ) for i in 1:8 ]
    agd.setnodatavalue!.(bands, valNoData)
    agd.fillraster!.(bands, valNoData)
    band_mats = agd.read.(bands)

    for r in 1:rows, c in 1:cols
        if !ismissing(mat[r, c]) && !isempty(mat[r, c])
            for (i, j, val) in mat[r, c]
                band_mats[hash_adjacent(i, j)][r, c] = val
            end
        end
    end

    for i in 1:8
        agd.write!( target_ds, band_mats[i], i )
    end

    agd.destroy(target_ds)
end

using Plots

import ArchGDAL as agd

include("../Library/Functions.jl")

# Find the direction with the greatest difference in height from (r, c)
function mindir( mat, r::Int64, c::Int64 )
    dir = 0
    min = Inf
    @inbounds for i in 1:8
        if mat[r, c, i] != -9999.0 && mat[r, c, i] < min
            dir = i
            min = mat[i][r, c]
        end
    end
    return (dir, min)
end
 #= MINDIR TEST VALUES
    arr = [
        [ 1  2  -9999.0; 1  2  3 ],
        [ 3  4  -9999.0; 4  5  6 ],
        [ 5  6  -9999.0; 7  8  9 ],
        [ 7  8  -9999.0; 10 11 12 ],
        [ 9  10 -9999.0; 13 14 15 ],
        [ 11 12 -9999.0; 16 17 18 ],
        [ 13 14 -9999.0; 19 20 21 ],
        [ 17 16 -9999.0; 22 23 24 ]
    ]
 =#

# Return the amount of water flowing along the path formed by the points in "flowpoints" at "instant" and update the value corresponding to the volume of water at each point
function moving_water!( flowpoints::Vector, instant::Int64, rain::Matrix, permeability::Matrix )
    v_moving = 0 # Water flowing from already visited cells following the path
    for (rₚ, cₚ, Δhminₚ, vₚ) in flowpoints
        vₚ = (vₚ + rain[instant][rₚ, cₚ] + v_moving) * (1 - permeability[rₚ, cₚ])
        if Δhmin <= 0 # If the following cell in the path is below the current one
            # All the water on the cell flows to the next one
            v_moving = vₚ
            vₚ = 0
        else # If the cell is below all the surrounding cells
            # The water accumulates on the cell, if the volume of the water added to the height of the cell exceeds the height of the lowest adjacent cell
             # the water overflows
            if vₚ > Δhminₚ
                v_moving = vₚ - Δhminₚ
                vₚ -= v_moving
            else
                v_moving = 0
            end
        end
    end
    return v_moving
end

function flow( flowpoints::Vector, r::Int64, c::Int64, rows::Int64, cols::Int64, instants::Int64, band::Matrix{Float32}, rev_rain_band::Matrix, permeability_band::Matrix )
    while 1 < r < rows && 1 < c < cols
        # While its still raining
        v_moving = 0 # Volume of water flowing FROM the current cell to the next in the current instant
         # "moving_water!" computes the water flowing TO the current cell and so is an instant behind accounting for "v_moving"
        while instants > 0
            # Ad ogni istante ricalcoliamo il volume d'acqua che sta scorrendo lungo il percorso del flusso, in questo modo il volume d'acqua
             # di ogni nuova cella tiene conto del fatto che la pioggia influisce su tutte le celle appertenenti al percorso
            # Volume of water on the current cell
            v = ( sum(rain_band[instants:end][r, c]) + v_moving + moving_water!(flowpoints, instant, rev_rain_band, permeability_band) ) * (1 - permeability_band[r, c])
            dir, Δhmin = mindir(band, r, c)
            if Δhmin <= 0 # If there is a lower adjacent cell all the water flows there
                v_moving = v
                v = 0
                instants -= 1
            else # If the cell is a local minimum the water keeps accumulating untill it overflows
                while instants > 0 && Δhmin > v
                    v = ( v + rain_band[instants][r, c] + moving_water!(flowpoints, instant, rev_rain_band, permeability_band) ) * (1 - permeability[r, c])
                    instants -= 1
                end
                # If the rain stops before the water overflows the stream alts on the cell
                v_moving = instants == 0 ? 0 : v - Δhmin
                v -= v_moving
            end

            push!(flowpoints, (r, c, Δhmin, v))

            # If there is no water flowing the flow stops
            if v_moving == 0 && instants == 0
                break
            end

            # Flow to the next cell
            Δr, Δc = hash_adjacent(dir)
            r += Δr
            c += Δc
        end
        # After the rain stops but there is still water flowing
        while v_moving > 0
        end
    end
end

function mindir( mat, r::Int64, c::Int64, noDataValue::Real )
    dir = [0]
    min = Inf
    @inbounds for i in 1:8
        #   println( "$i) $(mat[i][r, c])" )
        if mat[r, c, i] != noDataValue
            if mat[r, c, i] < min
                dir = [i]
                min = mat[r, c, i]
                continue
            end
            if mat[r, c, i] == min
                push!(dir, i)
            end
        end
    end
    return (dir, min)
end

function flow( band::Matrix{Float32}, rain, permeability, noDataValue::Real )
    rows, cols = size(band)
    water = zeros(Float64, rows, cols, 2) # water[:,:,1]: water volume, water[:,:,2]: volume of incoming water
    flow = Array{ Union{ Vector{Float64}, Float64 } }(noDataValue, rows, cols)

    for rain_epoch in rain
        # Compute the amount of water for every cell
        water[:,:,1] = (water[:,:,1] + rain_epoch + water[:,:,2]) .* (1 .- permeability)

        # For every cell compute the water flowing from that cell to another and update the volume of water in the cell
        @inbounds for r in 1:rows, c in 1:cols
            if band[r,c] == noDataValue
                water[r, c, 1] = water[r, c, 2] = noDataValue
                flow[r, c] = noDataValue
                continue
            end
            # Get the minimum difference in height with the neighbouring cells and the directions to the cells at that height
            dir, Δhmin = mindir(band, r, c)
            # Get the displacements to reach said cells from the current one
            flow[r, c] = hash_adjacent.(dir)

            # If the volume of water is greater than the difference in height between the current cell and the lowest adjacent, the water flows
            if water[r, c, 1] > Δhmin
                ld = length(flow[r,c])
                # The volume of flowing water depends on the difference of height between the cell and the adjacent one
                v = Δhmin <= 0 ? water[r, c, 1] : water[r, c, 1] - Δhmin
                # To the incoming water of cell "(r+Δr, c+Δc)" is added a portion of the water of cell "(r, c)" depending on the number of paths
                for delta in flow[r,c]
                    water[r+delta[1], c+delta[2], 2] += v / ld
                end
                # The remaining water is equal to the original volume less the amount that passed to the adjacent cell
                water[r, c, 1] -= v
            end
        end
    end

    return flow, water
end

#=
Problema:
Trovare il percorso seguito dall'acqua tenendo conto che se il percorso trova un'avvallamento e si ferma li, l'acqua si accumula e con il
    passare del tempo potrebbe riuscire a fuoriuscire e proseguire la corsa.

Idea:
Trovare il percorso fino al primo stop (primo avvallamento) calcolando il cammino minimo discendente (il cammino in cui da ongi punto al successivo si scende),
    trovato il primo stop, verificare se l'acqua ad un certo punto (temporale) può straripare (se l'avvallamento è profondo come il fosso delle marianne non succederà, o,
    qunatomeno, non in tempi abbastanza brevi da avere senso per lo studio), se ciò può avvenire (perchè la profondità dell'avvallamento lo consente), riapplicare il calcolo
    del cammino minimo da questo punto fino al prossimo, ripetere finchè non si trova uno stop forte (= buco molto profondo).
Calcolato questo percorso ideale utlizzando i dati della piovosità ecc calcolare l'effettivo flusso d'acqua (il flusso seguirà per forza questo percorso essendo il
    minimo, quindi la problematica principale sarà individuare il punto di stop del flusso).
=#

function flowpath!( mat::Matrix{Float32}, band::Matrix{Float32}, noDataValue::Real, row::Int64, col::Int64, )
    if band[row, col] != noDataValue
        r, c = row, col
        mat[r, c] = 1.f0
        directions, Δhmin = mindir(band, row, col)
        while Δhmin < 0
            if lenght(directions) > 1
                for dir in directions[2:end]
                    Δr, Δc = hash_adjacent(dir)
                    flowpath!( mat, band, noDataValue, r + Δr, c + Δc )
                end
            end
            r, c .+= hash_adjacent(directions[1])
            directions, Δhmin = mindir(band, r, c)
        end
    end
end

GC.gc()
# Per porre a 1 tutte le celle di un raster che stanno dentro and un poligono delimitato da uno shapefile
using Rasters
using Shapefile

dtm_file = split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\DTM_wgs84.tiff"
# sat_file = split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\sat WGS84\\sette_sorelle.shp"
csoil_file = "D:\\Z_Tirocinio_Dati\\Classi suolo WGS84\\Classi suolo.shp"
perm_file = "D:\\Z_Tirocinio_Dati\\Permeabilità suolo WGS84\\Permeabilità suolo.shp"
ccs_file = "D:\\Z_Tirocinio_Dati\\ccs WGS84\\ccs.shp"
ccs_file = "C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\ccs WGS84\\ccs.shp"

dtm = Raster(dtm_file)
# sat_shp = Shapefile.Handle(sat_file).shapes[1]
csoil_shp = Shapefile.Table(csoil_file)
perm_shp = Shapefile.Table(perm_file)
ccs_shp = Shapefile.Table(ccs_file)

# Creation of the final raster
originX = dtm.dims[1][1]
originY = dtm.dims[2][1]
stepX = dtm.dims[1][2] - originX
stepY = dtm.dims[2][2] - originY
dims =  X( Projected( originX:stepX:dtm.dims[1][end]+stepX; order=Rasters.ForwardOrdered(), span=Rasters.Regular(stepX), sampling=Rasters.Intervals(Rasters.Start()), crs=EPSG(4326) ) ),
        Y( Projected( originY:stepY:dtm.dims[2][end]+stepY; order=Rasters.ReverseOrdered(), span=Rasters.Regular(stepY), sampling=Rasters.Intervals(Rasters.Start()), crs=EPSG(4326) ) )

#   csoil = Raster( fill(dtm.missingval, dims); missingval=dtm.missingval )
#   perm = Raster( fill(dtm.missingval, dims); missingval=dtm.missingval )
ccs = Raster( fill(dtm.missingval, dims); missingval=dtm.missingval )

# CAMBIARE IL VALORE ASSOCIATO A "X" IN missingval O SIMILI, METTERE VALORI SENSATI E SISTEMARE I VALORI DELLE COPPIE
#   cdict = Dict( val => Float32(i) for (i, val) in enumerate(unique(csoil_shp.gr_idrolog)) )
#   for (shape, value) in zip(csoil_shp.geometry, csoil_shp.gr_idrolog)
#       rasterize!( csoil, shape; fill=cdict[value], order=(X, Y) )
#   end
#
#   Rasters.write( "D:\\Z_Tirocinio_Dati\\Classi suolo WGS84\\Classi suolo.tiff", csoil )

#   pdict = Dict( val => isnothing(tryparse(Float32, val)) ? 0.f0 : tryparse(Float32, val) for (i, val) in enumerate(unique(perm_shp.perm_clas)) )
#   for (shape, value) in zip(perm_shp.geometry, perm_shp.perm_clas)
#       rasterize!( perm, shape; fill=pdict[value], order=(X, Y) )
#   end
#
#   Rasters.write( "D:\\Z_Tirocinio_Dati\\Permeabilità suolo WGS84\\Permeabilità suolo.tiff", perm )

ccdict = Dict( desc => Float32(i) for (i, desc) in enumerate(unique(ccs_shp.legenda)) )
for (shape, value) in zip(ccs_shp.geometry[1:2], ccs_shp.legenda[1:2])
    rasterize!( ccs, shape; fill=ccdict[value], order=(X, Y) )
end

rasterize!( ccs, ccs_shp.geometry[1]; fill=1.f0, order=(X, Y) )
rasterize!( ccs, ccs_shp.geometry[394066]; fill=3.f0, order=(X, Y) )
rasterize!( ccs, ccs_shp.geometry[2]; fill=2.f0, order=(X, Y) )

Rasters.write( "D:\\Z_Tirocinio_Dati\\test.tiff", ccs )

# I poligoni 1 e 2 sono molto vicini ma non sono esattamente adiacenti (i loro mbr si intersecano) in mezzo a loro due c'è il poligono 394066.
#   Dopo la rasterizzazione i due poligoni si toccano coprendo il poligono che dovrebbe stare tra i due.
#   Rasterizzado tutti e tre non sembrano esserci sovrapposizioni nelle celle, ma i poligoni sono di piccole dimensioni e non sono rasterizzati con precisione.

using ArchGDAL
const agd = ArchGDAL

# DA "unique(ccs)" VENGONO 58 VALORI DIVERSI, AL POSTO DI 86 (un valore per ogni poligono più uno per il "-9999.0") QUINDI DEI POLIGONI VENGONO COPERTI DA ALTRI
res_ids = tryparse.( Int64, split( read("C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\res_ids.txt", String), "\n", keepempty=false ) )
#   ccs_shp = agd.read("D:\\Z_Tirocinio_Dati\\ccs WGS84\\ccs.shp")
#= IL SORTING VIENE EFFETUATO PRIMA DI SALVARE GLI ID
features =  collect( agd.getlayer( agd.read("C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\ccs WGS84\\ccs.shp"), 0 ) )
sort!(res_ids, lt=(x, y) -> isless( agd.geomarea.( agd.getgeom.(features[[x, y]]) )... ) )
=#
for (i, id) in enumerate(resids)
    rasterize!( ccs, ccs_shp.geometry[id], fill=Float32(i), order=(X,Y) )
end
Rasters.write( "C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\sat_polys.tiff", ccs )
Rasters.write( "D:\\Z_Tirocinio_Dati\\sat_polys.tiff", ccs )

using Rasters
using SpatialGraphs

dtm = Raster(split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\DTM_wgs84.tiff")
csoil = Raster("C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\Classi suolo WGS84\\Classi suolo.tiff")
perm = Raster("C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\Permeabilità suolo WGS84\\Permeabilità suolo.tiff")
perm = Raster("D:\\Z_Tirocinio_Dati\\Permeabilità suolo WGS84\\Permeabilità suolo.tiff")

perm2 = Raster(
    replace(
        perm.data,
        0.f0 => -9999.f0,
        2.f0 => 0.036f0,
        3.f0 => 0.36f0,
        4.f0 => 3.6f0,
        5.f0 => 36.f0,
        6.f0 => 360.f0
    ),
    perm.dims,
    missingval = perm.missingval
)

# "weights" sarà il raster delle resistenze di ongi cella potrebbe essere semplicemente il raster della permeabilità, oppure una combinazione di permeabilità e altezza
wrg = weightedrastergraph(
    perm2,
    directed = true,
    condition_raster = dtm,
    condition = ( hs, hd ) -> hs != dtm.missingval && hd != dtm.missingval && hs >= hd
)

function run_runoff( dem_file::AbstractString, ccs_file::AbstractString, source_file::AbstractString, target_file::AbstractString, resolution::Int64, folder::AbstractString=".\\" )

    source = agd.read(source_file)

    if agd.geomdim(source) != 0
        throw(DomainError(source_file, "The source shapefile must contain a single point."))
    end

    target = agd.read(target_file)
    ccs = agd.read(ccs_file)
    target_layer, landcover_layer = agd.getlayer([target, ccs], 0)

    if agd.geomdim(target_layer) != 2
        throw(DomainError(target_file, " The target shapefile must contain a polygon."))
    end

    if agd.geomdim(landcover_layer) != 2
        throw(DomainError(ccs_file, "Not a valid landcover geometry."))
    end

    refsys = agd.getspatialref(source)

    dem = agd.read(dem_file)

    if agd.getspatialref(target_layer) != agd.getspatialref(source) || agd.getspatialref(landcover) != agd.getspatialref(source) || agd.getspatialref(dem) != agd.getspatialref(source)
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

    path_temp_landcover = folder * "\\temp_lc.tiff"
    path_temp_soil = folder * "\\temp_soil.tiff"

    clc_df = Functions.cn_list_extract()

    soil_control = 0

 """ PRINT DI COSE
    messaggio+='ALGORITMO UTILIZZATO: calcolo della separazione delle componenti infiltrazione e ruscellamento tramite metodo SCS-CN; US Department of Agriculture Soil Conservation Service, 1972. National Engineering Handbook, Section 4, Hydrology. US Government Printing Office, Washington, DC, 544pp.\n\n'
 """

    output_path = folder * "\\runoff.tiff"
    intervallo = max(getCellDims(dem))

    # NON SO SE FUNZIONA
    bbox_src = agd.boundingbox(source)
    bbox_trgt = agd.boundingbox(target)
    # Geometry containing both the source and the target
    # area = agd.union( bbox_src, bbox_trgt )
    # Se union non ritorna una geometria unica ma un'unica geometria composta di due elementi disgiunti
    area = agd.boundingbox(agd.union( bbox_src, bbox_trgt ))

 """ PRENDE LA PORZIONE DI `landcover` RAPPRESENTATA DA `area`
    lc_clip_proc = processing.run('qgis:clip', {'INPUT':self.lc, 'OVERLAY':self.areastudio, 'OUTPUT':self.path_working+'/clip.gpkg'})
    lc_clip=QgsVectorLayer(lc_clip_proc['OUTPUT'], 'lc_clip', 'ogr')
    lc_clip.setCrs(self.source.crs())

    #path__layer_lc=lc_clip['OUTPUT'].dataProvider().dataSourceUri()
    path__layer_lc=lc_clip.dataProvider().dataSourceUri()
    path_lc=path__layer_lc.split("|")
    source_ds_lc = ogr.Open(path_lc[0])
    lc_layer = source_ds_lc.GetLayer()
 """
    # NON SONO CERTO PRESERVI LE INFORMAZIONI DI landcover
    landcover_clip = agd.intersection(landcover, area)
    landcover_layer = agd.getlayer(landcover_clip, 0)

    rows = agd.height(landcover_layer)
    cols = agd.width(landcover_layer)
    valNoData = -9999.0
    gtiff_driver = agd.getdriver("GTiff")

    landcover_ds = agd.create( path_temp_landcover, driver=gtiff_driver, width=rows, height=cols, nbands=1, dtype=Float32 )
    agd.setgeotransform!( landcover_ds, agd.getgeotransform(dem) )
    agd.setproj!(landcover_ds, refsys)
    band_lc = agd.getband(landcover_ds, 1)
    agd.setnodatavalue!(band_lc, valNoData)
    agd.fillraster!(band_lc, valNoData)
    bandlc = agd.read(band_lc)

 """ TRASFORMA IN RASTER LA PORZIONE DI `landcover` PRESA PRIMA """
    gdal.RasterizeLayer(lc_ds, [1], lc_layer,options=["ATTRIBUTE="+self.text_lcfield])
    lc_ds=None

    lc_layer=QgsRasterLayer(self.path_temp_lc,"lc_layer")
 """"""
    landcover_layer = agd.gdalrasterize( x -> x, landcover_ds )

    if soil_text == "Valore campo"
        soil_control = 1

        agd.getlayer(path_lc, 0)

        soil_ds = agd.create( path_temp_soil, driver=gtiff_driver, width=round(Int64, x_res), height=round(Int64, y_res), nbands=1, dtype=Float32 )
        agd.setgeotransform!( soil_ds, agd.getgeotransform(dem) )
        agd.setproj!(soil_ds, refsys)
        band_sl = agd.getband(soil_ds, 1)
        agd.setnodatavalue!( band_sl, Float64(valNoData) )
        bandsl = agd.read(band_sl)
        agd.fillraster!(bandsl, valNoData)

     """ RASTERIZZA IL  VETTORIALE DEL TIPO DI SUOLO """
        gdal.RasterizeLayer(lc_ds, [1], soil_layer,options=["ATTRIBUTE="+self.text_soilfield])
        soil_ds=None

        soil_layer=QgsRasterLayer(self.path_temp_soil,"soil_layer")
     """"""
        soil_layer = agd.gdalrasterize( x -> x, landcover_ds )
    end

 # VANNO PASSATI A flow I RANGE DELL'AREA IN CUI PIOVE

 # ASSUMENDO CHE target VENGA PASSATO COME FILE VETTORIALE
    target_layer = collect(agd.getlayer(target, 0))
    target_geom = agd.getgeom(agd.getgeom(target_layer[1], 0), 0)
 # SE target VIENE PASSATO COME GEOMETRIA ELIMINARE LE DUE RIGHE SOPRA
    length = agd.ngeom(target_geom)
    target_points = [ agd.getpoint(target_geom, i)[1:2] for i in 1:length ]

    x_source = agd.getx(source, 0)
    y_source = agd.gety(source, 0)
    mat = agd.read(".\\connectivity.tiff")
    res = flow(mat, x_source, y_source)

    volumes = []
    for point in res
        lc = landcover_layer[point...]
        if soil_control == 1
            soil = soil_layer[point...]
        else
            soil = text_soil
        end
        # IL BLOCCO TRY CATCH NON LO VEDO PARTICOLARMENTE NECESSARIO
        try
         # DA MODIFICARE TENENDO CONTO DEI VETTORIALI DEL SUOLO
         # PER LA CLASSE DEL SUOLO SI PUO' FAR RIFERIMENTO AL DIZIONARIO DEL FILE
         # CREDO STIA PRENDENDO LA RESISTENZA DEL SUOLO NEL PUNTO CORRENTE
            # Se passiamo come parameto un raster delle resisteze questa riga va sostituita
            cn = clc_df[lc, soil] # Classisoi ottenuto da "substance.db"
            S = 254.0((100 / cn) - 1)
        catch
            S = 0
        end
        pcheck = (p0 - 0.2S)^2 / (p0 - 0.2S + S)
        push!(volumes, pcheck > 0.2S ? pcheck : 0 )
    end

end

end # module

#=
function run_runoff( dem, source, target, landcover, soil_text::AbstractString, resolution::Integer, folder::AbstractString=".\\" )

    """ NON SO QUALE SIA L'EQUIVALENTE
       if not self.dem.isValid():
           QMessageBox.warning(self,"Warning", "The dem file is not valid" )
           return
    """

       if agd.geomdim(source) != 0
           throw(DomainError(source, "`source` must be a point"))
       end

       if agd.geomdim(targer) != 2
           throw(DomainError(source, "`target` must be a polygon"))
       end

       if agd.geomdim(landcover) != 2
           throw(DomainError(source, "Not a valid `landcover` geometry"))
       end

       if agd.getspatialref(target) != agd.getspatialref(source) || agd.getspatialref(landcover) != agd.getspatialref(source) || agd.getspatialref(dem) != agd.getspatialref(source)
           throw(DomainError("The reference systems are not uniform. Aborting analysis."))
       end

       refsys = agd.importEPSG(agd.fromWKT(agd.getspatialref(source)))

       path_temp_landcover = folder * "\\temp_lc.tiff"
       path_temp_soil = folder * "\\temp_soil.tiff"

    """ NON SO COSA SIANO QUESTE VARIABILI
       self.text_p = str(self.combo_fieldp.currentText())
    """

       clc_list = Functions.cn_list_extract()

       soil_control = 0

    """ PRINT DI COSE
       messaggio+='ALGORITMO UTILIZZATO: calcolo della separazione delle componenti infiltrazione e ruscellamento tramite metodo SCS-CN; US Department of Agriculture Soil Conservation Service, 1972. National Engineering Handbook, Section 4, Hydrology. US Government Printing Office, Washington, DC, 544pp.\n\n'
    """

       output_path = folder * "\\runoff.tiff"

       valNoData = -9999

       gtiff_driver = agd.getdriver("GTiff")
       target_ds = agd.create( output_path, gtiff_driver, round(Int64, ), round(Int64, ), 1, agd.GDAL.GDT_Float32 )
       agd.setgeotransform!(target_ds, [ , resolution, 0.0, , 0.0, -resolution ])
       agd.setproj!(target_ds, refsys)
    """ NON SO QUALE SIA IL COMANDO PER SETTARE I METADATI CON `ArchGDAL`
       target_ds.SetMetadata(
           Dict(
               "credits" => "Envifate - Francesco Geri, Oscar Cainelli, Paolo Zatelli, Gianluca Salogni, Marco Ciolli - DICAM Università degli Studi di Trento - Regione Veneto",
               "modulo" => "Analisi ruscellamento",
               "descrizione" => "Analisi di ruscellamento di un inquinante attraverso il metodo della separazione delle componenti",
               "srs" => refsys,
               "data" => today()
           )
       )
    """
       band1 = agd.getband(target_ds, 1)
       agd.setnodatavalue!( band1, Float64(valNoData) )
       band = agd.read(band1)
       agd.fillraster!(band, valNoData)
       xsize = agd.width(band)
       ysize = agd.height(band)

       intervallo = getCellDims(dem)

    """ PRENDE LA PORZIONE DI `landcover` RAPPRESENTATA DA `area` """
       lc_clip_proc = processing.run('qgis:clip', {'INPUT':self.lc, 'OVERLAY':self.areastudio, 'OUTPUT':self.path_working+'/clip.gpkg'})
       lc_clip=QgsVectorLayer(lc_clip_proc['OUTPUT'], 'lc_clip', 'ogr')
       lc_clip.setCrs(self.source.crs())

       #path__layer_lc=lc_clip['OUTPUT'].dataProvider().dataSourceUri()
       path__layer_lc=lc_clip.dataProvider().dataSourceUri()
       path_lc=path__layer_lc.split("|")
       source_ds_lc = ogr.Open(path_lc[0])
       lc_layer = source_ds_lc.GetLayer()
    """"""

       landcover_ds = agd.create( path_temp_landcover, gtiff_driver, round(Int64, x_res), round(Int64, y_res), 1, agd.GDAL.GDT_Float32 )
       agd.setgeotransform!(landcover_ds, [ x_min, 25.0, 0.0, y_max, 0.0, -25.0 ])
       agd.setproj!(landcover_ds, refsys)
       band_lc = agd.getband(landcover_ds, 1)
       agd.setnodatavalue!( band_lc, Float64(valNoData) )
       bandlc = agd.read(band_lc)
       agd.fillraster!(bandlc, valNoData)
       xsize = agd.width(band)
       ysize = agd.height(band)

    """ TRASFORMA IN RASTER LA PORZIONE DI `landcover` PRESA PRIMA """
       gdal.RasterizeLayer(lc_ds, [1], lc_layer,options=["ATTRIBUTE="+self.text_lcfield])
       lc_ds=None

       lc_layer=QgsRasterLayer(self.path_temp_lc,"lc_layer")
    """"""

       if soil_text == "Valore campo"
           soil_control = 1

        # NON SO SE SIA EQUIVALENTE
           #   source_ds_soil = ogr.Open(path_lc[0])
           #   soil_layer = source_ds_soil.GetLayer()
           agd.getlayer(path_lc, 0)

           soil_ds = agd.create( path_temp_soil, gtiff_driver, round(Int64, x_res), round(Int64, y_res), 1, agd.GDAL.GDT_Float32 )
           agd.setgeotransform!(soil_ds, [ x_min, 25.0, 0.0, y_max, 0.0, -25.0 ])
           agd.setproj!(soil_ds, refsys)
           band_sl = agd.getband(soil_ds, 1)
           agd.setnodatavalue!( band_sl, Float64(valNoData) )
           bandsl = agd.read(band_sl)
           agd.fillraster!(bandsl, valNoData)

        """ RASTERIZZA IL  VETTORIALE DEL TIPO DI SUOLO """
           gdal.RasterizeLayer(lc_ds, [1], soil_layer,options=["ATTRIBUTE="+self.text_soilfield])
           soil_ds=None

           soil_layer=QgsRasterLayer(self.path_temp_soil,"soil_layer")
        """"""
       end

    """ CALCOLA UN VETTORILE TRAMITE `r.drain` """
       grass_area=str(x_min)+','+str(x_max)+','+str(y_min)+','+str(y_max)+' ['+str(self.areastudio.crs().authid())+']'
       # grass_coord=str(x_source)+','+str(y_source)+' ['+str(self.source.crs().authid())+']'

       namewatershed=self.path_working+'/watershed'
       namedrain=self.path_working+'/wshed.shp'

       params = { 'GRASS_RASTER_FORMAT_OPT' : '','GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' :grass_area,
                  'GRASS_VECTOR_EXPORT_NOCAT' : False, '-a' : False, 'start_coordinates' : None,  '-n' : False,
                  'input' : self.dem.dataProvider().dataSourceUri(),'-c' : True, 'drain' : namedrain, 'GRASS_MIN_AREA_PARAMETER' : 0.0001,
                  'start_points' : self.source.dataProvider().dataSourceUri(), 'output' : namewatershed }

       waterwshed_proc = processing.run('grass7:r.drain', params)

       #aggiungo per controllo la viewshed alla toc
       #iface.addVectorLayer(namedrain,'watershed','ogr')
       #watershed=QgsProject.instance().mapLayersByName('watershed.shp')

       vdrain = QgsVectorLayer(namedrain, 'vdrain', 'ogr')

       idxlevel = self.source.fields().indexFromName(self.text_p)
       idxcat = vdrain.fields().indexFromName('cat')

       idxtargetname = self.target.fields().indexFromName(self.text_targetfield)
    """"""

       features = agd.getgeom.(collect(agd.features(vdrain)))
       nfeat = 0
       polygons_t = collect(agd.getfeature(target))

       #   start_time = time.time()

       for f in features
           length = agd.geomlength(f)
           currentdistance = intervallo
           nfeat += 1
           featlines = []

           #   firstpoint=geom.interpolate(0)
           firstpoint = agd.getpoint(geom, 0)
           old_x = agd.getx(firstpoint, 0)
           old_y = agd.gety(firstpoint, 0)

           fileoutput = folder * "\\drain$nfeat.shp"

        """ NON SO COSA FACCIA STA ROBA """
           vline = QgsVectorLayer("LineString?crs=EPSG:"+self.refsys, "drain"+str(nfeat), "memory")

           prline = vline.dataProvider()
           prlfield=prline.addAttributes( [ QgsField("concentrazione", QVariant.Double) ] )

           idf=f.attributes()[idxcat]
           feat_drain = next(self.source.getFeatures(QgsFeatureRequest().setFilterFid(idf-1)))
           p00=feat_drain.attributes()[idxlevel]
        """"""

           index_progress = 0
           while currentdistance < length
               if index_progress == 0
                   p0 = p00
               else
                   p0 = pe
               end

               #   point = geom.interpolate(currentdistance)
               point = agd.getpoint(geom, currentdistance)
               x = agd.getx(point, 0)
               y = agd.gety(point, 1)
               clc = lc_layer[x, y]
               if soil_control == 1
                   r, c = toIndexes(soil_layer, x, y)
                   soil = soil_layer[r, c]
               else
                   soil = text_soil
               end
               try
                   cn = listaclc[clc][classisoil[soil]]
                   S = 254.0((100 / cn) - 1)
               catch
                   S = 0
               end

               pcheck = (p0 - 0.2S)^2 / (p0 - 0.2S + S)
               if pcheck > 0.2S
                   pe = pcheck

                """ NON SO COSA FACCIA STA ROBA """
                  fetline = QgsFeature()
                  fetline.setGeometry( QgsGeometry.fromPolyline( [QgsPoint(old_x,old_y),QgsPoint(x,y)] ))
                  fetline.initAttributes(1)
                  fetline.setAttribute(0,pe)
                  vline.updateFeature(fetline)

                  featlines.append(fetline)
                """"""

                   index_progress += 1

                   for polygon in polygons_t
                       p_geom = agd.getgeom(polygon)
                       if agd.within(point, p_geom)
                           nometarget = pol_t.attributes()[idxtargetname]
                        """ MESSAGIO
                           messaggio = "\nIl vettore drain$nfeat ha raggiunto l'area bersaglio denominata $nometarget con un volume pari a: $(round(pe,3))mm\n"
                           self.console.appendPlainText(messaggio)
                        """
                           currentdistance = length + 1
                       end
                   end
               else
                   pe = 0
                   currentdistance = length + 1
               end
                   old_x = x
                   old_y = y
                   currentdistance += intervallo
               end

            """ NON SO COSA FACCIA STA ROBA """
               prline.addFeatures(featlines)
               vline.updateFields()
               QgsProject.instance().addMapLayer(vline)
            """"""

       end
   end
=#

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
## ------------------------------------------------------------------
_crossref_meta_dir() = joinpath(depotdir(), ".crossref")

function _download_crossref_meta(doi; dir = pwd(), force = false)
    name = string(_format_doi_for_filename(doi), ".json")
    outpath = joinpath(dir, name)
    force && rm(outpath; force)
    if !isfile(outpath)
        doi = _doi_to_url(doi)
        isempty(doi) && error("We don't like your DOI: '$(doi)'")
        # TODO: Understand what is going on here!
        mkpath(dir)
        @info("Downloading", doi)
        # cmdstr = """curl -L -H "Accept: application/vnd.crossref.unixsd+json" "$(doi)" > "$(outpath)" """
        cmdstr = """curl -G "https://api.crossref.org/works/$(doi)" > "$(outpath)" """
        run(`bash -c $(cmdstr)`; wait = true)
        sleep(0.2) # Avoid loading the line
        @info("Done")
    end
    return outpath
end

## ------------------------------------------------------------------
function _crossref_references(doi::String; force = false)
    outdir = _crossref_meta_dir()
    jsonfile = _download_crossref_meta(doi; force, dir = outdir)
    try
        jsonstr = read(jsonfile, String)
        jsondict = JSON.parse(jsonstr)
        val = _find_key(jsondict, ["message", "reference"])
        return isnothing(val) ? [] : val
    catch err
        rm(jsonfile; force = true)
        rethrow(err)
    end
    return []
end

## ------------------------------------------------------------------
# Access
# Interface with crossref api json

# const _DOI_REGEX = Regex("\\d{1,10}\\.\\d{3,20}/[-._;()/:a-zA-Z0-9]+")
# function _crossref_entry_doi(ref::AbstractDict)
#     !haskey(ref, "DOI") && return ""
#     doistr = string(ref["DOI"])
#     m = match(_DOI_REGEX, doistr)
#     isnothing(m) ? "" : string(m.match)
# end

_crossref_entry_doi(ref::AbstractDict) = get(ref, "DOI", "")
_crossref_entry_author(ref::AbstractDict) = get(ref, "author", "")
_crossref_entry_year(ref::AbstractDict) = get(ref, "year", "")
_crossref_entry_keywords(ref::AbstractDict) = get(ref, "keywords", "")

function _crossref_entry_title(ref::AbstractDict)
    hint = "-title"
    for key in keys(ref)
        key == "journal-title" && continue
        endswith(key, "-title") && return ref[key]
    end
    return ""
end

function _format_crossref_data!(ref::Dict)

    # clear content
    # doi
    if haskey(ref, "DOI")
        doistr = _doi_to_url(ref["DOI"])
        !isempty(doistr) && (ref["DOI"] = doistr)
    end

    for (key, val) in ref
        ref[key] = _clear_bibtex_entry(val)
    end
end

function _crossref_to_rbref(dict::AbstractDict)
    _format_crossref_data!(dict)

    ref = RBRef()

    # add firlds
    set_bibkey!(ref, "")
    set_author!(ref, _crossref_entry_author(dict))
    set_year!(ref, _crossref_entry_year(dict))
    set_title!(ref, _crossref_entry_title(dict))
    set_doi!(ref, _crossref_entry_doi(dict))
    add_tag!(ref, _crossref_entry_keywords(dict))
    refdict!(ref, dict)

    return ref
end

## ------------------------------------------------------------------
function crossrefs(doi::String; force = false)
    ref_dicts = _crossref_references(doi; force)
    refs = RBRef[]
    for dict in ref_dicts
        push!(refs, _crossref_to_rbref(dict))
    end
    return reflist(refs)
end


text:
function test_df_index_integer()
    ts = TS(df_integer_index, 1)
    @test typeof(ts) == TSx.TS
    @test ts.coredata == df_integer_index
end

function test_df_index_timetype()
    ts = TS(df_timetype_index, 1)
    @test typeof(ts) == TSx.TS
    @test ts.coredata == df_timetype_index
end

function test_df_index_symbol()
    ts = TS(df_integer_index, :Index)
    @test typeof(ts) == TSx.TS
    @test ts.coredata == df_integer_index
end

function test_df_index_string()
    ts = TS(df_integer_index, "Index")
    @test typeof(ts) == TSx.TS
    @test ts.coredata == df_integer_index
end

function test_df_index_range()
    ts = TS(df_vector, index_range)
    @test typeof(ts) == TSx.TS
    @test ts.coredata[!, :data] == df_vector[!, :data]
end

function test_vector_index_vector_integer()
    ts = TS(data_vector, index_integer)
    @test typeof(ts) == TSx.TS
    @test ts.coredata[!, :Index] == index_integer
    @test ts.coredata[!, 2] == data_vector
end

function test_vector_index_vector_timetype()
    ts = TS(data_vector, index_timetype)
    @test typeof(ts) == TSx.TS
    @test ts.coredata[!, :Index] == index_timetype
    @test ts.coredata[!, 2] == data_vector
end

function test_vector()
    ts = TS(data_vector)
    @test typeof(ts) == TSx.TS
    @test ts.coredata[!, :Index] == collect(1:length(data_vector))
    @test ts.coredata[!, 2] == data_vector
end

function test_array()
    ts = TS(data_array)
    @test typeof(ts) == TSx.TS
    @test typeof(ts.coredata) == DataFrames.DataFrame
    @test ts.coredata[!, :Index] == collect(1:size(data_vector)[1])
    @test Matrix(ts.coredata[!, Not(:Index)]) == data_array
end

# Run each test
# NOTE: Do not forget to add any new test-function created above
# otherwise that test won't run.
test_df_index_integer()
test_df_index_timetype()
test_df_index_symbol()
test_df_index_string()
test_df_index_range()
test_vector_index_vector_integer()
test_vector_index_vector_timetype()
test_vector()
test_array()


text:
using Test

@testset "Unit test" begin
    include("unittest.jl")
end

@testset "milp example" begin
    include("../examples/milp.jl")
    @test tree.best_incumbent == -3.0
end

@testset "minlp example" begin
    include("../examples/minlp.jl")
    @test isapprox(tree.best_incumbent, -1)
end


text:
module DemonstratingLinearInterpolatorsInterpolations

using Printf
using LazyAlgebra
using LazyAlgebra: Direct, Adjoint
using LinearInterpolators
using InterpolationKernels
using LinearInterpolators.Interpolations

using PyCall
# pygui(:gtk); # can be: :wx, :gtk, :qt
using LaTeXStrings
import PyPlot
const plt = PyPlot

function rundemos(::Type{T} = Float64) where {T<:AbstractFloat}

    z = T[0.5, 0.3, 0.1, 0.0, -0.2, -0.7, -0.7, 0.0, 1.7, 1.9, 2.1]

    # 1-D example
    t = range(-3, 14, length=2000);
    I0 = SparseInterpolator(RectangularSpline(T, Flat), t, length(z))
    I1 = SparseInterpolator(LinearSpline(T, Flat), t, length(z))
    I2 = SparseInterpolator(CatmullRomSpline(T, SafeFlat), t, length(z))
    I3 = SparseInterpolator(LanczosKernel(T, 8), t, length(z))
    plt.figure(1)
    plt.clf()
    plt.plot(t, I0(z), color="darkgreen",
             linewidth=1.5, linestyle="-", label="Rectangular spline");
    plt.plot(t, I1(z), color="darkred",
                 linewidth=1.5, linestyle="-", label="Linear spline");
    plt.plot(t, I2(z), color="orange",
             linewidth=1.5, linestyle="-", label="Catmull-Rom");
    plt.plot(t, I3(z), color="violet",
             linewidth=1.5, linestyle="-", label="Lanczos 8");
    plt.plot(1:length(z), z, "o", color="darkblue", label="Points")
    plt.legend()
    plt.title("Interpolations with Cardinal Kernels")

    # Interpolation + smoothing
    t = range(-3, 14, length=2000);
    I5 = SparseInterpolator(QuadraticSpline(T), t, length(z))
    I6 = SparseInterpolator(CubicSpline(T, Flat), t, length(z))
    I7 = SparseInterpolator(MitchellNetravaliSpline(T, Flat), t, length(z))
    plt.figure(2)
    plt.clf()
    plt.plot(t, I6(z), color="darkred",
                 linewidth=1.5, linestyle="-", label="Cubic B-spline");
    plt.plot(t, I5(z), color="darkgreen",
             linewidth=1.5, linestyle="-", label="Quadratic B-spline");
    #plt.plot(t, I7(z), color="orange",
    #         linewidth=1.5, linestyle="-", label="Mitchell & Netravali");
    plt.plot(1:length(z), z, "o", color="darkblue", label="Points")
    plt.legend()
    plt.title("Interpolations with Smoothing Kernels")

    # Test conversion to a sparse matrix.
    S1 = sparse(I1)
    S2 = sparse(I2)
    println("max. error 1: ", maximum(abs.(I1(z) - S1*z)))
    println("max. error 2: ", maximum(abs.(I2(z) - S2*z)))

    # 2-D example
    t = reshape(range(-3, 14, length=20*21), (20,21));
    I = SparseInterpolator(CatmullRomSpline(T, Flat), t, length(z))
    plt.figure(3)
    plt.clf()
    plt.imshow(I(z));
end
#rundemos()
#------------------------------------------------------------------------------

const CPU = 2.94e9

function printtest(prefix::String, count::Integer, tup)
    elt, mem, gct = tup[2:4]
    @printf("%s   %6.0f cycles (Mem. = %8d bytes   GC = %3.0f%%)\n",
            prefix, CPU*elt/count, mem, 1e2*gct/elt)
end

function testdirect!(dst::Vector{T}, ker::Kernel{T,S,B},
                     x::AbstractVector{T}, src::Vector{T},
                     cnt::Integer) where {T,S,B}
    for k in 1:cnt
        apply!(dst, Direct, ker, x, src)
    end
    return dst
end

function testadjoint!(dst::Vector{T}, ker::Kernel{T,S,B},
                      x::AbstractVector{T}, src::Vector{T},
                      cnt::Integer) where {T,S,B}
    for k in 1:cnt
        apply!(dst, Adjoint, ker, x, src)
    end
    return dst
end

function runtests(::Type{T} = Float64,
                  len::Integer = 1000, cnt::Integer = 100) where T<:AbstractFloat
    dim = 500
    w = 10
    K01 = RectangularSpline(T, Flat)
    K02 = RectangularSpline(T, SafeFlat)
    K03 = LinearSpline(T, Flat)
    K04 = LinearSpline(T, SafeFlat)
    K05 = QuadraticSpline(T, Flat)
    K06 = QuadraticSpline(T, SafeFlat)
    K07 = CardinalCubicSpline(T, -1, Flat)
    K08 = CardinalCubicSpline(T, -1, SafeFlat)
    K09 = CatmullRomSpline(T, Flat)
    K10 = CatmullRomSpline(T, SafeFlat)
    x = rand(T(1 - w):T(dim + w), len)
    y = randn(T, dim)
    z = Array{T}(undef, len)
    # compile methods
    testdirect!(z,  K01, x, y, 3)
    testadjoint!(y, K01, x, z, 3)
    testdirect!(z,  K02, x, y, 3)
    testadjoint!(y, K02, x, z, 3)
    testdirect!(z,  K03, x, y, 3)
    testadjoint!(y, K03, x, z, 3)
    testdirect!(z,  K04, x, y, 3)
    testadjoint!(y, K04, x, z, 3)
    testdirect!(z,  K05, x, y, 3)
    testadjoint!(y, K05, x, z, 3)
    testdirect!(z,  K06, x, y, 3)
    testadjoint!(y, K08, x, z, 3)
    testdirect!(z,  K07, x, y, 3)
    testadjoint!(y, K07, x, z, 3)
    testdirect!(z,  K08, x, y, 3)
    testadjoint!(y, K08, x, z, 3)
    testdirect!(z,  K09, x, y, 3)
    testadjoint!(y, K09, x, z, 3)
    testdirect!(z,  K10, x, y, 3)
    testadjoint!(y, K10, x, z, 3)
    # warm-up
    tmp = randn(T, 1_000_000)
    # run tests
    n = cnt*len
    printtest("RectangularSpline   Flat     direct ", n, @timed(testdirect!(z,  K01, x, y, cnt)))
    printtest("RectangularSpline   Flat     adjoint", n, @timed(testadjoint!(y, K01, x, z, cnt)))
    printtest("RectangularSpline   SafeFlat direct ", n, @timed(testdirect!(z,  K02, x, y, cnt)))
    printtest("RectangularSpline   SafeFlat adjoint", n, @timed(testadjoint!(y, K02, x, z, cnt)))
    println()
    printtest("LinearSpline        Flat     direct ", n, @timed(testdirect!(z,  K03, x, y, cnt)))
    printtest("LinearSpline        Flat     adjoint", n, @timed(testadjoint!(y, K03, x, z, cnt)))
    printtest("LinearSpline        SafeFlat direct ", n, @timed(testdirect!(z,  K04, x, y, cnt)))
    printtest("LinearSpline        SafeFlat adjoint", n, @timed(testadjoint!(y, K04, x, z, cnt)))
    println()
    printtest("QuadraticSpline     Flat     direct ", n, @timed(testdirect!(z,  K05, x, y, cnt)))
    printtest("QuadraticSpline     Flat     adjoint", n, @timed(testadjoint!(y, K05, x, z, cnt)))
    printtest("QuadraticSpline     SafeFlat direct ", n, @timed(testdirect!(z,  K06, x, y, cnt)))
    printtest("QuadraticSpline     SafeFlat adjoint", n, @timed(testadjoint!(y, K06, x, z, cnt)))
    println()
    printtest("CardinalCubicSpline Flat     direct ", n, @timed(testdirect!(z,  K07, x, y, cnt)))
    printtest("CardinalCubicSpline Flat     adjoint", n, @timed(testadjoint!(y, K07, x, z, cnt)))
    printtest("CardinalCubicSpline SafeFlat direct ", n, @timed(testdirect!(z,  K08, x, y, cnt)))
    printtest("CardinalCubicSpline SafeFlat adjoint", n, @timed(testadjoint!(y, K08, x, z, cnt)))
    println()
    printtest("CatmullRomSpline    Flat     direct ", n, @timed(testdirect!(z,  K09, x, y, cnt)))
    printtest("CatmullRomSpline    Flat     adjoint", n, @timed(testadjoint!(y, K09, x, z, cnt)))
    printtest("CatmullRomSpline    SafeFlat direct ", n, @timed(testdirect!(z,  K10, x, y, cnt)))
    printtest("CatmullRomSpline    SafeFlat adjoint", n, @timed(testadjoint!(y, K10, x, z, cnt)))

end

end # module


text:
import .Winston

Image(height, width, a::Winston.FramedPlot) =
    Image(height, width, ImageFileData(a))

function ImageFileData(a::Winston.FramedPlot)
    filename = tempname()*".pdf"
    Winston.savefig(a, filename)
    r = read(filename)
    rm(filename)
    ImageFileData(r, :pdf)
end


text:
module TestPkg

using Checkpoints: register, checkpoint, with_checkpoint_tags, Session

# We aren't using `@__MODULE__` because that would return TestPkg on 0.6 and Main.TestPkg on 0.7
const MODULE = "TestPkg"

__init__() = register(MODULE, ["foo", "bar", "baz", "qux_a", "qux_b", "deprecated"])

function foo(x::Matrix, y::Matrix)
    # Save multiple variables to 1 foo.jlso file by passing in pairs of variables
    checkpoint(MODULE, "foo", :x => x, :y => y)
    return x * y
end

function bar(a::Vector)
    # Save a single value for bar.jlso. The object name in that file defaults to :data.
    with_checkpoint_tags(:date=>"2017-01-01") do
        checkpoint(MODULE, "bar", a)
        return a * a'
    end
end

function baz(data::Dict)
    # Check that saving multiple values to a Session works.
    Session(MODULE, "baz") do s
        for (k, v) in data
            checkpoint(s, k => v)
        end
    end
end

function qux(a::Dict, b::Vector)
    # Check that saving multiple values to multiple Sessions also works.
    Session(MODULE, ["qux_a", "qux_b"]) do sa, sb
        for (k, v) in a
            checkpoint(sa, k => v)
        end

        checkpoint(sb, b)
    end
end

function deprecated_checkpoint_syntax()
    checkpoint(MODULE, "deprecated", [1, 2, 3]; date="2017-01-01")
end

end


text:
using ForwardDiffOverMeasurements
using Measurements
using ForwardDiff
using Test

f1(x,y) = x+y
f2(x,y) = x-y
f3(x,y) = x*y
f4(x,y) = x/y

@testset "ForwardDiffOverMeasurements.jl" begin
    x1 = 1.0 ± 0.1
    y1 = 2.0 ± 0.001
    for op in (:f1,:f2,:f3,:f4)
        @eval begin
            @test ForwardDiff.derivative(x->$(op)(x,$y1),$x1) isa Measurement
            @test ForwardDiff.derivative(y->$(op)($x1,y),$y1) isa Measurement
        end
    end
end


text:
using Test
include(joinpath(@__DIR__, "../src/FLS.jl"))
include(joinpath(@__DIR__, "../data/sketches/fg.jl"))

f2 = fls_from_json(to_json(fls))
@test fls_from_json(to_json(fls)) == fls

xmod,xrel = fls.cset(), fls.crel()
[add_part!(xmod, x) for x in [:A,:B,:C]]
[set_subpart!(xmod, 1, x, 1) for x in [:f, :g, :c]]

xmodrel = initrel(fls, xmod)
@test crel_to_cset(fls, xmodrel) == xmod

@test isempty(ChaseStepData())


text:
# Autogenerated wrapper script for FastTransforms_jll for i686-w64-mingw32
export libfasttransforms

using CompilerSupportLibraries_jll
using FFTW_jll
using MPFR_jll
using OpenBLAS_jll
JLLWrappers.@generate_wrapper_header("FastTransforms")
JLLWrappers.@declare_library_product(libfasttransforms, "libfasttransforms.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FFTW_jll, MPFR_jll, OpenBLAS_jll)
    JLLWrappers.@init_library_product(
        libfasttransforms,
        "bin\\libfasttransforms.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export ItemMean

"""
    ItemMean(data::DataAccessor)

Recommend based on global item mean rating.
"""
struct ItemMean <: Recommender
    data::DataAccessor
    scores::AbstractVector

    function ItemMean(data::DataAccessor)
        n_item = size(data.R, 2)
        new(data, vector(n_item))
    end
end

isbuilt(recommender::ItemMean) = isfilled(recommender.scores)

function build!(recommender::ItemMean)
    n_item = size(recommender.data.R, 2)

    for i in 1:n_item
        v = recommender.data.R[:, i]
        recommender.scores[i] = sum(v) / count(!iszero, v)
    end
end

function predict(recommender::ItemMean, u::Int, i::Int)
    check_build_status(recommender)
    recommender.scores[i]
end


text:
export dde_stability

"""
    dde_stability(A, B, τ; N=10)

Calculate (a finite number of) the characteristic roots of the linear delay differential
equation

```math
x'(t) = Ax(t) + Bx(t-τ)
```

where `A` and `B` can be (square) matrices. `N` controls the order of the Chebyshev
polynomial approximation to the infinitesimal generator; higher values of `N` result in a
more accurate approximation. A property of the algorithm used is that characteristic roots
closest to the origin are approximated best.

## Reference

Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential
Equations; D. Breda, S. Maset, and R. Vermiglio, SIAM Journal on Scientific Computing 2005
27:2, 482-495
"""
function dde_stability(A, B, τ; N=10)
    @assert size(A) == size(B)
    @assert size(A, 1) == size(A, 2)
    T = eltype(A) == Int ? Float64 : eltype(A)
    n = size(A, 1)
    (x, D) = cheb_diff(T, N)
    M = kron(-D, Matrix(I / convert(T, τ), n, n))
    M[1:n, 1:n] .= A
    M[1:n, (n + 1):(end - n)] .= 0
    M[1:n, (end - n + 1):end] .= B
    return eigen(M)
end

# Tests
# 1. x'(t) = 2x(t) - exp(1)*x(t - 1) has a double root at λ = 1
# 2. x'(t) = -x(t - τ) goes unstable at τ=π/2
# 3. x'(t) = α*A*x(t) + (1-α)*A*x(t-τ)
#    A = [
#          50  284   41   23   50   32
#        -280  -46  -19  -37  -10  -28
#          35   -1   26  143   35   17
#          5   -31 -139  -22    5  -13
#          20  -16   11   -7 -115  137
#         -10  -46  -19  -37 -145 -163
#    ]
#    with α=0.9 is stable for 3π/(270*(2*0.9 - 1)) ≤ τ ≤ 4π/270
#    that is 0.04363323129985824 ≤ τ ≤ 0.046542113386515455


text:
using Documenter, TextAnalysis

makedocs(
    modules = [TextAnalysis],
    sitename = "TextAnalysis",
    format = Documenter.HTML(
        canonical = "https://juliatext.github.io/TextAnalysis.jl/stable/",
        #prettyurls = false,
    ),
    assets = ["assets/custom.css", "assets/custom.js", "assets/favicon.ico"],
    pages = [
        "Home" => "index.md",
        "Documents" => "documents.md",
        "Corpus" => "corpus.md",
        "Features" => "features.md",
        "Semantic Analysis" => "semantic.md",
        "Classifier" => "classify.md",
        "Extended Example" => "example.md",
        "Evaluation Metrics" => "evaluation_metrics.md",
        "Conditional Random Fields" => "crf.md",
        "Named Entity Recognition" => "ner.md"
    ],
)

deploydocs(repo = "github.com/JuliaText/TextAnalysis.jl.git")


text:
@userplot SpectralPortrait

"""
    spectralportrait(A; npts=100)

Plot the spectral portrait of the square matrix `A`.

# Keyword arguments
* `npts::Integer = 100` The size of the grid used for computing the levels of the pseudospectrum
# Examples
```jldoctest
julia> A = Pseudospectra.grcar(40)
julia> spectralportrait(A)
```
"""
function spectralportrait end

RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
    if length(p.args) != 1 || !(typeof(p.args[1]) <: AbstractMatrix)
        error("spectralportrait must be given an AbstractMatrix. Got $(typeof(p.args))")
    end

    # Below is the computation of of the spectral portrait taken from the base implementation
    A₀ = p.args[1]

    local ps_data
    try
        ps_data = new_matrix(A₀)
    catch JE
        @warn "spectralportrait only works for simple cases."
        rethrow(JE)
    end

    A = ps_data.matrix
    ps_dict = ps_data.ps_dict
    B = get(ps_dict, :matrix2, I)
    eigA = ps_dict[:ews]
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    if isempty(zoom.ax)
        zoom.ax = vec2ax(eigA)
    end
    psa_opts = _basic_psa_opts(zoom, ps_dict)

    Z, xs, ys, t_levels, err, Tproj, eigAproj, algo = psa_compute(A, npts, zoom.ax, eigA, psa_opts, B)

    # Recipe code below

    @series begin
        seriestype := :contour
        #linecolor --> eigtool_cgrad
        xs, ys, log10.(Z)
    end

    @series begin
        seriestype := :scatter
        markercolor := :black
        markersize := 2
        label --> "eigvals"
        real.(eigA), imag.(eigA)
    end

end