text:
using Distances
using Statistics
using StatsBase
using Random
using ProgressMeter
using CalculatedABC
using JLD

function impurity(Y::Vector{Float64})
    if length(Y)<=1
        return 0
    end
    return var(Y)
end

function impurity(Y::Vector{Int})
    return gini_coeff(Y)
end

function findall3(f, a::AbstractArray{T, N}) where {T, N}
    j = 1
    b = Vector{Int}(undef, length(a))
    @inbounds for i in eachindex(a)
        @inbounds if f(a[i])
            b[j] = i
            j += 1
        end
    end
    resize!(b, j-1)
    sizehint!(b, length(b))
    return b
end

function impurity_split(Y::Vector{T},sp::AbstractVector{Float64}) where T
    l::Int64=length(Y)
    gauche::Vector{Int64} = findall3(x->x==1,sp)
    droite::Vector{Int64} = findall3(x->x==2,sp)
    if length(gauche)==0 || length(droite)==0
        return Inf
    end
    return impurity(Y[gauche])*length(gauche)/l + impurity(Y[droite])*length(droite)/l
end

### Maintenant la fonction de découpage aléatoire :::

function ERvar_split(Y::Vector{T},Im::Array{Float64,3}, ntry::Int64, dim::Tuple{Int,Int,Int}, dist) where T

    l::Int64=dim[2]

    if l==2

        D::Vector{Float64} = zeros(dim[3])
        @inbounds for p in 1:dim[3]
            D[p]= @views mean(sqrt.((Im[:,1,p]-Im[:,2,p]).^2))
        end
        return findmax(D)[2], [1, 2], [1 ,2], 0

    else
        impur = fill(Inf,dim[3])
        Splits_prime = fill(1,dim[3],dim[2])
        centers = fill(0,2, dim[3])
        whichsplit = 0
        n_dec = min(ntry,l)
        centers_prime = fill(0,2,n_dec)

        for p in 1:dim[3]

            ## plus de tests que d'images ??
            splits = ones(n_dec, l)
            impur_prime = fill(Inf, n_dec)

            for k in 1:n_dec
                centers_prime[:,k] = sample(1:dim[2], 2, replace=false)
                d_g = @views Distances.colwise(dist, Im[:,:,p], Im[:,centers_prime[1,k],p]) ## distance à gauche
                d_d = @views Distances.colwise(dist, Im[:,:,p], Im[:,centers_prime[2,k],p]) ## le fameux dd
                qui = findall3(x->x>=0, d_g-d_d)::Vector{Int}
                splits[k,qui].=  2
                impur_prime[k]= @views impurity_split(Y,splits[k,:])::Float64
            end

            impur[p], whichsplit = findmin(impur_prime)
            Splits_prime[p,:] = @view splits[whichsplit,:]
            centers[:,p]= @view centers_prime[:,whichsplit]
        end

        Imp, whichsplit = findmin(impur)
        return whichsplit, centers[:,whichsplit], Splits_prime[whichsplit,:], Imp
    end

end

### Maintenant on construit les arbres maximaux ::

function ERtmax(Im::Array{Float64,3},Y::Vector{Float64},mtry::Int64,ntry::Int64, minElem::Int,dim::Tuple{Int,Int,Int}, dist)
     ## On va éviter de calculer ça tout le temps !!!! A optimiser !!
    p::Int=dim[3]
    n::Int64= dim[2]
    boot::Vector{Int64} = unique(sample(1:n,n))
    V_split::Array{Float64,2} = zeros(4,length(boot)-1)
    Pred::Array{Float64,2}= ones(2,2*length(boot)-1)

    ### Il faut ensuite construire les données d'apprentissage :
    X_boot::Array{Float64,3} = @view Im[:,boot,:]

    feuilles::Vector{Float64} = ones(length(boot))
    feuilles_prime::Vector{Float64} = feuilles
    ### On passe ensuite à la boucle ::
    decoupe::Int=0
    courant::Int=1
    courant_pred::Int=2

    Pred[2,1]= mean(Y)

    @inbounds for i in 1:length(boot)-1
        ## On va tirer un sous-échantillon de variables possibles ::
        feuilles_prime = feuilles
        if decoupe == length(unique(feuilles))
            break

        else
            @inbounds for j in unique(feuilles)

                qui= collect(findall3(x->x==j, feuilles))
                V::Vector{Int}= sort(sample(1:p, mtry, replace=false))

                ### Maintenant il faut se restreindre aux éléments dans la feuilles
                if length(qui)>minElem
                    Split = ERvar_split(Y[boot][qui],X_boot[:,qui,V], ntry,(dim[1],length(qui),length(V)), dist)
                    gauche::Vector{Int}=findall3(x->x==1, Split[3])
                    droite::Vector{Int}= findall3(x->x==2, Split[3])

                    if length(gauche)>0 && length(droite)>0
                        V_split[1,courant] = j
                        V_split[2,courant] = V[Split[1]]
                        V_split[3,courant] = boot[qui[Split[2][1]]]
                        V_split[4,courant] = boot[qui[Split[2][2]]]

                        feuilles_prime[qui[gauche]].= 2*j#*ones(length(gauche))
                        feuilles_prime[qui[droite]].= (2*j + 1)#*ones(length(droite))

                        Pred[1,courant_pred]=2*j
                        Pred[1,courant_pred+1]=2*j+1

                        Pred[2,courant_pred]= mean(@view Y[boot][qui[gauche]])
                        Pred[2,courant_pred+1] = mean(@view Y[boot][qui[droite]])

                        courant += 1
                        courant_pred +=2
                    else
                        decoupe +=1
                    end

                else
                    decoupe = decoupe+1
                end
            end

            if decoupe < length(unique(feuilles))
                decoupe=0
            end
            feuilles = feuilles_prime
        end
    end
    return V_split[:,1:courant-1], Pred[:,1:courant_pred-1], boot
end

function ERtmax(Im::Array{Float64,3},Y::Vector{Int64},mtry::Int64,ntry::Int64, minElem::Int,dim::Tuple{Int,Int,Int}, dist)
    ## On va éviter de calculer ça tout le temps !!!! A optimiser !!
   p=dim[3]
   n= dim[2]
   boot = unique(sample(1:n,n))
   V_split = zeros(4,length(boot)-1)
   Pred= ones(Int64, 2,2*length(boot)-1)

   ### Il faut ensuite construire les données d'apprentissage :
   X_boot = @view Im[:,boot,:]

   feuilles = ones(length(boot))
   feuilles_prime = feuilles
   ### On passe ensuite à la boucle ::
   decoupe=0
   courant=1
   courant_pred=2

   Pred[2,1]= findmax(countmap(Y))[2]

   @inbounds for i in 1:length(boot)-1
       ## On va tirer un sous-échantillon de variables possibles ::
       feuilles_prime = feuilles
       if decoupe == length(unique(feuilles))
           break

       else
           @inbounds for j in unique(feuilles)

               qui= collect(findall3(x->x==j, feuilles))

               impur = gini_coeff(Y[boot][qui])

               if impur > 0

                    V= sort(sample(1:p, mtry, replace=false))
                    if length(qui)>minElem
                        Split = ERvar_split(Y[boot][qui],X_boot[:,qui,V], ntry,(dim[1],length(qui),length(V)), dist)
                        gauche::Vector{Int}=findall3(x->x==1, Split[3])
                        droite::Vector{Int}= findall3(x->x==2, Split[3])
                        if length(gauche)>0 && length(droite)>0
                            V_split[1,courant] = j
                            V_split[2,courant] = V[Split[1]]
                            V_split[3,courant] = boot[qui[Split[2][1]]]
                            V_split[4,courant] = boot[qui[Split[2][2]]]
                            feuilles_prime[qui[gauche]].= 2*j#*ones(length(gauche))
                            feuilles_prime[qui[droite]].= (2*j + 1)#*ones(length(droite))

                            Pred[1,courant_pred]=2*j
                            Pred[1,courant_pred+1]=2*j+1

                            Pred[2,courant_pred]= findmax(countmap(@view Y[boot][qui[gauche]]))[2]
                            Pred[2,courant_pred+1] = findmax(countmap(@view Y[boot][qui[droite]]))[2]

                            courant += 1
                            courant_pred +=2
                        else
                            decoupe +=1
                        end
                    else
                        decoupe = decoupe+1
                    end
                else
                    decoupe += 1
                end
            end

            if decoupe < length(unique(feuilles))
                decoupe=0
            end

           feuilles = feuilles_prime
       end
   end
   return V_split[:,1:courant-1], Pred[:,1:courant_pred-1], boot
end

### On va passer aux forêts aléatoires ::::

function ERFRF(Im::Array{Float64,3}, Y::Vector{Float64}, mtry::Int, ntree::Int, ntry::Int, minElem::Int, dist)
    ### Allez c'est parti !!
    ### on transforme les images en données spatio-Temporelles::
    #### On passe à la suite :::

    dim::Tuple{Int,Int,Int}=size(Im)
    boot::Array{Float64} = zeros(ntree, dim[2])
    trees::Array{Float64,3} = zeros(4,2*dim[2]-1,ntree)
    P::Array{Float64,3} = zeros(2,2*dim[2]-1,ntree)

    ### On va essayer la barre de progression ::
    p=Progress(ntree);
    update!(p,0)
    jj = Threads.Atomic{Int}(0)
    l = Threads.SpinLock()

    Threads.@threads for k in 1:ntree
        trees_temp , pred_temp, boot_temp = ERtmax(Im,Y,mtry,ntry,minElem, dim, dist)
        trees[:,1:size(trees_temp,2),k] = trees_temp
        P[:,1:size(pred_temp,2),k]= pred_temp
        boot[k,1:length(boot_temp)]= boot_temp
        ## on affiche la barre de progression ::
        Threads.atomic_add!(jj, 1)
        Threads.lock(l)
        update!(p, jj[])
        Threads.unlock(l)
    end

    return trees, P, boot
end

function ERFRF(Im::Array{Float64,3}, Y::Vector{Int64}, mtry::Int, ntree::Int, ntry::Int, minElem::Int, dist)
    ### Allez c'est parti !!
    ### on transforme les images en données spatio-Temporelles::
    #### On passe à la suite :::

    dim::Tuple{Int,Int,Int}=size(Im)
    boot::Array{Float64} = zeros(ntree, dim[2])
    trees::Array{Float64,3} = zeros(4,2*dim[2]-1,ntree)
    P::Array{Int64,3} = zeros(Int64,2,2*dim[2]-1,ntree)

    ### On va essayer la barre de progression ::
    p=Progress(ntree);
    update!(p,0)
    jj = Threads.Atomic{Int}(0)
    l = Threads.SpinLock()

    Threads.@threads for k in 1:ntree
        trees_temp , pred_temp, boot_temp = ERtmax(Im,Y,mtry,ntry,minElem, dim, dist)
        trees[:,1:size(trees_temp,2),k] = trees_temp
        P[:,1:size(pred_temp,2),k]= pred_temp
        boot[k,1:length(boot_temp)]= boot_temp
        ## on affiche la barre de progression ::
        Threads.atomic_add!(jj, 1)
        Threads.lock(l)
        update!(p, jj[])
        Threads.unlock(l)
    end

    return trees, P, boot
end

### On passe à la fonction de prédiction sur les arbres :::

function pred_tree(tree::AbstractArray{Float64,2},Pred::AbstractArray{T,2},X::AbstractArray{Float64,3},X_init::AbstractArray{Float64,3}, dist) where T

    pred::Vector{T} = zeros(T,size(X,2))
    nodes::Vector{Int64} = ones(Int64,size(X,2))
    feuilles::Vector{Int64}= zeros(Int64,size(X,2))

    #tree = convert(Array{Int64,3},tree) ## A tester pour voir si on peut pas couper dans le gras ?

    while sum(feuilles)<size(X,2)
        @inbounds for i in unique(nodes)
            qui::Vector{Int64} = findall3(x->x==i,nodes)
            col::Vector{Int64} = findall3(x->x==i,@view tree[1,:])
            if length(col)>0
                d_g::Vector{Float64}= @views Distances.colwise(dist, X[:,qui,convert(Int64,tree[2,col][1])], X_init[:,convert(Int64,tree[3,col][1]),convert(Int64,tree[2,col][1])])
                d_d::Vector{Float64}= @views Distances.colwise(dist, X[:,qui,convert(Int64,tree[2,col][1])], X_init[:,convert(Int64,tree[4,col][1]),convert(Int64,tree[2,col][1])])
                gauche::Vector{Int64} = findall3(x->x<=0, d_g-d_d)
                droite::Vector{Int64} = findall3(x->x>0, d_g-d_d)
                nodes[qui[gauche]] = 2*nodes[qui[gauche]]
                nodes[qui[droite]] = 2*nodes[qui[droite]] .+ 1

            elseif unique(feuilles[qui])[1]==0.0
                pred[qui] .=  @view Pred[2,findall3(x->x==i,@view Pred[1,:])]
                feuilles[qui] .= 1
            end

        end
    end
    return pred
end

function Importance(frf::Array{Float64,3},X::Array{Float64,3}, Y::Vector{Float64},P::Array{Float64,3}, id::Array{Float64,2}, dist)
    dim=size(X)::Tuple{Int,Int,Int}
    variables::Int = dim[3]
    imp::Vector{Float64}= zeros(dim[3])
    ntree::Int=size(frf,3)
    err_courante::Vector{Float64}=zeros(ntree)

    prog =Progress(variables);
    update!(prog ,0)

    ID::Array{Int64,2} = convert(Array{Int64,2},id)

    @inbounds for p in 1:variables
        err_courante.=0
        Threads.@threads for i in 1:ntree
            boot::Vector{Int64}= ID[i,findall3(x->x>0.0,ID[i,:])]
            ### Trouver les éléments qui ne sont pas dans l'échantillon OOB :
            OOB::Vector{Int64} = setdiff(1:dim[2],boot)
            X_permute::Array{Float64,3} =  X[:,OOB,:]
            X_permute[:,:,p] = @view X[:,OOB[randperm(length(OOB))],p]
            ## Il faut maintenant regarder la différence en erreur de prédiction ::
            err_courante[i] = @views mean((Y[OOB].-pred_tree(frf[:,:,i], P[:,:,i],X_permute,X, dist)).^2)-mean((Y[OOB].-pred_tree(frf[:,:,i],P[:,:,i],X[:,OOB,:],X, dist)).^2)
        end
        imp[p]= mean(err_courante)
        update!(prog, p)
    end
    return imp
end

function Importance(frf::Array{Float64,3},X::Array{Float64,3}, Y::Vector{Int64},P::Array{Int64,3}, id::Array{Float64,2}, dist)
    dim=size(X)::Tuple{Int,Int,Int}
    variables::Int = dim[3]
    imp::Vector{Float64}= zeros(dim[3])
    ntree::Int=size(frf,3)
    err_courante::Vector{Float64}=zeros(ntree)

    prog =Progress(variables);
    update!(prog ,0)

    ID::Array{Int64,2} = convert(Array{Int64,2},id)

    @inbounds for p in 1:variables
        err_courante.=0
        Threads.@threads for i in 1:ntree
            boot::Vector{Int64}= ID[i,findall3(x->x>0.0,ID[i,:])]
            ### Trouver les éléments qui ne sont pas dans l'échantillon OOB :
            OOB::Vector{Int64} = setdiff(1:dim[2],boot)
            X_permute::Array{Float64,3} =  X[:,OOB,:]
            X_permute[:,:,p] = @view X[:,OOB[randperm(length(OOB))],p]
            ## Il faut maintenant regarder la différence en erreur de prédiction ::
            err_courante[i] = @views mean((Y[OOB].!=pred_tree(frf[:,:,i], P[:,:,i],X_permute,X, dist)))-mean((Y[OOB].!=pred_tree(frf[:,:,i],P[:,:,i],X[:,OOB,:],X, dist)))
        end
        imp[p]= mean(err_courante)
        update!(prog, p)
    end
    return imp
end

### On fait la fonction du calcul du ‰ de variance expliquée ::
### On va aussi rajouter les prediction OOB pour la forêt complète ::

function FRFERR(frf::Array{Float64,3},X::Array{Float64,3}, Y::Vector{Float64},P::Array{Float64,3}, id::Array{Float64,2}, dist)
    dim=size(X)::Tuple{Int,Int,Int}
    #variables::Int = dim[3]
    ntree::Int=size(frf,3)

    pred_OOB::Vector{Float64}=zeros(dim[2])
    ZZ::Array{Float64,3} = zeros(dim[1],2,dim[3])

    @inbounds for i in 1:dim[2]
        Pred_courante = fill(Inf,2,ntree)
        ZZ[:,1,:], ZZ[:,2,:] = X[:,i,:], X[:,i,:]
        Threads.@threads for k in 1:ntree
            if length(findall3(x->x==i,id[k,:]))==0
                Pred_courante[1,k] = 1.0
                Pred_courante[2,k] = @views pred_tree(frf[:,:,k],P[:,:,k],ZZ,X, dist)[1]
            end

            ## Il faut maintenant regarder la différence en erreur de prédiction ::
        end
        pred_OOB[i] = mean(@views Pred_courante[2,findall3(x->x<Inf,Pred_courante[1,:])])
    end

    mse = mean((Y.-pred_OOB).^2)

    varex = 1 - (mse/impurity(Y))
    return pred_OOB, mse, varex
end

function FRFERR(frf::Array{Float64,3},X::Array{Float64,3}, Y::Vector{Int64},P::Array{Int64,3}, id::Array{Float64,2}, dist)
    dim=size(X)
    #variables::Int = dim[3]
    ntree=size(frf,3)

    pred_OOB=zeros(Int64,dim[2])
    ZZ = zeros(dim[1],2,dim[3])

    @inbounds for i in 1:dim[2]
        Pred_courante = fill(-1,2,ntree)
        ZZ[:,1,:], ZZ[:,2,:] = X[:,i,:], X[:,i,:]
        Threads.@threads for k in 1:ntree
            if length(findall3(x->x==i,id[k,:]))==0
                Pred_courante[1,k] = 0
                Pred_courante[2,k] = @views pred_tree(frf[:,:,k],P[:,:,k],ZZ,X, dist)[1]
            end

            ## Il faut maintenant regarder la différence en erreur de prédiction ::
        end
        pred_OOB[i] = findmax(countmap(@views Pred_courante[2,findall3(x->x==0,Pred_courante[1,:])]))[2]
    end

    mse = mean(Y.!=pred_OOB)

    init_imp = mean(Y .!= findmax(countmap(Y))[2])

    varex = 1 - (mse/init_imp)
    return pred_OOB, mse, varex
end

### On regarde maintenant les predictions OOB de la forêt de Fréchet pour calculer le % de var expliquée:::

function Importance_Unique(frf::Array{Char},X::Array{Float64,3}, Y::Vector{Float64},variables::UnitRange{Int64}, dist)
    dim=size(X)::Tuple{Int,Int,Int}
    imp::Array{Float64,2}= zeros(2,length(variables))

    trees=readdir(frf; join = true)
    ntree::Int=length(trees)

    err_courante::Vector{Float64}=zeros(ntree)
    imp[1,:] = variables
    ID::Array{Int64,2} = convert(Array{Int64,2},id)

    @inbounds for p in 1:length(variables)
        err_courante.=0
        Threads.@threads for i in 1:ntree

            infos = load(trees[k])
            tree = infos["tree"]
            ID = infos["boot"]
            P = infos["P"]

            boot::Vector{Int64}= @view ID[i,findall3(x->x>0.0,ID[i,:])]
            ### Trouver les éléments qui ne sont pas dans l'échantillon OOB :
            OOB::Vector{Int64} = setdiff(1:dim[2],boot)

            X_permute::Array{Float64,3} = @views X[:,OOB,:]
            X_permute[:,:,variables[p]] = @views X[:,OOB[randperm(length(OOB))],variables[p]]

            ## Il faut maintenant regarder la différence en erreur de prédiction ::
            err_courante[i] = @views mean((Y[OOB].-pred_tree(tree, P,X_permute,X, dist)).^2)-mean((Y[OOB].-pred_tree(tree,P,X[:,OOB,:],X, dist)).^2)
        end
        imp[2,p]= mean(err_courante)
    end
    return imp
end

function Importance_Unique(frf::Array{Char},X::Array{Float64,3}, Y::Vector{Int64}, variables::UnitRange{Int64}, dist)
    dim=size(X)::Tuple{Int,Int,Int}
    imp::Array{Float64,2}= zeros(2,length(variables))

    trees=readdir(frf; join = true)
    ntree::Int=length(trees)

    err_courante::Vector{Float64}=zeros(ntree)
    imp[1,:] = variables

    @inbounds for p in 1:length(variables)
        err_courante.=0
        Threads.@threads for i in 1:ntree

            infos = load(trees[k])
            tree = infos["tree"]
            ID = infos["boot"]
            P = infos["P"]

            boot::Vector{Int64}= @view ID[i,findall3(x->x>0.0,ID[i,:])]
            ### Trouver les éléments qui ne sont pas dans l'échantillon OOB :
            OOB::Vector{Int64} = setdiff(1:dim[2],boot)
            X_permute::Array{Float64,3} = @views X[:,OOB,:]
            X_permute[:,:,variables[p]] = @views X[:,OOB[randperm(length(OOB))],variables[p]]
            ## Il faut maintenant regarder la différence en erreur de prédiction ::
            err_courante[i] = @views mean((Y[OOB].!=pred_tree(tree, P,X_permute,X, dist)))-mean((Y[OOB].!=pred_tree(tree,P,X[:,OOB,:],X, dist)))
        end
        imp[2,p]= mean(err_courante)
    end
    return imp
end

function ExtraFrechetRF(X::Array{Float64,3}, Y::Vector{Float64}, mtry::Int, ntree::Int, ntry::Int, imp::Bool,minElem::Int, dist)
    println("Building the Extra Fréchet Forest:")
    frf, P, boot  = ERFRF(X,Y, mtry, ntree, ntry,minElem, dist)
    println("OOB errors and % of explained variance:")
    pred_OOB, mse, varex = FRFERR(frf,X, Y,P, boot, dist)

    if imp==true
        println("Variables Importance Scores:")
        Imp = Importance(frf,X,Y,P, boot, dist)
        return Dict("FrechetRF"=> frf, "Imp"=>Imp, "varex"=> varex, "mse"=>mse, "OOB pred"=>pred_OOB, "P"=>P, "mtry"=>mtry, "ntry"=>ntry, "ntree"=>ntree, "dist"=>dist)
    end

    return Dict("FrechetRF"=> frf, "Imp"=>false, "varex"=> varex, "mse"=>mse, "OOB pred"=>pred_OOB, "P"=>P, "mtry"=>mtry, "ntry"=>ntry, "ntree"=>ntree, "dist"=>dist)
end

function ExtraFrechetRF(X::Array{Float64,3}, Y::Vector{Int64}, mtry::Int, ntree::Int, ntry::Int, imp::Bool,minElem::Int, dist)
    println("Building the Extra Fréchet Forest:")
    frf, P, boot  = ERFRF(X,Y, mtry, ntree, ntry,minElem, dist)
    println("OOB errors and % of explained variance:")
    pred_OOB, mse, varex = FRFERR(frf,X, Y,P, boot, dist)

    if imp==true
        println("Variables Importance Scores:")
        Imp = Importance(frf,X,Y,P, boot, dist)
        return Dict("FrechetRF"=> frf, "Imp"=>Imp, "varex"=> varex, "Class error"=>mse, "OOB pred"=>pred_OOB, "P"=>P, "mtry"=>mtry, "ntry"=>ntry, "ntree"=>ntree, "dist"=>dist)
    end

    return Dict("FrechetRF"=> frf, "Imp"=>false, "varex"=> varex, "Class error"=>mse, "OOB pred"=>pred_OOB, "P"=>P, "mtry"=>mtry, "ntry"=>ntry, "ntree"=>ntree, "dist"=>dist)
end

function prediction_rf(pred::Array{Float64,2})
    return mean(pred,dims=2)
end

function prediction_rf(pred::Array{Int64,2})

    sortie = fill(0, size(pred,2))
    for i in 1:size(pred,2)
        sortie[i] = findmax(countmap(pred[:,i]))[2]
    end
    return sortie
end

function pred_rf(frf::Dict{String, Any}, X::Array{Float64,3}, X_init::Array{Float64,3})
    ntree::Int = get(frf, "ntree",1)
    pred::Array{Float64,2}= zeros(size(X,2), ntree)

    forest = get(frf, "FrechetRF", 1)
    Pred = get(frf, "P", 1)
    dist = get(frf, "dist",1)

    p=Progress(ntree);
    update!(p,0)
    jj = Threads.Atomic{Int}(0)
    l = Threads.SpinLock()

    Threads.@threads for i in 1:ntree
        pred[:,i]= pred_tree(forest[:,:,i], Pred[:,:,i], X,X_init, dist)

        Threads.atomic_add!(jj, 1)
        Threads.lock(l)
        update!(p, jj[])
        Threads.unlock(l)
    end

    return prediction_rf(pred)
end

function OOB_unique(frf::Dict{String, Array}, X::Array{Float64,3}, indiv::Vector{Int64}, dist)

    dim=size(X)

    ntree = size(frf["trees"],3)
    type = eltype(frf["P"][1,1,1])

    pred_OOB=zeros(type,length(indiv))
    proba = zeros(length(indiv))
    ZZ = zeros(dim[1],2,dim[3])

    @inbounds for i in 1:length(indiv)
        Pred_courante = fill(-1,2,ntree)
        ZZ[:,1,:], ZZ[:,2,:] = X[:,indiv[i],:], X[:,indiv[i],:]
        Threads.@threads for k in 1:ntree

            if length(findall3(x->x==indiv[i],frf["id"][k,:])) == 0
                Pred_courante[1,k] = 0
                Pred_courante[2,k] = @views pred_tree(frf["trees"][:,:,k],frf["P"][:,:,k],ZZ,X, dist)[1]
            end

            ## Il faut maintenant regarder la différence en erreur de prédiction ::
        end

        if type==Int

            elem  = findmax(countmap(@views Pred_courante[2,findall3(x->x==0,Pred_courante[1,:])]))
            pred_OOB[i] = elem[2]
            proba[i] = elem[1]/length(findall3(x->x==0,Pred_courante[1,:]))
        else
            pred_OOB[i] = mean(@views Pred_courante[2,findall3(x->x==0,Pred_courante[1,:])])
        end
    end
    return pred_OOB, proba
end

function Combine(lien::String,X::Array{Float64,3})

    dim = size(X)

    Trees=readdir(lien, join = true)
    ntree=length(Trees)

    type = eltype(load(Trees[1])["P"][1,1])

    boot = zeros(ntree, dim[2])
    trees = zeros(4,2*dim[2]-1,ntree)
    P = zeros(type,(2,2*dim[2]-1,ntree))

    for k in 1:ntree

        infos = load(Trees[k])
        tree_temp = infos["tree"]
        id_temp = infos["boot"]
        P_temp = infos["P"]

        trees[:,1:size(tree_temp,2),k] = tree_temp
        P[:,1:size(P_temp,2),k]= P_temp
        boot[k,1:length(id_temp)]= id_temp

    end

    return Dict("trees"=>trees, "P"=>P, "id"=>boot)
end

function Combine_OOB(racine::String,Y)

    cd(string(racine,"/OOB"))
    files = readdir(pwd(), join=true)
    pred_OOB = []
    proba_OOB = []
    var_init = mean(Y.!=findmax(countmap(Y))[2])

    for i in 1:length(files)
        infos = load(string("OOB_",i,".jld"))
        pred_OOB = vcat(pred_OOB,infos["pred"])
        proba_OOB = vcat(proba_OOB, infos["proba"])
    end

    mse = mean(Y.!= pred_OOB)
    varex = 1 - mse/var_init

    return pred_OOB, proba_OOB, mse, varex
end


text:
module Foo
import Base.show
export MyType, foo

struct MyType
  x
end

bar(x) = 2x
foo(a::MyType) = bar(a.x) + 1
show(io::IO, a::MyType) = print(io, "MyType $(a.x)")

end

ext: jl


text:
################################################################################
#
#  Contained in the alternating group
#
################################################################################

function is_discriminant_square(IdG::GAP.GapObj)
  G = GAP.Globals.SmallGroup(IdG)
  mp = GAP.Globals.RegularActionHomomorphism(G)
  S = GAP.Globals.ImagesSource(mp)
  lg = GAP.Globals.GeneratorsOfGroup(S)
  for i = 1:length(lg)
    s = GAP.Globals.SignPerm(lg[i])
    if !isone(s)
      return false
    end
  end
  return true
end

###############################################################################
#
#  Conductors
#
###############################################################################

function _conductors_using_cocycles(F::FieldsTower, st::Vector{Int}, l_cond::Vector, E)
  lp = ramified_primes(F)
  auts = automorphisms(F.field, copy = false)
  projections = F.projections_for_conductors
  G = GAP.Globals.ImagesSource(projections[1])
  E = GAP.Globals.Source(projections[1])
  D = F.isomorphism
  n = prod(st)
  O = maximal_order(F)
  for p in lp
    lP = prime_decomposition(O, p)
    Hp = inertia_subgroup(lP[1][1])
    gHp = small_generating_set(Hp)
    els = [D[g] for g in gHp]
    sub = GAP.Globals.Subgroup(G, GAP.julia_to_gap(els))
    ord = GAP.Globals.Size(sub)
    sizes_preimages = Int[]
    for proj in projections
      subgs = Vector{Vector{GAP.GapObj}}(undef, n)
      for i = 1:n
        subgs[i] = Vector{GAP.GapObj}(undef, length(els))
      end
      for j = 1:length(els)
        pels = GAP.Globals.List(GAP.Globals.PreImages(proj, els[j]))
        for i = 1:length(pels)
          subgs[i][j] = pels[i]
        end
      end
      for lelem in subgs
        sub = GAP.Globals.Subgroup(E, GAP.julia_to_gap(lelem))
        onew = GAP.Globals.Size(sub)
        push!(sizes_preimages, onew)
      end
    end
    if minimum(sizes_preimages) != ord
      #The prime must ramify!
      l1 = Vector{Tuple{Int, Dict{NfOrdIdl, Int}}}()
      if !divisible(fmpz(n), p)
        for x in l_cond
          if divisible(fmpz(x[1]), p)
            push!(l1, x)
          end
        end
      else
        for x in l_cond
          if !isempty(x[2])
            lI = keys(x[2])
            found = false
            for k in lI
              if minimum(k) == p
                found = true
                break
              end
            end
            if found
              push!(l1, x)
            end
          end
        end
      end
      l_cond = l1
    elseif maximum(sizes_preimages) == ord && !divisible(fmpz(n), p)
      #The prime must be unramified!
      l1 = Vector{Tuple{Int, Dict{NfOrdIdl, Int}}}()
      if !divisible(fmpz(n), p)
        for x in l_cond
          if !divisible(fmpz(x[1]), p)
            push!(l1, x)
          end
        end
      else
        for x in l_cond
          if !isempty(x[2])
            lI = keys(x[2])
            found = false
            for k in lI
              if minimum(k) == p
                found = true
                break
              end
            end
            if !found
              push!(l1, x)
            end
          end
        end
      end
      l_cond = l1
    end
  end
  return l_cond

end

function conductors_with_restrictions(F::FieldsTower, st::Vector{Int}, IdG::GAP.GapObj, bound::fmpz)

  O = maximal_order(F)
  l_cond = Hecke.conductors(O, st, bound)
  G = GAP.Globals.SmallGroup(IdG)
  new_conds = _conductors_using_cocycles(F, st, l_cond, G)
  if length(st) != 1 || !isprime(st[1]) || isempty(new_conds)
    return new_conds
  end
  #If the extension is cyclic, I take care of the discriminant being a square or not for the wild ramification
  issquare = is_discriminant_square(IdG)
  p = st[1]
  v = valuation(discriminant(O), p)
  is_square_disc_base_field = iszero(mod(v*p, 2))
  td = prime_decomposition_type(O, p)
  if iseven(length(td)) || iseven(td[1][1]) || isodd(p)
    #Regardless of the exponents, the norm of the discriminant will be a square
    if issquare && is_square_disc_base_field
      return new_conds
    elseif issquare
      return typeof(new_conds)()
    else
      return new_conds
    end
  end
  #Now, p must be 2.
  if issquare && is_square_disc_base_field
    #Only the even exponents are allowed!
    newer_conds = typeof(new_conds)()
    for i = 1:length(new_conds)
      if isempty(new_conds[i][2])
        push!(newer_conds, new_conds[i])
        continue
      end
      if iszero(mod(first(values(new_conds[i][2])), 2))
        push!(newer_conds, new_conds[i])
      end
    end
  elseif issquare
    #Only the odd exponents are allowed!
    newer_conds = typeof(new_conds)()
    for i = 1:length(new_conds)
      if !isempty(new_conds[i][2]) && !iszero(mod(first(values(new_conds[i][2])), 2))
        push!(newer_conds, new_conds[i])
      end
    end
  else
    newer_conds = new_conds
  end

  #Now, tame ramification.
  list_tame = Int[x[1] for x in newer_conds]
  list_tame = coprime_base(list_tame)
  l = length(list_tame)
  for i = 1:length(list_tame)
    x = list_tame[i]
    if !isone(x) && !isprime(x)
      append!(list_tame, Hecke.divisors(x))
    end
  end
  list_tame = coprime_base(list_tame)
  for q in list_tame
    q == 1 && continue
    #@assert isprime(q)
    v = valuation(discriminant(O), q)
    is_square_disc_base_field = iszero(mod(v*p, 2))
    td = prime_decomposition_type(O, q)
    if iszero(mod(length(td) * td[1][1] * (p-1), 2))
      #Regardless of the exponents, the norm of the discriminant will be a square
      if issquare && is_square_disc_base_field
        continue
      elseif issquare || is_square_disc_base_field
        return typeof(new_conds)()
      else
        continue
      end
    end
    #Now, p must be 2.
    if issquare && is_square_disc_base_field
      #Only the even exponents are allowed!
      #Therefore the prime can't ramify
      newest_conds = typeof(new_conds)()
      for i = 1:length(newer_conds)
        if !iszero(mod(newer_conds[i][1], q))
          push!(newer_conds, newer_conds[i])
        end
      end
    elseif issquare
      #Only the odd exponents are allowed!
      #Therefore the prime must ramify
      newest_conds = typeof(new_conds)()
      for i = 1:length(newer_conds)
        if iszero(mod(newer_conds[i][1], q))
          push!(newest_conds, newer_conds[i])
        end
      end
    else
      newest_conds = newer_conds
    end
    newer_conds = newest_conds
  end
  return newer_conds
end


text:
using Documenter, DiffResults

makedocs(
    modules=[DiffResults],
    doctest = false,
    sitename = "DiffResults",
    pages = ["Documentation" => "index.md"])

deploydocs(
    repo = "github.com/JuliaDiff/DiffResults.jl.git")


text:
"""
    Bernoulli(p)

A *Bernoulli distribution* is parameterized by a success rate `p`, which takes value 1
with probability `p` and 0 with probability `1-p`.

```math
P(X = k) = \\begin{cases}
1 - p & \\quad \\text{for } k = 0, \\\\
p & \\quad \\text{for } k = 1.
\\end{cases}
```

```julia
Bernoulli()    # Bernoulli distribution with p = 0.5
Bernoulli(p)   # Bernoulli distribution with success rate p

params(d)      # Get the parameters, i.e. (p,)
succprob(d)    # Get the success rate, i.e. p
failprob(d)    # Get the failure rate, i.e. 1 - p
```

External links:

* [Bernoulli distribution on Wikipedia](http://en.wikipedia.org/wiki/Bernoulli_distribution)
"""
struct Bernoulli{T<:Real} <: DiscreteUnivariateDistribution
    p::T

    function Bernoulli{T}(p::T) where T
        @check_args(Bernoulli, zero(p) <= p <= one(p))
        new{T}(p)
    end

end

Bernoulli(p::T) where {T<:Real} = Bernoulli{T}(p)
Bernoulli(p::Integer) = Bernoulli(Float64(p))
Bernoulli() = Bernoulli(0.5)

@distr_support Bernoulli 0 1

#### Conversions
convert(::Type{Bernoulli{T}}, p::Real) where {T<:Real} = Bernoulli(T(p))
convert(::Type{Bernoulli{T}}, d::Bernoulli{S}) where {T <: Real, S <: Real} = Bernoulli(T(d.p))

#### Parameters

succprob(d::Bernoulli) = d.p
failprob(d::Bernoulli) = 1 - d.p

params(d::Bernoulli) = (d.p,)
@inline partype(d::Bernoulli{T}) where {T<:Real} = T

#### Properties

mean(d::Bernoulli) = succprob(d)
var(d::Bernoulli) =  succprob(d) * failprob(d)
skewness(d::Bernoulli) = (p0 = failprob(d); p1 = succprob(d); (p0 - p1) / sqrt(p0 * p1))
kurtosis(d::Bernoulli) = 1 / var(d) - 6

mode(d::Bernoulli) = ifelse(succprob(d) > 1/2, 1, 0)

function modes(d::Bernoulli)
    p = succprob(d)
    p < 1/2 ? [0] :
    p > 1/2 ? [1] : [0, 1]
end

median(d::Bernoulli) = ifelse(succprob(d) <= 1/2, 0, 1)

function entropy(d::Bernoulli)
    p0 = failprob(d)
    p1 = succprob(d)
    (p0 == 0 || p0 == 1) ? zero(d.p) : -(p0 * log(p0) + p1 * log(p1))
end

#### Evaluation

pdf(d::Bernoulli, x::Bool) = x ? succprob(d) : failprob(d)
pdf(d::Bernoulli, x::Int) = x == 0 ? failprob(d) :
                            x == 1 ? succprob(d) : zero(d.p)

cdf(d::Bernoulli, x::Bool) = x ? failprob(d) : one(d.p)
cdf(d::Bernoulli, x::Int) = x < 0 ? zero(d.p) :
                            x < 1 ? failprob(d) : one(d.p)

ccdf(d::Bernoulli, x::Bool) = x ? succprob(d) : one(d.p)
ccdf(d::Bernoulli, x::Int) = x < 0 ? one(d.p) :
                             x < 1 ? succprob(d) : zero(d.p)

function quantile(d::Bernoulli{T}, p::Real) where T<:Real
    0 <= p <= 1 ? (p <= failprob(d) ? zero(T) : one(T)) : T(NaN)
end
function cquantile(d::Bernoulli{T}, p::Real) where T<:Real
    0 <= p <= 1 ? (p >= succprob(d) ? zero(T) : one(T)) : T(NaN)
end

mgf(d::Bernoulli, t::Real) = failprob(d) + succprob(d) * exp(t)
cf(d::Bernoulli, t::Real) = failprob(d) + succprob(d) * cis(t)

#### Sampling

rand(rng::AbstractRNG, d::Bernoulli) = rand(rng) <= succprob(d)

#### MLE fitting

struct BernoulliStats <: SufficientStats
    cnt0::Float64
    cnt1::Float64

    BernoulliStats(c0::Real, c1::Real) = new(Float64(c0), Float64(c1))
end

fit_mle(::Type{Bernoulli}, ss::BernoulliStats) = Bernoulli(ss.cnt1 / (ss.cnt0 + ss.cnt1))

function suffstats(::Type{Bernoulli}, x::AbstractArray{T}) where T<:Integer
    n = length(x)
    c0 = c1 = 0
    for i = 1:n
        @inbounds xi = x[i]
        if xi == 0
            c0 += 1
        elseif xi == 1
            c1 += 1
        else
            throw(DomainError())
        end
    end
    BernoulliStats(c0, c1)
end

function suffstats(::Type{Bernoulli}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    n = length(x)
    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions."))
    c0 = c1 = 0
    for i = 1:n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        if xi == 0
            c0 += wi
        elseif xi == 1
            c1 += wi
        else
            throw(DomainError())
        end
    end
    BernoulliStats(c0, c1)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 10058780-939b-11eb-17b4-275d951b8958
using ImageDraw,Images,ImageCore,ImageView

# ╔═╡ 49c72e40-93aa-11eb-273f-33087023fe1e
begin
abstract type AbstractPolyFillAlgorithm end
	function draw(img::AbstractArray{T,2},verts::Vector{CartesianIndex{2}},f::AbstractPolyFillAlgorithm) where T <: Colorant
    f(res,verts,f.x,f.y,f.fill_color,f.boundarycolor)
end
function draw(img::AbstractArray{T,2},
    verts::Vector{CartesianIndex{2}},
    f::AbstractPolyFillAlgorithm) where {T<:Colorant}
   res = copy(img)
   f(res,verts,f.x,f.y,f.fill_color,f.boundarycolor)
end
end

# ╔═╡ 7413c77e-939a-11eb-0988-5503daa6d125
# begin
# """
#     AbstractPolyFillAlgorithm

# The root of polygon filling algorithms type system
# """

# abstract type AbstractPolyFillAlgorithm end

# struct boundaryfill{T<:Colorant} <: AbstractPolyFillAlgorithm
#     boundarycolor::T
#     # function boundaryfill(boundarycolor::T) where {T<:Colorant}
#     #     boundarycolor < 0 && throw(ArgumentError("window_size should be non-negative."))
#     #     println("Test 1")
#     #     new{T}(boundarycolor)
#     # end
# end

# function boundaryfill(;boundarycolor::Colorant=RGB(1))
#     println("Object Creation")
#     boundaryfill(boundarycolor)
# end
# function boundaryfill(x::Int, y::Int, fill_color::T, boundary_color::T) where T<:Colorant
# 	println("Enter Boundary Fill")
#     if (res[y, x] != boundary_color && res[y, x] !=fill_color)
#             if checkbounds(Bool, res, y, x) res[y, x] = fill_color end
#              boundaryfill(x + 1, y, fill_color, boundary_color)
#              boundaryfill(x, y + 1, fill_color, boundary_color)
#              boundaryfill(x - 1, y, fill_color, boundary_color)
#              boundaryfill(x, y - 1, fill_color, boundary_color)
#     end
#     img
# end

# function (f::boundaryfill)(img::AbstractArray{T,2},verts::Vector) where {T<:Colorant}
#     println("Enter boundary fill API")
# 	y=4
# 	x=4
# 	# println(verts)
# 	# println(img)
# 	fill_color = RGB(1)
# 	boundary_color = RGB(1)
# 	for i in 1:length(verts)-1
# 		draw!(res, LineSegment(verts[i], verts[i+1]), fill_color)
# 		imshow(res)
#     end

# 	boundaryfill(x::Int, y::Int, fill_color::T, boundary_color::T)
#    res
# end
# function drawnewtype(img::AbstractArray{T,2},verts::Vector{CartesianIndex{2}},f::AbstractPolyFillAlgorithm) where T <: Colorant
#     println("Test 2")
#     f(img,verts);
# end
# function drawnewtype(img::AbstractArray{T,2},
#     verts::Vector{CartesianIndex{2}},
#     f::AbstractPolyFillAlgorithm) where {T<:Colorant}
#    println("Start")
#    f(res,verts)
#    # println(verts)
#    res
# end
# end

# ╔═╡ e60baf34-93a9-11eb-00ac-09ab117e4923
begin

img = zeros(RGB, 7, 7)
expected = copy(img)
expected[2:6, 2] .= RGB{N0f8}(1)
expected[2:6, 6] .= RGB{N0f8}(1)
expected[2, 2:6] .= RGB{N0f8}(1)
expected[6, 2:6] .= RGB{N0f8}(1)
verts = [CartesianIndex(2, 2), CartesianIndex(2, 6),CartesianIndex(6, 6), CartesianIndex(6, 2),CartesianIndex(2,2)]

struct boundaryfill{T<:Colorant} <: AbstractPolyFillAlgorithm
	x::Int
	y::Int
	fill_color::T
    boundarycolor::T
    function boundaryfill(x::Int,y::Int,fill_color::T,boundarycolor::T) where {T<:Colorant}
        println("Test 1")
        new{T}(x,y,fill_color,boundarycolor)
    end
end

function boundaryfill(;x::Int=0,y::Int=0,fill_color::Colorant=RGB(1),boundarycolor::Colorant=RGB(1))
		boundaryfill(x,y,fill_color,boundarycolor)
end
function boundaryfill(res::AbstractArray{T,2},x::Int, y::Int, fill_color::T, boundary_color::T) where T<:Colorant
    if (res[y, x] != boundary_color && res[y, x] !=fill_color)
            if checkbounds(Bool, res, y, x) res[y, x] = fill_color end
             boundaryfill(res,x + 1, y, fill_color, boundary_color)
             boundaryfill(res,x, y + 1, fill_color, boundary_color)
             boundaryfill(res,x - 1, y, fill_color, boundary_color)
             boundaryfill(res,x, y - 1, fill_color, boundary_color)
    end
end

function (f::boundaryfill)(res::AbstractArray{T,2},verts::Vector,x::Int,y::Int,fill_color::T,boundarycolor::T) where {T<:Colorant}
	for i in 1:length(verts)-1
		draw!(res, LineSegment(verts[i], verts[i+1]), fill_color)
    end
	boundaryfill(res,f.x, f.y, f.fill_color, f.boundarycolor)
   res
end
end

# ╔═╡ 15b0f17c-93a5-11eb-17af-57808dba7004
img

# ╔═╡ a271c9c4-939b-11eb-10b1-2140a3edc2c7
draw(img,verts,boundaryfill(x=4,y=4,fill_color=RGB(1),boundarycolor=RGB(1)))

# ╔═╡ dc8384a4-942c-11eb-1590-3177edbb685e
# begin
# struct floodfill{T<:Colorant} <: AbstractPolyFillAlgorithm
# 	x::Int
# 	y::Int
# 	fill_color::T
#     boundarycolor::T
#     # function boundaryfill(boundarycolor::T) where {T<:Colorant}
#     #     boundarycolor < 0 && throw(ArgumentError("window_size should be non-negative."))
#     #     println("Test 1")
#     #     new{T}(boundarycolor)
#     # end
# end

# function floodfill(;boundarycolor::Colorant=RGB(1))
#     println("Object Creation")
#     floodfill(x,y,fill_color,boundarycolor)
# end
# # function floodfill(x::Int, y::Int, fill_color::T, boundary_color::T) where T<:Colorant
# # 	println("Enter Boundary Fill")
# #     if (res[y, x] != boundary_color && res[y, x] !=fill_color)
# #             if checkbounds(Bool, res, y, x) res[y, x] = fill_color end
# #              boundaryfill(x + 1, y, fill_color, boundary_color)
# #              boundaryfill(x, y + 1, fill_color, boundary_color)
# #              boundaryfill(x - 1, y, fill_color, boundary_color)
# #              boundaryfill(x, y - 1, fill_color, boundary_color)
# #     end
# #     img
# # end
# function floodfillcolor(res, x, y, current_color, fill_color)
#            if (res[y, x] != current_color || res[y, x] == fill_color) return end
#            if checkbounds(Bool, res, y, x) res[y, x] = fill_color end
#            floodfillcolor(res , x + 1, y, current_color, fill_color)
#            floodfillcolor(res , x - 1, y, current_color, fill_color)
#            floodfillcolor(res , x, y + 1, current_color, fill_color)
#            floodfillcolor(res , x, y - 1, current_color, fill_color)
# end

# function floodfill(res, x, y, fill_color)
#            current_color = res[y,x]
# 		   current_color= RGB(0)
#            floodfillcolor(res, x, y, current_color, fill_color)
# end

# function (f::floodfill)(res::AbstractArray{T,2},verts::Vector) where {T<:Colorant}
#     println("Enter boundary fill API")
# 	y=4
# 	x=4
# 	# println(verts)
# 	# println(img)
# 	imshow(res)
# 	fill_color = RGB(1)
# 	boundary_color = RGB(1)
# 	for i in 1:length(verts)-1
# 		draw!(res, LineSegment(verts[i], verts[i+1]), fill_color)
# 		imshow(res)
#     end
# 	floodfill(x::Int, y::Int, fill_color::T, boundary_color::T)
#    res
# end
# end

# ╔═╡ e6df086e-93a1-11eb-147a-a7d188f4059d
expected

# ╔═╡ 05e7f2a6-93da-11eb-3a50-31f63d244130
# begin

# vert=CartesianIndex{2}[]
# 	push!(vert, CartesianIndex(2,2))
# 	push!(vert, CartesianIndex(3,2))
#   push!(vert, CartesianIndex(4,2))
# 	push!(vert, CartesianIndex(5,2))
#   push!(vert, CartesianIndex(6,3))
# 	push!(vert, CartesianIndex(7,4))
# 	push!(vert, CartesianIndex(8,4))
# 	push!(vert, CartesianIndex(9,3))
# 	push!(vert, CartesianIndex(10,2))
# 	push!(vert, CartesianIndex(11,3))
# 	push!(vert, CartesianIndex(12,4))
# 	push!(vert, CartesianIndex(13,5))
# 	push!(vert, CartesianIndex(12,6))
# 	push!(vert, CartesianIndex(11,7))
# 	push!(vert, CartesianIndex(10,6))#u
# 	push!(vert, CartesianIndex(9,7))
# 	push!(vert, CartesianIndex(8,8))
# 	push!(vert, CartesianIndex(8,9))
# 	push!(vert, CartesianIndex(8,10))
# 	push!(vert, CartesianIndex(8,11))
# 	push!(vert, CartesianIndex(9,11))
# 	push!(vert, CartesianIndex(10,11))
# 	push!(vert, CartesianIndex(11,11))
# 	push!(vert, CartesianIndex(12,11))
# 	push!(vert, CartesianIndex(13,11))
# 	push!(vert, CartesianIndex(13,12))
# 	push!(vert, CartesianIndex(13,13))
# 	push!(vert, CartesianIndex(12,13))
# 	push!(vert, CartesianIndex(11,13))
# 	push!(vert, CartesianIndex(10,13))
# 	push!(vert, CartesianIndex(9,13))
# 	push!(vert, CartesianIndex(8,13))
# 	push!(vert, CartesianIndex(7,13))
# 	push!(vert, CartesianIndex(6,13))
# 	push!(vert, CartesianIndex(6,12))
# 	push!(vert, CartesianIndex(6,11))#u
# 	push!(vert, CartesianIndex(5,11))
# 	push!(vert, CartesianIndex(4,11))
# 	push!(vert, CartesianIndex(3,11))
# 	push!(vert, CartesianIndex(2,11))#u
# 	push!(vert, CartesianIndex(3,10))
# 	push!(vert, CartesianIndex(4,9))
# 	push!(vert, CartesianIndex(5,8))
# 	push!(vert, CartesianIndex(4,7))
# 	push!(vert, CartesianIndex(4,6))
# 	push!(vert, CartesianIndex(5,5))
# 	push!(vert, CartesianIndex(4,4))
#   push!(vert, CartesianIndex(3,3))
# 	push!(vert, CartesianIndex(2,2))

# end

# ╔═╡ a7810aa4-93d9-11eb-3dce-955d8f4e15b8
# begin
# struct edgetabletuple
# 	initial::CartesianIndex
# 	final::CartesianIndex
# end

# img = draw(zeros(Gray{Bool},14,14), Polygon(vert));
# function createedgetable()
# edgetable= []
# push!(edgetable, edgetabletuple(CartesianIndex(2,2),CartesianIndex(5,2)))
# push!(edgetable, edgetabletuple(CartesianIndex(5,2),CartesianIndex(7,4)))
# push!(edgetable, edgetabletuple(CartesianIndex(7,4),CartesianIndex(8,4)))
# push!(edgetable, edgetabletuple(CartesianIndex(8,4),CartesianIndex(10,2)))
# push!(edgetable, edgetabletuple(CartesianIndex(10,2),CartesianIndex(13,5)))
# push!(edgetable, edgetabletuple(CartesianIndex(13,5),CartesianIndex(11,7)))
# push!(edgetable, edgetabletuple(CartesianIndex(11,7),CartesianIndex(10,6)))
# push!(edgetable, edgetabletuple(CartesianIndex(10,6),CartesianIndex(8,8)))
# push!(edgetable, edgetabletuple(CartesianIndex(8,8),CartesianIndex(8,11)))
# push!(edgetable, edgetabletuple(CartesianIndex(8,11),CartesianIndex(13,11)))
# push!(edgetable, edgetabletuple(CartesianIndex(13,11),CartesianIndex(13,13)))
# push!(edgetable, edgetabletuple(CartesianIndex(13,13),CartesianIndex(6,13)))
# push!(edgetable, edgetabletuple(CartesianIndex(6,13),CartesianIndex(6,11)))
# push!(edgetable, edgetabletuple(CartesianIndex(6,11),CartesianIndex(2,11)))
# push!(edgetable, edgetabletuple(CartesianIndex(2,11),CartesianIndex(5,8)))
# push!(edgetable, edgetabletuple(CartesianIndex(5,8),CartesianIndex(4,8)))
# push!(edgetable, edgetabletuple(CartesianIndex(4,8),CartesianIndex(4,7)))
# push!(edgetable, edgetabletuple(CartesianIndex(4,7),CartesianIndex(4,6)))
# push!(edgetable, edgetabletuple(CartesianIndex(4,6),CartesianIndex(5,5)))
# push!(edgetable, edgetabletuple(CartesianIndex(5,5),CartesianIndex(2,2)))
# return edgetable
# end
# #find ymin and ymax
# function yminmax(vert)
# 	ymax = vert[1][1];
# 	ymin = vert[1][1];
# 	for i in 1:length(vert)
# 	   if(vert[i][2] > ymax)
# 		 ymax = vert[i][2]
# 	   elseif(vert[i][2] < ymin)
# 		 ymax = vert[i][2]
# 	   end
# 	end
# 	return ymin, ymax
# end

# function findintersections(edgetable, yvalue)
# 	points =[]
#     for i in 1:length(edgetable)
# 		if(edgetable[i].final[2] > yvalue && edgetable[i].initial[2] > yvalue)
# 			continue
# 		end
# 		if(edgetable[i].final[2] <= yvalue && edgetable[i].initial[2] <= yvalue)
# 			continue
# 		end
# 		x=1
# 		deltay = edgetable[i].final[1]-edgetable[i].initial[1]
# 		deltax = edgetable[i].final[2]-edgetable[i].initial[2]
# 		alpha = deltay/deltax
# 		constant = edgetable[i].initial[1]
# 		x = alpha * (yvalue - edgetable[i].initial[2]) + constant
# 		if (x == -Inf || isnan(x) || x == Inf) continue end
# 		push!(points,CartesianIndex(ceil(Int,x),yvalue))
# 	end
# 	return points
# end
# function timetocolor(points, fill_color::T)where T<:Colorant
#    for i in 1:2:length(points)-1
# 	  draw!(img, LineSegment(points[i], points[i+1]),fill_color)
#    end
# end
# function scanline(vert, fill_color::T) where T<:Colorant
# 	ymin, ymax = yminmax(vert)
# 	edgetable = createedgetable()
# 	pyramid = []
# 	for i in ymin:ymax
# 		points = findintersections(edgetable, i)
# 		points = sort!(points, by = x -> x[1])
# 		timetocolor(points, fill_color)
# 		push!(pyramid,img[:,:])
# 	end
# 	pyramid
# end
# fill_color = Gray{Bool}(1.0)

# end

# ╔═╡ e747c024-93d9-11eb-2c75-a92984725097
# pyramid = scanline(vert, fill_color)

# ╔═╡ Cell order:
# ╠═10058780-939b-11eb-17b4-275d951b8958
# ╠═49c72e40-93aa-11eb-273f-33087023fe1e
# ╟─7413c77e-939a-11eb-0988-5503daa6d125
# ╠═e60baf34-93a9-11eb-00ac-09ab117e4923
# ╠═15b0f17c-93a5-11eb-17af-57808dba7004
# ╠═a271c9c4-939b-11eb-10b1-2140a3edc2c7
# ╠═dc8384a4-942c-11eb-1590-3177edbb685e
# ╠═e6df086e-93a1-11eb-147a-a7d188f4059d
# ╠═05e7f2a6-93da-11eb-3a50-31f63d244130
# ╠═a7810aa4-93d9-11eb-3dce-955d8f4e15b8
# ╠═e747c024-93d9-11eb-2c75-a92984725097

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
@doc doc"""
    Nvim([address,] port)
    Nvim(path)

Represents a particular NeoVim instance. Can be instantiated with any method
signature which would return a `TCPSocket` or `Pipe` when passed to `connect`.
""" ->
type Nvim{T}
    recstream::T
    sendstream::T

    reqid::Int
end

function Nvim{T <: Base.AsyncStream}(ins::T, outs::T)
    n = Nvim(ins, outs, 1)
    finalizer(n, close)
    n
end

function Nvim(::Type{Val{:Embedded}})
    ins, outs = STDIN, STDOUT

    # TODO: use proper logging (Logging.jl, LumberJack.jl) and a more correct file path
    logfile = open(joinpath(homedir(), ".neovim.jl.log"), "a")
    redirect_stdout(logfile)
    redirect_stderr(logfile)
    redirect_stdin()
    # use now() when 0.4 is ubiquitous
    #println("Log starting $(now())\n===")
    println("Starting Log\n============")

    Nvim(ins, outs)
end

function Nvim(::Type{Val{:Spawn}})
    # TODO: have a way to allow users to specify program path
    ins, outs, proc = readandwrite(`nvim --embed`)
    Nvim(ins, outs), proc
end

# TODO: would it be better to split this up?
function Nvim(args...)
    conn = connect(args...)
    Nvim(conn, conn)
end

function Base.close(n::Nvim)
    close(n.recstream)
    close(n.sendstream)
end

function getindex(n::Nvim, name::Symbol)
    if name === :buffers
        return vim_get_buffers(n)
    elseif name === :windows
        return vim_get_windows(n)
    elseif name === :tabpages
        return vim_get_tabpages(n)
    else
        throw(KeyError(name))
    end
end

MsgPack.pack(n::Nvim, v) = MsgPack.pack(n.sendstream, v)
MsgPack.unpack(n::Nvim) = MsgPack.unpack(n.recstream)


text:
# activation functions must be applied to vectors
#  and return a vector, as well as their derivatives

export tanhprime, softmax, softmaxprime, σ

#============================#
#            TANH            #
#============================#
function tanh(l::FFNNLayer)
    if l.bias
        return vcat([1.0], map(tanh, l[1:end]))
    else
        return map(tanh, l[1:end])
    end
end

function tanhprime(l::FFNNLayer)
    diagm(map(xi -> sech(xi)^2, l.neurons))
end

#============================#
#          SOFTMAX           #
#============================#
function softmax(l::FFNNLayer)
    Z = sum(exp(l.neurons)) # Normalization
    return [exp(xi)/Z for xi in l.neurons]
end

function softmax(l::FFNNLayer, i::Int)
    Z = sum(exp(l.neurons)) # Normalization
    return exp(l.neurons[i])/Z
end

function softmaxprime(l::FFNNLayer)
    delta(i,j) = (i==j) ? 1 : 0
    return Float64[softmax(l, i)*(delta(i,j) - softmax(l,j)) for i in 1:length(l.neurons), j in 1:length(l.neurons)]
end


text:
using Knet

# Load MNIST

if !isdefined(:MNIST)
    include(Pkg.dir("Knet/examples/mnist.jl"))
end

# Prepare data

batchsize = 100
dtrn = Any[]
for i=1:batchsize:size(MNIST.xtrn,2)
    j=i+batchsize-1
    push!(dtrn, (MNIST.xtrn[:,i:j],MNIST.ytrn[:,i:j]))
end
dtst = Any[]
for i=1:batchsize:size(MNIST.xtst,2)
    j=i+batchsize-1
    push!(dtst, (MNIST.xtst[:,i:j],MNIST.ytst[:,i:j]))
end

# Prepare model

@knet function mnistsimple(x1)
    w1 = par(;dims=(64,0))
    z1 = dot(w1,x1)
    b1 = par(;dims=(0,))
    y1 = add(b1,z1)
    x2 = relu(y1)
    w2 = par(;dims=(10,0))
    z2 = dot(w2,x2)
    b2 = par(;dims=(0,))
    y2 = add(b2,z2)
    ou = soft(y2)
end

model = FNN(mnistsimple)

# Train and evaluate

setopt!(model, lr=0.5)
for epoch=1:10
    (l,w,g) = train(model, dtrn, softloss)
    (l1,a1) = (test(model, dtrn, softloss), 1-test(model, dtrn, zeroone))
    (l2,a2) = (test(model, dtst, softloss), 1-test(model, dtst, zeroone))
    println("epoch=$epoch tstacc=$a2 trnacc=$a1 tstloss=$l2 trnloss=$l1 wnorm=$w gnorm=$g")
end

### SAMPLE RUN

# epoch=1 tstacc=0.9362 trnacc=0.9387500000000001 tstloss=0.19834305 trnloss=0.19569032 wnorm=18.56881 gnorm=3.0631502
# epoch=2 tstacc=0.9589000000000001 trnacc=0.96365 tstloss=0.1326303 trnloss=0.12078557 wnorm=22.30527 gnorm=4.1069307
# epoch=3 tstacc=0.9663 trnacc=0.9712666666666667 tstloss=0.11089823 trnloss=0.09421845 wnorm=24.958153 gnorm=3.7882617
# epoch=4 tstacc=0.9686 trnacc=0.9760666666666669 tstloss=0.100210935 trnloss=0.078416765 wnorm=27.135939 gnorm=3.2001014
# epoch=5 tstacc=0.972 trnacc=0.9791333333333335 tstloss=0.092702836 trnloss=0.06718973 wnorm=28.992907 gnorm=2.4996655
# epoch=6 tstacc=0.9716 trnacc=0.9807333333333335 tstloss=0.09260971 trnloss=0.06175735 wnorm=30.643934 gnorm=2.2204309
# epoch=7 tstacc=0.9716 trnacc=0.9817333333333335 tstloss=0.09355682 trnloss=0.056200992 wnorm=32.127914 gnorm=2.1083603
# epoch=8 tstacc=0.9729 trnacc=0.9836000000000001 tstloss=0.09248012 trnloss=0.049878325 wnorm=33.49334 gnorm=2.1731966
# epoch=9 tstacc=0.9719 trnacc=0.9832333333333334 tstloss=0.09715143 trnloss=0.04935321 wnorm=34.74469 gnorm=1.9201015
# epoch=10 tstacc=0.9732000000000001 trnacc=0.9856166666666668 tstloss=0.09411778 trnloss=0.04231144 wnorm=35.89667 gnorm=1.8556117


text:
# test solver initialization
let
    DefaultAStarSC()
    AStarSC{Int}()
    ISPS()
    CBSSolver()
    TaskGraphsMILPSolver()
    NBSSolver()
    solver = NBSSolver()
    best_cost(solver)
    NBSSolver(TaskGraphsMILPSolver(),CBSSolver())
end
let
    solver = NBSSolver()
    increment_iteration_count!(solver)
    increment_iteration_count!(route_planner(solver))
    @test max_iterations(solver) == 1
    @test iterations(solver) == 1
    reset_solver!(solver)
    @test max_iterations(solver) == 1
    @test iterations(solver) == 0
    hard_reset_solver!(solver)
    @test max_iterations(solver) == 0
    @test max_iterations(route_planner(solver)) == 0
    @test iterations(solver) == 0
end
# in-depth testing of solver stages on a single problem
let
    # init search env
    f = pctapf_problem_4
    cost_model = SumOfMakeSpans()
    solver = NBSSolver(assignment_model = TaskGraphsMILPSolver(GreedyAssignment()))
    pc_tapf = f(solver;cost_function=cost_model,verbose=false)
    base_env = pc_tapf.env
    prob = formulate_assignment_problem(solver.assignment_model,pc_tapf)
    sched, cost = solve_assignment_problem!(solver.assignment_model,prob,pc_tapf)
    @test validate(sched)
    # Test AStarSC
    let
        search_env = construct_search_env(solver,deepcopy(sched),base_env)
        pc_mapf = PC_MAPF(search_env)
        path_planner = DefaultAStarSC()
        node = initialize_root_node(search_env)
        # plan for Robot Start
        node_id = RobotID(1)
        schedule_node = get_node_from_id(search_env.schedule, node_id)
        v = get_vtx(search_env.schedule, node_id)
        @test plan_path!(path_planner, pc_mapf,search_env, node, schedule_node, v)
        # plan for GO
        v2 = outneighbors(search_env.schedule,v)[1]
        schedule_node = get_node_from_vtx(search_env.schedule, v2)
        @test plan_path!(path_planner, pc_mapf,search_env, node, schedule_node, v2)
        # @show convert_to_vertex_lists(node.solution)
    end
    let
        search_env = construct_search_env(solver,deepcopy(sched),base_env)
        pc_mapf = PC_MAPF(search_env)
        node = initialize_root_node(search_env)
        path_planner = ISPS()
        for i in 1:nv(search_env.schedule)
            @test plan_next_path!(path_planner,pc_mapf,node.solution,node)
        end
    end
    # Test ISPS
    let
        search_env = construct_search_env(solver,deepcopy(sched),base_env)
        path_planner = ISPS()
        node = initialize_root_node(search_env)
        @test compute_route_plan!(path_planner,PC_MAPF(search_env),node)#,node)
        # @show convert_to_vertex_lists(node.solution)
        # @show convert_to_vertex_lists(get_route_plan(search_env))
    end
    let
        search_env = construct_search_env(solver,deepcopy(sched),base_env)
        path_planner = ISPS()
        solution, cost = solve!(path_planner,PC_MAPF(search_env))
        # @show convert_to_vertex_lists(solution.route_plan)
    end
    let
        search_env = construct_search_env(solver,deepcopy(sched),base_env)
        path_planner = ISPS()
        node = initialize_root_node(search_env)
        @test low_level_search!(path_planner,PC_MAPF(search_env),node)
        # @show convert_to_vertex_lists(node.solution)
    end
    # Test CBS
    let
        search_env = construct_search_env(solver,deepcopy(sched),base_env)
        pc_mapf = PC_MAPF(search_env)
        path_planner = CBSSolver(ISPS())
        solution, cost = solve!(path_planner,pc_mapf)
        # solution, cost = CRCBS.cbs!(path_planner,pc_mapf)
        # @show convert_to_vertex_lists(solution.route_plan)
    end
    let
        # NOTE Freezing here
        path_planner = CBSSolver(ISPS())
        env, cost = plan_route!(path_planner,deepcopy(sched),pc_tapf)
        # @show convert_to_vertex_lists(get_route_plan(env))
    end
    # test NBS
    let
        solver = NBSSolver(assignment_model = TaskGraphsMILPSolver(GreedyAssignment()))
        # solver = NBSSolver(assignment_model = TaskGraphsMILPSolver())
        set_iteration_limit!(solver,3)
        env, cost = solve!(solver,base_env)
        # @show convert_to_vertex_lists(get_route_plan(env))
        # @show get_logger(solver)
        # @show optimality_gap(solver) > 0
        @test validate(env.schedule)
        @test validate(env.schedule,convert_to_vertex_lists(get_route_plan(env)))
    end
end
# test task assignment
let
    # init search env
    for (i, f) in enumerate(pctapf_test_problems())
        for cost_model in [SumOfMakeSpans(), MakeSpan()]
            let
                costs = Float64[]
                for solver in [
                        NBSSolver(assignment_model = TaskGraphsMILPSolver(AssignmentMILP())),
                        NBSSolver(assignment_model = TaskGraphsMILPSolver(SparseAdjacencyMILP())),
                        NBSSolver(assignment_model = TaskGraphsMILPSolver(GreedyAssignment())),
                        ]
                    # set_verbosity!(assignment_solver(solver),3)
                    # @show i, f, solver
                    pc_tapf = f(solver;cost_function=cost_model,verbose=false);
                    search_env = pc_tapf.env
                    prob = formulate_assignment_problem(solver.assignment_model,pc_tapf)
                    sched, cost = solve_assignment_problem!(solver.assignment_model,prob,pc_tapf)
                    if !isa(solver.assignment_model.milp,GreedyAssignment)
                        push!(costs, cost)
                    end
                    @test validate(sched)
                    @test cost != typemax(Int)
                    @test cost != Inf
                end
                # @show costs
                @test all(costs .== costs[1])
            end
        end
    end
end
# test full solver
let
    # init search env
    for (i, f) in enumerate([
        pctapf_problem_1,
        pctapf_problem_2,
        pctapf_problem_3,
        pctapf_problem_4,
        pctapf_problem_5,
        pctapf_problem_6,
        pctapf_problem_7,
        pctapf_problem_8,
        ])
        for cost_model in [SumOfMakeSpans(), MakeSpan()]
            let
                costs = Float64[]
                for solver in [
                        NBSSolver(assignment_model = TaskGraphsMILPSolver(AssignmentMILP())),
                        NBSSolver(assignment_model = TaskGraphsMILPSolver(SparseAdjacencyMILP())),
                        NBSSolver(assignment_model = TaskGraphsMILPSolver(GreedyAssignment())),
                        NBSSolver(
                            assignment_model = TaskGraphsMILPSolver(GreedyAssignment()),
                            path_planner = PIBTPlanner{NTuple{3,Float64}}()
                            ),
                        ]
                    # set_verbosity!(solver,3)
                    set_iteration_limit!(solver,1)
                    set_iteration_limit!(route_planner(solver),100)
                    pc_tapf = f(solver;cost_function=cost_model,verbose=false);
                    env, cost = solve!(solver,pc_tapf)
                    # @show convert_to_vertex_lists(get_route_plan(env))
                    if !isa(solver.assignment_model.milp,GreedyAssignment)
                        push!(costs, cost[1])
                    end
                    @test validate(env.schedule)
                    @test validate(env.schedule,convert_to_vertex_lists(get_route_plan(env)))
                    @test cost != typemax(Int)
                    @test cost != Inf
                end
                # @show costs
                @test all(costs .== costs[1])
            end
        end
    end
end
# Test that optimal collision-free planners obtain the correct costs
let
    cost_model = MakeSpan()
    for (i, (f,expected_cost)) in enumerate([
        (pctapf_problem_1,5),
        (pctapf_problem_2,8),
        (pctapf_problem_3,10),
        (pctapf_problem_4,3),
        (pctapf_problem_5,4),
        ])
        for solver in [
                NBSSolver(assignment_model = TaskGraphsMILPSolver(AssignmentMILP())),
                NBSSolver(assignment_model = TaskGraphsMILPSolver(SparseAdjacencyMILP())),
                ]
            let
                pc_tapf = f(solver;cost_function=cost_model,verbose=false);
                env, cost = solve!(solver,pc_tapf)
                @test expected_cost == cost
            end
        end
    end
    cost_model = SumOfMakeSpans()
    for (i, (f,expected_cost)) in enumerate([
        (pctapf_problem_1,5),
        (pctapf_problem_2,8),
        (pctapf_problem_3,10),
        (pctapf_problem_4,3),
        (pctapf_problem_5,4),
        (pctapf_problem_8,16),
        ])
        for solver in [
                NBSSolver(assignment_model = TaskGraphsMILPSolver(AssignmentMILP())),
                NBSSolver(assignment_model = TaskGraphsMILPSolver(SparseAdjacencyMILP())),
                ]
            let
                pc_tapf = f(solver;cost_function=cost_model,verbose=false);
                env, cost = solve!(solver,pc_tapf)
                @test expected_cost == cost
            end
        end
    end
end
# verify that ISPS+A*sc avoids a collision by exploiting slack
let
    solver = NBSSolver()
    cost_model = MakeSpan()
    pc_tapf = pctapf_problem_3(solver;cost_function=cost_model,verbose=false);
    base_env = pc_tapf.env

    prob = formulate_assignment_problem(assignment_solver(solver),pc_tapf)
    project_schedule, cost = solve_assignment_problem!(assignment_solver(solver),prob,pc_tapf)
    @test termination_status(prob) == MOI.OPTIMAL

    env = construct_search_env(solver,project_schedule,base_env)
    pc_mapf = PC_MAPF(env)
    node = initialize_root_node(solver,pc_mapf)
    low_level_search!(route_planner(solver),pc_mapf,node)

    path1 = convert_to_vertex_lists(node.solution.route_plan.paths[1])
    path2 = convert_to_vertex_lists(node.solution.route_plan.paths[2])
    @test path2[2] == 5 # test that robot 2 will indeed wait for robot 1
    @test path1 == [2, 6, 10, 14, 18, 22, 26, 30, 31, 32, 32]
    @test path2 == [5, 5, 6,  7,  8,  12, 12, 12, 12, 12, 16]
    # @show path1
    # @show path2
end
# Non-zero process time
let
    """
    Verify that the algorithm correctly handles scenarios where process time
    is non-zero
    """
    println("Test non-zero wait time:")
    for Δt in [0,1,4]
        let
            solver = NBSSolver()
            cost_model=MakeSpan()
            pc_tapf = pctapf_problem_6(solver;
                cost_function=cost_model,
                verbose=false,
                Δt_op=Δt
                );
            env, cost = solve!(solver,pc_tapf)
            paths = convert_to_vertex_lists(get_route_plan(env))
            # @show paths[1]
            # @show paths[2]
            @test paths[1][8+Δt] == 13
            @test paths[1][9+Δt] == 17
        end
    end
end
# Non-zero collection time
let
    """
    Verify that the algorithm correctly handles scenarios where collection time
    and deposit time are non-zero.
    """
    println("Test non-zero collection time:")
    Δt_deposit_1 = 1
    for (Δt_collect_2, true_cost) in zip([0,1,2,3,4],[6,7,8,8,8])
        solver = NBSSolver()
        cost_model = MakeSpan()
        pc_tapf = pctapf_problem_7(solver;
            cost_function=cost_model,
            verbose=false,
            Δt_op=0,
            Δt_collect=[0,Δt_collect_2,0],
            Δt_deposit=[Δt_deposit_1,0,0]
            );
        env, cost = solve!(solver,pc_tapf)
        paths = convert_to_vertex_lists(get_route_plan(env))
        @test cost == true_cost
    end
end
# Collaborative transport problems
let
    cost_model = MakeSpan()
    for (i, (f,expected_cost,expected_paths)) in enumerate([
        (pctapf_problem_11, 6, [
            [1,2,6,10,14,13,9],
            [4,3,7,11,15],
            [15,11,10,9,5,1],
        ]),
        ])
        for solver in [
                NBSSolver(assignment_model = TaskGraphsMILPSolver(SparseAdjacencyMILP())),
                ]
            let
                pc_tapf = f(solver;cost_function=cost_model);
                c_pc_tapf = C_PC_TAPF(pc_tapf.env)
                # set_verbosity!(solver,4)
                env, cost = solve!(solver,c_pc_tapf)
                paths = convert_to_vertex_lists(get_route_plan(env))
                @log_info(-1,verbosity(solver),sprint_route_plan(get_route_plan(env)))
                # @show paths
                @test cost == expected_cost
                for (expected_path, path) in zip(expected_paths,paths)
                    @test length(path) <= maximum(map(length,expected_paths))
                    @test path[1:length(expected_path)] == expected_path
                end
            end
        end
    end
end
# Gap repair problems
let
    solver = NBSSolver()
    set_iteration_limit!(low_level(route_planner(solver)),1)
    pc_tapf = pctapf_problem_12(solver)
    env, cost = solve!(solver,pc_tapf)
    cost
    # @show get_route_plan(env)
end
# Collaborative gap repair problems
let
    solver = NBSSolver()
    set_iteration_limit!(low_level(route_planner(solver)),1)
    pc_tapf = C_PC_TAPF(pctapf_problem_13(solver).env)
    env, cost = solve!(solver,pc_tapf)
    cost
    # @show get_route_plan(env)
end


text:
using Test
using Latexify
using DiffEqBiological

### Tips for generating latex tests:
### Latexify has an unexported macro:
###
### Latexify.@generate_test
###
### which generates a test using a given latexify function.
### For example:
###
### Latexify.@generate_test latexify([1, 2, 3], [4, 5, 6]; env=:mdtable)
###
### This puts a ready-made test in your clipboard which you can paste into the
### test file.
###
### Just be sure to remove all such macros before you commit a change since it
### will cause issues with Travis.

r = @reaction_network begin
    hillR(X2,v1,K1,n1)*hill(X4,v1,K1,n1), ∅ → X1
    hill(X5,v2,K2,n2), ∅ → X2
    hill(X3,v3,K3,n3), ∅ → X3
    hillR(X1,v4,K4,n4), ∅ → X4
    hill(X2,v5,K5,n5), ∅ → X5
    (k1,k2), X2 ⟷ X1 + 2X4
    (k3,k4), X4 ⟷ X3
    (k5,k6), 3X5 + X1 ⟷ X2
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d5

@test latexify(r; noise=true, cdot=false) ==
raw"\begin{align}
\mathrm{dX1}\left( t \right) =& \left( \frac{v1 K1^{n1}}{K1^{n1} + X2^{n1}} \frac{v1 X4^{n1}}{K1^{n1} + X4^{n1}} + k1 X2 - \frac{k2}{2} X1 X4^{2} - \frac{k5}{6} X5^{3} X1 + k6 X2 - d1 X1 \right) dt + \sqrt{\left\|\frac{v1 K1^{n1}}{K1^{n1} + X2^{n1}} \frac{v1 X4^{n1}}{K1^{n1} + X4^{n1}}\right\|} \mathrm{dW_1}\left( t \right) + \sqrt{\left\|k1 X2\right\|} \mathrm{dW_6}\left( t \right) - \sqrt{\left\|\frac{k2}{2} X1 X4^{2}\right\|} \mathrm{dW_7}\left( t \right) - \sqrt{\left\|\frac{k5}{6} X5^{3} X1\right\|} \mathrm{dW_10}\left( t \right) + \sqrt{\left\|k6 X2\right\|} \mathrm{dW_11}\left( t \right) - \sqrt{\left\|d1 X1\right\|} \mathrm{dW_12}\left( t \right) \\
\mathrm{dX2}\left( t \right) =& \left( \frac{v2 X5^{n2}}{K2^{n2} + X5^{n2}} - k1 X2 + \frac{k2}{2} X1 X4^{2} + \frac{k5}{6} X5^{3} X1 - k6 X2 - d2 X2 \right) dt + \sqrt{\left\|\frac{v2 X5^{n2}}{K2^{n2} + X5^{n2}}\right\|} \mathrm{dW_2}\left( t \right) - \sqrt{\left\|k1 X2\right\|} \mathrm{dW_6}\left( t \right) + \sqrt{\left\|\frac{k2}{2} X1 X4^{2}\right\|} \mathrm{dW_7}\left( t \right) + \sqrt{\left\|\frac{k5}{6} X5^{3} X1\right\|} \mathrm{dW_10}\left( t \right) - \sqrt{\left\|k6 X2\right\|} \mathrm{dW_11}\left( t \right) - \sqrt{\left\|d2 X2\right\|} \mathrm{dW_13}\left( t \right) \\
\mathrm{dX3}\left( t \right) =& \left( \frac{v3 X3^{n3}}{K3^{n3} + X3^{n3}} + k3 X4 - k4 X3 - d3 X3 \right) dt + \sqrt{\left\|\frac{v3 X3^{n3}}{K3^{n3} + X3^{n3}}\right\|} \mathrm{dW_3}\left( t \right) + \sqrt{\left\|k3 X4\right\|} \mathrm{dW_8}\left( t \right) - \sqrt{\left\|k4 X3\right\|} \mathrm{dW_9}\left( t \right) - \sqrt{\left\|d3 X3\right\|} \mathrm{dW_14}\left( t \right) \\
\mathrm{dX4}\left( t \right) =& \left( \frac{v4 K4^{n4}}{K4^{n4} + X1^{n4}} + 2 k1 X2 -2 \frac{k2}{2} X1 X4^{2} - k3 X4 + k4 X3 - d4 X4 \right) dt + \sqrt{\left\|\frac{v4 K4^{n4}}{K4^{n4} + X1^{n4}}\right\|} \mathrm{dW_4}\left( t \right) + 2 \sqrt{\left\|k1 X2\right\|} \mathrm{dW_6}\left( t \right) -2 \sqrt{\left\|\frac{k2}{2} X1 X4^{2}\right\|} \mathrm{dW_7}\left( t \right) - \sqrt{\left\|k3 X4\right\|} \mathrm{dW_8}\left( t \right) + \sqrt{\left\|k4 X3\right\|} \mathrm{dW_9}\left( t \right) - \sqrt{\left\|d4 X4\right\|} \mathrm{dW_15}\left( t \right) \\
\mathrm{dX5}\left( t \right) =& \left( \frac{v5 X2^{n5}}{K5^{n5} + X2^{n5}} -3 \frac{k5}{6} X5^{3} X1 + 3 k6 X2 - d5 X5 \right) dt + \sqrt{\left\|\frac{v5 X2^{n5}}{K5^{n5} + X2^{n5}}\right\|} \mathrm{dW_5}\left( t \right) -3 \sqrt{\left\|\frac{k5}{6} X5^{3} X1\right\|} \mathrm{dW_10}\left( t \right) + 3 \sqrt{\left\|k6 X2\right\|} \mathrm{dW_11}\left( t \right) - \sqrt{\left\|d5 X5\right\|} \mathrm{dW_16}\left( t \right)
\end{align}
"

r = @reaction_network begin
    (hill(B, p_a, k, n), d_a), 0 ↔ A
    (p_b, d_b), 0 ↔ B
    (r_a, r_b), 3B ↔ A
end p_a k n d_a p_b d_b r_a r_b

@test latexify(r; noise=true, bracket=true) ==
raw"\begin{align}
\mathrm{\mathrm{d}\left[A\right]}\left( t \right) =& \left( \frac{p_{a} \cdot \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}} - d_{a} \cdot \left[ A \right] + \frac{r_{a}}{6} \cdot \left[ B \right]^{3} - r_{b} \cdot \left[ A \right] \right) \cdot dt + \sqrt{\left\|\frac{p_{a} \cdot \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}}\right\|} \cdot \mathrm{dW_1}\left( t \right) - \sqrt{\left\|d_{a} \cdot \left[ A \right]\right\|} \cdot \mathrm{dW_2}\left( t \right) + \sqrt{\left\|\frac{r_{a}}{6} \cdot \left[ B \right]^{3}\right\|} \cdot \mathrm{dW_5}\left( t \right) - \sqrt{\left\|r_{b} \cdot \left[ A \right]\right\|} \cdot \mathrm{dW_6}\left( t \right) \\
\mathrm{\mathrm{d}\left[B\right]}\left( t \right) =& \left( p_{b} - d_{b} \cdot \left[ B \right] -3 \cdot \frac{r_{a}}{6} \cdot \left[ B \right]^{3} + 3 \cdot r_{b} \cdot \left[ A \right] \right) \cdot dt + \sqrt{\left\|p_{b}\right\|} \cdot \mathrm{dW_3}\left( t \right) - \sqrt{\left\|d_{b} \cdot \left[ B \right]\right\|} \cdot \mathrm{dW_4}\left( t \right) -3 \cdot \sqrt{\left\|\frac{r_{a}}{6} \cdot \left[ B \right]^{3}\right\|} \cdot \mathrm{dW_5}\left( t \right) + 3 \cdot \sqrt{\left\|r_{b} \cdot \left[ A \right]\right\|} \cdot \mathrm{dW_6}\left( t \right)
\end{align}
"

@test latexify(r; noise_only=true, bracket=true) ==
raw"\begin{align}
\mathrm{\mathrm{d}\left[A\right]}\left( t \right) ∝& \sqrt{\left\|\frac{p_{a} \cdot \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}}\right\|} \cdot \mathrm{dW_1}\left( t \right) - \sqrt{\left\|d_{a} \cdot \left[ A \right]\right\|} \cdot \mathrm{dW_2}\left( t \right) + \sqrt{\left\|\frac{r_{a}}{6} \cdot \left[ B \right]^{3}\right\|} \cdot \mathrm{dW_5}\left( t \right) - \sqrt{\left\|r_{b} \cdot \left[ A \right]\right\|} \cdot \mathrm{dW_6}\left( t \right) \\
\mathrm{\mathrm{d}\left[B\right]}\left( t \right) ∝& \sqrt{\left\|p_{b}\right\|} \cdot \mathrm{dW_3}\left( t \right) - \sqrt{\left\|d_{b} \cdot \left[ B \right]\right\|} \cdot \mathrm{dW_4}\left( t \right) -3 \cdot \sqrt{\left\|\frac{r_{a}}{6} \cdot \left[ B \right]^{3}\right\|} \cdot \mathrm{dW_5}\left( t \right) + 3 \cdot \sqrt{\left\|r_{b} \cdot \left[ A \right]\right\|} \cdot \mathrm{dW_6}\left( t \right)
\end{align}
"

@test latexify(r; noise_only=true, bracket=true, cdot=false) ==
raw"\begin{align}
\mathrm{\mathrm{d}\left[A\right]}\left( t \right) ∝& \sqrt{\left\|\frac{p_{a} \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}}\right\|} \mathrm{dW_1}\left( t \right) - \sqrt{\left\|d_{a} \left[ A \right]\right\|} \mathrm{dW_2}\left( t \right) + \sqrt{\left\|\frac{r_{a}}{6} \left[ B \right]^{3}\right\|} \mathrm{dW_5}\left( t \right) - \sqrt{\left\|r_{b} \left[ A \right]\right\|} \mathrm{dW_6}\left( t \right) \\
\mathrm{\mathrm{d}\left[B\right]}\left( t \right) ∝& \sqrt{\left\|p_{b}\right\|} \mathrm{dW_3}\left( t \right) - \sqrt{\left\|d_{b} \left[ B \right]\right\|} \mathrm{dW_4}\left( t \right) -3 \sqrt{\left\|\frac{r_{a}}{6} \left[ B \right]^{3}\right\|} \mathrm{dW_5}\left( t \right) + 3 \sqrt{\left\|r_{b} \left[ A \right]\right\|} \mathrm{dW_6}\left( t \right)
\end{align}
"

@test latexify(r; noise=true, noise_var=:Noise, bracket=true) ==
raw"\begin{align}
\mathrm{\mathrm{d}\left[A\right]}\left( t \right) =& \left( \frac{p_{a} \cdot \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}} - d_{a} \cdot \left[ A \right] + \frac{r_{a}}{6} \cdot \left[ B \right]^{3} - r_{b} \cdot \left[ A \right] \right) \cdot dt + \sqrt{\left\|\frac{p_{a} \cdot \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}}\right\|} \cdot \mathrm{dNoise_1}\left( t \right) - \sqrt{\left\|d_{a} \cdot \left[ A \right]\right\|} \cdot \mathrm{dNoise_2}\left( t \right) + \sqrt{\left\|\frac{r_{a}}{6} \cdot \left[ B \right]^{3}\right\|} \cdot \mathrm{dNoise_5}\left( t \right) - \sqrt{\left\|r_{b} \cdot \left[ A \right]\right\|} \cdot \mathrm{dNoise_6}\left( t \right) \\
\mathrm{\mathrm{d}\left[B\right]}\left( t \right) =& \left( p_{b} - d_{b} \cdot \left[ B \right] -3 \cdot \frac{r_{a}}{6} \cdot \left[ B \right]^{3} + 3 \cdot r_{b} \cdot \left[ A \right] \right) \cdot dt + \sqrt{\left\|p_{b}\right\|} \cdot \mathrm{dNoise_3}\left( t \right) - \sqrt{\left\|d_{b} \cdot \left[ B \right]\right\|} \cdot \mathrm{dNoise_4}\left( t \right) -3 \cdot \sqrt{\left\|\frac{r_{a}}{6} \cdot \left[ B \right]^{3}\right\|} \cdot \mathrm{dNoise_5}\left( t \right) + 3 \cdot \sqrt{\left\|r_{b} \cdot \left[ A \right]\right\|} \cdot \mathrm{dNoise_6}\left( t \right)
\end{align}
"

@test latexify(r; bracket=true) ==
raw"\begin{align}
\frac{d\left[ A \right](t)}{dt} =& \frac{p_{a} \cdot \left[ B \right]^{n}}{k^{n} + \left[ B \right]^{n}} - d_{a} \cdot \left[ A \right] + \frac{r_{a}}{6} \cdot \left[ B \right]^{3} - r_{b} \cdot \left[ A \right] \\
\frac{d\left[ B \right](t)}{dt} =& p_{b} - d_{b} \cdot \left[ B \right] -3 \cdot \frac{r_{a}}{6} \cdot \left[ B \right]^{3} + 3 \cdot r_{b} \cdot \left[ A \right]
\end{align}
"

@test latexify(r; bracket=false) ==
raw"\begin{align}
\frac{dA(t)}{dt} =& \frac{p_{a} \cdot B^{n}}{k^{n} + B^{n}} - d_{a} \cdot A + \frac{r_{a}}{6} \cdot B^{3} - r_{b} \cdot A \\
\frac{dB(t)}{dt} =& p_{b} - d_{b} \cdot B -3 \cdot \frac{r_{a}}{6} \cdot B^{3} + 3 \cdot r_{b} \cdot A
\end{align}
"


text:
"""
    test_init_field(N, L, ftf)

Test that the field initialized by the field type function `ftf` on the grid g
has the correct size.
"""
function correct_field_size(arch, grid, fieldtype)
    f = fieldtype(arch, grid)
    return size(f) == size(grid)
end

"""
    test_set_field(N, L, ftf, val)

Test that the field initialized by the field type function `ftf` on the grid g
can be correctly filled with the value `val` using the `set!(f::AbstractField, v)`
function.
"""
function correct_field_value_was_set(arch, grid, fieldtype, val::Number)
    f = fieldtype(arch, grid)
    set!(f, val)
    return data(f) ≈ val * ones(size(f))
end

function correct_field_addition(arch, grid, fieldtype, val1::Number, val2::Number)
    f1 = fieldtype(arch, grid)
    f2 = fieldtype(arch, grid)

    set!(f1, val1)
    set!(f2, val2)
    f3 = f1 + f2

    val3 = convert(mm.float_type, val1) + convert(mm.float_type, val2)
    f_ans = val3 * ones(size(f1))
    return f3.data ≈ f_ans
end

@testset "Fields" begin
    println("Testing fields...")

    N = (4, 6, 8)
    L = (2π, 3π, 5π)

    fieldtypes = (CellField, FaceFieldX, FaceFieldY, FaceFieldZ)

    @testset "Field initialization" begin
        println("  Testing field initialization...")
        for arch in archs, FT in float_types
            grid = RegularCartesianGrid(FT, N, L)

            for fieldtype in fieldtypes
                @test correct_field_size(arch, grid, fieldtype)
            end
        end
    end

    int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4), Int128(-5)]
    uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10), UInt128(11)]
    float_vals = Any[0.0, -0.0, 6e-34, 1.0f10]
    rational_vals = Any[1//11, -23//7]
    other_vals = Any[π]
    vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals)

    @testset "Setting fields" begin
        println("  Testing field setting...")

        for arch in archs, FT in float_types
            grid = RegularCartesianGrid(FT, N, L)

            for fieldtype in fieldtypes, val in vals
                @test correct_field_value_was_set(arch, grid, fieldtype, val)
            end
        end
    end
end


text:
using Gadfly

module ComposeVideo

	using Gadfly

	export
	VideoBuilder,
	add_frame,
	display_video

	include("to_video.jl")

end


text:
using Tirocinio
using Test

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
# DualPol_simulator.jl
#
# A dual-polarization radar simulator
#
# By Michael M. Bell
# Copyright 2015
#
# This program will simulate the backscatter from
# weather targets using a single, double, or fully
# explicit bin dropsize distribution. It is currently
# configured to work with WRF output, and was used in
# Brown, B. R., M. M. Bell, and A. J. Frambach, 2016:
# "Validation of Simulated Hurricane Drop Size Distributions
# using Polarimetric Radar", Geophys. Res. Lett., 42,
# doi:10.1002/2015GL067278.
#
# Inputs: Type of DSD (single, double, or full)
#         Name of WRF output file
# Outputs: NetCDF file with polarimetric radar variables
#
# Current output is radar reflectivity factor at horizontal
# and vertical polarization, and differential reflectivity
# for raindrops only. Additional variables and hydrometeor
# species will be added in due course
#

using ArgParse
using NetCDF
using DataStructures
#using Debug

# Global constants
xam_r = pi*997.0/6.0
xbm_r = 3.0
xmu_r = 0.0
xobmr = 1.0/xbm_r
xcre = 1. + 2.*xbm_r + xmu_r

# Angular moments from Jung et al (2010,JAMC) Eqns (4)
# Canting angle distribution width (sigma) equal 10 degrees only for rain & oblate crystals
sig = deg2rad(10.0)
r = exp(-2.0*sig*sig)
A = ( 0.375 + (0.5)*r + (0.125)*r^4 )
B = ( 0.375 - (0.5)*r + (0.125)*r^4 )
C = (0.125)*( 1 - r^4 )

# Parse the command line arguments
function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--dsdtype","-d"
    help = "single, double, full"
    arg_type = String
    required = true
    "--output","-o"
    help = "path to output file"
    arg_type = String
    default = "dualpol_radar.nc"
    "input"
    help = "path to input file"
    arg_type = String
    required = true
  end
  return parse_args(s)
end

# Calculate the radar variables using a prescribed gamma distribution
function calc_radar_gamma_dsd(N0,lambda)
  # Define radar constants

  # Fixed complex index of refraction
  eps = 8.88 + 0.63im

  # Radar wavelength (mm)
  wavelength = 100

  # Dielectric constant
  k = (abs2(eps) - 1)/(abs2(eps) +2)

  # Define the scattering amplitudes
  # Pre-calculated using T-matrix scattering code from Mishchenko (2000)
  s_amp = Array{Complex64}(16,2)
  s_amp = [ [5.93590e-05+(3.07940e-07)im] [-5.94150e-05+(-3.08530e-07)im]
            [4.70300e-04+(2.29610e-06)im] [-4.76460e-04+(-2.35930e-06)im]
            [1.55640e-03+(6.77020e-06)im] [-1.61730e-03+(-7.37320e-06)im]
            [3.58870e-03+(1.28470e-05)im] [-3.86500e-03+(-1.54380e-05)im]
            [6.77460e-03+(1.72310e-05)im] [-7.62290e-03+(-2.45720e-05)im]
            [1.12570e-02+(1.34310e-05)im] [-1.33110e-02+(-2.91580e-05)im]
            [1.71230e-02+(-9.24210e-06)im] [-2.13580e-02+(-1.72640e-05)im]
            [2.44090e-02+(-6.71130e-05)im] [-3.21760e-02+(3.40790e-05)im]
            [3.31110e-02+(-1.84310e-04)im] [-4.61260e-02+(1.67320e-04)im]
            [4.31850e-02+(-3.96140e-04)im] [-6.34650e-02+(4.59880e-04)im]
            [5.45400e-02+(-7.54470e-04)im] [-8.42680e-02+(1.05350e-03)im]
            [6.70130e-02+(-1.33600e-03)im] [-1.08320e-01+(2.21210e-03)im]
            [8.03320e-02+(-2.25460e-03)im] [-1.34920e-01+(4.44170e-03)im]
            [9.40660e-02+(-3.67760e-03)im] [-1.62650e-01+(8.75570e-03)im]
            [1.07580e-01+(-5.84290e-03)im] [-1.88770e-01+(1.73100e-02)im]
            [1.20010e-01+(-9.05960e-03)im] [-2.08160e-01+(3.50790e-02)im]
            [0.0+(0.0)im] [0.0+(0.0)im] ];

  # Calculate the reflectivity using 0.5 mm bins
  zh = 0.0
  zv = 0.0
  ql = 0.0
  mu = 0.0
  h = 8.0/48.0 # Simpson's composite rule (8.0 - 0.0)/ (16 * 3)

  for i in 1:17
    if (i == 1) || (i == 17)
      intcoeff = 1.0
    elseif (i%2 == 0)
      intcoeff = 4.0
    else
      intcoeff = 2.0
    end
      D = i*0.5
      N = N0*(D^mu)*exp(-lambda*D)
      ql += xam_r*N0*(D^(3.0+mu))*exp(-lambda*D)*intcoeff

      # Jung et al (2010,JAMC) Eqns (4) and (3)
      zv += ( B*abs2(s_amp[i,2]) + A*abs2(s_amp[i,1]) + 2*C*real( conj(s_amp[i,1])*s_amp[i,2] ) )*N*intcoeff
      zh += ( A*abs2(s_amp[i,2]) + B*abs2(s_amp[i,1]) + 2*C*real( conj(s_amp[i,1])*s_amp[i,2] ) )*N*intcoeff
   end

  ql *= h*1.0e-9
  zv *= h*(4*wavelength^4)/(pi^4*k^2)
  zh *= h*(4*wavelength^4)/(pi^4*k^2)
  za = N0*gamma(xcre)/(lambda^xcre)

  if (zh > 10.0^(-3.5)) && (zv > 10.0^(-3.5))
   Zdr = 10*log10(zh/zv)
   Zv = 10*log10(zv)
   Zh = 10*log10(zh)
  else
   Zdr = 0.0
   Zv = -35.0
   Zh = -35.0
  end

  if (za > 10.0^(-3.5))
   Za = 10*log10(za)
  else
   Za = -35.0
  end

  return Zdr, Zv, Zh, Za, ql
end

# Calculate the radar variables using a full bin DSD
# This DSD follows the HUJI bin model output
function calc_radar_fulldsd(dsd, flag, qrv, N0,lambda)
  # Fixed complex index of refraction
  eps = 8.88 + 0.63im

  # Radar wavelength (mm)
  wavelength = 100

  # Dielectric constant
  k = (abs2(eps) - 1)/(abs2(eps) +2)

  # Define the scattering amplitudes
  s_amp = Array{Complex64}(33,2)
  s_amp = [ [3.03150e-11+(1.59470e-13)im] [-3.04840e-11+(-1.61260e-13)im] #0.00400
          [6.06310e-11+(3.18960e-13)im] [-6.09670e-11+(-3.22510e-13)im] #0.00500
          [1.21270e-10+(6.37940e-13)im] [-1.21930e-10+(-6.44990e-13)im] #0.00640
          [2.42540e-10+(1.27600e-12)im] [-2.43860e-10+(-1.28990e-12)im] #0.00800
          [4.85100e-10+(2.55210e-12)im] [-4.87720e-10+(-2.57980e-12)im] #0.01000
          [9.70240e-10+(5.10470e-12)im] [-9.75410e-10+(-5.15930e-12)im] #0.01260
          [1.94060e-09+(1.02110e-11)im] [-1.95080e-09+(-1.03180e-11)im] #0.01600
          [3.88150e-09+(2.04240e-11)im] [-3.90140e-09+(-2.06340e-11)im] #0.02020
          [7.76370e-09+(4.08550e-11)im] [-7.80230e-09+(-4.12630e-11)im] #0.02540
          [1.55290e-08+(8.17280e-11)im] [-1.56040e-08+(-8.25140e-11)im] #0.03200
          [3.10630e-08+(1.63500e-10)im] [-3.12050e-08+(-1.65000e-10)im] #0.04040
          [6.21370e-08+(3.27090e-10)im] [-6.24050e-08+(-3.29910e-10)im] #0.05080
          [1.24300e-07+(6.54380e-10)im] [-1.24800e-07+(-6.59620e-10)im] #0.06400
          [2.48660e-07+(1.30920e-09)im] [-2.49560e-07+(-1.31870e-09)im] #0.08060
          [4.97470e-07+(2.61930e-09)im] [-4.99040e-07+(-2.63580e-09)im] #0.10160
          [9.95270e-07+(5.23990e-09)im] [-9.97890e-07+(-5.26760e-09)im] #0.12800
          [1.99120e-06+(1.04810e-08)im] [-1.99540e-06+(-1.05240e-08)im] #0.16120
          [3.98400e-06+(2.09560e-08)im] [-3.98980e-06+(-2.10170e-08)im] #0.20320
          [7.97050e-06+(4.18720e-08)im] [-7.97760e-06+(-4.19470e-08)im] #0.25600
          [1.59440e-05+(8.35630e-08)im] [-1.59520e-05+(-8.36450e-08)im] #0.32260
          [3.18850e-05+(1.66430e-07)im] [-3.18990e-05+(-1.66570e-07)im] #0.40640
          [6.37300e-05+(3.30320e-07)im] [-6.37970e-05+(-3.31020e-07)im] #0.51200
          [1.27260e-04+(6.51900e-07)im] [-1.27630e-04+(-6.55740e-07)im] #0.64500
          [2.53730e-04+(1.27450e-06)im] [-2.55450e-04+(-1.29240e-06)im] #0.81280
          [5.04610e-04+(2.45330e-06)im] [-5.11700e-04+(-2.52600e-06)im] #1.02400
          [9.99600e-04+(4.59930e-06)im] [-1.02620e-03+(-4.86780e-06)im] #1.29020
          [1.96830e-03+(8.22720e-06)im] [-2.06200e-03+(-9.14360e-06)im] #1.62540
          [3.84180e-03+(1.34180e-05)im] [-4.15370e-03+(-1.63240e-05)im] #2.04800
          [7.40440e-03+(1.73950e-05)im] [-8.39530e-03+(-2.58320e-05)im] #2.58040
          [1.40260e-02+(5.27710e-06)im] [-1.70300e-02+(-2.63510e-05)im] #3.25100
          [2.59700e-02+(-8.41070e-05)im] [-3.46010e-02+(5.15370e-05)im] #4.09600
          [4.67000e-02+(-4.92210e-04)im] [-6.97740e-02+(6.08330e-04)im] #5.16060
          [8.03860e-02+(-2.25910e-03)im] [-1.35030e-01+(4.45380e-03)im] #6.50200
          ];

  # Calculate the reflectivity using equal mass bins
  zh = 0.0
  zv = 0.0
  ql = 0.0
  NT = 0.0
  mu = 0.0
  h = 0.5
  r3 = 8.0e-9
  Dprevious = 0.0
  zv1 = 0.0
  zh1 = 0.0
  for i in 1:33
    intcoeff = 1.0
    radius = r3^(1.0 / 3.0)
    mass = (4.0/3.0) * pi * r3 * 1.0e-6
    D = radius*2.0
    binbottom = (D + Dprevious)/2.0
    r3 = r3 * 2.0
    bintop = (D + 2.0 * r3^(1.0 / 3.0))/2.0
    intcoeff = 0.5 * (bintop - binbottom)
    q_bin = dsd[i]
    N = dsd[i]/(mass * (bintop - binbottom))
    Nexp = N0*(D^mu)*exp(-lambda*D)
    ql += q_bin
    NT += N*intcoeff
    if (flag)
      println("$D, $N, $Nexp, $q_bin")
    end
    # Jung et al (2010,JAMC) Eqns (4) and (3)
    zv2 = ( B*abs2(s_amp[i,2]) + A*abs2(s_amp[i,1]) + 2*C*real( conj(s_amp[i,1])*s_amp[i,2] ) ) * N
    zh2 = ( A*abs2(s_amp[i,2]) + B*abs2(s_amp[i,1]) + 2*C*real( conj(s_amp[i,1])*s_amp[i,2] ) ) * N
    zv += (zv1 + zv2)*intcoeff
    zh += (zh1 + zh2)*intcoeff
    Dprevious = D
    zv1 = zv2
    zh1 = zh2
   end
  #ql *= 1.0e-9
  #println("\n$NT, $nrv, $ql, $qrv, $rho")
  zv *= (4*wavelength^4)/(pi^4*k^2)
  zh *= (4*wavelength^4)/(pi^4*k^2)
  za = 10.0^(-3.5)

  if (zh > 10.0^(-3.5)) && (zv > 10.0^(-3.5))
   Zdr = 10*log10(zh/zv)
   Zv = 10*log10(zv)
   Zh = 10*log10(zh)
  else
   Zdr = 0.0
   Zv = -35.0
   Zh = -35.0
  end

  if (za > 10.0^(-3.5))
   Za = 10*log10(za)
  else
   Za = -35.0
  end

  return Zdr, Zv, Zh, Za, ql
end

# Read in the WRF output
function read_nc_var(filename,varnames::Array)
  ###Read in Cartesian nc-file
  ###Only works as intended if you read in more than 1 variable
  #println("Read in nc file ...")
  data = OrderedDict()

  for varname in varnames
    data[string(varname)] = ncread(filename, string(varname))
  end

  return collect(values(data))
end

# Write out the radar variables
function write_ncfile(filename,lat,lon,lev,times,varnames)
  ###Write to nc-file
  println("Write to nc file ...")

  ncvars = NcVar[]
  xatts = Dict("long_name" => "x (longitude)", "units" => "deg", "missing_value" => -999.0, "_FillValue" => -999.0)
  yatts = Dict("long_name" => "y (latitude)",  "units" => "deg", "missing_value" => -999.0, "_FillValue" => -999.0)
  zatts = Dict("long_name" => "z (eta)",  "units" => "unitless", "missing_value" => -999.0, "_FillValue" => -999.0)
  tatts = Dict("long_name" => "time (minutes)",  "units" => "min", "missing_value" => -999.0, "_FillValue" => -999.0)
  x_dim = NcDim("east_west",collect(1:length(lon[:,1,1])),xatts)
  y_dim = NcDim("south_north",collect(1:length(lat[1,:,1])),yatts)
  z_dim = NcDim("bottom_top",collect(1:length(lev[:,1])),zatts)
  t_dim = NcDim("Time",collect(1:length(times)),tatts)

  for varname in varnames
    atts  = Dict("long_name" => varname, "units" => "???", "missing_value" => -999.0, "_FillValue" => -999.0)
    push!(ncvars,NcVar(varname,[x_dim,y_dim,z_dim,t_dim],atts,Float64))
  end
  atts  = Dict("long_name" => "Latitude", "units" => "deg", "missing_value" => -999.0, "_FillValue" => -999.0)
  push!(ncvars,NcVar("XLAT",[x_dim,y_dim,t_dim],atts,Float64))
  atts  = Dict("long_name" => "Longitude", "units" => "deg", "missing_value" => -999.0, "_FillValue" => -999.0)
  push!(ncvars,NcVar("XLON",[x_dim,y_dim,t_dim],atts,Float64))
  atts  = Dict("long_name" => "Time", "units" => "deg", "missing_value" => -999.0, "_FillValue" => -999.0)
  push!(ncvars,NcVar("XTIME",[t_dim],atts,Float64))
  atts  = Dict("long_name" => "Eta Levels", "units" => "deg", "missing_value" => -999.0, "_FillValue" => -999.0)
  push!(ncvars,NcVar("ZNU",[z_dim,t_dim],atts,Float64))
  nc = NetCDF.create(filename,ncvars)

  NetCDF.putvar(nc,"REFL_10CM",refl)
  NetCDF.putvar(nc,"ZDR",ZDR)
  NetCDF.putvar(nc,"ZV",ZV)
  NetCDF.putvar(nc,"ZH",ZH)
  NetCDF.putvar(nc,"Rayleigh",DBZ)
  NetCDF.putvar(nc,"Zdiff",Zdiff)
  NetCDF.putvar(nc,"QRAIN",qrain)
  NetCDF.putvar(nc,"QNRAIN",qnrain)
  NetCDF.putvar(nc,"XLAT",lat)
  NetCDF.putvar(nc,"XLON",lon)
  NetCDF.putvar(nc,"XTIME",times)
  NetCDF.putvar(nc,"ZNU",lev)

  NetCDF.close(nc)
  return 0
end

# Main program
args = parse_commandline()
filein = args["input"]
fileout = args["output"]
dsdtype = args["dsdtype"]

### Read in Cartesian analysis
println("Reading ",filein)
fileinfo = ncinfo(filein)

dims = ["XLAT","XLONG","ZNU","XTIME"]
lat,lon,lev,times = read_nc_var(filein,dims)

if (dsdtype == "single")
  # Single moment
  vars = ["QRAIN","QVAPOR","PB","P","T","REFL_10CM"]
  qrain,qv,pb,pp,theta,refl = read_nc_var(filein,vars)
elseif (dsdtype == "double")
  # Double moment
  vars = ["QRAIN","QNRAIN","QVAPOR","PB","P","T","REFL_10CM"]
  qrain,qnrain,qv,pb,pp,theta,refl = read_nc_var(filein,vars)
elseif (dsdtype == "full")
  # Spectral bin
  vars = ["QRAIN","QNRAIN","QVAPOR","PB","P","T","REFL_10CM",
  "ff1i01","ff1i02","ff1i03","ff1i04","ff1i05","ff1i06","ff1i07","ff1i08",
  "ff1i09","ff1i10","ff1i11","ff1i12","ff1i13","ff1i14","ff1i15","ff1i16",
  "ff1i17","ff1i18","ff1i19","ff1i20","ff1i21","ff1i22","ff1i23","ff1i24",
  "ff1i25","ff1i26","ff1i27","ff1i28","ff1i29","ff1i30","ff1i31","ff1i32",
  "ff1i33"]
  dummyvar =  ncread(filein,"QRAIN")
  ff1 = Array(Array{Float32,4},33)
  for d in 1:33
    ff1[d] = similar(dummyvar); fill!(ff1[d], -999.0)
  end
  ( qrain,qnrain,qv,pb,pp,theta,refl,
  ff1[1],ff1[2],ff1[3],ff1[4],ff1[5],ff1[6],ff1[7],ff1[8],
  ff1[9],ff1[10],ff1[11],ff1[12],ff1[13],ff1[14],ff1[15],ff1[16],
  ff1[17],ff1[18],ff1[19],ff1[20],ff1[21],ff1[22],ff1[23],ff1[24],
  ff1[25],ff1[26],ff1[27],ff1[28],ff1[29],ff1[30],ff1[31],ff1[32],ff1[33] ) = read_nc_var(filein,vars)
else
  println(dsdtype, " not recognized")
  exit
end

d1,d2,d3,d4 = size(qrain)

# Initialize new radar variables
ZDR      = similar(qrain); fill!(ZDR, -999.0)
ZV       = similar(qrain); fill!(ZV, -999.0)
ZH       = similar(qrain); fill!(ZH, -999.0)
DBZ      = similar(qrain); fill!(DBZ, -999.0)
Zdiff    = similar(qrain); fill!(Zdiff, -999.0)
ql       = similar(qrain); fill!(ql, -999.0)
qdiff    = similar(qrain); fill!(qdiff, -999.0)
if (dsdtype == "single")
  qnrain = similar(qrain); fill!(qnrain, -999.0)
elseif (dsdtype == "full")
  dsd = Array(Float32,33)
end

# Calculate the radar variables
println("Calculating radar variables...")
gamma1 = gamma(1. + xbm_r + xmu_r)
gamma2 = gamma(1. + xmu_r)
xorg2 = 1.0/gamma2
for t in 1:d4
  for k in 1:d3
    for j in 1:d2
      for i in 1:d1
        qvapor = max(1.0e-10,qv[i,j,k,t])
        pressure = pb[i,j,k,t]+pp[i,j,k,t]
        tempk = (theta[i,j,k,t]+300.0)*(pressure/100000.0)^(2.0/7.0)
        rho = 0.622*pressure/(287.15*tempk*(qvapor+0.622))
        if (qrain[i,j,k,t] > 1.0e-9)
          qrv = qrain[i,j,k,t]*rho
          if (dsdtype == "single")
            N0 = 8000
            lambda = ((xam_r*gamma1*N0*1000.0/qrv)^(1./(1. + xbm_r)))/1000.0
            qnrain[i,j,k,t] = N0*gamma2/(lambda^(1. + xmu_r))
          else
            if (qnrain[i,j,k,t] > 0.0)
              nrv = qnrain[i,j,k,t]*rho
              lambda = ((xam_r*gamma1*xorg2*nrv/qrv)^xobmr)/1000.0
            else
              nrv = 0.0
              lambda = 1./20.E-3
            end
          end
          if (lambda < 1./2800.E-3)
            lambda = 1./2800.E-3
          elseif (lambda > 1./20.E-3)
            lambda = 1./20.E-3
          end
          if (dsdtype == "single")
            ZDR[i,j,k,t], ZV[i,j,k,t], ZH[i,j,k,t], DBZ[i,j,k,t], ql[i,j,k,t] = calc_radar_gamma_dsd(N0,lambda)
          elseif (dsdtype == "double")
            N0 = nrv*xorg2*lambda^(1. + xmu_r)
            ZDR[i,j,k,t], ZV[i,j,k,t], ZH[i,j,k,t], DBZ[i,j,k,t], ql[i,j,k,t] = calc_radar_gamma_dsd(N0,lambda)
          else
            if (i == 215 && j == 200 && k == 10 && t == 3)
              flag = true
              println("Found it!")
            else
              flag = false
            end
            for d in [1:33]
              dsd[d] = ff1[d][i,j,k,t]*rho
            end
            N0 = nrv*xorg2*lambda^(1. + xmu_r)
            ZDR[i,j,k,t], ZV[i,j,k,t], ZH[i,j,k,t], DBZ[i,j,k,t], ql[i,j,k,t] = calc_radar_fulldsd(dsd, flag, qrv, N0, lambda)
          end

          Zdiff[i,j,k,t] = refl[i,j,k,t] - ZH[i,j,k,t]
          ql[i,j,k,t] = 1.e-9*xam_r*N0*gamma1/(rho*lambda^(1. + xbm_r + xmu_r))
          if (DBZ[i,j,k,t] > 100.0)
            println("Large Z :",DBZ[i,j,k,t])
            println(lat[i,j,t],",",lon[i,j,t])
            println(i,",",j,",",k,":",nrv, " ", qrv, " ", rho, " ", refl[i,j,k,t])
            println(N0,",",lambda)
          elseif (DBZ[i,j,k,t] < -35.0)
            println("Small Z :",DBZ[i,j,k,t])
            println(lat[i,j,t],",",lon[i,j,t])
            println(i,",",j,",",k,":",nrv, " ", qrv, " ", rho, " ", refl[i,j,k,t])
            println(N0,",",lambda)
          end
        else
          qrv = 1.0e-12
          nrv = 1.0e-12
          ql[i,j,k,t] = 1.0e-9
          qrain[i,j,k,t] = 1.0e-9
          ZDR[i,j,k,t] = 0.0
          ZV[i,j,k,t] = ZH[i,j,k,t] = DBZ[i,j,k,t] = -35.0
       end
       qdiff[i,j,k,t] = ql[i,j,k,t] / qrain[i,j,k,t]
     end
   end
 end
end

vars_out = ["REFL_10CM","ZDR","ZV","ZH","Rayleigh","Zdiff","QRAIN","QNRAIN"]
println("Writing ", fileout)
write_ncfile(fileout,lat,lon,lev,times,vars_out)


text:
function InterMap(L::Lattice{D}, cr::CartesianRange{CartesianIndex{D}}, filename::AbstractString, cutoff::Int64, origin::Vector{Float64}) where D
    vals = Vector{Float64}()
    distances = Vector{Float64}()
    indices = Vector{Tuple{NTuple{D,Int64},Int64}}()
    vectors = Vector{Vector{Float64}}()
    for ci in cr
        for c = 1:L.cell_size
            push!(indices, (ci.I, c))
            vector = L.basis_vectors*collect(ci.I) + L.cell_vectors[c]-origin
            push!(distances, norm(vector))
            push!(vals, L.Jfunc(vector))
            push!(vectors, vector)
        end
    end
    p = sortperm(vals; by = abs, rev = true)
    permute!(vals, p)
    permute!(distances, p)
    permute!(indices, p)
    permute!(vectors, p)
    io = open(filename, "w")
    for i = 1:min(cutoff, length(vals))
        print(io, vals[i], "\t\t", distances[i], "\t\t", vectors[i], "\t\t", indices[i][1], "\t", indices[i][2], "\n")
    end
    close(io)
end

function get_num_eff(L::Lattice{D}, cr::CartesianRange{CartesianIndex{D}}, origin::Vector{Float64}) where D
    jj = 0.0
    j = 0.0
    for ci in cr
        for c = 1:L.cell_size
            vector = L.basis_vectors*collect(ci.I) + L.cell_vectors[c]-origin
            if norm(vector) != 0.0
                val = L.Jfunc(vector)^2
                jj += val
                j += val^2
            end
        end
    end
    return jj^2/j
end


text:
Base.promote_rule(::Type{C3}, ::Type{Cgray}) where {C3<:Color3,Cgray<:AbstractGray} = base_colorant_type(C3){promote_type(eltype(C3), eltype(Cgray))}
Base.promote_rule(::Type{C3}, ::Type{Cagray}) where {C3<:Color3,Cagray<:AbstractAGray} = alphacolor(base_colorant_type(C3)){promote_type(eltype(C3), eltype(Cagray))}
Base.promote_rule(::Type{C3}, ::Type{Cgraya}) where {C3<:Color3,Cgraya<:AbstractGrayA} = coloralpha(base_colorant_type(C3)){promote_type(eltype(C3), eltype(Cgraya))}

Base.promote_rule(::Type{C3}, ::Type{Cgray}) where {C3<:Transparent3,Cgray<:AbstractGray} = base_colorant_type(C3){promote_type(eltype(C3), eltype(Cgray))}
Base.promote_rule(::Type{C3}, ::Type{Cgray}) where {C3<:Transparent3,Cgray<:TransparentGray} = base_colorant_type(C3){promote_type(eltype(C3), eltype(Cgray))}


text:
using Base.Test
using QuantumOptics

b = FockBasis(3)

u = Ket[fockstate(b, 1), fockstate(b, 2)]
v = Ket[fockstate(b, 2), fockstate(b, 1)]

bu = SubspaceBasis(u)
bv = SubspaceBasis(v)

T1 = projector(bu, b)
T2 = projector(bv, b)
T12 = projector(bu, bv)

state = fockstate(b, 2)
state_u = Ket(bu, [0, 1])
state_v = Ket(bv, [1., 0])

@test T1*state == state_u
@test T2*state == state_v

state_v = Ket(bv, [1, -1])
state_u = Ket(bu, [-1, 1])

@test T12*state_v == state_u

u2 = Ket[1.5*fockstate(b, 1), fockstate(b, 1) + fockstate(b, 2)]
bu2_orth = subspace.orthonormalize(SubspaceBasis(u))
@test bu2_orth.basisstates == bu.basisstates


text:
abstract type SteadyStateDiffEqAlgorithm <: DiffEqBase.AbstractSteadyStateAlgorithm end

struct SSRootfind{F} <: SteadyStateDiffEqAlgorithm
  nlsolve::F
end
SSRootfind(;nlsolve=(f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,ftol = abstol);res.zero)) = SSRootfind(nlsolve)

struct DynamicSS{A,AT,RT,TS} <: SteadyStateDiffEqAlgorithm
  alg::A
  abstol::AT
  reltol::RT
  tspan::TS
end
DynamicSS(alg; abstol = 1e-8, reltol = 1e-6, tspan = Inf) =
  DynamicSS(alg, abstol, reltol, tspan)

# Backward compatibility:
DynamicSS(alg, abstol, reltol) = DynamicSS(alg; abstol=abstol, reltol=reltol)


text:
pll   = PLL()
ƒ     = -0.1
ϕ     = 3.0
n     = 300
t     = 0:n-1
x     = exp( im*( ƒ*t + ϕ) )
(y,e) = exec( pll, x )

using Winston
plotTable = Table(3,1)
plotTable[1,1] = plot( t, real(x), t, real(y), "r", title = "Real" )
plotTable[2,1] = plot( t, imag(x), t, imag(y), "r", title = "Imaginary" )
plotTable[3,1] = plot( t, e, title = "ϕ error" )
display( plotTable )


text:
using SnoopCompile
invalidations = @snoopr using EAGO

@show length(invalidations)
trees = invalidation_trees(invalidations)
#methinvs = trees[end]
#root = methinvs.backedges[end]
#ascend(root)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function [w,x] = GaussLobattoQuad(n)
w=zeros(n+1,1);
x=zeros(n+1,1);
switch n
   case 0
     w(1)=1;
     x(1)=0;
   case 1
      w(1)=1;
      w(2)=1;
      x(1)=-1;
      x(2)=1;

   case 2
      w(1)=1/3;
      w(2)=4/3;
      w(3)=1/3;

      x(1)=-1;
      x(2)=0;
      x(3)=1;

   case 3
      w(1)=1/6;
      w(2)=5/6;
      w(3)=5/6;
      w(4)=1/6;

      x(1)=-1;
      x(2)=-1/sqrt(5);
      x(3)=1/sqrt(5);
      x(4)=1;

   case 4
      w(1)=1/10;
      w(2)=49/90;
      w(3)=32/45;
      w(4)=49/90;
      w(5)=1/10;

      x(1)=-1;
      x(2)=-sqrt(3/7);
      x(3)=0;
      x(4)=sqrt(3/7);
      x(5)=1;

   case 5
      w(1)=1/15;
      w(2)=(14-sqrt(7))/30;
      w(3)=(14+sqrt(7))/30;
      w(4)=(14+sqrt(7))/30;
      w(5)=(14-sqrt(7))/30;
      w(6)=1/15;

      x(1)=-1;
      x(2)=-sqrt(1/3+2*sqrt(7)/21);
      x(3)=-sqrt(1/3-2*sqrt(7)/21);
      x(4)=sqrt(1/3-2*sqrt(7)/21);
      x(5)=sqrt(1/3+2*sqrt(7)/21);
      x(6)=1;

   case 6
      w(1)=1/21;
      w(2)=(124-7*sqrt(15))/350;
      w(3)=(124+7*sqrt(15))/350;
      w(4)=256/525;
      w(5)=(124+7*sqrt(15))/350;
      w(6)=(124-7*sqrt(15))/350;
      w(7)=1/21;

      x(1)=-1;
      x(2)=-sqrt(5/11+2/11*sqrt(5/3));
      x(3)=-sqrt(5/11-2/11*sqrt(5/3));
      x(4)=0;
      x(5)=sqrt(5/11-2/11*sqrt(5/3));
      x(6)=sqrt(5/11+2/11*sqrt(5/3));
      x(7)=1;

   otherwise
      error('ord1+ord2 zu gro� oder zu klein')

end
end