text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # Quadratically constrained programs

# A simple quadratically constrained program based on an [example from Gurobi](https://www.gurobi.com/documentation/9.0/examples/qcp_c_c.html).

using JuMP
import Ipopt
import Test

function example_qcp(; verbose = true)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, x)
    @variable(model, y >= 0)
    @variable(model, z >= 0)
    @objective(model, Max, x)
    @constraint(model, x + y + z == 1)
    @constraint(model, x * x + y * y - z * z <= 0)
    @constraint(model, x * x - y * z <= 0)
    optimize!(model)
    if verbose
        print(model)
        println("Objective value: ", objective_value(model))
        println("x = ", value(x))
        println("y = ", value(y))
    end
    Test.@test termination_status(model) == MOI.LOCALLY_SOLVED
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT
    Test.@test objective_value(model) ≈ 0.32699 atol = 1e-5
    Test.@test value(x) ≈ 0.32699 atol = 1e-5
    Test.@test value(y) ≈ 0.25707 atol = 1e-5
    return
end

example_qcp()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module FlattenM
    using LoopVectorization

    mutable struct Flatten
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple

        output::Array{Float32}
        δ::Array{Float32}

        function Flatten(;input_shape::Tuple)
            s = 1
            for i in input_shape
                s *= i
            end
            output_shape = (s,)

            new(save_Flatten, load_Flatten, activate_Flatten, init_Flatten, update_Flatten, input_shape, output_shape)
        end
    end

    function init_Flatten(layer::Flatten, mini_batch::Int64)
        layer.output = zeros(Float32, layer.output_shape..., mini_batch)
        layer.δ = zeros(Float32, layer.input_shape..., mini_batch)
    end

    function activate_Flatten(layer::Flatten, input::Array{Float32})
        @avxt for i in eachindex(input)
            layer.output[i] = input[i]
        end
    end

    function update_Flatten(layer::Flatten, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=1)
        @avxt for i in eachindex(δₗ₊₁)
            layer.δ[i] = δₗ₊₁[i]
        end
    end

    function save_Flatten(layer::Flatten, file::Any, id::String)
        write(file, id, "Flatten")
        write(file, id*"input_shape", collect(layer.input_shape))
    end

    function load_Flatten(layer::Flatten, file::Any, id::String)
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        return (input_shape=input_shape,)
    end
end

ext: jl


text:
#=
##
## TODO:	mtie has great performance for phase data but extremely worse performance for frequency data
##			Is this due to branch misprediction?
##
=#
"""
mtie(data, rate; [frequency=false], [overlapping=true], [taus=Octave])
Calculates the maximal time interval error

# parameters:
* `<data>`:			The data array to calculate the deviation from either as as phases or frequencies.
* `<rate>`:			The rate of the data given.
* `[frequency]`:	True if `data` contains frequency data otherwise (default) phase data is assumed.
* `[overlapping]`:	True (default) to calculate overlapping deviation, false otherwise.
* `[taus]`:			Taus to calculate the deviation at. This can either be an AllanTauDescriptor type `AllTaus, Decadade, HalfDecade, Octave, HalfOctave, QuarterOctave`, an array of taus to calculate at, a float number to build a custom log-scale on or an integer to build a specific number of log spaced points.

# returns: named tupple (tau, deviation, error, count)
* `tau`:		Taus which where used.
* `deviation`:	Deviations calculated.
* `error`:		Respective errors.
* `count`:		Number of contributing terms for each deviation.
"""
function mtie(
		data::Array{T, 1},
		rate::AbstractFloat;
		frequency::Bool = false,
		overlapping::Bool = true,
		taus::Union{Type{U}, Integer, AbstractFloat, Array{Float64}} = 192) where {T, U <: AllanTauDescriptor}

	#frequency to phase calculation
	if frequency
		data = frequencytophase(data, rate)
	end

	n = length(data)
	if n < 2
		error("Length for `data` in mtie must be at least 2 or greater")
	end

	if !overlapping #warn for consecutive execution
		@warn "It is highly unusual to use the mtie in the non overlapping form. Do not use this for definite interpretation or publication."
	end

	#tau calculations
	m = taudescription_to_m(taus, rate, n)

	dev = zeros(T, length(m)) #mtie
	deverr = zeros(T, length(m)) #mtie error
	devcount = zeros(Int, length(m)) #sum term count

	mStride = 1 #overlapping - can be overwritten in loop for consecutive
	@inbounds for (index, τ) in enumerate(m)

		if !overlapping #overwrite stride for consecutive operation
			mStride = τ
		end

		#mtie: https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nistspecialpublication1065.pdf
		terms = n - τ
		if terms < 2
			break
		end

		submin = data[1]
		submax = data[1]
		for j = 1:(1 + τ)
			if data[j] < submin
				submin = data[j]
			elseif data[j] > submax
				submax = data[j]
			end
		end
		delta = submax - submin
		maximumv = delta
		for i = (1 + mStride):(mStride):(n - τ)

			#max pipe
			if data[i - mStride] == submax #rolling max-pipe is obsolete
				submax = data[i]
				for j = i:(i + τ)
					if data[j] > submax
						submax = data[j]
					end
				end
				delta = submax - submin
			elseif data[i + τ] > submax #if new element is bigger than the old one
				submax = data[i + τ]
				delta = submax - submin
			end

			#min pipe
			if data[i - mStride] == submin #rolling min-pipe is obsolete
				submin = data[i]
				for j = i:(i + τ)
					if data[j] < submin
						submin = data[j]
					end
				end
				delta = submax - submin
			elseif data[i + τ] < submin #if new element is smaller than the old one
				submin = data[i + τ]
				delta = submax - submin
			end

			#comparer
			if delta > maximumv
				maximumv = delta
			end

		end

		dev[index] = maximumv
		deverr[index] = dev[index] / sqrt(terms)
		devcount[index] = terms
	end

	selector = devcount .> 1 #select only entries, where 2 or more terms contributed to the deviation
	(tau = m[selector] ./ rate, deviation = dev[selector], error = deverr[selector], count = devcount[selector])
end


text:
function phantom(N::Int64, n::Int64) #N:이미지 크기, n:팬텀의 타입
    #Shepp-Logan phantom n=1
    #            A     a       b        x0     y0        phi
    #            k1    k2      k3       k4     k5        k6
    # ---------------------------------------------------------
    sheplog = [ 1.0   0.69    0.92     0.0    0.0       0.0   ;
               -0.98  0.6624  0.8740   0.0   -0.0184    0.0   ;
                0.01  0.2100  0.2500   0.0    0.35      0.0   ;
                0.01  0.0460  0.0460   0.0    0.1       0.0   ;
                0.01  0.0460  0.0460   0.0   -0.1       0.0   ;
                0.01  0.0460  0.0230  -0.08  -0.605     0.0   ;
                0.01  0.0230  0.0230   0.0   -0.606     0.0   ;
                0.01  0.0230  0.0460   0.06  -0.605     0.0   ;
               -0.02  0.1100  0.3100   0.22   0.0      -0.314 ;
               -0.02  0.1600  0.4100  -0.22   0.0       0.314 ]

    #modified-Shepp_Logan phantom n=2
    #              A    a       b        x0     y0       phi
    #              k1   k2      k3       k4     k5       k6
    #----------------------------------------------------------
    modishep = [  1.0  0.69    0.92     0.0    0.0      0.0   ;
                 -0.8  0.6624  0.8740   0.0   -0.0184   0.0   ;
                  0.1  0.2100  0.2500   0.0    0.35     0.0   ;
                  0.1  0.0460  0.0460   0.0    0.1      0.0   ;
                  0.1  0.0460  0.0460   0.0   -0.1      0.0   ;
                  0.1  0.0460  0.0230  -0.08  -0.605    0.0   ;
                  0.1  0.0230  0.0230   0.0   -0.606    0.0   ;
                  0.1  0.0230  0.0460   0.06  -0.605    0.0   ;
                 -0.2  0.1100  0.3100   0.22   0.0     -0.314 ;
                 -0.2  0.1600  0.4100  -0.22   0.0      0.314 ]

    if n==1
      sl=sheplog
    elseif n==2
      sl=modishep
    end

    P = zeros(Float64, N, N)
    x = LinRange(-1,1,N)
    y = - x

    for k=1:10

      if k < 9 #회전이 없는 타원
        for i=1:N, j=1:N #i, j 헷갈리지 않게 주의
          if ( (x[j]-sl[k,4])/sl[k,2] )^2 + ( (y[i]-sl[k,5])/sl[k,3] )^2 < 1
            P[i, j]+= sl[k,1]
          end
        end

      else #회전이 있는 타원
        for i=1:N, j=1:N
          if ( (cos(sl[k,6])*(x[j]-sl[k,4])+sin(sl[k,6])*(y[i]-sl[k,5]))/sl[k,2] )^2+( (sin(sl[k,6])*(x[j]-sl[k,4])-cos(sl[k,6])*(y[i]-sl[k,5]))/sl[k,3] )^2< 1
            P[i, j]+= sl[k,1]
          end
        end
      end

    end

    return P

end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module PilotResponse

export
    AbstractPilotResponse,

    initialize,
    update,

    updatePilotResponse,

    SimplePilotResponse,
    SimplePRResolutionAdvisory,
    SimplePRCommand,

    StochasticLinearPR,
    StochasticLinearPRCommand,
    StochasticLinearPRRA,

    LLDetPR,
    LLDetPRAircraft,
    LLDetPRCommand,
    LLDetPRRA,
    LLDetPROutput

import CommonInterfaces: initialize, update

using AbstractPilotResponseImpl

using SimplePilotResponseImpl

using StochasticLinearPRImpl

using LLDetPRImpl

end


text:
# Case: src[i]
function _srctoken!(index, ex::Expr)::Symbol
    @assert ex.head == :ref
    src, token = ex.args[1], ex.args[2]
    @assert isa(src, Symbol) & isa(token, Symbol)
    # TODO: What happens in case of src[i], src[j]?
    haskey(index, token) && throw(ArgumentError("Non-distinct row tokens."))
    index[token] = src
    index[src] = token
    return src
end

# Case: src
function _srctoken!(index, src::Symbol)::Symbol
    token = gensym()
    in(src, values(index)) && throw(ArgumentError("Non-distinct sources."))
    index[token] = src
    index[src] = token
    return src
end

macro with(srcs, verbs::Expr)
    # We use the same object as a map from tokens to sources and from sources
    # to tokens
    index = Dict{Symbol, Symbol}()
    # extract source/token from each source in srcs
    primary = Nullable{Symbol}()
    if isa(srcs, Expr)
        if srcs.head == :ref
            _srctoken!(index, srcs)
        elseif srcs.head == :tuple
            # foreach(x -> _srctoken!(index, x), srcs.args)
            # the primary source is the first tokenless source
            tokenless = false
            for arg in srcs.args
                src = _srctoken!(index, srcs)
                if !tokenless
                    primary, tokenless = Nullable(src), true
                end
            end
        else
            # TODO: informative error message
            error()
        end
    elseif isa(srcs, Symbol)
        primary = Nullable(_srctoken!(index, srcs))
    end

    # generate graph expression from `verbs`
    if verbs.head == :call
        graph_ex = _graph_ex([verbs], index, primary)
    elseif (verbs.head == :tuple) | (verbs.head == :block)
        graph_ex = _graph_ex(verbs.args, index, primary)
    else
        # TODO: informative error message
        error()
    end
    return quote
        q = $graph_ex
        # @show q
        _with(q, source(q))
    end
end

_with(q::Node, src) = Cursor(q)


text:
"""
    sub_tss(h, metaV, metaE; opt=false)

Compute the target set exploiting a subtractive approach.
"""
function sub_tss(h, metaV, metaE; opt=false, printme=false)
    Vsub = Dict{Int,Tuple{Int,Int}}() # node => (degree, threshold)
    Esub = Dict{Int,Tuple{Int,Int}}() # edge => (size, threshold)

    nCount=0
    eCount=0
    nDegreeSum=0
    nThSum=0
    eSizeSum=0
    eThSum=0

    avg_degree = 0
    avg_size = 0

	#init structures
    for v=1:nhv(h)
        push!(Vsub, v => (length(gethyperedges(h,v)), metaV[v]))
        avg_degree += length(gethyperedges(h,v))
        nDegreeSum += length(gethyperedges(h,v))
        nThSum += metaV[v]
        nCount += 1
    end

    for e=1:nhe(h)
        push!(Esub, e => (length(getvertices(h,e)), metaE[e]))
        avg_size += length(getvertices(h,e))
        eSizeSum += length(getvertices(h,e))
        eThSum += metaE[e]
        eCount += 1
    end

    avg_degree /= nhv(h)
    avg_size /= nhe(h)

    printme && println("avg degree ", avg_degree, " avg size ", avg_size)
    printme && println("Degree sum ", nDegreeSum, " nCount ", nCount, " nThSum ", nThSum)
    printme && println("Size sum ", eSizeSum, " eCount ", eCount, " eThSum ", eThSum)

    S = Int[] #the seed set
    LU = Set{Int}() #limbo nodes
    LHU = Set{Int}() #limbo edges
    U = deepcopy(Vsub) #clone vertices
    HU = deepcopy(Esub) #clone edges

	upsidedown = Set{Int}()

	empty_th_nodes  = Set{Int}()
	empty_th_edges  = Set{Int}()

	case1a = case1b = case2 = case3a = case3b = 0

	#while the graph is not empty (no nodes or edges)
    while length(U) != 0
        #CASE  1
		#if exits a node with threshold 0 it is self-activated and we remove from the current graph
        if length(empty_th_nodes) > 0
			#CASE 1 A

			while length(empty_th_nodes) != 0
				case1a += 1
				node_with_th_zero = pop!(empty_th_nodes)

				printme &&  println("Case1A ", node_with_th_zero)
                #remove minny from U and update U

				#if remove the node we have to reduce the threshold and the size of the containing edges
	            #updateHU!(h, node_with_th_zero, HU, empty_th_edges)
                eThSum -= local_updateHUThOnly!(h, node_with_th_zero, HU, empty_th_edges, LHU)
                if !(node_with_th_zero in LU)
                    nCount-=1
                    nDegreeSum-= U[node_with_th_zero][1]
                    eSizeSum -= local_updateHUSizeOnly!(h, node_with_th_zero, HU, LHU)
                else
                    delete!(LU, node_with_th_zero)
                end

                delete!(U, node_with_th_zero)

                if (node_with_th_zero in upsidedown)
                    delete!(upsidedown, node_with_th_zero)
                end
            end
        else
			#CASE 1 B
			#here we find an edge with threshold 0
            if length(empty_th_edges) > 0

				while length(empty_th_edges) != 0
					case1b += 1
					candidatedge = pop!(empty_th_edges)

					printme &&  println("Case1B ", candidatedge)
                    #remove the edge from the current graph

					#update the nodes inside the selected edge by decreasing threshold and degree
                    #updateU!(h, candidatedge, U, empty_th_nodes)
                    nThSum -= local_updateUThOnly!(h, candidatedge, U, empty_th_nodes,LU,upsidedown)
                    if !(candidatedge in LHU)
                        eCount-=1
                        eSizeSum-= HU[candidatedge][1]
                        nDegreeSum -=local_updateUDegreeOnly!(h, candidatedge, U, upsidedown, LU)
                    else
                        delete!(LHU, candidatedge)
                    end

                    delete!(HU, candidatedge)
				end
            else
				# CASE 2
				# NO EDGES AND NODES HAVE threshold equal to 0
				# We look for a node with the threshold value greater than its degree
                #upsidedown = filter(x -> x[2][2] > x[2][1], collect(U)) #returns a list with this feature
				if length(upsidedown) > 0
					# CASE 2 A
                	# One node has t(v) > d(v), for this reason in order to be activated
					# we have to insert in S because it cannot be activated by its neighbors
					# while length(upsidedown) != 0
						case2 += 1
						node = pop!(upsidedown)

						printme &&  println("Case2 ", node)

						push!(S, node)
						#update the current graph by removing this node, same case of CASE 1 A

                        nCount-=1
                        nDegreeSum-= U[node][1]
                        nThSum-= U[node][2]

                        delete!(U, node)
                        #updateHU!(h, node, HU, empty_th_edges)
                        eThSum-= local_updateHUThOnly!(h, node, HU, empty_th_edges,LHU)
                        if !(node in LU)
                            eSizeSum -=local_updateHUSizeOnly!(h, node, HU, LHU)
                        else
                            printme &&  printf("maybe there is a problem: removing node from LU!")
                            delete!(LU, node)
                        end
                else

                    avg_degree = nDegreeSum/ nCount
                    avg_nth = nThSum/nCount
                    #nscale= (avg_nth)/(avg_degree*(avg_degree+1))
                    nscale = op(avg_nth, avg_degree)
                    # nscale=(avg_nth)/(avg_degree)
                    avg_size   = eSizeSum/eCount
                    avg_eth = eThSum/eCount
                    escale = op(avg_eth, avg_size)
                    # escale= (avg_eth)/(avg_size*(avg_size+1))
                    # escale= (avg_eth)/(avg_size)

                    printme && println("avg nscale ",nscale," avg escale ",escale)

					#CASE 3
					#max node for a formula value
                    U_selected = filter(x -> !(x[1] in LU), collect(U))
					candidate_nodes = sort!(collect(U_selected), by= x-> op(x[2][2], x[2][1]), rev = true)

                    candidate_node = nothing
					for u in candidate_nodes
						ok = true
						for he in gethyperedges(h, u.first)
							if haskey(HU, he.first) && HU[he.first][1] == HU[he.first][2]
								ok = false
								break
							end
						end
						if ok
							candidate_node = u
							break
						end
					end
					#max edge for a formula value
                    HU_selected = filter(x -> !(x[1] in LHU), collect(HU))
                    candidate_edge = sort!(collect(HU_selected), by=x-> op(x[2][2], x[2][1]), rev = true)[1]

                 #   if ( candidate_node != nothing)
                if !isnothing(candidate_node) &&
                    op(candidate_node[2][2], candidate_node[2][1] * nscale) > op(candidate_edge[2][2], candidate_edge[2][1] * escale)

                # if ( candidate_node !== nothing &&
				# 	((candidate_node[2][2] / (candidate_node[2][1] * (candidate_node[2][1]+1))) * (nscale)  )  >
                #    ((candidate_edge[2][2] / (candidate_edge[2][1] * (candidate_edge[2][1]+1)))  * (escale)  )
                #    )
                # if ( candidate_node != nothing &&
                #         ((candidate_node[2][2] / (candidate_node[2][1] )) *  (nscale)  )  >
                #        ((candidate_edge[2][2] / (candidate_edge[2][1] ))  * (escale)  )
                #       )

						case3a += 1
						#println("Case 3A ", candidate_node[1])
                        #delete!(U, candidate_node[1])

                        nCount-=1
                        nDegreeSum-= U[candidate_node[1]][1]
                        nThSum-= U[candidate_node[1]][2]

                        push!(LU, candidate_node[1])
                        eSizeSum -= local_updateHUSizeOnly!(h, candidate_node[1], HU, LHU)
					else
						case3b += 1
						printme &&  println("Case 3B ", candidate_edge[1])
                        #delete!(HU, candidate_edge[1])

                        eCount-=1
                        eSizeSum-= HU[candidate_edge[1]][1]
                        eThSum-= HU[candidate_edge[1]][2]

                        push!(LHU, candidate_edge[1])
						nDegreeSum-=local_updateUDegreeOnly!(h, candidate_edge[1], U, upsidedown, LU)
					end
                end
            end
        end
    end

    actE = zeros(Bool, nhe(h))
    actV = zeros(Bool, nhv(h))

    for s in S
         actV[s] = true
    end
    simres = simulate!(h, actV, actE, metaV, metaE; printme = false)

    if simres.actvs != nhv(h)
        println("There's a problem here ", simres.actvs, " ", nhv(h))
		i = 1
		for v in simres.actvsnodes
			v == false && println(i," ",v, " ", length(gethyperedges(h,i)), " ", metaV[i]," ", gethyperedges(h,i))
			i+=1
		end

    end

    if opt
        new_seed = optimize_seed_set(h, Set(keys(S)), metaV, metaE)
        return length(new_seed), case1a, case1b, case2, case3a, case3b
    end

    length(S), case1a, case1b, case2, case3a, case3b
end

#
# HELPER FUNCTIONS
#

function op(a, b)
    return a / (b * (b+1))
end

function local_updateUDegreeOnly!(h, e, U, list, LU)
    s = 0
    for w in getvertices(h,e)
        if haskey(U, w.first) && U[w.first][1] !=0
            if !(w.first in LU)
                s+=1
            end
            push!(U, w.first => (max(0, U[w.first][1]-1), U[w.first][2]))
            if U[w.first][1] <  U[w.first][2]
                if !(w.first in LU)
                     push!(list, w.first)
                 end
            end
        end
    end
    return s
end

function local_updateHUSizeOnly!(h, v, Esub, LHU)
    s = 0
    for he in gethyperedges(h, v)
        if haskey(Esub, he.first) &&  Esub[he.first][1] != 0
            if !(he.first in LHU)
                s += 1
            end
            push!(Esub, he.first => (max(0, Esub[he.first][1]-1), Esub[he.first][2]))
			if Esub[he.first][1] <  Esub[he.first][2]
				println("Oh my God!")
			end
        end
    end
    return s
end

#Update the graph by decreasing the size and threshold of edges, due to a new node is removed from the current graph
function local_updateHUThOnly!(h, v, Esub, list, LHU)
    t = 0
    for he in gethyperedges(h, v)
		#if the edge is in the current graph
        if haskey(Esub, he.first) && Esub[he.first][2] != 0
            if !(he.first in LHU)
                t += 1
            end
			#we decrese of 1 both size and threshold
            push!(Esub, he.first =>
						(Esub[he.first][1], max(0, Esub[he.first][2]-1)))
			if  Esub[he.first][2] == 0
				push!(list, he.first)
			end
        end
    end
    return t
end

#Update the graph by decreasing the degree and threshold of nodes, due to a removing edge
function local_updateUThOnly!(h, e, U, list, LU, list2)
    t = 0
    for w in getvertices(h,e)
        if haskey(U, w.first) && U[w.first][2] != 0
            if !(w.first in LU)
                t+=1
            end
            push!(U, w.first =>
					(U[w.first][1], max(0, U[w.first][2]-1)))
			if U[w.first][2] == 0
				push!(list, w.first)
            else
                if w.first in list2 && U[w.first][1] >=  U[w.first][2]
                    delete!(list2, w.first)
                end
            end
        end
    end
    return t
end


text:
function writefullfile{T<:Geometry}(geoms::Vector{T},lstfile::String=tempname())
    # Write full list file
    open(f->println(f,listfull(geoms)),lstfile,"w")
    return lstfile
end

function writefilestodir{T<:Geometry}(geoms::Vector{T},dir::String=mktempdir())
    # Create folder if it does not exist
    !ispath(dir)&&mkdir(dir)
    # Write one shape file per unique shape
    map(unique(shapes(geoms))) do shape
        open(f->println(f,getfile(shape)),joinpath(dir,shape.name),"w")
    end
    # Write list file
    lstfile = joinpath(dir,"list.lst")
    open(f->println(f,list(geoms)),lstfile,"w")
    return lstfile
end

function runandparse_ffs(cmd)
    # Run command
    stdout, process = open(cmd,"r")
    wait(process)
    res = split(strip(readstring(stdout)),NL)
    # Abort and show log if error occured
    if process.exitcode!=0
        error("FastCap2 error $(process.exitcode)",NL,replace(replace(join(res,"§"),"§§",NL),"§",""))
    end
    # Parse output (FastCap2 or FasterCap)
    nc    = parse(Int,res[end]) # number of conductors
    names = Vector{String}(res[end-nc:end-1]) # conductor names
    mem   = Int64(reinterpret(UInt32,parse(Int32,res[end-nc-1]))) # memory used by solver, correct the overflow bug
    time  = parse(Float64,res[end-nc-2]) # time used by solver
    capm  = readcsv(IOBuffer(join(res[end-2nc-2:end-nc-3],NL)),Float64,dims=(nc,nc))
    log   = replace(replace(join(res[1:end-2nc-3],"§"),"§§",NL),"§","")
    # Check if matrix is strictly diagonally dominant
    rows_bad = sum(capm,2).<0
    if any(rows_bad)
        rows = find(rows_bad)
        warn("Capacitance matrix is not stricly diagonally dominant due to row$(length(rows)>1?"s":"") $(join(rows,", "," and ")).")
    end
    # Return result
    return FastCapResult(capm,names,mem,time,log)
end

function runandparse_wr(cmd)
    # Run command while timing
    time = @elapsed log = strip(readstring(cmd))
    # Parse output
    res    = split(log,NL)
    matidx = find(x->contains(x,"CAPACITANCE MATRIX"),res)[1]
    resmat = readdlm(IOBuffer(join(res[matidx+2:end],NL)))
    names  = Vector{String}(resmat[:,1])
    capm   = Matrix{Float64}(resmat[:,3:end])
    # Return result
    return FastCapResult(capm,names,0,time,log)
end

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
# Attributes

"""
    AbstractRecognizeAttribute

Abstract supertype for attribute objects that will be used in the recognize function
in order to specify which features of the provided model to look into.

"""

export ConstantObjective, DifferenceConstraints, AllIntegerVariables, AllIntegerConstraintBounds

abstract type AbstractRecognizeAttribute end

# TODO: Describe each attribute

struct DifferenceConstraints <: AbstractRecognizeAttribute end
struct AllIntegerVariables <: AbstractRecognizeAttribute end
struct ConstantObjective <: AbstractRecognizeAttribute end
struct AllIntegerConstraintBounds <: AbstractRecognizeAttribute end

"""
    AbstractSolveAttribute

Abstract supertype for attribute objects that will be used in the solve function
in order to specify which ways to solve the model.

"""

export ShortestPath

abstract type AbstractSolveAttribute end

struct ShortestPath <: AbstractSolveAttribute end

"""
    AbstractOptimizeAttribute

Abstract supertype for attribute objects that will be used in order to specify what to do when optimizing the model.

"""

abstract type AbstractOptimizeAttribute end

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
using ITensors
using Printf
using Profile
using StatProfilerHTML

function main()
  N = 100
  sites = spinOneSites(N)
  H = Heisenberg(sites)
  psi = randomMPS(sites)
  sw = Sweeps(5)
  maxdim!(sw,10,20,100,100,200)
  cutoff!(sw,1E-11)

  smallsw = Sweeps(1)
  maxdim!(smallsw,10)
  cutoff!(smallsw,1E-4)
  @profile dmrg(H,psi,smallsw,maxiter=2)

  Profile.clear()
  @profile dmrg(H,psi,sw,maxiter=2)
  statprofilehtml()

  #energy,psi = @time dmrg(H,psi,sw,maxiter=2)
  #@printf "Final energy = %.12f\n" energy
  #@show inner(psi,H,psi)
end
main()


text:
module ShootingStar

using LinearAlgebra
using DifferentialEquations
using FiniteDiff

include("twostage_shooting.jl")

export twostage_shooting

end # module


text:
#################### Mathematical Operators ####################
function cummean(x::AbstractArray)
    return mapslices(cummean, x, dims = 1)
end

function cummean(x::AbstractVector)
    y = similar(x, Float64)
    xs = 0.0
    fill!(y, xs)

    c = 0
    for i in 1:length(x)
        if !ismissing(x[i])
            c += 1
            xs += x[i]
        end
        y[i] = xs / c
    end
    return y
end

## Csorgo S and Faraway JJ. The exact and asymptotic distributions of the
## Cramer-von Mises statistic. Journal of the Royal Statistical Society,
## Series B, 58: 221-234, 1996.
function pcramer(q::Real)
    p = 0.0
    for k in 0:3
        c1 = 4.0 * k + 1.0
        c2 = c1^2 / (16.0 * q)
        p += gamma(k + 0.5) / factorial(k) * sqrt(c1) * exp(-c2) * besselk(0.25, c2)
    end
    return p / (pi^1.5 * sqrt(q))
end

# This sorting stack was sourced from https://rosettacode.org/wiki/Natural_sorting#Julia
splitbynum(x) = split(x, r"(?<=\D)(?=\d)|(?<=\d)(?=\D)")
numstringtonum(arr) = [(n = tryparse(Float32, e)) != nothing ? n : e for e in arr]
function natural(x, y)
    xarr = numstringtonum(splitbynum(x))
    yarr = numstringtonum(splitbynum(y))
    for i in 1:min(length(xarr), length(yarr))
        if typeof(xarr[i]) != typeof(yarr[i])
            a = string(xarr[i]); b = string(yarr[i])
        else
             a = xarr[i]; b = yarr[i]
        end
        if a == b
            continue
        else
            return a < b
        end
    end
    return length(xarr) < length(yarr)
end

function _dict2namedtuple(d::Dict)
    t_keys = ntuple(x -> Symbol(collect(keys(d))[x]), length(keys(d)))
    t_vals = ntuple(x -> collect(values(d))[x], length(values(d)))
    return NamedTuple{t_keys}(t_vals)
end

function _namedtuple2dict(n::NamedTuple)
    return Dict(pairs(n))
end

function _ntdictmerge(n::NamedTuple, args::NamedTuple...)
    ndict = _namedtuple2dict(n)
    ndict2 = map(x -> _namedtuple2dict(x), args)
    return _dict2namedtuple(merge(ndict, ndict2...))
end

"""
    cskip(x)

Wrapper for `collect(skipmissing(x))`.
"""
cskip(x) = collect(skipmissing(x))


text:
@testset "pass" begin

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(mod, "SomeFunction", ft)

    bb = BasicBlock(fn, "SomeBasicBlock")
    position!(builder, bb)

    ret!(builder)

    verify(mod)

    # module pass

    function runOnModule(cur_mod::LLVM.Module)
        @test cur_mod == mod
        return false
    end

    let pass = ModulePass("SomeModulePass", runOnModule)

        ModulePassManager() do mpm
            add!(mpm, pass)
            run!(mpm, mod)
        end
    end

    # function pass

    function runOnFunction(cur_fn::LLVM.Function)
        @test cur_fn == fn
        return false
    end

    let pass = FunctionPass("SomeFunctionPass", runOnFunction)

        FunctionPassManager(mod) do fpm
            add!(fpm, pass)
            run!(fpm, fn)
        end
    end

    # basic block pass

    function runOnBasicBlock(cur_bb::BasicBlock)
        @test cur_bb == bb
        return false
    end

    if VERSION < v"1.4.0-DEV.589"
        let pass = BasicBlockPass("SomeBasicBlockPass", runOnBasicBlock)

            FunctionPassManager(mod) do fpm
                add!(fpm, pass)
                run!(fpm, fn)
            end
        end
    end
end
end
end

end


text:
import ..StaticNets

import ..StaticNets:NetModeErgmPml, ErgmDirBin0Rec0

Base.@kwdef struct  SdErgmPml <: SdErgm
    staticModel::NetModeErgmPml
    indTvPar :: BitArray{1}
    scoreScalingType::String = "FISH_D" # String that specifies the rescaling of the score. For a list of possible choices see function scalingMatGas
    options::SortedDict{Any, Any} = SortedDict("integrated"=>false, "initMeth"=>"estimateFirstObs")
end
export SdErgmPml

SdErgmPml(staticModel::NetModeErgmPml) = SdErgmPml(staticModel=staticModel, indTvPar = trues(staticModel.nErgmPar))

SdErgmPml(ergmTermsString::String, isDirected::Bool) = SdErgmPml(NetModeErgmPml(ergmTermsString, isDirected))

function Base.getproperty(x::SdErgmPml, p::Symbol)
    if p in fieldnames(typeof(x))
        return Base.getfield(x, p)
    else
        return Base.getproperty(x.staticModel, p)
    end
end

function name(x::SdErgmPml)
    if isempty(x.options)
        optString = ""
    else
        optString = ", " * reduce(*,["$k = $v, " for (k,v) in x.options])[1:end-2]
    end

    if x.staticModel.ergmTermsString == ergm_term_string(ErgmDirBin0Rec0())
        return  "SdErgmDirBin0Rec0_pmle($(x.indTvPar), scal = $(x.scoreScalingType)$optString)"
    else
        return  "SdErgmPML($(x.indTvPar), scal = $(x.scoreScalingType)$optString)"
    end
end

stats_from_mat(model::SdErgmPml, A ::Matrix{<:Real}) = StaticNets.change_stats(model.staticModel, A)

function static_estimate(model::SdErgmPml, A_T)
    staticPars = ErgmRcall.get_one_mple(A_T, model.staticModel.ergmTermsString)
    return staticPars
end

function target_function_t(model::SdErgmPml, obs_t, N, par)

    changeStat, response, weights = ErgmRcall.decomposeMPLEmatrix(obs_t)

    pll = StaticNets.obj_fun( par, changeStat, response, weights)

    return pll
end

function target_function_t_grad(model::SdErgmPml, obs_t, N, f_t)

    changeStat, response, weights = ErgmRcall.decomposeMPLEmatrix(obs_t)

    return StaticNets.grad_obj_fun(f_t, changeStat, response, weights)
end

function target_function_t_hess(model::SdErgmPml, obs_t, N, f_t)

    changeStat, response, weights = ErgmRcall.decomposeMPLEmatrix(obs_t)

    return StaticNets.hessian_obj_fun(f_t, changeStat, weights)

end

function target_function_t_fisher(model::SdErgmPml, obs_t, N,  f_t)

    changeStat, response, weights = ErgmRcall.decomposeMPLEmatrix(obs_t)

    return  StaticNets.fisher_info_obj_fun(f_t, changeStat, weights)

end

function reference_model(model::SdErgmPml)
    if model.staticModel.ergmTermsString == ergm_term_string(ErgmDirBin0Rec0())
        return SdErgmDirBin0Rec0_mle(scoreScalingType=model.scoreScalingType, indTvPar = model.indTvPar)
    else
        return model
    end
end

#region list example models
model_edge_gwd(decay_par) = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges + gwidegree(decay = $decay_par, fixed = TRUE, cutoff=10) + gwodegree(decay = $decay_par, fixed = TRUE, cutoff=10)", true), indTvPar = [true, true, true], scoreScalingType="FISH_D")

model_edge_mutual = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges + mutual", true), indTvPar = [true, true], scoreScalingType="FISH_D")

model_edge_gwesp = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges + gwesp(decay = 0.25, fixed = TRUE, cutoff=10)", true), indTvPar = [true, true], scoreScalingType="FISH_D")

model_edge_mutual_gwesp = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges + mutual + gwesp(decay = 0.25, fixed = TRUE, cutoff=10)", true), indTvPar = [true, true, true], scoreScalingType="FISH_D")

model_edge_mutual_gwd(decay_par) = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges + mutual + gwidegree(decay = $decay_par, fixed = TRUE, cutoff=10) + gwodegree(decay = $decay_par, fixed = TRUE, cutoff=10)", true), indTvPar = [true, true, true, true], scoreScalingType="FISH_D")

model_rec_p_star = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges + mutual ", true), indTvPar = [true, true], scoreScalingType="FISH_D")
#endregion

function sample_time_var_par_from_dgp(model::SdErgmPml, dgpType, N, T; minVals = -Inf * ones(1), maxVals = Inf * ones(1), meanVals = nothing, sigma = [0.01], B = [0.95], A=[0.01], indTvPar=trues(number_ergm_par(model)), maxAttempts = 25000, plotFlag = false)

    @debug "[sample_time_var_par_from_dgp][init][$( (;model, dgpType, N, T, minVals, maxVals)) ]"

    nErgmPar = model.staticModel.nErgmPar

    if isnothing(meanVals)
        meanVals = ( maxVals .+ minVals )./2
    end

    parDgpT = zeros(Real,nErgmPar,T)

    if dgpType=="AR"
        length(sigma) == 1 ? sigmaVec = sigma[1].*ones(nErgmPar) : sigmaVec = sigma
        length(B) == 1 ? BVec = B[1].*ones(nErgmPar) : BVec = B

        okSampleFlag = false
        for n=1:maxAttempts

            for p = 1:nErgmPar
                parDgpT[p, :] = dgpAR(meanVals[p], BVec[p], sigmaVec[p], T)
            end

            if all(parDgpT .< maxVals) & all(parDgpT .> minVals)
                okSampleFlag = true
                break
            end
        end

        okSampleFlag ? () : error("could not sample a path respecting boundaries last sample = $parDgpT")

    end

    if plotFlag
     fig, ax = subplots(nErgmPar,1)
        for p in 1:nErgmPar
            ax[1,p].plot(parDgpT[p,:], "k")
        end
    end

     @debug "[sample_time_var_par_from_dgp][end]"
    return parDgpT
end


text:
using AMDGPUnative
const AS = AMDGPUnative.AS
using LLVM, LLVM.Interop
using HSARuntime
using InteractiveUtils
using SpecialFunctions
using Test

agent_name = HSARuntime.get_name(get_default_agent())
agent_isa = get_first_isa(get_default_agent())
@info "Testing using device $agent_name with ISA $agent_isa"

@testset "AMDGPUnative" begin

include("util.jl")

include("base.jl")
include("pointer.jl")
# TODO: include("codegen.jl")

if AMDGPUnative.configured
    @test length(get_agents()) > 0
    if length(get_agents()) > 0
        # The original test :)
        @testset "Vector Addition Kernel" begin
            function vadd(a,b,c)
                i = threadIdx().x
                c[i] = a[i] + b[i]
                sync_workgroup()
                return nothing
            end

            dims = (8,)
            a = round.(rand(Float32, dims) * 100)
            b = round.(rand(Float32, dims) * 100)

            d_a = HSAArray(a)
            d_b = HSAArray(b)
            d_c = similar(d_a)
            len = prod(dims)

            @debug begin
                @show d_a.handle
                @show d_b.handle
                @show d_c.handle

                aspace = AMDGPUnative.AS.Global
                arrdims = ndims(a)
                arrT = ROCDeviceArray{Float32,arrdims,aspace}
                @debug "LLVM IR"
                AMDGPUnative.code_llvm(vadd, Tuple{arrT,arrT,arrT}; kernel=true)
                @debug "GCN Device Code"
                AMDGPUnative.code_gcn(vadd, Tuple{arrT,arrT,arrT}; kernel=true)
                ""
            end

            @roc groupsize=len gridsize=len vadd(d_a, d_b, d_c)
            @debug begin
                @show d_a
                @show d_b
                @show d_c
                ""
            end
            c = Array(d_c)
            @test a+b ≈ c
        end

        include("device/synchronization.jl")
        if Base.libllvm_version >= v"7.0"
            include("device/math.jl")
        else
            @warn "Testing with LLVM 6; some tests will be disabled!"
            @test_skip "Math Intrinsics"
        end
        #include("device/codegen.jl")
        #include("device/execution.jl")
        #include("device/pointer.jl")
        #include("device/array.jl")
        #include("device/intrinsics.jl")
    end
else
    @warn("AMDGPUnative.jl has not been configured; skipping on-device tests.")
end

end


text:
"""
Internal function that wraps HTTP.get
"""
function _get_request(uri::String)
    resp = HTTP.get(uri)
    _check_api_limit(resp)
    _check_status_code(resp)
    return resp
end

"""
Internal function to check whether the response contains
    'Note' that indicates the API limit was reached
"""
function _check_api_limit(resp::HTTP.Messages.Response)
    content_type = HTTP.Messages.header(resp.headers, "Content-Type")
    if content_type == "application/json"
        body = _parse_response(resp, "json")
        if haskey(body, "Note")
            error("API limit exceeded")
        end
    end
end

"""
Internal function to check response status code
"""
function _check_status_code(resp::HTTP.Messages.Response)
    if resp.status != 200
        error("Expected status code 200 but received $resp.status")
    end
end

"""
Internal function that parses the response
"""
function _parse_response(data, datatype::Union{String, Nothing})
    if datatype == "csv"
        return readdlm(data.body, ',', header=true)
    elseif datatype == "json"
        body = copy(data.body)  # TODO: re-write to avoid copying
        return JSON.Parser.parse(String(body))
    else
        raw = readdlm(data.body, ',', header=true)
        return AlphaVantageResponse(raw)
    end
end

"""
Internal function to build the request URI
"""
function _build_uri(scheme::String, host::String, path::String, params::Dict)
    @argcheck in(scheme, ["https", "http"])
    params = filter(p -> !isnothing(p.second), params)
    params = collect(pairs(params))
    query = join(map(p -> "$(p.first)=$(p.second)", params), "&")
    return "$(scheme)://$(host)/$(path)?$(query)"
end

_parser(p, datatype) = p
function _parser(p::AbstractString, datatype)
    if p == "default"
        return x -> _parse_response(x, datatype)
    else
        # if parser is unknown, then return raw data
        return identity
    end
end


text:
"""
wrapper to clearly indicate that something should be treated as a Callable
"""
struct Callable{F}
  f::F
end

(callable::Callable)(args...; kwargs...) = callable.f(args...;kwargs...)


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

using Markdown
using InteractiveUtils

# ╔═╡ 9b58babd-410d-41ef-b702-d2e0a9eda454
# Activate the Julia environment of the project
begin
	using Pkg
	Pkg.activate("..\\")
end

# ╔═╡ 47e86052-d08f-4c34-8b65-298844ffa8ba
# Load the necessary packages
begin
	using ArchGDAL
	using GeoStats
	using Plots
end

# ╔═╡ 9b6db6b0-9a07-11ec-112e-5904bc97f187
# Load the necessary modules
begin
	include("..\\src\\Data Gathering\\GroundData.jl")
	include("..\\src\\Data Gathering\\SatelliteData.jl")

	include("..\\src\\Analysis\\Aquifers.jl")
	include("..\\src\\Analysis\\Lakes.jl")
	include("..\\src\\Analysis\\Noises.jl")
	include("..\\src\\Analysis\\Plumes.jl")
	include("..\\src\\Analysis\\Sediments.jl")
end

# ╔═╡ 1769fb6b-c6a6-4ab4-91eb-8a1f245b7d21
# Loading of modules and packages

# ╔═╡ a12a30b3-8db3-4eaf-803f-b86aa4946d2a
# Aliasing of the modules to shorten the code
begin
	const agd = ArchGDAL

 # Data gathering
	const grdt = GroundData
	const stdt = SatelliteData
end

# ╔═╡ 810f6349-b471-44df-a60c-cfc0602cb3f9
# Code testing

# ╔═╡ cee41885-9767-4bbc-aeca-fdcdc72add39
# Data Gathering testing

# ╔═╡ 3996736f-6216-4c5c-aa7e-844133995648
# Data from the measurement stations in various regions (Alto Adige, Friuli Venezia Giulia, Lombardia, Trentino and Veneto)
resmt, mresmt = grdt.getGroundData(:METEO, grdt.AA, grdt.FVG, grdt.L, grdt.T, grdt.V)

# ╔═╡ 81190d46-872b-490f-95b7-17fcd02dd2e1
unique(resmt.parameter)

# ╔═╡ 0447a923-fcc3-4425-b17d-6a5bd4c38d83
begin
 # All relative humidity measurements
	measurements = resmt[ resmt.parameter .== "Umidità relativa", : ]

 # All the measurement's values
	values = ( val = measurements.value, )

 # All the coordinates of the measurements
	coord = Tuple{Float64, Float64}[ (c[1], c[2]) for c in eachrow(measurements[:, [:longitude, :latitude]]) ]
end

# ╔═╡ 28ddd2f5-f53b-46c0-993c-d05b10fca309

#= Based on:
	https://juliaearth.github.io/GeoStats.jl/stable/
=#
begin
 # Georeference data
	D = georef(values, coord)

 # Estimation domain
	G = CartesianGrid(100, 100)

 # Estimation problem
	problem = EstimationProblem(D, G, :val)

 # Solver from the list of solvers
	solver = Kriging( :val => ( variogram = GaussianVariogram(range=35.0), ) )

 # Solving problem
	solution = solve(problem, solver)

 # Solution plot
	contourf(solution, clabels=true)
end

# ╔═╡ d2ca35b9-1912-4329-8ddd-2020b5852ae9
# Data Analysis testing

# ╔═╡ 9e8d41ce-4582-44ff-890c-7d89cdce9984
# Initial Setup

# ╔═╡ 4450aaa2-f049-44df-9bd7-43c13fdaaa3e
# Digital Terrain Model, the main raster for the analysis, reppresents Veneto
begin
	dtm_file = "..\\resources\\Analysis data\\DTM_32.tiff"
	dtm = agd.readraster(dtm_file)
end

# ╔═╡ b7699c02-99a1-4010-afc8-816da12e2c25
# Create a source point to test the analysis functions
begin
 # Coordinates of the source
    # lat, lon = (11.930065824163105, 45.425861311724816) # WGS84
	lat, lon = (726454.9302346368, 5.025993899219433e6)
 # Path to the source shapefile
	source_dir = "..\\resources\\Analysis data\\source_shapefile"
 # Directory holding the `source` files
	!isdir(source_dir) && mkdir(source_dir)
 # Creation of a shapefile containing the point if not already existing
	if "source_32.shp" ∉ readdir(source_dir)
		agd.create(
			source_dir*"\\source_32.shp",
			driver = agd.getdriver("ESRI Shapefile")
		) do ds
			agd.createlayer(
				geom = agd.wkbPoint,
				# spatialref = agd.importEPSG(4326) # WGS84
            	spatialref = agd.importEPSG(32632)
			) do layer
        		agd.createfeature(layer) do feature
            		agd.setgeom!(feature, agd.createpoint(lat, lon))
        		end
    			agd.copy(layer, dataset=ds)
			end
		end
	end
	source_file = source_dir*"\\source_32.shp"
end

# ╔═╡ e960080e-05a5-480a-b27f-4149153300e4
# Analysis functions execution

# ╔═╡ 4b7f821c-c019-4b41-a58a-96edee0566c1
# Pollutants diffusion in an aquifer

# ╔═╡ 75eaae0c-55e9-4de5-9d31-f405f60f5863
#=
contaminants: (vettore di String)
	nomi della sostanze da usare per estrarre i valori "c_henry", "koc_kd" dal database.
concentrations: (vettore di Float64)
	concentrazioni delle sostanze.
aquifer_depth: (Float64)
	profondità della falda.
acquifer_flow_direction: (Int64)
	angolo del flusso dell'aqua nella falda in gradi (credo?)
mean_rainfall: (Float64)
	pioggia media
texture: (String)
	credo sia la tipologia di terreno in cui filtrano gli inquinanti, da usare per estrarre i valori "tot_por", "c_water_avg", "effective_infiltration",
	"por_eff", "grain" dal database.
orthogonal_extension: (Float64)
	estensione ortogonale?
soil_density: (Float64)
	densità del suolo
source_thickness: (Float64)
	spessore fonte?
darcy_velocity: (Float64)
	velocità di Darcy?
mixed_zone_depth: (Float64)
	profondità zoma mista?
decay_coeff: (Float64)
	coefficiente di decadimento?
=#
Aquifers.run_leaching(
	dem_file = dtm_file,
	source_file = source_file,
	contaminantCASNum="108-88-3",
	concentration = 100.0,
	tollerance = 2,
	aquifer_depth = 1000.0,
	acquifer_flow_direction = 180,
	mean_rainfall = 20.0,
	texture = "sand",
	resolution = 25.0,
	time = 10,
	orthogonal_width = 10.0,
	darcy_velocity = 0.000025,
	mixed_zone_depth = 1580.0,
	algorithm = :domenico,
	output_path = "..\\resources\\Analysis data\\Analysis results\\test_aquifer.tiff")

# ╔═╡ d8f5f781-6d0d-48e5-b7be-63c6c41fe94a
# Polutants dispersion in lakes

# ╔═╡ f62e94e9-9c0c-4ed7-b87d-c0d958257d59
#=
wind_direction: (Int64)
	direzione del vento come angolo in gradi.
pollutant_mass: (Float64)
	massa dell'inquinante.
flow_mean_speed: (Float64)
	velocità del flusso.
hours: (Int64)
	durata dell'analisi in ore.
fickian_x: (Float64)
	?
fickian_y: (Float64)
	?
λk: (Float64)
	?
=#
Lakes.run_lake(
	dem_file = dtm_file,
	source_file = source_file,
	wind_direction = 180,
	pollutant_mass = 2000.0,
	tollerance = 2,
	mean_flow_speed = 0.03,
	resolution = 25.0,
	hours = 10.0,
	fickian_x = 4.0,
	fickian_y = 3.0,
	output_path = "..\\resources\\Analysis data\\Analysis results\\test_lake.tiff"
)

# ╔═╡ 537c45ac-88c6-4a4b-96ec-a5021f69ee5a
# Noise pollution

# ╔═╡ 725222e8-960b-4618-8b3c-aa493caf7c16
#=
"A" sostituisce il raster delle impedenze (nel codice uso una matrice contenente solo 0)
temperatura atmosferica in kelvin
umidità relativa (credo come percentuale in decimale)
intensità del suono
frequenza del suono
=#
Noises.run_noise(
    dem_file = dtm_file,
    source_file = source_file,
    temperature_K = 293.15,
    relative_humidity = 0.2,
    intensity_dB = 110.0,
    frequency = 400.0
)

# ╔═╡ 83a9cd5e-041c-4997-ae8d-d1027f297b96
# Airborne pollutants dispersion from a stack

# ╔═╡ 7c9a6143-5331-42dd-9541-052d411bc284
#=
stability: (String)
	Informazioni sul tempo meteorologico?
	Usato per ottenere i campi "sf_ing", "sf_inal", "iur", "rfd_ing", "rfd_inal", "rfc" a questi corrispondono i messaggi "Slope Factor per ingestione",
	"Slope Factor per inalazione","Inhalation Unit Risk","Reference Dose per ingestione","Reference Dose per inalazione","Reference Concentration" dal database.
outdoor: (String)
	Usato insieme a "stability" per ottenere valori dal database.
wind_direction: (Int64)
	direzione del vento come angolo.
concentration: (Float64)
	concentrazione della sostanza.
wind_speed: (Float64)
	velocità del vento.
stack_height: (Float64)
	altezza del camino.
gas_speed: (Float64)
	velocità dei fumi?
stack_diameter: (Float64)
	diametro del camino.
smoke_temperature: (Float64)
	temperatura dei fumi.
temperature: (Float64)
	temperatura dell'ambiente.
=#
Plumes.run_plume(
	dem_file = dtm_file,
	source_file = source_file,
	stability = "a",
	outdoor = "c",
	concentration = 10000.0,
	tollerance = 2,
	resolution = 25.0,
	wind_direction = 135,
	wind_speed = 0.1,
	stack_height = 80.0,
	stack_diameter = 1.0,
	gas_velocity = 0.1,
	gas_temperature = 150.0,
	temperature = 18.0,
	output_path = "..\\resources\\Analysis test data\\Analysis results\\test_plumes.tiff"
)

# ╔═╡ 3173f593-c0af-4a14-86a6-93d91109edca
# Pollutants sedimentation

# ╔═╡ 928d1a74-86a8-488c-92fd-b660794db328
#=
mean_flow_speed: (Float64)
	velocità media della corrente.
mean_depth: (Float64)
	profondità media.
x_dispersion_coeff: (Float64)
	coefficiente di dispersione su x?
y_dispersion_coeff: (Float64)
	coefficiente di dispersione su y?
dredged_mass: (Float64)
	massa trasportata.
flow_direction: (Int64)
	direzione del flusso come angolo in gradi.
mean_sedimentation_velocity: (Float64)
	velocità di sedimentazione.
time: (Int64)
	tempo di inizio del modello.
time_intreval: (Int64)
	lunghezza di un'epoca.
current_oscillatory_amplitude: (Float64)
	ampiezza di oscillazione della corrente?
tide: (Int64)
	ciclo di marea in ore?.
=#
Sediments.run_sediment(
	dem_file = dtm_file,
	source_file = source_file,
	resolution = 25.0,
	mean_flow_speed = 0.03,
	mean_depth = 13.0,
	x_dispersion_coeff = 1.0,
	y_dispersion_coeff = 10.0,
	dredged_mass = 4.0,
	tollerance = 2,
	flow_direction = 180,
	mean_sedimentation_velocity = 0.0359,
	time = 100,
	time_intreval = 1,
	current_oscillatory_amplitude = 0.0,
	tide = 0,
	output_path = "..\\resources\\Analysis data\\Analysis results\\test_sediments.tiff"
)

# ╔═╡ Cell order:
# ╠═1769fb6b-c6a6-4ab4-91eb-8a1f245b7d21
# ╠═9b58babd-410d-41ef-b702-d2e0a9eda454
# ╠═47e86052-d08f-4c34-8b65-298844ffa8ba
# ╠═9b6db6b0-9a07-11ec-112e-5904bc97f187
# ╠═a12a30b3-8db3-4eaf-803f-b86aa4946d2a
# ╠═810f6349-b471-44df-a60c-cfc0602cb3f9
# ╠═cee41885-9767-4bbc-aeca-fdcdc72add39
# ╠═3996736f-6216-4c5c-aa7e-844133995648
# ╠═81190d46-872b-490f-95b7-17fcd02dd2e1
# ╠═0447a923-fcc3-4425-b17d-6a5bd4c38d83
# ╠═28ddd2f5-f53b-46c0-993c-d05b10fca309
# ╠═d2ca35b9-1912-4329-8ddd-2020b5852ae9
# ╠═9e8d41ce-4582-44ff-890c-7d89cdce9984
# ╠═4450aaa2-f049-44df-9bd7-43c13fdaaa3e
# ╠═b7699c02-99a1-4010-afc8-816da12e2c25
# ╠═e960080e-05a5-480a-b27f-4149153300e4
# ╠═4b7f821c-c019-4b41-a58a-96edee0566c1
# ╠═75eaae0c-55e9-4de5-9d31-f405f60f5863
# ╠═d8f5f781-6d0d-48e5-b7be-63c6c41fe94a
# ╠═f62e94e9-9c0c-4ed7-b87d-c0d958257d59
# ╠═537c45ac-88c6-4a4b-96ec-a5021f69ee5a
# ╠═725222e8-960b-4618-8b3c-aa493caf7c16
# ╠═83a9cd5e-041c-4997-ae8d-d1027f297b96
# ╠═7c9a6143-5331-42dd-9541-052d411bc284
# ╠═3173f593-c0af-4a14-86a6-93d91109edca
# ╠═928d1a74-86a8-488c-92fd-b660794db328

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020
# Copyright: Georgia Institute of Technology, 2020

module InvertibleNetworks

import Base.size, Base.getindex, Flux.glorot_uniform, Base.reverse, Base.reverse!
using LinearAlgebra, Random, NNlib, Flux, Statistics, Wavelets, Zygote

# For function signatur definition
import CUDA: CuArray

export clear_grad!, glorot_uniform, get_params, put_params!

# Getters for DenseConvDims fields
# (need to redefine here as they are not public methods in NNlib)
input_size(c::DenseConvDims) = c.I
kernel_size(c::DenseConvDims{N,K,C_in,C_out,S,P,D,F}) where {N,K,C_in,C_out,S,P,D,F} = K
channels_in(c::DenseConvDims{N,K,C_in,C_out,S,P,D,F}) where {N,K,C_in,C_out,S,P,D,F} = C_in
channels_out(c::DenseConvDims{N,K,C_in,C_out,S,P,D,F}) where {N,K,C_in,C_out,S,P,D,F} = C_out

# Utils
include("utils/parameter.jl")
include("utils/objective_functions.jl")
include("utils/dimensionality_operations.jl")
include("utils/activation_functions.jl")
include("utils/test_distributions.jl")
include("utils/neuralnet.jl")

# Single network layers (invertible and non-invertible)
include("conditional_layers/conditional_layer_residual_block.jl")
include("layers/layer_residual_block.jl")
include("layers/layer_flux_block.jl")
include("layers/layer_affine.jl")
include("layers/invertible_layer_actnorm.jl")
include("layers/invertible_layer_conv1x1.jl")
include("layers/invertible_layer_basic.jl")
include("layers/invertible_layer_irim.jl")
include("layers/invertible_layer_glow.jl")
include("layers/invertible_layer_hyperbolic.jl")
include("layers/invertible_layer_hint.jl")
include("layers/invertible_layer_slim_additive.jl")
include("layers/invertible_layer_slim_affine.jl")
include("layers/invertible_layer_slim_learned.jl")

# Invertible network architectures
include("networks/invertible_network_irim.jl")  # i-RIM: Putzky and Welling (2019)
include("networks/invertible_network_glow.jl")  # Glow: Dinh et al. (2017), Kingma and Dhariwal (2018)
include("networks/invertible_network_hyperbolic.jl")    # Hyperbolic: Lensink et al. (2019)

# Conditional layers and nets
include("conditional_layers/conditional_layer_hint.jl")
include("conditional_layers/conditional_layer_slim.jl")
include("networks/invertible_network_conditional_hint.jl")
include("networks/invertible_network_conditional_hint_multiscale.jl")

# gpu
include("utils/cuda.jl")

end


text:
using PanelMatrices
using Test
using StaticArrays
using StaticNumbers
using LinearAlgebra

@testset "set/get" begin
    y = reshape(collect(1:100), 10, 10)
    x = PanelMatrix(y)
    @test x[1,1] == y[1,1]
    @test all(x .== y)

    v = PanelMatrices.panel_view(x, 1, 1)
    @test all(v .== [1 11 21 31; 2 12 22 32; 3 13 23 33; 4 14 24 34])
    @test all(PanelMatrices.panel_view(x, 3, 1) .== [9 19 29 39; 10 20 30 40])
    @test all(PanelMatrices.get_panel(x, 3, 1) .== [9 19 29 39; 10 20 30 40])
    @test all(PanelMatrices.panel_view(x, 1, 3) .== [81 91; 82 92; 83 93; 84 94])
    @test all(PanelMatrices.get_panel(x, 1, 3) .== [81 91; 82 92; 83 93; 84 94])
    @test all(PanelMatrices.panel_view(x, 3, 3) .== [89 99; 90 100])
    @test all(PanelMatrices.get_panel(x, 3, 3) .== [89 99; 90 100])
    @test size(PanelMatrices.full_panel_view(x, 3, 3)) == (4, 4)

    p = PanelMatrices.get_full_panel(x, 1, 1)
    @test p === @SMatrix [1 11 21 31; 2 12 22 32; 3 13 23 33; 4 14 24 34]

    z = reshape(101:116, 4, 4)
    PanelMatrices.set_full_panel!(x, z, 2, 2)
    @test x[5,5] == z[1,1]
    @test all([x[i,j] for i ∈ 5:8, j ∈ 5:8] .== z)

    # TODO: Write tests that ensure that there is no allocation in critical routines.
end

@testset "matmul" begin
    for M = (5,8), N = (6,8), K = (7,8)
        A0 = reshape(collect(1.0:M*K), M, K)
        B0 = reshape(collect(1.0:K*N), K, N)
        C0 = A0*B0

        A = PanelMatrix(A0, static.((4,4)))
        B = PanelMatrix(B0, static.((4,4)))
        C = PanelMatrix{Float64}(undef, (M,N), static.((4,4)))

        mul!(C, A, B)
        @test C ≈ C0
    end
end


text:
#= MIT License

Copyright (c) 2020, 2021 Uwe Fechner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. =#

# data structures for the flight state and the flight log
# functions for creating a demo flight state, demo flight log, loading and saving flight logs
# function se() for reading the settings
# the parameter P is the number of points of the tether, equal to segments+1
# in addition helper functions for working with rotations

"""
    rot3d(ax, ay, az, bx, by, bz)

Calculate the rotation matrix that needs to be applied on the reference frame (ax, ay, az) to match
the reference frame (bx, by, bz).
All parameters must be 3-element vectors. Both refrence frames must be orthogonal,
all vectors must already be normalized.

Source: [TRIAD_Algorithm](http://en.wikipedia.org/wiki/User:Snietfeld/TRIAD_Algorithm)
"""
function rot3d(ax, ay, az, bx, by, bz)
    R_ai = hcat(ax, az, ay)
    R_bi = hcat(bx, bz, by)
    return R_bi * R_ai'
end

"""
    rot(pos_kite, pos_before, v_app)

Calculate the rotation matrix of the kite based on the position of the
last two tether particles and the apparent wind speed vector.
"""
function rot(pos_kite, pos_before, v_app)
    delta = pos_kite - pos_before
    @assert norm(delta) > 0.0 "Error in function rot() ! pos_kite must be not equal to pos_before. "
    c = -delta
    z = normalize(c)
    y = normalize(cross(-v_app, c))
    x = normalize(cross(y, c))
    rot = rot3d([0,-1.0,0], [1.0,0,0], [0,0,-1.0], z, y, x)
end

"""
    ground_dist(vec)

Calculate the ground distance of the kite from the groundstation based on the kite position (x,y,z, z up).
"""
function ground_dist(vec)
    sqrt(vec[1]^2 + vec[2]^2)
end

"""
    calc_elevation(vec)

Calculate the elevation angle in radian from the kite position.
"""
function calc_elevation(vec)
    atan(vec[3] / ground_dist(vec))
end

"""
    azimuth_east(vec)

Calculate the azimuth angle in radian from the kite position in ENU reference frame.
Zero east. Positive direction clockwise seen from above.
Valid range: -π .. π.
"""
function azimuth_east(vec)
    return -atan(vec[2], vec[1])
end

"""
    acos2(arg)

Calculate the acos of arg, but allow values slightly above one and below
minus one to avoid exceptions in case of rounding errors. Returns an
angle in radian.
"""
@inline function acos2(arg)
   arg2 = min(max(arg, -1.0), 1.0)
   acos(arg2)
end


text:
# # [Rising Thermal Bubble](@id EX-RTB-docs)
#
# In this example, we demonstrate the usage of the `ClimateMachine`
# [AtmosModel](@ref AtmosModel-docs) machinery to solve the fluid
# dynamics of a thermal perturbation in a neutrally stratified background state
# defined by its uniform potential temperature. We solve a flow in a
# [`FlatOrientation`](@ref LESConfig) (Box) configuration - this is
# representative of a large-eddy simulation. Several versions of the problem
# setup may be found in literature, but the general idea is to examine the
# vertical ascent of a thermal _bubble_ (we can interpret these as simple
# representation of convective updrafts).
#
# ## Description of experiment
# 1) Dry Rising Bubble (circular potential temperature perturbation)
# 2) Boundaries
#    - `Impenetrable(FreeSlip())` - no momentum flux, no mass flux through
#      walls.
#    - `Impermeable()` - non-porous walls, i.e. no diffusive fluxes through
#       walls.
#    - Laterally periodic
# 3) Domain - 2500m (horizontal) x 2500m (horizontal) x 2500m (vertical)
# 4) Resolution - 50m effective resolution
# 5) Total simulation time - 1000s
# 6) Mesh Aspect Ratio (Effective resolution) 1:1
# 7) Overrides defaults for
#    - CPU Initialisation
#    - Time integrator
#    - Sources
#    - Smagorinsky Coefficient
# 8) Default settings can be found in `src/Driver/<files>.jl`

#md # !!! note
#md #     This experiment setup assumes that you have installed the
#md #     `ClimateMachine` according to the instructions on the landing page.
#md #     We assume the users' familiarity with the conservative form of the
#md #     equations of motion for a compressible fluid (see the
#md #     [`AtmosModel`](@ref AtmosModel-docs) page).
#md #
#md #     The following topics are covered in this example
#md #     - Package requirements
#md #     - Defining a `model` subtype for the set of conservation equations
#md #     - Defining the initial conditions
#md #     - Applying boundary conditions
#md #     - Applying source terms
#md #     - Choosing a turbulence model
#md #     - Adding tracers to the model
#md #     - Choosing a time-integrator
#md #     - Choosing diagnostics (output) configurations
#md #
#md #     The following topics are not covered in this example
#md #     - Defining new boundary conditions
#md #     - Defining new turbulence models
#md #     - Building new time-integrators
#md #     - Adding diagnostic variables (beyond a standard pre-defined list of
#md #       variables)
#
# ## Boilerplate (Using Modules)
#
# #### [Skip Section](@ref init)
#
# Before setting up our experiment, we recognize that we need to import some
# pre-defined functions from other packages. Julia allows us to use existing
# modules (variable workspaces), or write our own to do so.  Complete
# documentation for the Julia module system can be found
# [here](https://docs.julialang.org/en/v1/manual/modules/#).

# We need to use the `ClimateMachine` module! This imports all functions
# specific to atmospheric and ocean flow modelling.  While we do not cover the
# ins-and-outs of the contents of each of these we provide brief descriptions
# of the utility of each of the loaded packages.
using ClimateMachine
ClimateMachine.init()

using ClimateMachine.Atmos
# - Required so that we inherit the appropriate model types for the large-eddy
#   simulation (LES) and global-circulation-model (GCM) configurations.
using ClimateMachine.ConfigTypes
# - Required so that we may define diagnostics configurations, e.g. choice of
#   file-writer, choice of output variable sets, output-frequency and directory,
using ClimateMachine.Diagnostics
# - Required so that we may define (or utilise existing functions) functions
#   that are `called-back` or executed at frequencies of either timesteps,
#   simulation-time, or wall-clock time.
using ClimateMachine.GenericCallbacks
# - Required so we load the appropriate functions for the time-integration
#   component. Contains ODESolver methods.
using ClimateMachine.ODESolvers
# - Required for utility of spatial filtering functions (e.g. positivity
#   preservation)
using ClimateMachine.Mesh.Filters
# - Required so functions for computation of temperature profiles.
using ClimateMachine.TemperatureProfiles
# - Required so functions for computation of moist thermodynamic quantities is
#   enabled.
using ClimateMachine.MoistThermodynamics
# - Required so we may access our variable arrays by a sensible naming
#   convention rather than by numerical array indices.
using ClimateMachine.VariableTemplates
# - Required so we may access planet parameters
#   ([CLIMAParameters](https://CliMA.github.io/CLIMAParameters.jl/latest/)
#   specific to this problem include the gas constant, specific heats,
#   mean-sea-level pressure, gravity and the Smagorinsky coefficient)

# In ClimateMachine we use `StaticArrays` for our variable arrays.
using StaticArrays
# We also use the `Test` package to help with unit tests and continuous
# integration systems to design sensible tests for our experiment to ensure new
# / modified blocks of code don't damage the fidelity of the physics. The test
# defined within this experiment is not a unit test for a specific
# subcomponent, but ensures time-integration of the defined problem conditions
# within a reasonable tolerance. Immediately useful macros and functions from
# this include `@test` and `@testset` which will allow us to define the testing
# parameter sets.
using Test

using CLIMAParameters
using CLIMAParameters.Atmos.SubgridScale: C_smag
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# ## [Initial Conditions](@id init)
# This example of a rising thermal bubble can be classified as an initial value
# problem. We must (at the very least) assign values for the initial variables
# in a sensible manner. This example demonstrates the use of functions defined
# in the [`MoistThermodynamics`](@ref MoistThermodynamics-docs) package to
# generate the appropriate initial state for our problem.

#md # !!! note
#md #     The following variables are assigned in the initial condition
#md #     - `state.ρ` = Scalar quantity for initial density profile
#md #     - `state.ρu`= 3-component vector for initial momentum profile
#md #     - `state.ρe`= Scalar quantity for initial total-energy profile
#md #       humidity
#md #     - `state.tracers.ρχ` = Vector of four tracers (here, for demonstration
#md #       only; we can interpret these as dye injections for visualisation
#md #       purposes)
function init_risingbubble!(bl, state, aux, (x, y, z), t)
    # Problem float-type
    FT = eltype(state)

    # Unpack constant parameters
    R_gas::FT = R_d(bl.param_set)
    c_p::FT = cp_d(bl.param_set)
    c_v::FT = cv_d(bl.param_set)
    p0::FT = MSLP(bl.param_set)
    _grav::FT = grav(bl.param_set)
    γ::FT = c_p / c_v

    # Define bubble center and background potential temperature
    xc::FT = 1250
    yc::FT = 1250
    zc::FT = 1000
    r = sqrt((x - xc)^2 + (y - yc)^2 + (z - zc)^2)
    rc::FT = 500
    # TODO: clean this up, or add convenience function:
    # This is configured in the reference hydrostatic state
    θ_ref::FT = bl.ref_state.virtual_temperature_profile.T_surface
    Δθ::FT = 0

    # Compute temperature difference over bubble region
    if r <= rc
        Δθ = FT(5) * cospi(r / rc / 2)
    end

    # Compute perturbed thermodynamic state:
    θ = θ_ref + Δθ                                      # potential temperature
    π_exner = FT(1) - _grav / (c_p * θ) * z             # exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas)      # density
    T = θ * π_exner
    e_int = internal_energy(bl.param_set, T)
    ts = PhaseDry(bl.param_set, e_int, ρ)
    ρu = SVector(FT(0), FT(0), FT(0))                   # momentum
    #State (prognostic) variable assignment
    e_kin = FT(0)                                       # kinetic energy
    e_pot = gravitational_potential(bl.orientation, aux)# potential energy
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)         # total energy

    ρχ = FT(0)                                          # tracer

    # We inject tracers at the initial condition at some specified z coordinates
    if 500 < z <= 550
        ρχ += FT(0.05)
    end

    # We want 4 tracers
    ntracers = 4

    # Define 4 tracers, (arbitrary scaling for this demo problem)
    ρχ = SVector{ntracers, FT}(ρχ, ρχ / 2, ρχ / 3, ρχ / 4)

    # Assign State Variables
    state.ρ = ρ
    state.ρu = ρu
    state.ρe = ρe_tot
    state.tracers.ρχ = ρχ
end

# ## [Model Configuration](@id config-helper)
# We define a configuration function to assist in prescribing the physical
# model. The purpose of this is to populate the
# [`ClimateMachine.AtmosLESConfiguration`](@ref LESConfig) with arguments
# appropriate to the problem being considered.
function config_risingbubble(FT, N, resolution, xmax, ymax, zmax)

    # Choose an Explicit Multi-rate Solver from the existing [ODESolvers](@ref
    # ODESolvers-docs) options Apply the outer constructor to define the
    # `ode_solver`. Here `AtmosAcousticGravityLinearModel` splits the
    # acoustic-gravity wave components from the advection-diffusion dynamics.
    # The 1D-IMEX method is less appropriate for the problem given the current
    # mesh aspect ratio (1:1)
    ode_solver = ClimateMachine.MultirateSolverType(
        linear_model = AtmosAcousticGravityLinearModel,
        slow_method = LSRK144NiegemannDiehlBusch,
        fast_method = LSRK144NiegemannDiehlBusch,
        timestep_ratio = 10,
    )

    # Since we want four tracers, we specify this and include the appropriate
    # diffusivity scaling coefficients (normally these would be physically
    # informed but for this demonstration we use integers corresponding to the
    # tracer index identifier)
    ntracers = 4
    δ_χ = SVector{ntracers, FT}(1, 2, 3, 4)

    # The model coefficient for the turbulence closure is defined via the
    # [CLIMAParameters
    # package](https://CliMA.github.io/CLIMAParameters.jl/latest/) A reference
    # state for the linearisation step is also defined.
    T_surface = FT(300)
    T_min_ref = FT(0)
    T_profile = DryAdiabaticProfile{FT}(param_set, T_surface, T_min_ref)
    ref_state = HydrostaticState(T_profile)

    # The fun part! Here we assemble the `AtmosModel`.
    #md # !!! note
    #md #     Docs on model subcomponent options can be found here:
    #md #     - [`param_set`](https://CliMA.github.io/CLIMAParameters.jl/latest/)
    #md #     - [`turbulence`](@ref Turbulence-Closures-docs)
    #md #     - [`hyperdiffusion`](@ref Hyperdiffusion-docs)
    #md #     - [`source`](@ref atmos-sources)
    #md #     - [`tracers`](@ref Tracers-docs)
    #md #     - [`init_state`](@ref init)

    _C_smag = FT(C_smag(param_set))
    model = AtmosModel{FT}(
        AtmosLESConfigType,                           # Flow in a box, requires the AtmosLESConfigType
        param_set;                                    # Parameter set corresponding to earth parameters
        turbulence = SmagorinskyLilly(_C_smag),       # Turbulence closure model
        moisture = DryModel(),                        # Exclude moisture variables
        hyperdiffusion = StandardHyperDiffusion(60),  # Hyperdiffusion (4th order) model
        source = (Gravity(),),                        # Gravity is the only source term here
        tracers = NTracers{ntracers, FT}(δ_χ),        # Tracer model with diffusivity coefficients
        ref_state = ref_state,                        # Reference state
        init_state_conservative = init_risingbubble!, # Apply the initial condition
    )

    # Finally, we pass a `Problem Name` string, the mesh information, and the
    # model type to  the [`AtmosLESConfiguration`] object.
    config = ClimateMachine.AtmosLESConfiguration(
        "DryRisingBubble",       # Problem title [String]
        N,                       # Polynomial order [Int]
        resolution,              # (Δx, Δy, Δz) effective resolution [m]
        xmax,                    # Domain maximum size [m]
        ymax,                    # Domain maximum size [m]
        zmax,                    # Domain maximum size [m]
        param_set,               # Parameter set.
        init_risingbubble!,      # Function specifying initial condition
        solver_type = ode_solver,# Time-integrator type
        model = model,           # Model type
    )
    return config
end

#md # !!! note
#md #     `Keywords` are used to specify some arguments (see appropriate source
#md #     files).

# ## [Diagnostics](@id config_diagnostics)
# Here we define the diagnostic configuration specific to this problem.
function config_diagnostics(driver_config)
    interval = "10000steps"
    dgngrp = setup_atmos_default_diagnostics(interval, driver_config.name)
    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

function main()
    # These are essentially arguments passed to the
    # [`config_risingbubble`](@ref config-helper) function.  For type
    # consistency we explicitly define the problem floating-precision.
    FT = Float64
    # We need to specify the polynomial order for the DG discretization,
    # effective resolution, simulation end-time, the domain bounds, and the
    # courant-number for the time-integrator. Note how the time-integration
    # components `solver_config` are distinct from the spatial / model
    # components in `driver_config`. `init_on_cpu` is a helper keyword argument
    # that forces problem initialisation on CPU (thereby allowing the use of
    # random seeds, spline interpolants and other special functions at the
    # initialisation step.)
    N = 4
    Δh = FT(50)
    Δv = FT(50)
    resolution = (Δh, Δh, Δv)
    xmax = FT(2500)
    ymax = FT(2500)
    zmax = FT(2500)
    t0 = FT(0)
    timeend = FT(1000)
    CFL = FT(20)

    # Assign configurations so they can be passed to the `invoke!` function
    driver_config = config_risingbubble(FT, N, resolution, xmax, ymax, zmax)
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFL,
    )
    dgn_config = config_diagnostics(driver_config)

    # User defined filter (TMAR positivity preserving filter)
    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do (init = false)
        Filters.apply!(solver_config.Q, 6, solver_config.dg.grid, TMARFilter())
        nothing
    end

    # Invoke solver (calls `solve!` function for time-integrator), pass the driver, solver and diagnostic config
    # information.
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbtmarfilter,),
        check_euclidean_distance = true,
    )

    # Check that the solution norm is reasonable.
    @test isapprox(result, FT(1); atol = 1.5e-3)
end

# The experiment definition is now complete. Time to run it.

# ## Running the Experiment
# `julia --project /experiments/AtmosLES/risingbubble.jl` will run the
# experiment from the main ClimateMachine.jl directory, with diagnostics output
# at the intervals specified in [`config_diagnostics`](@ref
# config_diagnostics).  You can also prescribe command line arguments (docs
# pending, `Driver.jl`) for simulation update and output specifications.  For
# rapid turnaround, we recommend that you run this experiment on a GPU.

# ## [Output Visualisation](@id output-viz)
#md # See the `ClimateMachine` [command line arguments](@ref ClimateMachine-args)
#md # for generating output.
#md #
#md # Given VTK output,
#md # - [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit/)
#md # - [Paraview](https://wci.llnl.gov/simulation/computer-codes/visit/)
#md # are two commonly used programs for `.vtu` files.
#md #
#md # For NetCDF or JLD2 diagnostics you may use Julia's
#md # [`NCDatasets`](https://github.com/Alexander-Barth/NCDatasets.jl) and
#md # [`JLD2`](https://github.com/JuliaIO/JLD2.jl) packages with a suitable
#md # plotting program.

#

main()


text:
include("obc.jl")
include("mpm.jl")


text:
"""
    gdalinfo(dataset::AbstractDataset, options = String[])

List various information about a GDAL supported raster dataset.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdalinfo utility.

### Returns
String corresponding to the information about the raster dataset.
"""
function gdalinfo(dataset::AbstractDataset, options = String[])::String
    gdal_info_options = GDAL.gdalinfooptionsnew(options, C_NULL)
    return try
        GDAL.gdalinfo(dataset.ptr, gdal_info_options)
    finally
        GDAL.gdalinfooptionsfree(gdal_info_options)
    end
end

"""
    unsafe_gdaltranslate(
        dataset::AbstractDataset,
        options = String[];
        dest = "/vsimem/tmp")

Convert raster data between different formats.

### Parameters
* **dataset**: The dataset to be translated.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdal_translate utility.

### Returns
The output dataset.
"""
function unsafe_gdaltranslate(
    dataset::AbstractDataset,
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdaltranslateoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdaltranslate(dest, dataset.ptr, options, usage_error)
    GDAL.gdaltranslateoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdalwarp(
        datasets::Vector{<:AbstractDataset},
        options = String[];
        dest = "/vsimem/tmp")

Image reprojection and warping function.

### Parameters
* **datasets**: The list of input datasets.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdalwarp utility.

### Returns
The output dataset.
"""
function unsafe_gdalwarp(
    datasets::Vector{<:AbstractDataset},
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdalwarpappoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdalwarp(
        dest,
        C_NULL,
        length(datasets),
        [ds.ptr for ds in datasets],
        options,
        usage_error,
    )
    GDAL.gdalwarpappoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdalvectortranslate(
        datasets::Vector{<:AbstractDataset},
        options = String[];
        dest = "/vsimem/tmp")

Convert vector data between file formats.

### Parameters
* **datasets**: The list of input datasets (only 1 supported currently).
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the ogr2ogr utility.

### Returns
The output dataset.
"""
function unsafe_gdalvectortranslate(
    datasets::Vector{<:AbstractDataset},
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdalvectortranslateoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdalvectortranslate(
        dest,
        C_NULL,
        length(datasets),
        [ds.ptr for ds in datasets],
        options,
        usage_error,
    )
    GDAL.gdalvectortranslateoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdaldem(
        dataset::AbstractDataset,
        processing::String,
        options = String[];
        dest = "/vsimem/tmp",
        colorfile)

Tools to analyze and visualize DEMs.

### Parameters
* **dataset**: The source dataset.
* **pszProcessing**: the processing to apply (one of "hillshade", "slope",
    "aspect", "color-relief", "TRI", "TPI", "Roughness").
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdaldem utility.

# Keyword Arguments
* **colorfile**: color file (mandatory for "color-relief" processing,
    should be empty otherwise).

### Returns
The output dataset.
"""
function unsafe_gdaldem(
    dataset::AbstractDataset,
    processing::String,
    options = String[];
    dest = "/vsimem/tmp",
    colorfile = C_NULL,
)::Dataset
    if processing == "color-relief"
        @assert colorfile != C_NULL
    end
    options = GDAL.gdaldemprocessingoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdaldemprocessing(
        dest,
        dataset.ptr,
        processing,
        colorfile,
        options,
        usage_error,
    )
    GDAL.gdaldemprocessingoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdalnearblack(
        dataset::AbstractDataset,
        options = String[];
        dest = "/vsimem/tmp")

Convert nearly black/white borders to exact value.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the nearblack utility.

### Returns
The output dataset.
"""
function unsafe_gdalnearblack(
    dataset::AbstractDataset,
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdalnearblackoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdalnearblack(dest, C_NULL, dataset.ptr, options, usage_error)
    GDAL.gdalnearblackoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdalgrid(
        dataset::AbstractDataset,
        options = String[];
        dest = "/vsimem/tmp")

Create a raster from the scattered data.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdal_grid utility.

### Returns
The output dataset.
"""
function unsafe_gdalgrid(
    dataset::AbstractDataset,
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdalgridoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdalgrid(dest, dataset.ptr, options, usage_error)
    GDAL.gdalgridoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdalrasterize(
        dataset::AbstractDataset,
        options = String[];
        dest = "/vsimem/tmp")

Burn vector geometries into a raster.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdal_rasterize utility.

### Returns
The output dataset.
"""
function unsafe_gdalrasterize(
    dataset::AbstractDataset,
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdalrasterizeoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdalrasterize(dest, C_NULL, dataset.ptr, options, usage_error)
    GDAL.gdalrasterizeoptionsfree(options)
    return Dataset(result)
end

"""
    unsafe_gdalbuildvrt(
        datasets::Vector{<:AbstractDataset},
        options = String[];
        dest = "/vsimem/tmp")

Build a VRT from a list of datasets.

### Parameters
* **datasets**: The list of input datasets.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdalbuildvrt utility.

### Returns
The output dataset.
"""
function unsafe_gdalbuildvrt(
    datasets::Vector{<:AbstractDataset},
    options = String[];
    dest = "/vsimem/tmp",
)::Dataset
    options = GDAL.gdalbuildvrtoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.gdalbuildvrt(
        dest,
        length(datasets),
        [ds.ptr for ds in datasets],
        C_NULL,
        options,
        usage_error,
    )
    GDAL.gdalbuildvrtoptionsfree(options)
    return Dataset(result)
end


text:
function get_xyz_vector_from_rφz_field_vector_at_rφz(field,r,φ,z,ir,iφ,iz)::Vector
    startpoint_vector = get_xyz_vector_from_rφz_vector([r,φ,z])
    endpoint_vector = get_xyz_vector_from_rφz_vector([r,φ,z]+field[ir,iφ,iz])
    xyz_vector = endpoint_vector-startpoint_vector
    for ic in 1:size(xyz_vector,1)
        isapprox(xyz_vector[ic],0.0) ? xyz_vector[ic] = 0.0 : nothing
    end
    xyz_vector
end

function get_xy_magnitude(xyz_vector::AbstractArray)
    sqrt(xyz_vector[1]^2+xyz_vector[2]^2)
end

@recipe function f(electrical_field::Array{SVector{3,T},3}, grid::CylindricalGrid{T}; view=:Components, plane=:rz, i_fixed=3, spacing = 8, vectorscale = 0.0018, SI_factor=1/1000.) where{T <: SSDFloat}

    vectorfield = electrical_field.*SI_factor
    units = Dict(1e-3=>"mm",1e-2=>"cm",1e-1=>"dm",1=>"m",1.0=>"m")
    if view == :Components
        if plane == :rz
            layout := (2,2)
            size := (800,900)
            vectorfield_r = get_component_field(vectorfield);
            vectorfield_φ = get_component_field(vectorfield,:phi);
            vectorfield_z = get_component_field(vectorfield,:z);
            vectorfield_magn = get_magnitude_of_rφz_vector.(vectorfield);

            st := :heatmap
            # colorbar_title := "Field Strength / V m\$\^\{-1\}\$"
            @series begin
                subplot := 1
                title := "r_component"
                ylabel --> "z ["*units[SI_factor]*"]"
                grid[:r] ./ SI_factor, grid[:z] ./ SI_factor, vectorfield_r[:,i_fixed,:]'
            end
            @series begin
                subplot := 2
                colorbar_title --> "Field Strength [V / "*units[SI_factor]*"]"
                title := "φ_component"

                grid[:r] ./SI_factor, grid[:z] ./SI_factor, vectorfield_φ[:,i_fixed,:]'
            end
            @series begin
                subplot := 3
                title := "z_component"
                ylabel --> "z ["*units[SI_factor]*"]"
                xlabel --> "r ["*units[SI_factor]*"]"
                grid[:r] ./SI_factor, grid[:z] ./SI_factor, vectorfield_z[:,i_fixed,:]'
            end
            @series begin
                subplot := 4
                colorbar_title --> "Field Strength [V / "*units[SI_factor]*"]"
                xlabel --> "r ["*units[SI_factor]*"]"
                title:= "magnitude"
                grid[:r] ./SI_factor, grid[:z] ./SI_factor, vectorfield_magn[:,i_fixed,:]'
            end
        end
    elseif view == :ef
        if plane == :rφ
            vectorfield_xyz = Array{Vector{Float32}}(undef,size(vectorfield,1),size(vectorfield,2),size(vectorfield,3));
            for (iz,z) in enumerate(grid[:z])
                for (iφ,φ) in enumerate(grid[:φ])
                    for (ir,r) in enumerate(grid[:r])
                        vectorfield_xyz[ir,iφ,iz]=get_xyz_vector_from_rφz_field_vector_at_rφz(vectorfield,r,φ,z,ir,iφ,iz)
                    end
                end
            end
        elseif plane == :xy
            vectorfield_xyz = electrical_field.field
        end
        vectorfield_xy_magn = map(x->get_xy_magnitude(x),vectorfield_xyz[:,:,i_fixed])
        max_magn = maximum(vectorfield_xy_magn)
        diff_magn = max_magn-minimum(vectorfield_xy_magn)

        size := (800,800)
        line := (:arrow,:blue,2)
        label := ""
        ylabel := "y "
        xlabel := "x "
        title := "z = $(round(grid[:z][i_fixed]/SI_factor,digits=2)) / mm"
        xlims := (-1.2/SI_factor*maximum(grid[:r]),1.2/SI_factor*maximum(grid[:r]))
        ylims := (-1.2/SI_factor*maximum(grid[:r]),1.2/SI_factor*maximum(grid[:r]))
        for (ir,r) in enumerate(grid[:r][1:spacing:end])
        for (iφ,φ) in enumerate(grid.φ)
                x= r*cos(φ)
            y= r*sin(φ)
            ir_actual=findfirst(x->x==r,grid[:r])
            iφ_actual=findfirst(x->x==φ,grid.φ)
                xy_magn = vectorfield_xy_magn[ir_actual,iφ_actual]
                vector=vectorfield_xyz[ir_actual,iφ_actual,i_fixed]/xy_magn
                vector*=((vectorfield_xy_magn[ir_actual,iφ_actual]-0.8*minimum(vectorfield_xy_magn))/diff_magn)
                vector*=vectorscale
                @series begin
                    [x-0.5*vector[1],x+0.5*vector[1]]/SI_factor, [y-0.5*vector[2],y+0.5*vector[2]]/SI_factor
                end
            end
        end
    end
end
#
# @recipe function f(electrical_field::Array{SVector{3,T},3}, d::SolidStateDetector, grid::CylindricalGrid; φ_value=deg2rad(0), spacing=0.003, steps=1000, myscale=1, potential=true) where{T <: SSDFloat}
#     size --> (900,1100)
#     d.bulk_type == :ptype ? interpolated_efield = setup_interpolated_vectorfield(electrical_field, grid) : nothing
#     d.bulk_type == :ntype ? interpolated_efield = setup_interpolated_vectorfield(map(x->-1*x,electrical_field), grid) : nothing
#     if potential==true
#         @series begin
#             φ --> φ_value
#             grid
#         end
#     end
#
#     corner_offset = 5e-5
#     spawn_positions::Array{Array{T,1},1}=[]
#
#     for (i,tuple) in enumerate(d.segmentation_r_ranges)
#         if tuple[1]==tuple[2]
#             for z in corner_offset + d.segmentation_z_ranges[i][1]:spacing:d.segmentation_z_ranges[i][2] - corner_offset
#                 push!(spawn_positions,[tuple[1],φ_value,z])
#             end
#         end
#     end
#
#     for (i,tuple) in enumerate(d.segmentation_z_ranges)
#         if tuple[1]==tuple[2]
#             for r in corner_offset+d.segmentation_r_ranges[i][1]:spacing:d.segmentation_r_ranges[i][2] - corner_offset
#                 push!(spawn_positions,[r,φ_value,tuple[1]])
#             end
#         end
#     end
#
#     for (i,orientation) in enumerate(d.segmentation_types)
#         if orientation != "Tubs"
#             if orientation[1]=='c' o=-1 else o=1 end
#             for z in d.segmentation_z_ranges[i][1]:spacing:d.segmentation_z_ranges[i][2]
#                 push!(spawn_positions,[o * corner_offset+analytical_taper_r_from_φz(φ_value, z, orientation, d.segmentation_r_ranges[i],
#                                                                                d.segmentation_phi_ranges[i],
#                                                                                d.segmentation_z_ranges[i]),
#                                                                                φ_value, z])
#             end
#         end
#     end
#
#     for i in eachindex(spawn_positions[1:end-1])
#
#         xpath, ypath, zpath = driftonecharge(d, get_xyz_vector_from_rφz_vector(spawn_positions[i]), steps, myscale*1e-9, :e, interpolated_efield, interpolated_efield)
#         rpath=[]
#
#         for ir in eachindex(xpath)
#             push!(rpath, sqrt(xpath[ir]^2+ypath[ir]^2))
#         end
#
#         @series begin
#             c --> :white
#             label --> ""
#             rpath, zpath
#         end
#     end
# end
using Base.Math
@recipe function f( electric_field::Array{ <:StaticVector{3, T}, 3}, det::SolidStateDetector; φ=missing, z = missing, scaling_factor=1.0, spacings = [3,20]) where T
    ismissing(φ) && !ismissing(z) ? iz = searchsortednearest(det.zs,T(z)) : nothing
    !ismissing(φ) && ismissing(z) ? iφ = searchsortednearest(det.φs,T(φ)) : nothing
    guidefontsize --> 15
    tickfontsize -->15
    # labelfontsize --> 13
    if ismissing(φ)
        xy = T[0, 0]
        values = T[0, 0]
        magnitudes_xy = []
        for (iφ ,φ)  in enumerate(det.φs)
            if (iφ+spacings[1]-1)%spacings[1]==0
                for (ir, r) in enumerate(det.rs)
                    if (ir+spacings[2]-1)%spacings[2] == 0
                        xy = hcat(xy,[ r * cos(φ), r * sin(φ)])
                        values = hcat(values, electric_field[ir, iφ, iz][1:2])
                        push!(magnitudes_xy,sqrt(electric_field[ir, iφ, iz][1]^2 + electric_field[ir, iφ, iz][2]^2))
                    end
                end
            end
        end

        num = minimum([mean(diff(det.φs)*spacings[1]),mean(diff(det.rs))*spacings[2]])
        nom = maximum(magnitudes_xy)
        scaling = num/nom * scaling_factor
        for i in 1:size(xy,2)
            @series begin
                title --> "z = $(round(det.zs[iz] .* 1000,digits=2)) mm"
                xlabel --> "x / mm"
                ylabel --> "y / mm"
                aspect_ratio --> 1
                label --> ""
                color --> :red
                arrow --> true
                [xy[1,i] - 0.5*scaling*values[1,i], xy[1,i] + 0.5*scaling*values[1,i]] .* 1000, [xy[2,i] - 0.5*scaling*values[2,i], xy[2,i] + 0.5*scaling*values[2,i]] .* 1000
            end
        end

    elseif ismissing(z)
        rz = T[0, 0]
        values = T[0, 0]
        magnitudes_rz = []
        for (iz ,z)  in enumerate(det.zs)
            if (iz+spacings[1])%spacings[1]==0
                for (ir, r) in enumerate(det.rs)
                    if (ir+spacings[2]-1)%spacings[2] == 0
                        rz = hcat(rz,[ r, z])
                        values = hcat(values, electric_field[ir, iφ, iz][1:2:3])
                        push!(magnitudes_rz,sqrt(electric_field[ir, iφ, iz][1]^2 + electric_field[ir, iφ, iz][3]^2))
                    end
                end
            end
        end
        num = minimum([mean(diff(det.zs)*spacings[1]),mean(diff(det.rs))*spacings[2]])
        nom = maximum(magnitudes_rz)
        scaling = num/nom * scaling_factor
        for i in 1:size(rz,2)
            @series begin
                title --> "φ = $(round(rad2deg(det.φs[iφ]) , digits=2)) deg"
                xlabel --> "r [mm]"
                ylabel --> "z [mm]"
                aspect_ratio --> 1
                label --> ""
                color --> :red
                arrow --> true
                [rz[1,i] - 0.5*scaling*values[1,i], rz[1,i] + 0.5*scaling*values[1,i]] .* 1000, [rz[2,i] - 0.5*scaling*values[2,i], rz[2,i] + 0.5*scaling*values[2,i]] .* 1000
            end
        end
    end
end

@userplot MyQuiver
@recipe function f(gdd::MyQuiver; scaling=1.0)
    xy::Matrix = gdd.args[1]
    values = gdd.args[2]
    for i in 1:size(xy,2)
        @series begin
            arrow --> true
            [xy[1,i] - 0.5*scaling*values[1,i], xy[1,i] + 0.5*scaling*values[1,i]], [xy[2,i] - 0.5*scaling*values[2,i], xy[2,i] + 0.5*scaling*values[2,i]]
        end
    end
end

# @recipe function f( electric_field::Array{ <:StaticVector{3, T}, 3}, det::SolidStateDetector, ep::ElectricPotential{T};
#             φ=missing, spacing=T(0.003), n_steps=3000, potential=true, contours_equal_potential=true, offset = T(5e-5)) where {T <: SSDFloat}
#     size --> (700,900)
#     det.bulk_type == :ptype ? interpolated_efield = setup_interpolated_vectorfield(electric_field, ep.grid) : nothing
#     det.bulk_type == :ntype ? interpolated_efield = setup_interpolated_vectorfield(map(x -> -1*x, electric_field), ep.grid) : nothing
#     if ismissing(φ)
#         φ = 0
#     end
#     φ_rad = deg2rad(φ)
#     aspect_ratio --> 1
#     title --> "Electric Field Lines @φ=$(φ)°"
#     xlabel --> L"$r$ / m"
#     ylabel --> L"$z$ / m"
#
#     if potential==true
#         @series begin
#             contours_equal_potential --> contours_equal_potential
#             φ --> φ
#             ep
#         end
#     end
#
#     spawn_positions_cyl::Vector{CylindricalPoint} = []
#
#     for contact in det.contacts
#         if (typeof(contact) == SSD.Contact{T,:N} && det.bulk_type == :ntype) || (typeof(contact) == SSD.Contact{T,:P} && det.bulk_type == :ptype)
#             nothing
#         else
#             for g in contact.geometry_positive
#                 if typeof(g) == SSD.Tube{T}
#                     rStart,rStop = get_r(g)
#                     zStart,zStop = get_z(g)
#                     for r in rStart:spacing:rStop
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(r,φ,zStart+offset))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(r,φ,zStart-offset))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(r,φ,zStop+offset))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(r,φ,zStop-offset))
#                     end
#                     for z in zStart:spacing:zStop
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(rStart+offset,φ,z))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(rStart-offset,φ,z))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(rStop+offset,φ,z))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(rStop-offset,φ,z))
#                     end
#                 elseif typeof(g) == SSD.ConeMantle{T}
#                     zStart,zStop = get_z(g)
#                     for z in zStart:spacing:zStop
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(get_diagonal_r_from_z(g.cone, z)+offset,φ,z))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(get_diagonal_r_from_z(g.cone, z)-offset,φ,z))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(get_diagonal_r_from_z(g.cone, z)+offset,φ,z))
#                         push!(spawn_positions_cyl,CylindricalPoint{T}(get_diagonal_r_from_z(g.cone, z)-offset,φ,z))
#                     end
#                 end
#             end
#         end
#     end
#     filter!(x -> x in det && !in(x, det.contacts),spawn_positions_cyl)
#     spawn_positions_xyz::Vector{CartesianPoint{T}} = map(x -> CartesianPoint(CylindricalPoint{T}(x[1],x[2],x[3])), spawn_positions_cyl)
#     for i in eachindex(spawn_positions_cyl[1:end])
#         path = [@SVector zeros(T,3) for i in 1:n_steps]
#         drift_charge!(path, det, spawn_positions_xyz[i], T(1f-9), interpolated_efield)
#         @series begin
#             c --> :white
#             label --> ""
#             map(x->sqrt(x[1]^2+x[2]^2),path), map(x->x[3],path)
#         end
#     end
# end

@userplot Plot_electric_field
@recipe function f(gdd::Plot_electric_field; φ = missing, spacing = 10, grid_spacing=[0.0005, deg2rad(1.0), 0.0005], n_steps=3000, potential=true, contours_equal_potential=true, offset = (5e-5))
    setup = gdd.args[1]
    T = Float32
    φ = ismissing(φ) ? T(0) : T(φ)
    φ_rad = deg2rad(φ)
    aspect_ratio --> 1
    title --> "Electric Field Lines @φ=$(round(φ, digits=2))°"
    xlabel --> L"$r$ / m"
    ylabel --> L"$z$ / m"

    if potential==true
        @series begin
            contours_equal_potential --> contours_equal_potential
            φ --> φ
            setup.electric_potential
        end
    end

    contacts_to_spawn_charges_for = filter!(x -> x.id !=1, SSD.Contact{T}[c for c in setup.detector.contacts])
    spawn_positions = CylindricalPoint{T}[]
    grid = Grid(setup.detector, init_grid_spacing = grid_spacing, full_2π = true)
    pt_offset = T[offset,0.0,offset]

    for c in contacts_to_spawn_charges_for
        ongrid_positions= map(x-> CylindricalPoint{T}(grid[x...]),SSD.paint_object(c, grid, φ_rad))
        for position in ongrid_positions
            push!(spawn_positions, CylindricalPoint{T}((position + pt_offset)...))
            push!(spawn_positions, CylindricalPoint{T}((position - pt_offset)...))
        end
    end

    filter!(x -> x in setup.detector && !in(x, setup.detector.contacts), spawn_positions)

    el_field_itp = SSD.get_interpolated_drift_field(setup.electric_field.data, setup.electric_field.grid)
    el_field_itp_inv = SSD.get_interpolated_drift_field(setup.electric_field.data .* -1, setup.electric_field.grid)
    for (ipos, pos) in enumerate(spawn_positions)
        if ((spacing-1)+ipos)%spacing == 0
            path = CartesianPoint{T}[CartesianPoint{T}(0.0,0.0,0.0) for i in 1:n_steps]
            SSD.drift_charge!(path, setup.detector, setup.electric_potential.grid, CartesianPoint(pos), T(2e-9), el_field_itp, verbose = false )
            @series begin
                c --> :white
                label --> ""
                map(x->sqrt(x[1]^2+x[2]^2),path), map(x->x[3],path)
            end

            path = CartesianPoint{T}[CartesianPoint{T}(0.0,0.0,0.0) for i in 1:n_steps]
            SSD.drift_charge!(path, setup.detector, setup.electric_potential.grid, CartesianPoint(pos), T(2e-9), el_field_itp_inv, verbose = false )
            @series begin
                c --> :white
                label --> ""
                map(x->sqrt(x[1]^2+x[2]^2),path), map(x->x[3],path)
            end

        end
    end
end