text:
##
##!------------------------------------------------------------------------------------------------!
##  digraph.jl
##
##  Digraphs and related functionality.
##
##  Creation date: 18/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##

abstract Digraph

abstract BareDigraph

type DigraphException <: Exception
end

type DigraphDynamic <: Digraph
    order::UInt32
    size::UInt64
    adjacencies_out::Vector{Vector{UInt32}}
    adjacencies_in::Vector{Vector{UInt32}}
    vertex_info::Vector{Int64}
    edge_info::Vector{Vector{Int64}}
    vertex_indicator::Vector{Bool} # vertex_indicator[u] is true if u is in the graph.
end

type DigraphStatic <: Digraph
    order::UInt32
    size::UInt64
    adjacencies_out::Vector{Vector{UInt32}}
    adjacencies_in::Vector{Vector{UInt32}}
    vertex_info::Vector{Int64}
    edge_info::Vector{Vector{Int64}}
end

type BareDigraphDynamic <: BareDigraph
    order::UInt32
    size::UInt64
    adjacencies_out::Vector{Vector{UInt32}}
    adjacencies_in::Vector{Vector{UInt32}}
    vertex_indicator::Vector{Bool} # vertex_indicator[u] is true if u is in the graph.
end

type BareDigraphStatic <: BareDigraph
    order::UInt32
    size::UInt64
    adjacencies_out::Vector{Vector{UInt32}}
    adjacencies_in::Vector{Vector{UInt32}}
end

function DigraphDynamic(order::UInt32)
    adjacencies_out = Vector{Vector{UInt32}}(order)
    adjacencies_in = Vector{Vector{UInt32}}(order)
    vertex_info = Vector{Int64}(order)
    edge_info = Vector{Vector{Int64}}(order)
    vertex_indicator = Vector{Bool}(order)

    for i in 1:order
        adjacencies_out[i] = Vector{UInt32}()
        adjacencies_in[i] = Vector{UInt32}()
        edge_info[i] = Vector{Int64}()
        vertex_indicator[i] = true
    end
    return DigraphDynamic(order, UInt64(0), adjacencies_out, adjacencies_in,
        vertex_info, edge_info, vertex_indicator)
end

function DigraphStatic(order::UInt32)
    adjacencies_out = Vector{Vector{UInt32}}(order)
    adjacencies_in = Vector{Vector{UInt32}}(order)
    vertex_info = Vector{Int64}(order)
    edge_info = Vector{Vector{Int64}}(order)

    for i in 1:order
        adjacencies_out[i] = Vector{UInt32}()
        adjacencies_in[i] = Vector{UInt32}()
        edge_info[i] = Vector{Int64}()
    end
    return DigraphStatic(order, UInt64(0), adjacencies_out, adjacencies_in,
        vertex_info, edge_info)
end

function BareDigraphDynamic(order::UInt32)
    adjacencies_out = Vector{Vector{UInt32}}(order)
    adjacencies_in = Vector{Vector{UInt32}}(order)
    vertex_indicator = Vector{Bool}(order)

    for i in 1:order
        adjacencies_out[i] = Vector{UInt32}()
        adjacencies_in[i] = Vector{UInt32}()
        vertex_indicator[i] = true
    end
    return BareDigraphDynamic(order, UInt64(0), adjacencies_out, adjacencies_in,
        vertex_indicator)
end

function BareDigraphStatic(order::UInt32)
    adjacencies_out = Vector{Vector{UInt32}}(order)
    adjacencies_in = Vector{Vector{UInt32}}(order)

    for i in 1:order
        adjacencies_out[i] = Vector{UInt32}()
        adjacencies_in[i] = Vector{UInt32}()
    end
    return BareDigraphStatic(order, UInt64(0), adjacencies_out, adjacencies_in)
end

## Properties
order(g::Union{Digraph, BareDigraph}) = g.order
size(g::Union{Digraph, BareDigraph}) = g.size
degree(g::Union{Digraph, BareDigraph}, u::UInt32) = degree_out(g, u) + degree_in(g, u)
degree_out(g::Union{Digraph, BareDigraph}, u::UInt32) = length(g.adjacencies_out[u])
degree_in(g::Union{Digraph, BareDigraph}, u::UInt32) = length(g.adjacencies_in[u])

## Access
function vertices(g::Union{DigraphDynamic, BareDigraphDynamic})
    vertices = Vector{UInt32}(g.order)
    # order keeps changing; use length of vertex_indicator.
    i, j = length(g.vertex_indicator), g.order
    while i >= 1
        if g.vertex_indicator[i]
            vertices[j] = i
            j -= 1
        end
        i -= 1
    end
    return vertices
end

function vertices(g::Union{DigraphStatic, BareDigraphStatic})
    return Vector{UInt32}(1:g.order)
end

function edges(g::Union{DigraphDynamic, BareDigraphDynamic})
    edges = Vector{Tuple{UInt32, UInt32}}(g.size)
    j = 1
    l = length(g.vertex_indicator)
    for i in 1:l
        if !g.vertex_indicator[i]
            continue
        end
        for t in g.adjacencies_out[i]
            edges[j] = (i, t)
            j += 1
        end
    end
    return edges
end

function edges(g::Union{DigraphStatic, BareDigraphStatic})
    edges = Vector{Tuple{UInt32, UInt32}}(g.size)
    j = 1
    for i in 1:g.order
        for t in g.adjacencies_out[i]
            edges[j] = (i, t)
            j += 1
        end
    end
    return edges
end

function exists(g::Union{DigraphDynamic, BareDigraphDynamic}, u::UInt32)
    return g.vertex_indicator[u]
end

function exists(g::Union{DigraphStatic, BareDigraphStatic}, u::UInt32)
    return 0 < u <= g.order
end

# Finds the index of u in l or returns 0, if u is not in l.
function _find_index(l::Vector{UInt32}, u::UInt32)
    for i in 1:length(l)
        if u == l[i]
            return i
        end
    end
    return 0
end

# Does a binary search on l[i:j] and returns the index if u is found, or 0 if not.
function _binary_search(l::Vector{UInt32}, i::UInt32, j::UInt32, u::UInt32)
    # l contains tuples of the form (target, edgeinfo).
    if i == j
        return l[i] == u ? i : 0
    end

    k = uint32((i + j) / 2)
    v = l[k]
    if u < v
        return _binary_search(l, i, k, u)
    elseif u > v
        return _binary_search(l, k + 1, j, u)
    else
        return k
    end
end

function adjacent(g::Union{DigraphDynamic, BareDigraphDynamic}, u::UInt32, v::UInt32)
    return _find_index(g.adjacencies_out[u], v) != 0
end

function adjacent(g::Union{DigraphStatic, BareDigraphStatic}, u::UInt32, v::UInt32)
    # For static graphs, the adjacency list is sorted.
    return _binary_search(g.adjacencies_out[u], 1, length(g.adjacencies_out[u]), v) != 0
end

function neighbours_out(g::Union{Digraph, BareDigraph}, u::UInt32)
    return g.adjacencies_out[u]
end

function neighbours_in(g::Union{Digraph, BareDigraph}, u::UInt32)
    return g.adjacencies_in[u]
end

function vertex_info(g::Union{DigraphDynamic, DigraphStatic}, u::UInt32)
    return g.vertex_info[u]
end

function edge_info(g::DigraphDynamic, u::UInt32, v::UInt32)
    l = g.edge_info[u]
    return l[_find_index(g.adjacencies_out[u], v)]
end

function edge_info(g::DigraphStatic, u::UInt32, v::UInt32)
    l = g.edge_info[u]
    return l[_binary_search(g.adjacencies_out[u], 1, length(g.adjacencies_out[u]), v)]
end

## Mutation
function add_vertex!(g::DigraphDynamic)
    # The vertex is added at the end.
    push!(g.adjacencies_out, Vector{UInt32}())
    push!(g.adjacencies_in, Vector{UInt32}())
    push!(g.vertex_info, Int64(0))
    push!(g.edge_info, Vector{Int64}())
    push!(g.vertex_indicator, true)
    g.order += 1
end

function add_vertex!(g::DigraphDynamic, u::UInt32)
    # Assume u was previously deleted.
    g.vertex_indicator[u] = true
    g.order += 1
end

function add_vertex!(g::BareDigraphDynamic)
    # The vertex is added at the end.
    push!(g.adjacencies_out, Vector{UInt32}())
    push!(g.adjacencies_in, Vector{UInt32}())
    push!(g.vertex_indicator, true)
    g.order += 1
end

function add_vertex!(g::BareDigraphDynamic, u::UInt32)
    # Assume u was previously deleted.
    g.vertex_indicator[u] = true
    g.order += 1
end

function add_edge!(g::DigraphDynamic, u::UInt32, v::UInt32)
    push!(g.adjacencies_out[u], v)
    push!(g.adjacencies_in[v], u)
    push!(g.edge_info[u], v)
    g.size += 1
end

function add_edge!(g::BareDigraphDynamic, u::UInt32, v::UInt32)
    push!(g.adjacencies_out[u], v)
    push!(g.adjacencies_out[v], u)
    g.size += 1
end

function remove_vertex!(g::Union{DigraphDynamic, BareDigraphDynamic}, u::UInt32)
    # Do not actually remove the vertex, only set vertex_indicator[u] to false.
    g.vertex_indicator[u] = false
    # Remove all edges incident with u.
    n = copy(g.adjacencies_out[u])
    for v in n
        remove_edge!(g, u, v)
    end
    g.order -= 1
end

function remove_edge!(g::DigraphDynamic, u::UInt32, v::UInt32)
    no = g.adjacencies_out[u]
    ni = g.adjacencies_out[v]
    l = g.edge_info[u]
    # Find the index of v in the adjacency list.
    i = _find_index(no, v)
    j = _find_index(ni, u)
    no[i] = no[end]
    pop!(no)
    ni[j] = ni[end]
    pop!(ni)
    l[i] = l[end]
    pop!(l)

    g.size -= 1
end

function remove_edge!(g::BareDigraphDynamic, u::UInt32, v::UInt32)
    no = g.adjacencies_out[u]
    ni = g.adjacencies_out[u]
    # Find the index of v in the adjacency list.
    i = _find_index(no, v)
    j = _find_index(ni, u)
    no[i] = no[end]
    pop!(no)
    ni[j] = ni[end]
    pop!(ni)

    g.size -= 1
end

function vertex_info!(g::Union{DigraphDynamic, DigraphStatic}, u::UInt32, info::Int64)
    g.vertex_info[u] = info
end

function edge_info!(g::Union{DigraphDynamic, DigraphStatic}, u::UInt32, v::UInt32, info::Int64)
    i = _find_index(g.adjacencies_out[u], v)
    g.edge_info[u][i] = info
end


text:
"Supertype of trait types which define a kind of duality."
abstract type DualType end

"Orthogonal dual type: duality with respect to the continuous measure of the dictionary."
struct OrthogonalDual <: DualType end
"Discrete orthogonal dual type: duality with respect to the discrete measure of the dictionary."
struct DiscreteOrthogonalDual <: DualType end
"Duality defined by the inverse of the Gram matrix with a given measure."
struct GramDual <: DualType end

"Return the dictionary dual to `dict` in the sense defined by the dual type."
function dual(::DualType, dict; options...) end

# Use orthogonal dual as the default
dual(dict::Dictionary, args...; options...) = dual(OrthogonalDual(), dict, args...; options...)

dual(::OrthogonalDual, dict; options...) = gramdual(dict, measure(dict); options...)
dual(::DiscreteOrthogonalDual, dict; options...) = gramdual(dict, discretemeasure(dict); options...)
# For GramDual, the measure has to be specified as a keyword argument
dual(::GramDual, dict; measure, options...) = gramdual(dict, measure; options...)

gramdual(dict::Dictionary, measure::Measure; options...) =
    default_gramdual(dict, measure; options...)

function default_gramdual(dict::Dictionary, measure::Measure; options...)
    try
        conj(inv(gram(dict, measure; options...))) * dict
    catch DimensionMismatch
        # @warn "Convert DictionaryOperator to dense ArrayOperator"
        # ArrayOperator(inv(conj(Matrix(gram(dict, measure; options...)))),dict,dict) * dict
        error("When does this happen? Follow the stack trace and fix the underlying issue.")
    end
end

ext: jl


text:
mutable struct Substation <: Bus
    id::Int # for internal use only, don't mess with this.
    coords::LatLon # latitude and longitude
    voltages::Vector{Real} # kV
    load::Real # MW
    generation::Real # MW
    population::Integer
    connected_to::Set{Substation}
    grouping::Vector{Bus}

    function Substation(id, coords, volts, load, gen, pop, con = Set(), group = [])
        return new(id, coords, volts, load, gen, pop, con, group)
    end
end

function Substation(
    id,
    coords::Tuple{T,T} where T <: Real,
    volts,
    load,
    gen,
    pop,
    con=Set(),
    group=[]
)
    return Substation(
        id,
        LatLon(coords[1], coords[2]),
        volts,
        load,
        gen,
        pop,
        con,
        group
    )
end

function Substation(coords::LatLon, volts, load, gen, pop, con = Set(), group = [])
    return Substation(-1, coords, volts, load, gen, pop, con, group)
end

function Substation(
    coords::Tuple{T,T} where T <: Real,
    volts,
    load,
    gen,
    pop,
    con=Set(),
    group=[]
)
    return Substation(
        -1,
        LatLon(coords[1], coords[2]),
        volts,
        load,
        gen,
        pop,
        con,
        group
    )
end

function show(io::IO, bus::Substation)
    println(io, "Substation(")
    println(io, "\tid=$(bus.id)")
    println(io, "\tcoords=$(bus.coords),")
    println(io, "\tvoltages=$(bus.voltages),")
    println(io, "\tload=$(bus.load), ")
    println(io, "\tgeneration=$(bus.generation), ")
    println(io, "\tpopulation=$(bus.population), ")
    if length(bus.connected_to) > 0
        println(io, "\tconnected_to=Set{Substation}(...)")
    else
        println(io, "\tconnected_to=Set{Substation}()")
    end
    println(io, "\tgrouping=$(bus.grouping)")
    print(io, ")")
end

function hash(a::Substation, h::UInt)
    h = hash(a.id, h)
    h = hash(a.grouping, h)
    h = hash(a.generation, h)
    h = hash(a.population, h)
    h = hash(a.voltages, h)
    h = hash(a.load, h)
    return hash(a.coords, h)
end

function ==(a::Substation, b::Substation)
    compared = Set{Tuple{UInt, UInt}}()

    function subs_equal(a::Substation, b::Substation)
        a.id == b.id &&
        a.coords == b.coords &&
        a.load== b.load &&
        a.voltages == b.voltages &&
        a.population == b.population &&
        length(a.connected_to) == length(b.connected_to) &&
        a.grouping == b.grouping &&
        a.generation == b.generation || return false

        pair = (objectid(a), objectid(b))
        if !(pair in compared)
            push!(compared, pair)

            all(subs_equal(a, b) for (a, b) in zip(a.connected_to, b.connected_to)) ||
            return false
        end
        return true
    end

    return subs_equal(a, b)
end

voltage(s::Substation) = s.voltages
capacity(s::Substation) = s.generation
population(s::Substation) = s.population
buses(s::Substation) = s.grouping


text:
module SpringCollab2020FlagToggleStarTrackSpinner

using ..Ahorn, Maple

@pardef FlagToggleStarTrackSpinner(x1::Integer, y1::Integer, x2::Integer=x1 + 16, y2::Integer=y1, speed::String="Normal", startCenter::Bool=false, flag::String="flag_toggle_star_track_spinner", inverted::Bool=false) =
    Entity("SpringCollab2020/FlagToggleStarTrackSpinner", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], speed=speed, startCenter=startCenter, flag=flag, inverted=inverted)

const placements = Ahorn.PlacementDict(
    "Star (Track, $(uppercasefirst(speed)), Flag Toggle) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FlagToggleStarTrackSpinner,
        "line",
        Dict{String, Any}(
            "speed" => speed
        )
    ) for speed in Maple.track_spinner_speeds
)

Ahorn.editingOptions(entity::FlagToggleStarTrackSpinner) = Dict{String, Any}(
    "speed" => Maple.track_spinner_speeds
)

Ahorn.nodeLimits(entity::FlagToggleStarTrackSpinner) = 1, 1

function Ahorn.selection(entity::FlagToggleStarTrackSpinner)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    return [Ahorn.Rectangle(startX - 8, startY - 8, 16, 16), Ahorn.Rectangle(stopX - 8, stopY - 8, 16, 16)]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleStarTrackSpinner, room::Maple.Room)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = entity.data["nodes"][1]

    Ahorn.drawSprite(ctx, "danger/starfish13", stopX, stopY)
    Ahorn.drawArrow(ctx, startX, startY, stopX, stopY, Ahorn.colors.selection_selected_fc, headLength=10)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleStarTrackSpinner, room::Maple.Room)
    startX, startY = Ahorn.position(entity)

    Ahorn.drawSprite(ctx, "danger/starfish13", startX, startY)
end

end


text:
# https://github.com/JuliaLang/julia/issues/27613#issuecomment-400016111
function fargmax(f, itr)
    r = iterate(itr)
    r === nothing && error("empty collection")
    m, state = r
    f_m = f(m)
    while true
        r = iterate(itr, state)
        r === nothing && break
        x, state = r
        f_x = f(x)
        isless(f_m, f_x) || continue
        m, f_m = x, f_x
    end
    return m
end

function fargmin(f, itr)
    r = iterate(itr)
    r === nothing && error("empty collection")
    m, state = r
    f_m = f(m)
    while true
        r = iterate(itr, state)
        r === nothing && break
        x, state = r
        f_x = f(x)
        isless(f_x, f_m) || continue
        m, f_m = x, f_x
    end
    return m
end


text:
function AESDecryptBlock!(result::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key) where {B<:AESBlock}
	cache = gen_cache(key)
	AESDecryptBlock!(result, block, key, cache)
end

function AESDecryptBlock(block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key) where {B<:AESBlock}
	result = zeros(UInt8, 16)
	cache = gen_cache(key)
	AESDecryptBlock!(result, block, key, cache)
	result
end

function AESDecryptBlock(block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key, cache) where {B<:AESBlock}
	result = zeros(UInt8, 16)
	AESDecryptBlock!(result, block, key, cache)
	result
end

function AESDecryptBlock!(result::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key::AES128Key, cache::CipherCache) where {B<:AESBlock}
	current = result
	copyto!(current, block)
	K = cache.K
	tmp = cache.tmp

	for i in 1:4
		for j in 1:4
			K[i][j] = key[4i+j-4]
		end
	end
	# Generate last round key
	for i in 2:11
		copyto!(tmp, K[4])
		for j in 1:4
			tmp[j] = SBOX[tmp[j]+1]
		end
		tmp1 = tmp[1]
		tmp[1] = tmp[2]
		tmp[2] = tmp[3]
		tmp[3] = tmp[4]
		tmp[4] = tmp1
		K[1] .= K[1] .⊻ tmp
		K[1][1] = K[1][1] ⊻ RCON[i - 1]
		K[2] .= K[2] .⊻ K[1]
		K[3] .= K[3] .⊻ K[2]
		K[4] .= K[4] .⊻ K[3]
	end
	# Start Rounds
	for i in 1:11
		# AddRoundKey
		for j in 1:4
			for k in 1:4
				current[4j+k-4] = current[4j+k-4] ⊻ K[j][k]
			end
		end
		if i < 11
			# GenNextRoundKey
			K[4] .= K[4] .⊻ K[3]
			K[3] .= K[3] .⊻ K[2]
			K[2] .= K[2] .⊻ K[1]
			copyto!(tmp, K[4])
			for j in 1:4
				tmp[j] = SBOX[tmp[j]+1]
			end
			tmp1 = tmp[1]
			tmp[1] = tmp[2]
			tmp[2] = tmp[3]
			tmp[3] = tmp[4]
			tmp[4] = tmp1
			K[1] .= K[1] .⊻ tmp
			K[1][1] = K[1][1] ⊻ RCON[11 - i]
			# InvMixColumns
			if i != 1
				for j in 1:4
					tmp[1] = aes_mul(current[4j-3], 0xE) ⊻ aes_mul(current[4j-2], 0xB) ⊻ aes_mul(current[4j-1], 0xD) ⊻ aes_mul(current[4j], 0x9)
					tmp[2] = aes_mul(current[4j-3], 0x9) ⊻ aes_mul(current[4j-2], 0xE) ⊻ aes_mul(current[4j-1], 0xB) ⊻ aes_mul(current[4j], 0xD)
					tmp[3] = aes_mul(current[4j-3], 0xD) ⊻ aes_mul(current[4j-2], 0x9) ⊻ aes_mul(current[4j-1], 0xE) ⊻ aes_mul(current[4j], 0xB)
					tmp[4] = aes_mul(current[4j-3], 0xB) ⊻ aes_mul(current[4j-2], 0xD) ⊻ aes_mul(current[4j-1], 0x9) ⊻ aes_mul(current[4j], 0xE)
					current[4j-3] = tmp[1]
					current[4j-2] = tmp[2]
					current[4j-1] = tmp[3]
					current[4j] = tmp[4]
				end
			end
			# InvSubBytes
			for j in 1:16
				current[j] = INVSBOX[current[j]+1]
			end

			# InvShiftRows
			tmp1 = current[14]
			current[14] = current[10]
			current[10] = current[6]
			current[6] = current[2]
			current[2] = tmp1

			tmp1 = current[3]
			current[3] = current[11]
			current[11] = tmp1
			tmp1 = current[7]
			current[7] = current[15]
			current[15] = tmp1

			tmp1 = current[4]
			current[4] = current[8]
			current[8] = current[12]
			current[12] = current[16]
			current[16] = tmp1
		end
	end
	current
end

function AESDecryptBlock!(result::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key::AES192Key, cache::CipherCache) where {B<:AESBlock}
	current = result
	copyto!(current, block)
	K = cache.K
	tmp = cache.tmp

	for i in 1:6
		for j in 1:4
			K[i][j] = key[4i+j-4]
		end
	end
	# Generate last round key
	rid = 1
	for i in 2:9
			copyto!(tmp, K[6]) #DIFF
			for j in 1:4
				tmp[j] = SBOX[tmp[j]+1]
			end
			tmp1 = tmp[1]
			tmp[1] = tmp[2]
			tmp[2] = tmp[3]
			tmp[3] = tmp[4]
			tmp[4] = tmp1
			K[1] .= K[1] .⊻ tmp
			K[1][1] = K[1][1] ⊻ RCON[rid]
			rid += 1
			K[2] .= K[2] .⊻ K[1]
			K[3] .= K[3] .⊻ K[2]
			K[4] .= K[4] .⊻ K[3]
			K[5] .= K[5] .⊻ K[4] #DIFF
			K[6] .= K[6] .⊻ K[5] #DIFF
	end
	# Start Rounds
	rid = 8
	for i in 1:13
		# AddRoundKey
		if i % 3 == 1
			for j in 1:4
				for k in 1:4
					current[4j+k-4] = current[4j+k-4] ⊻ K[j][k]
				end
			end
		elseif i % 3 == 0
			for j in 1:4
				for k in 1:4
					current[4j+k-4] = current[4j+k-4] ⊻ K[(j + 3) % 6 + 1][k]
				end
			end
		else
			for j in 1:4
				for k in 1:4
					current[4j+k-4] = current[4j+k-4] ⊻ K[j + 2][k]
				end
			end
		end
		if i < 13
			# GenNextRoundKey
			if (i % 3 == 2) || (i % 3 == 1)
				K[6] .= K[6] .⊻ K[5]
				K[5] .= K[5] .⊻ K[4]
			end
			if (i % 3 == 0) || (i % 3 == 1)
				K[4] .= K[4] .⊻ K[3]
				K[3] .= K[3] .⊻ K[2]
				K[2] .= K[2] .⊻ K[1]
				copyto!(tmp, K[6])
				for j in 1:4
					tmp[j] = SBOX[tmp[j]+1]
				end
				tmp1 = tmp[1]
				tmp[1] = tmp[2]
				tmp[2] = tmp[3]
				tmp[3] = tmp[4]
				tmp[4] = tmp1
				K[1] .= K[1] .⊻ tmp
				K[1][1] = K[1][1] ⊻ RCON[rid]
				rid -= 1
			end
			# InvMixColumns
			if i != 1
				for j in 1:4
					tmp[1] = aes_mul(current[4j-3], 0xE) ⊻ aes_mul(current[4j-2], 0xB) ⊻ aes_mul(current[4j-1], 0xD) ⊻ aes_mul(current[4j], 0x9)
					tmp[2] = aes_mul(current[4j-3], 0x9) ⊻ aes_mul(current[4j-2], 0xE) ⊻ aes_mul(current[4j-1], 0xB) ⊻ aes_mul(current[4j], 0xD)
					tmp[3] = aes_mul(current[4j-3], 0xD) ⊻ aes_mul(current[4j-2], 0x9) ⊻ aes_mul(current[4j-1], 0xE) ⊻ aes_mul(current[4j], 0xB)
					tmp[4] = aes_mul(current[4j-3], 0xB) ⊻ aes_mul(current[4j-2], 0xD) ⊻ aes_mul(current[4j-1], 0x9) ⊻ aes_mul(current[4j], 0xE)
					current[4j-3] = tmp[1]
					current[4j-2] = tmp[2]
					current[4j-1] = tmp[3]
					current[4j] = tmp[4]
				end
			end
			# InvSubBytes
			for j in 1:16
				current[j] = INVSBOX[current[j]+1]
			end

			# InvShiftRows
			tmp1 = current[14]
			current[14] = current[10]
			current[10] = current[6]
			current[6] = current[2]
			current[2] = tmp1

			tmp1 = current[3]
			current[3] = current[11]
			current[11] = tmp1
			tmp1 = current[7]
			current[7] = current[15]
			current[15] = tmp1

			tmp1 = current[4]
			current[4] = current[8]
			current[8] = current[12]
			current[12] = current[16]
			current[16] = tmp1
		end
	end
	current
end


text:
# This code was taken from: https://github.com/ViralBShah/MatrixMarket.jl

module MatrixMarket

function mmread(filename::String, infoonly::Bool)
#      Reads the contents of the Matrix Market file 'filename'
#      into a matrix, which will be either sparse or dense,
#      depending on the Matrix Market format indicated by
#      'coordinate' (coordinate sparse storage), or
#      'array' (dense array storage).  The data will be duplicated
#      as appropriate if symmetry is indicated in the header. (Not yet
#      implemented).
#
#      If infoonly is true information on the size and structure is
#      returned.
  mmfile = open(filename,"r")
  tokens = split(chomp(readline(mmfile)))
  if length(tokens) != 5 error("Not enough words on header line") end
  if tokens[1] != "%%MatrixMarket" error("Not a valid MatrixMarket header.") end
  (head1, rep, field, symm) = map(lowercase, tokens[2:5])
  if head1 != "matrix"
      error("This seems to be a MatrixMarket $head1 file, not a MatrixMarket matrix file")
  end
  if field != "real" error("non-float fields not yet allowed") end

  ll   = readline(mmfile)         # Read through comments, ignoring them
  while length(ll) > 0 && ll[1] == '%' ll = readline(mmfile) end
  dd     = int(split(ll))         # Read dimensions
  rows   = dd[1]
  cols   = dd[2]
  entries = rep == "coordinate" ? dd[3] : rows * cols
  if infoonly return rows, cols, entries, rep, field, symm end
  if rep == "coordinate"
    rr = Array(Int, entries)
    cc = Array(Int, entries)
    xx = Array(Float64, entries)
    for i in 1:entries
        flds = split(readline(mmfile))
        rr[i] = int32(flds[1])
        cc[i] = int32(flds[2])
        xx[i] = float64(flds[3])
    end
    return sparse(rr, cc, xx, rows, cols)
  end
  reshape([float64(readline(mmfile)) for i in 1:entries], (rows,cols))
end

mmread(filename::String) = mmread(filename, false)

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ResourceMetricSource indicates how to scale on a resource metric known to Kubernetes, as specified in requests and limits, describing each pod in the current scale target (e.g. CPU or memory).  The values will be averaged together before being compared to the target.  Such metrics are built in to Kubernetes, and have special scaling options on top of those available to normal per-pod metrics using the \"pods\" source.  Only one \"target\" type should be set.

    IoK8sApiAutoscalingV2beta1ResourceMetricSource(;
        name=nothing,
        targetAverageUtilization=nothing,
        targetAverageValue=nothing,
    )

    - name::String : name is the name of the resource in question.
    - targetAverageUtilization::Int32 : targetAverageUtilization is the target value of the average of the resource metric across all relevant pods, represented as a percentage of the requested value of the resource for the pods.
    - targetAverageValue::IoK8sApimachineryPkgApiResourceQuantity : targetAverageValue is the target value of the average of the resource metric across all relevant pods, as a raw value (instead of as a percentage of the request), similar to the \"pods\" metric source type.
"""
mutable struct IoK8sApiAutoscalingV2beta1ResourceMetricSource <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    targetAverageUtilization::Any # spec type: Union{ Nothing, Int32 } # spec name: targetAverageUtilization
    targetAverageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: targetAverageValue

    function IoK8sApiAutoscalingV2beta1ResourceMetricSource(;name=nothing, targetAverageUtilization=nothing, targetAverageValue=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricSource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricSource, Symbol("targetAverageUtilization"), targetAverageUtilization)
        setfield!(o, Symbol("targetAverageUtilization"), targetAverageUtilization)
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricSource, Symbol("targetAverageValue"), targetAverageValue)
        setfield!(o, Symbol("targetAverageValue"), targetAverageValue)
        o
    end
end # type IoK8sApiAutoscalingV2beta1ResourceMetricSource

const _property_map_IoK8sApiAutoscalingV2beta1ResourceMetricSource = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("targetAverageUtilization")=>Symbol("targetAverageUtilization"), Symbol("targetAverageValue")=>Symbol("targetAverageValue"))
const _property_types_IoK8sApiAutoscalingV2beta1ResourceMetricSource = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("targetAverageUtilization")=>"Int32", Symbol("targetAverageValue")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1ResourceMetricSource))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1ResourceMetricSource[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1ResourceMetricSource[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1ResourceMetricSource)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }, name::Symbol, val)
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Andrey Copat, Estefane S. Horn, Marcos L. Alencastro (Revised Gerson B. Bicca)
#* $Id$
#*--------------------------------------------------------------------
type centrifugal_pump
	centrifugal_pump()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Molar Weight"
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pump Efficiency",
				:Default=>0.75,
				:Lower=>1E-3
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Mechanical efficiency",
				:Default=>0.95,
				:Lower=>1E-3
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"NPSH Options",
				:Valid=>["Default","Include Kinetic Head"],
				:Default=>"Default"
			)),
			acceleration(Dict{Symbol,Any}(
				:Brief=>"Gravity Acceleration",
				:Default=>9.81,
				:Hidden=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle Suction Area",
				:Default=>0.001
			)),
			flow_vol(Dict{Symbol,Any}(
				:Brief=>"Volumetric Flow Rate" ,
				:Protected=>true
			)),
			flow_mass(Dict{Symbol,Any}(
				:Brief=>"Inlet Mass Flow Rate" ,
				:Protected=>true
			)),
			flow_mass(Dict{Symbol,Any}(
				:Brief=>"Outlet Mass Flow Rate",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at inlet conditions",
				:Lower=>1E-6,
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at outlet conditions",
				:Lower=>1E-6,
				:Protected=>true
			)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Mixture Molar Weight" ,
				:Protected=>true
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Mixture Vapour Pressure" ,
				:Protected=>true
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Fluid Power"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Brake Power"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Eletrical Potency"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Symbol=>"P_{ratio}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Increase",
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"P_{incr}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Static Head"
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Actual Head",
				:Protected=>true
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Isentropic Head",
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Available Net Positive Suction Head" ,
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Velocity Head",
				:Protected=>true
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Velocity Inlet Nozzle",
				:Hidden=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.4025,
				:Symbol=>"_{in}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.20,
				:Symbol=>"_{out}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Work Inlet",
				:PosX=>0.5,
				:PosY=>1,
				:Protected=>true
			)),
			[
				:(Fvol = NozzleVelocity*SuctionArea),
				:(VelocityHead = 0.5*NozzleVelocity^2/g),
				:(Mwm = sum(Mw*Inlet.z)),
				:(rho_in = PP.LiquidDensity(Inlet.T, Inlet.P, Inlet.z)),
				:(rho_out= PP.LiquidDensity(Outlet.T, Outlet.P, Outlet.z)),
				:(Fw_in = Mwm*Inlet.F),
				:(Fw_out = Fw_in),
				:(Outlet.P = Inlet.P + Pincrease),
				:(Pvapor = PP.BubbleP(Inlet.T,Inlet.z)),
				:(Outlet.P = Inlet.P * Pratio),
				:(Outlet.P = Inlet.P - Pdrop),
				:(HeadIsentropic = -Pdrop/rho_in),
				:(Head = HeadIsentropic/PumpEfficiency),
				:(Head*Mwm = (Outlet.h-Inlet.h)),
				:(FluidPower = HeadIsentropic *Mwm* Inlet.F),
				:(BrakePower * PumpEfficiency = FluidPower),
				:(EletricPower = -WorkIn),
				:(BrakePower = EletricPower * MechanicalEff),
				:(Outlet.F = Inlet.F),
				:(Outlet.z = Inlet.z),
				:(Fvol = Fw_in/rho_in),
				:(NPSH_available = (Inlet.P - Pvapor)/(rho_in*g) + StaticHead),
				:(NPSH_available = (Inlet.P - Pvapor)/(rho_in*g)+VelocityHead+StaticHead),
			],
			[
				"Velocity Inlet Nozzle","Velocity Head","Average Molecular Weight","Mass Density at inlet conditions","Mass Density at outlet conditions","Inlet Flow Mass","Outlet Flow Mass","Pressure Increase","Mixture Vapour Pressure","Pressure Ratio","Pressure Drop","Isentropic Head","Pump Efficiency","Actual Head","Fluid Power","Brake Power","Eletric Power","Eletric Power","Molar Balance","Outlet Composition","Volumetric Flow Rate","Net Positive Suction Head Available - Without Velocity Head","Net Positive Suction Head Available - Included Velocity Head",
			],
			[:PP,:NComp,:Mw,:PumpEfficiency,:MechanicalEff,:NPSH_Options,:g,:SuctionArea,],
			[:Fvol,:Fw_in,:Fw_out,:rho_in,:rho_out,:Mwm,:Pvapor,:FluidPower,:BrakePower,:EletricPower,:Pratio,:Pdrop,:Pincrease,:StaticHead,:Head,:HeadIsentropic,:NPSH_available,:VelocityHead,:NozzleVelocity,:Inlet,:Outlet,:WorkIn,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Mw::Array{molweight}
	PumpEfficiency::positive
	MechanicalEff::positive
	NPSH_Options::DanaSwitcher
	g::acceleration
	SuctionArea::area
	Fvol::flow_vol
	Fw_in::flow_mass
	Fw_out::flow_mass
	rho_in::dens_mass
	rho_out::dens_mass
	Mwm::molweight
	Pvapor::pressure
	FluidPower::power
	BrakePower::power
	EletricPower::power
	Pratio::positive
	Pdrop::press_delta
	Pincrease::press_delta
	StaticHead::length
	Head::energy_mass
	HeadIsentropic::energy_mass
	NPSH_available::length
	VelocityHead::length
	NozzleVelocity::velocity
	Inlet::stream
	Outlet::streamPH
	WorkIn::power
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export centrifugal_pump
function set(in::centrifugal_pump)
	Mw = PP.MolecularWeight()
	 g = 9.81*"m/s^2"

end
function setEquationFlow(in::centrifugal_pump)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	let switch=NPSH_Options
		if switch=="Default"
			addEquation(22)
		elseif switch=="Include Kinetic Head"
			addEquation(23)
		end
	end
end
function atributes(in::centrifugal_pump,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Pump"
	fields[:Brief]="Model of a centrifugal pump."
	fields[:Info]="== Assumptions ==
* Steady State;
* Only Liquid;
* Adiabatic;
* Isentropic.
"
	drive!(fields,_)
	return fields
end
centrifugal_pump(_::Dict{Symbol,Any})=begin
	newModel=centrifugal_pump()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#####
##### Samples
#####

"""
    VALIDATE_SAMPLES_DEFAULT[]

Defaults to `true`.

When set to `true`, `Samples` objects will be validated upon construction for compliance
with the Onda specification.

Users may interactively set this reference to `false` in order to disable this extra layer
validation, which can be useful when working with malformed Onda datasets.

See also: [`Onda.validate`](@ref)
"""
const VALIDATE_SAMPLES_DEFAULT = Ref{Bool}(true)

"""
    Samples(data::AbstractMatrix, info::SamplesInfo, encoded::Bool;
            validate::Bool=Onda.VALIDATE_SAMPLES_DEFAULT[])

Return a `Samples` instance with the following fields:

- `data::AbstractMatrix`: A matrix of sample data. The `i` th row of the matrix
  corresponds to the `i`th channel in `info.channels`, while the `j`th
  column corresponds to the `j`th multichannel sample.

- `info::SamplesInfo`: The `SamplesInfo` object that describes the `Samples` instance.

- `encoded::Bool`: If `true`, the values in `data` are LPCM-encoded as prescribed
  by the `Samples` instance's `info`. If `false`, the values in `data` have
  been decoded into the `info`'s canonical units.

If `validate` is `true`, [`Onda.validate`](@ref) is called on the constructed `Samples`
instance before it is returned.

Note that `getindex` and `view` are defined on `Samples` to accept normal integer
indices, but also accept channel names or a regex to match channel names for row indices,
and `TimeSpan` values for column indices; see `Onda/examples/tour.jl` for a comprehensive
set of indexing examples.

See also: [`load`](@ref), [`store`](@ref), [`encode`](@ref), [`encode!`](@ref), [`decode`](@ref), [`decode!`](@ref)
"""
struct Samples{D<:AbstractMatrix,S<:SamplesInfo}
    data::D
    info::S
    encoded::Bool
    function Samples(data, info::SamplesInfo, encoded::Bool;
                     validate::Bool=VALIDATE_SAMPLES_DEFAULT[])
        samples = new{typeof(data),typeof(info)}(data, info, encoded)
        validate && Onda.validate(samples)
        return samples
    end
end

"""
    ==(a::Samples, b::Samples)

Returns `a.encoded == b.encoded && a.info == b.info && a.data == b.data`.
"""
Base.:(==)(a::Samples, b::Samples) = a.encoded == b.encoded && a.info == b.info && a.data == b.data

"""
    validate(samples::Samples)

Returns `nothing`, checking that the given `samples` are valid w.r.t. the
underlying `samples.info` and the Onda specification's canonical LPCM
representation. If a violation is found, an `ArgumentError` is thrown.

Properties that are validated by this function include:

- encoded element type matches `sample_type(samples.info)`
- the number of rows of `samples.data` matches the number of channels in `samples.info`
"""
function validate(samples::Samples)
    n_channels = channel_count(samples.info)
    n_rows = size(samples.data, 1)
    if n_channels != n_rows
        throw(ArgumentError("number of channels in info ($n_channels) " *
                            "does not match number of rows in data matrix " *
                            "($n_rows)"))
    end
    if samples.encoded && !(eltype(samples.data) === sample_type(samples.info))
        throw(ArgumentError("encoded `samples.data` matrix eltype does not match `sample_type(samples.info)`"))
    end
    return nothing
end

TimeSpans.istimespan(::Samples) = true
TimeSpans.start(::Samples) = Nanosecond(0)
TimeSpans.stop(samples::Samples) = TimeSpans.time_from_index(samples.info.sample_rate, size(samples.data, 2) + 1)

"""
    channel(samples::Samples, name)

Return `channel(samples.info, name)`.

This function is useful for indexing rows of `samples.data` by channel names.
"""
channel(samples::Samples, name) = channel(samples.info, name)

"""
    channel(samples::Samples, i::Integer)

Return `channel(samples.info, i)`.
"""
channel(samples::Samples, i::Integer) = channel(samples.info, i)

"""
    channel_count(samples::Samples)

Return `channel_count(samples.info)`.
"""
channel_count(samples::Samples) = channel_count(samples.info)

"""
    sample_count(samples::Samples)

Return the number of multichannel samples in `samples` (i.e. `size(samples.data, 2)`)
"""
sample_count(samples::Samples) = size(samples.data, 2)

#####
##### indexing
#####

for f in (:getindex, :view)
    @eval begin
        @inline function Base.$f(samples::Samples, rows, columns)
            rows = row_arguments(samples, rows)
            columns = column_arguments(samples, columns)
            channels = rows isa Colon ? samples.info.channels : samples.info.channels[rows]
            info = @compat SamplesInfo(Tables.rowmerge(samples.info; channels))
            return Samples($f(samples.data, rows, columns), info, samples.encoded; validate=false)
        end
    end
end

_rangify(i) = i
_rangify(i::Integer) = i:i

_indices_fallback(::Any, ::Samples, i::Union{Colon,AbstractRange,Integer}) = i
_indices_fallback(f, samples::Samples, x) = map(x -> f(samples, x), x)

row_arguments(samples::Samples, x) = _rangify(_row_arguments(samples, x))

_row_arguments(samples::Samples, x) = _indices_fallback(_row_arguments, samples, x)
_row_arguments(samples::Samples, name::AbstractString) = channel(samples, name)
_row_arguments(samples::Samples, name::Regex) = findall(c -> match(name, c) !== nothing, samples.info.channels)

column_arguments(samples::Samples, x) = _rangify(_column_arguments(samples, x))

function _column_arguments(samples::Samples, x)
    TimeSpans.istimespan(x) && return TimeSpans.index_from_time(samples.info.sample_rate, TimeSpan(x))
    return _indices_fallback(_column_arguments, samples, x)
end

#####
##### encoding utilities
#####

function encode_sample(::Type{S}, resolution_in_unit, offset_in_unit, sample_in_unit,
                       noise=zero(sample_in_unit)) where {S<:LPCM_SAMPLE_TYPE_UNION}
    sample_in_unit += noise
    isnan(sample_in_unit) && S <: Integer && return typemax(S)
    from_unit = clamp((sample_in_unit - offset_in_unit) / resolution_in_unit, typemin(S), typemax(S))
    return S <: Integer ? round(S, from_unit) : from_unit
end

function dither_noise!(rng::AbstractRNG, storage, step)
    rand!(rng, storage)
    broadcast!(_dither_noise, storage, storage, step + step, step)
    return storage
end

dither_noise!(storage, step) = dither_noise!(Random.GLOBAL_RNG, storage, step)

function _dither_noise(x, range, step)
    rs = range * x
    if rs < step
        return sqrt(rs * step) - step
    else
        return step - sqrt(range * (1 - x) * step)
    end
end

#####
##### `encode`/`encode!`
#####

"""
    encode(sample_type::DataType, sample_resolution_in_unit, sample_offset_in_unit,
           sample_data, dither_storage=nothing)

Return a copy of `sample_data` quantized according to `sample_type`, `sample_resolution_in_unit`,
and `sample_offset_in_unit`. `sample_type` must be a concrete subtype of `Onda.VALID_SAMPLE_TYPE_UNION`.
Quantization of an individual sample `s` is performed via:

    round(S, (s - sample_offset_in_unit) / sample_resolution_in_unit)

with additional special casing to clip values exceeding the encoding's dynamic range.

If `dither_storage isa Nothing`, no dithering is applied before quantization.

If `dither_storage isa Missing`, dither storage is allocated automatically and
triangular dithering is applied to the info prior to quantization.

Otherwise, `dither_storage` must be a container of similar shape and type to
`sample_data`. This container is then used to store the random noise needed for the
triangular dithering process, which is applied to the info prior to quantization.

If:

    sample_type === eltype(sample_data) &&
    sample_resolution_in_unit == 1 &&
    sample_offset_in_unit == 0

then this function will simply return `sample_data` directly without copying/dithering.
"""
function encode(::Type{S}, sample_resolution_in_unit, sample_offset_in_unit,
                sample_data, dither_storage=nothing) where {S<:LPCM_SAMPLE_TYPE_UNION}
    if (S === eltype(sample_data) &&
        sample_resolution_in_unit == 1 &&
        sample_offset_in_unit == 0)
        return sample_data
    end
    return encode!(similar(sample_data, S), S,
                   sample_resolution_in_unit, sample_offset_in_unit,
                   sample_data, dither_storage)
end

"""
    encode!(result_storage, sample_type::DataType, sample_resolution_in_unit,
            sample_offset_in_unit, sample_data, dither_storage=nothing)
    encode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit,
            sample_data, dither_storage=nothing)

Similar to `encode(sample_type, sample_resolution_in_unit, sample_offset_in_unit, sample_data, dither_storage)`,
but write encoded values to `result_storage` rather than allocating new storage.

`sample_type` defaults to `eltype(result_storage)` if it is not provided.

If:

    sample_type === eltype(sample_data) &&
    sample_resolution_in_unit == 1 &&
    sample_offset_in_unit == 0

then this function will simply copy `sample_data` directly into `result_storage` without dithering.
"""
function encode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit,
                 sample_data, dither_storage=nothing)
    return encode!(result_storage, eltype(result_storage), sample_resolution_in_unit,
                   sample_offset_in_unit, sample_data, dither_storage)
end

function encode!(result_storage, ::Type{S}, sample_resolution_in_unit, sample_offset_in_unit,
                 sample_data, dither_storage=nothing) where {S<:LPCM_SAMPLE_TYPE_UNION}
    if (S === eltype(sample_data) &&
        sample_resolution_in_unit == 1 &&
        sample_offset_in_unit == 0)
        copyto!(result_storage, sample_data)
    else
        if dither_storage isa Nothing
            broadcast!(encode_sample, result_storage, S,
                       sample_resolution_in_unit,
                       sample_offset_in_unit, sample_data)
        else
            if dither_storage isa Missing
                dither_storage = similar(sample_data)
            elseif size(dither_storage) != size(sample_data)
                throw(DimensionMismatch("dithering storage container does not match shape of sample_data"))
            end
            dither_noise!(dither_storage, sample_resolution_in_unit)
            broadcast!(encode_sample, result_storage, S,
                       sample_resolution_in_unit, sample_offset_in_unit,
                       sample_data, dither_storage)
        end
    end
    return result_storage
end

"""
    encode(samples::Samples, dither_storage=nothing)

If `samples.encoded` is `false`, return a `Samples` instance that wraps:

    encode(sample_type(samples.info),
           samples.info.sample_resolution_in_unit,
           samples.info.sample_offset_in_unit,
           samples.data, dither_storage)

If `samples.encoded` is `true`, this function is the identity.
"""
function encode(samples::Samples, dither_storage=nothing)
    samples.encoded && return samples
    return Samples(encode(sample_type(samples.info),
                          samples.info.sample_resolution_in_unit,
                          samples.info.sample_offset_in_unit,
                          samples.data, dither_storage),
                   samples.info, true; validate=false)
end

"""
    encode!(result_storage, samples::Samples, dither_storage=nothing)

If `samples.encoded` is `false`, return a `Samples` instance that wraps:

    encode!(result_storage,
            sample_type(samples.info),
            samples.info.sample_resolution_in_unit,
            samples.info.sample_offset_in_unit,
            samples.data, dither_storage)`.

If `samples.encoded` is `true`, return a `Samples` instance that wraps
`copyto!(result_storage, samples.data)`.
"""
function encode!(result_storage, samples::Samples, dither_storage=nothing)
    if samples.encoded
        copyto!(result_storage, samples.data)
    else
        encode!(result_storage, sample_type(samples.info),
                samples.info.sample_resolution_in_unit,
                samples.info.sample_offset_in_unit,
                samples.data, dither_storage)
    end
    return Samples(result_storage, samples.info, true; validate=false)
end

#####
##### `decode`/`decode!`
#####

"""
    decode(sample_resolution_in_unit, sample_offset_in_unit, sample_data)

Return `sample_resolution_in_unit .* sample_data .+ sample_offset_in_unit`.

If:

    sample_data isa AbstractArray &&
    sample_resolution_in_unit == 1 &&
    sample_offset_in_unit == 0

then this function is the identity and will return `sample_data` directly without copying.
"""
function decode(sample_resolution_in_unit, sample_offset_in_unit, sample_data)
    if sample_data isa AbstractArray
        sample_resolution_in_unit == 1 && sample_offset_in_unit == 0 && return sample_data
    end
    return sample_resolution_in_unit .* sample_data .+ sample_offset_in_unit
end

"""
    decode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit, sample_data)

Similar to `decode(sample_resolution_in_unit, sample_offset_in_unit, sample_data)`, but
write decoded values to `result_storage` rather than allocating new storage.
"""
function decode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit, sample_data)
    f = x -> sample_resolution_in_unit * x + sample_offset_in_unit
    return broadcast!(f, result_storage, sample_data)
end

"""
    decode(samples::Samples)

If `samples.encoded` is `true`, return a `Samples` instance that wraps

    decode(samples.info.sample_resolution_in_unit, samples.info.sample_offset_in_unit, samples.data)

If `samples.encoded` is `false`, this function is the identity.
"""
function decode(samples::Samples)
    samples.encoded || return samples
    return Samples(decode(samples.info.sample_resolution_in_unit,
                          samples.info.sample_offset_in_unit,
                          samples.data),
                   samples.info, false; validate=false)
end

"""
    decode!(result_storage, samples::Samples)

If `samples.encoded` is `true`, return a `Samples` instance that wraps

    decode!(result_storage, samples.info.sample_resolution_in_unit, samples.info.sample_offset_in_unit, samples.data)

If `samples.encoded` is `false`, return a `Samples` instance that wraps `copyto!(result_storage, samples.data)`.
"""
function decode!(result_storage, samples::Samples)
    if samples.encoded
        decode!(result_storage, samples.info.sample_resolution_in_unit,
                samples.info.sample_offset_in_unit, samples.data)
    else
        copyto!(result_storage, samples.data)
    end
    return Samples(result_storage, samples.info, false; validate=false)
end

#####
##### load/store
#####

"""
    load(signal[, span]; encoded::Bool=false)
    load(file_path, file_format::Union{AbstractString,AbstractLPCMFormat}, info::SamplesInfo[, span]; encoded::Bool=false)

Return the `Samples` object described by `signal`/`file_path`/`file_format`/`info`.

If `span` is present, return `load(...)[:, span]`, but attempt to avoid reading
unreturned intermediate sample data. Note that the effectiveness of this optimized method
versus the naive approach depends on the types of `file_path` (i.e. if there is a fast
method defined for `Onda.read_byte_range(::typeof(file_path), ...)`) and `file_format`
(i.e. does the corresponding format support random or chunked access).

If `encoded` is `true`, do not decode the `Samples` object before returning it.
"""
function load(signal, span...; encoded::Bool=false)
    return @compat load(signal.file_path, signal.file_format, extract_samples_info(signal), span...; encoded)
end

function load(file_path, file_format::AbstractString, info::SamplesInfo, span...; encoded::Bool=false)
    return @compat load(file_path, format(file_format, info), info, span...; encoded)
end

function load(file_path, file_format::AbstractLPCMFormat, info::SamplesInfo; encoded::Bool=false)
    samples = Samples(read_lpcm(file_path, file_format), info, true)
    return encoded ? samples : decode(samples)
end

function load(file_path, file_format::AbstractLPCMFormat, info::SamplesInfo, span; encoded::Bool=false)
    sample_range = TimeSpans.index_from_time(info.sample_rate, span)
    sample_offset, sample_count = first(sample_range) - 1, length(sample_range)
    sample_data = read_lpcm(file_path, file_format, sample_offset, sample_count)
    samples = Samples(sample_data, info, true)
    return encoded ? samples : decode(samples)
end

"""
    Onda.mmap(signal)

Return `Onda.mmap(signal.file_path, Onda.extract_samples_info(signal))`, throwing an `ArgumentError` if `signal.file_format != "lpcm"`.
"""
function mmap(signal)
    signal.file_format == "lpcm" || throw(ArgumentError("unsupported file_format for mmap: $(signal.file_format)"))
    return mmap(signal.file_path, extract_samples_info(signal))
end

"""
    Onda.mmap(mmappable, info::SamplesInfo)

Return `Samples(data, info, true)` where `data` is created via `Mmap.mmap(mmappable, ...)`.

`mmappable` is assumed to reference memory that is formatted according to the Onda Format's canonical
interleaved LPCM representation in accordance with `sample_type(info)` and `channel_count(info)`. No
explicit checks are performed to ensure that this is true.
"""
function mmap(mmappable, info::SamplesInfo)
    data = reshape(Mmap.mmap(mmappable, Vector{sample_type(info)}), (channel_count(info), :))
    return Samples(data, info, true)
end

"""
    store(file_path, file_format::Union{AbstractString,AbstractLPCMFormat}, samples::Samples)

Serialize the given `samples` to `file_format` and write the output to `file_path`.
"""
function store(file_path, file_format, samples::Samples)
    fmt = file_format isa AbstractLPCMFormat ? file_format : format(file_format, samples.info)
    return write_lpcm(file_path, fmt, encode(samples).data)
end

"""
    store(file_path, file_format::Union{AbstractString,AbstractLPCMFormat}, samples::Samples,
          recording::UUID, start::Period; custom...)

Serialize the given `samples` to `file_format` and write the output to `file_path`, returning
a `Signal` instance constructed from the provided arguments (any provided `custom` keyword
arguments are forwarded to an invocation of the `Signal` constructor).
"""
function store(file_path, file_format, samples::Samples, recording, start; custom...)
    store(file_path, file_format, samples)
    span = TimeSpan(start, Nanosecond(start) + TimeSpans.duration(samples))
    return @compat Signal(Tables.rowmerge(samples.info; recording, file_path, file_format, span, custom...))
end

#####
##### pretty printing
#####

function Base.show(io::IO, samples::Samples)
    if get(io, :compact, false)
        print(io, "Samples(", summary(samples.data), ')')
    else
        duration_in_seconds = size(samples.data, 2) / samples.info.sample_rate
        duration_in_nanoseconds = round(Int, duration_in_seconds * 1_000_000_000)
        println(io, "Samples (", TimeSpans.format_duration(duration_in_nanoseconds), "):")
        println(io, "  info.kind: ", repr(samples.info.kind))
        println(io, "  info.channels: ", string('[', join(map(repr, samples.info.channels), ", "), ']'))
        println(io, "  info.sample_unit: ", repr(samples.info.sample_unit))
        println(io, "  info.sample_resolution_in_unit: ", samples.info.sample_resolution_in_unit)
        println(io, "  info.sample_offset_in_unit: ", samples.info.sample_offset_in_unit)
        println(io, "  sample_type(info): ", sample_type(samples.info))
        println(io, "  info.sample_rate: ", samples.info.sample_rate, " Hz")
        println(io, "  encoded: ", samples.encoded)
        println(io, "  data:")
        show(io, "text/plain", samples.data)
    end
end

#####
##### Arrow conversion
#####

const SamplesArrowType{T,S} = NamedTuple{(:data, :info, :encoded),Tuple{Vector{T},S,Bool}}

const SAMPLES_ARROW_NAME = Symbol("JuliaLang.Samples")

Arrow.ArrowTypes.arrowname(::Type{<:Samples}) = SAMPLES_ARROW_NAME

Arrow.ArrowTypes.ArrowType(::Type{<:Samples{D,S}}) where {D,S} = SamplesArrowType{eltype(D),Arrow.ArrowTypes.ArrowType(S)}

function Arrow.ArrowTypes.toarrow(samples::Samples)
    return (data=vec(samples.data),
            info=Arrow.ArrowTypes.toarrow(samples.info),
            encoded=samples.encoded)
end

function Arrow.ArrowTypes.JuliaType(::Val{SAMPLES_ARROW_NAME}, ::Type{<:SamplesArrowType{T}}) where {T}
    return Samples{Matrix{T},<:SamplesInfo}
end

function Arrow.ArrowTypes.fromarrow(::Type{<:Samples}, arrow_data, arrow_info, arrow_encoded)
    info = Arrow.ArrowTypes.fromarrow(SamplesInfo, arrow_info...)
    data = reshape(arrow_data, (channel_count(info), :))
    return Samples(data, info, arrow_encoded)
end


text: past(h, m, s) = 1000(60(60h+m)+s)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Test
using MolecularConformation
using DelimitedFiles
using LinearAlgebra

include("testmainfunction.jl")


text:
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.LinearSolvers
using MUMPSjInv
using MAT
using LinearAlgebra
using SparseArrays

# filename= "2DmodelLShaped.mat"
filename= joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples","Sources3D")
file = matread(filename * ".mat")

# problem description
domain = file["domain"]
mf     = round.(Int64,file["m"])
v      = file["v"]
sig    = file["sig"]

# data and receiver location
dtrue  = file["dtrue"]
rec    = file["rec"]

Mfine = getRegularMesh(domain,mf);
x1c,x2c,x3c = getCellCenteredAxes(Mfine)
rec3D = [kron(ones(length(x3c)),rec) kron(x3c,ones(size(rec,1)))]

# set up experiments
Meshes       = [ getRegularMesh(domain,[96 48 48]);];
alphas       = exp10.(range(0,stop=-6,length=30))
regularizersName = ["wTVReg"];
regularizers = [wTVReg];
noiseLevel = 0.1;

for mk = 1:length(Meshes)
	M      = Meshes[mk]
	x1,x2,x3 = getNodalAxes(M)
	P = interpmat(x1,x2,x3, rec3D);
	for rk = 1:length(regularizers)
		regFun= regularizers[rk]
		regName= regularizersName[rk]

		# build linear interpolation matrix from nodes to receiver locations

		dnoise = randn(size(dtrue));
		dnoise /= norm(dnoise)
		dnoise *= noiseLevel*norm(dtrue)
		dobs = dtrue + dnoise

		pFor = getConvDiffFEMParam(M,v=v,sig=sig,P=P,Ainv=getMUMPSsolver());

		## configure misfit
		Wt         = ones(size(dobs))/sqrt(mf[3])
		sigback    = 0.0
		pMis       = getMisfitParam(pFor,Wt,dobs,SSDFun)

		## Configure regularization
		mref       = zeros(M.nc)
		alpha      = alphas[1]
		reg        = (m,mr,M,I=1.0) -> regFun(m,mr,M,eps=1e-8)
		# reg        = (m,mr,M,I=1.0) -> wTVReg(m,mr,M,eps=1e-8)

		## Configure optimization
		maxIter    = 20
		minUpdate  = 1e-3
		HesPrec    = getSSORRegularizationPreconditioner(1.0,1e-15,50)
		cgit       = 5
		pcgTol     = 1e-1
		modFun     = identityMod
		boundsLow  = 0*ones(M.nc)
		boundsHigh = 1*ones(M.nc)
		maxStep	   = 0.1*maximum(boundsHigh)

		## store the configuration
		pInv       = getInverseParam(M,modFun,reg,alpha,mref,
		                             boundsLow,boundsHigh,maxStep=maxStep,
		                            pcgMaxIter=cgit,pcgTol=pcgTol,minUpdate=minUpdate,maxIter=maxIter,
		                            HesPrec=HesPrec);

		mc = mref.+0.1

		Mis    = zeros(length(alphas))
		Reg    = zeros(length(alphas))
		Fields = zeros(size(P,2),length(alphas))
		Sources = zeros(length(mc),length(alphas))
	    PredictedData   = zeros(size(P,1),length(alphas))

		nPDEsolves = 0;

		runtimeLcurve = @elapsed begin
			for ak = 1:length(alphas)
				pInv.alpha = alphas[ak]
				println("\n\n--- Inversion on $(M.n) grid with $(regName) and alpha=$(pInv.alpha) --- ")
				# mc = max.(.01,min.(.99,mc))
				mc,Dc,flag,His = projGNCG(mc,pInv,pMis)
				uc = pFor.Fields;          # u variables
				wr = pInv.modelfun(mc)[1]; # w variables

				Fields[:,ak] = uc;
				Sources[:,ak] = wr
				PredictedData[:,ak] = Dc

				iter = findlast(His.F .> 0)
				Mis[ak] = His.F[iter]
				Reg[ak] = His.Rc[iter]/alphas[ak]
		    end
		end
		println("no. PDE solves:\t$(pFor.Ainv.nSolve)")
		println("time PDE factor:\t$(pFor.Ainv.facTime)")
		println("time PDE solves:\t$(pFor.Ainv.solveTime)")
		println("\n overall runtime: $runtimeLcurve")

		results = Dict("n"=>M.n, "Mis"=>Mis, "Reg"=>Reg, "alphas"=> alphas,
		"noiseLevel"=>noiseLevel,"Fields"=> Fields, "Sources" => Sources,
		"PredictedData" => PredictedData, "dobs" => dobs,"runtime"=>runtimeLcurve)
		resFile = "$(filename)-$regName-noise-$(noiseLevel)-$(M.n[1])x$(M.n[2])x$(M.n[3]).mat"
		matwrite(resFile,results)
	end
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test, SparseArrays, LinearAlgebra, Statistics

function get_stress_tensor(element, ip, time)
    haskey(element, "displacement") || return nothing
    gradu = element("displacement", ip, time, Val{:Grad})
    eps = 0.5*(gradu' + gradu)
    E = element("youngs modulus", ip, time)
    nu = element("poissons ratio", ip, time)
    mu = E/(2.0*(1.0+nu))
    la = E*nu/((1.0+nu)*(1.0-2.0*nu))
    S = la*tr(eps)*I + 2.0*mu*eps
    return S
end

""" Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """
function get_stress(element, ip, time)
    S = get_stress_tensor(element, ip, time)
    return [S[1,1], S[2,2], S[3,3], S[1,2], S[2,3], S[1,3]]
end

""" Return principal stresses. """
function get_stress_principal(element, ip, time)
    S = get_stress_tensor(element, ip, time)
    return sort(eigvals(S))
end

""" Make least squares fit for some field to nodes. """
function lsq_fit(elements, field, time)
    A = SparseMatrixCOO()
    b = SparseMatrixCOO()
    volume = 0.0
    for element in elements
        gdofs = get_connectivity(element)
        # increase integration order by 1 from default
        for ip in get_integration_points(element, 1)
            detJ = element(ip, time, Val{:detJ})
            w = ip.weight*detJ
            N = element(ip, time)
            f = field(element, ip, time)
            add!(A, gdofs, gdofs, w*kron(N', N))
            for i=1:length(f)
                add!(b, gdofs, w*f[i]*N, i)
            end
            volume += w
        end
    end
    @info("Mass matrix for least-squares fit is assembled. Total volume to fit: $volume")
    A = sparse(A)
    b = sparse(b)
    A = 1/2*(A + A')

    SparseArrays.droptol!(A, 1.0e-6)
    SparseArrays.dropzeros!(A)
    nz = get_nonzero_rows(A)
    F = ldlt(A[nz,nz])

    x = zeros(size(b)...)
    x[nz, :] = F \ b[nz, :]

    nodal_values = Dict(i => vec(x[i,:]) for i in nz)
    return nodal_values
end

#=
test subjects:
- surface pressure load in curved surface
- verification of elements wedge6 and wedge15
- stress interpolation from gauss points to nodes

from Code Aster, for linear model
 N38      -8.85861895037377E-01 -3.46944695195361E-18 -3.46944695195361E-18
coords of N38 = (1.0, 0.0, 0.0)
 N49      -4.50934684240566E-01 -4.46908405333021E-01 -5.92329377847994E-01
coords of N49 = (0.521002, 0.515482, 0.68032) # quite middle of surface

Analytical solution
uᵣ(r) = b³p/(2Er²(a³-b³)) * (a³(ν+1) + r³(-4ν+2)), where
a = inner surface radial distance, b = outer surface ...

if a=0.9, b=1.0, ν=1/3, E = 24580 and p = 7317 equation yields
-9/10 for radial displacement

http://mms2.ensmp.fr/emms_paris/plasticite3D/exercices/eSpherePress.pdf
=#
function test_wedge_sphere(model, u_CA, S_CA)
    mesh_file = @__DIR__() * "/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, model)

    body = Problem(Elasticity, "hollow sphere 1/8 model", 3)
    body_elements = create_elements(mesh, "HOLLOWSPHERE8")
    update!(body_elements, "youngs modulus", 24580.0)
    update!(body_elements, "poissons ratio", 1/3)
    add_elements!(body, body_elements)

    bc = Problem(Dirichlet, "symmetry bc", 3, "displacement")
    el1 = create_elements(mesh, "FACE1")
    update!(el1, "displacement 3", 0.0)
    el2 = create_elements(mesh, "FACE2")
    update!(el2, "displacement 2", 0.0)
    el3 = create_elements(mesh, "FACE3")
    update!(el3, "displacement 1", 0.0)
    add_elements!(bc, el1, el2, el3)

    load = Problem(Elasticity, "pressure load", 3)
    load_elements = create_elements(mesh, "OUTER")
    update!(load_elements, "surface pressure", 7317.0)
    add_elements!(load, load_elements)

    analysis = Analysis(Linear)
    add_problems!(analysis, body, load, bc)
    run!(analysis)

    X = load("geometry", 0.0)
    u = load("displacement", 0.0)
    nids = sort(collect(keys(X)))
    umag = Float64[norm(u[id]) for id in nids]
    um = mean(umag)
    us = std(umag)
    rtol = norm(um - 0.9) / max(norm(um), 0.9) * 100.0
    @debug("Displacement field statistics", mean=um, std=us, rtol=rtol)
    @test rtol < 1.5 # percents

    @debug("Verifying displacement against Code Aster solution.. ")
    pass = true
    for nid in keys(u_CA)
        rtol = norm(u[nid] - u_CA[nid]) / max(norm(u[nid]), norm(u_CA[nid])) * 100.0
        @debug("Node id $nid, rel diff to CA = $rtol %")
        pass &= isapprox(u[nid], u_CA[nid])
    end
    @test pass

    S = lsq_fit(body.elements, get_stress, 0.0)

    @debug("Verifying stress against Code Aster solution.. ")
    for nid in keys(S_CA)
        rtol = norm(S[nid] - S_CA[nid]) / max(norm(S[nid]), norm(S_CA[nid])) * 100.0
        @debug("Node id $nid, S=$(S[nid]), S_CA=$(S_CA[nid]), rel diff to CA = $rtol %")
        # @test isapprox(S[nid], S_CA[nid])
        # http://code-aster.org/doc/default/en/man_r/r3/r3.06.03.pdf
        pass &= (rtol < 10.0) # percents
    end
    @test pass

    # Calculate principal stresses in nodes
    Sp = Dict()
    for (nid, s) in S
        # order is: 11, 22, 33, 12, 23, 13
        stress_tensor = [
            s[1] s[4] s[6]
            s[4] s[2] s[5]
            s[6] s[5] s[3]]
        Sp[nid] = sort(eigvals(stress_tensor))
    end

    node_ids = sort(collect(keys(Sp)))
    for (i, nid) in enumerate(node_ids)
        @debug("$nid -> $(Sp[nid])")
        if i > 9
            @debug("...")
            break
        end
    end
    # Stress state in M141 element
    S_CA_gp = Dict()
    S_CA_gp[1] = [-2.31497155149131E+04, -2.30873958541450E+04, -2.59731750416165E+04, 1.16262216932559E+04, 1.05883267606309E+04, 1.04255059462671E+04]
    S_CA_gp[2] = [-2.20268425477595E+04, -2.59602179285979E+04, -2.45141753387488E+04, 1.08344827948196E+04, 1.02393080182646E+04, 1.18221520216293E+04]
    S_CA_gp[3] = [-2.12175032863505E+04, -2.38106746858742E+04, -2.75848287641255E+04, 1.23383626529382E+04, 9.45231798421561E+03, 1.04600859263009E+04]
    S_CA_gp[4] = [-2.88674347639302E+04, -2.90047144010521E+04, -3.19545639967296E+04, 1.24436733999204E+04, 1.14033189825355E+04, 1.09837505652098E+04]
    S_CA_gp[5] = [-2.76741820176108E+04, -3.20576001061179E+04, -3.04041166591832E+04, 1.16023096459047E+04, 1.10324243354755E+04, 1.24679360552984E+04]
    S_CA_gp[6] = [-2.68141147396280E+04, -2.97733271271234E+04, -3.36672334497799E+04, 1.32004501509553E+04, 1.01961070966382E+04, 1.10204979599285E+04]

    time = 0.0
    for element in body.elements
        element.id == 141 || continue
        for (i, ip) in enumerate(get_integration_points(element))
            S = get_stress(element, ip, time)
            rtol = norm(S - S_CA_gp[i]) / max(norm(S), norm(S_CA_gp[i]))
            @debug("$i $S, rtol=$rtol")
            pass &= isapprox(S, S_CA_gp[i])
        end
    end
    @test pass
end

# 1/8 hollow sphere with surface load
u_CA = Dict()
u_CA[38] = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
u_CA[49] = [-4.50934684240566E-01, -4.46908405333021E-01, -5.92329377847994E-01]
S_CA = Dict()
S_CA[38] = [-1.94504819940510E+03, -3.47014655438479E+04, -3.40876135857114E+04, 1.70379805227866E+03, 2.16707698388864E+03, 4.32009983342141E-12]
S_CA[49] = [-2.36067010168718E+04, -2.33820775692753E+04, -1.80606705820104E+04, 8.85783922092890E+03, 1.12122431934777E+04, 1.14035796957343E+04]
test_wedge_sphere("HOLLOWSPHERE8_WEDGE6", u_CA, S_CA)
#test_wedge_sphere("HOLLOWSPHERE8_WEDGE15")


text:
using Test, CatmullClark, Makie

include("test1.jl")
include("test2.jl")


text:
#------------------------------------------------------------------------------------------------------------------------#
#    This is central notebook which executes coupled heat and moisture transfer based on defined input (\"input.jl\")    #
#------------------------------------------------------------------------------------------------------------------------#

using PlotlyJS
using DelimitedFiles
using LinearAlgebra
using Distributed
using DistributedArrays
using SparseArrays
using Krylov
# using ArrayFire
# using CUDAdrv, CUDAnative, CuArrays

include("material_library.jl")
include("functions.jl")

###---1D problem analysed as 2D steady-state heat transfer---###

#  include("inputs/input_1D_ss.jl")
#
#  M=mesh(E,dl) #Mesh matrix -- function on matrix of geometry, materials and boundary conditions (E) included from input.jl
#  #M=readdlm("projects/HM_T/mesh.csv",',');	#Reading already defined mesh
#  Ns=nsort(E,M,dl)
#  Con=con1(E,M,dl,Ns)
#  K=Kmat(M,Con,dl)
#  @time begin
#  x=modconjgrad(K,M)
#  end
#
#  # @time begin
#  # K2=Kmat2(M,Con,dl)
#  # k2=K2[1];f=K2[2]
#  # kinv=inv(k2)
#  # x=kinv*f
#  M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#M2=[["x_coord" "y_coord"];M[:,3:4]]#M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#
#  Z0_vect=zeros(size(M2)[1]-1)#Z0_vect=zeros(size(M)[1])#Z0_vect=zeros(size(M2)[1]-1)#
#  Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]#
#  x1=["temp";x]#x1=["temp";x;M[size(x)[1]+1:end,7]]#x1=["temp";x]#
#  writedlm("results/paraview_1D_ss.csv",  [M2 Z0_vect x1], ',')#writedlm("results/paraview_results.csv",  [M2 Z0_vect x1], ',')
#  # end
#  xres=[x;M[size(x)[1]+1:end,7]]; #results with BC
#
# ###---1D problem analysed as 2D transient heat transfer---###
#
#  println("Calculation started!")
#  @time begin
#  include("inputs/input_1D_transient.jl")
#  M=mesh(E,dl)
#  M2=copy(M)
#  Ns=nsort(E,M,dl)
#  Con=con1(E,M,dl,Ns)
#  K=Kmat(M,Con,dl,dt)
#  x_in=10*ones(size(K)[1])
#  Temp_e=readdlm("inputs/Cos_1_240.csv",',')
#  t_e_ind1=floor(Int64,Ns[4][1,1])	#global IDs of left - exterior BC
#  t_e_ind2=floor(Int64,Ns[4][end,1])	#
#  Temp_i=22*ones(size(K)[1])
#  t_i_ind1=floor(Int64,Ns[6][1,1])	#global IDs of right - interior BC
#  t_i_ind2=floor(Int64,Ns[6][end,1])	#
#  for i=1:floor(Int64,round(T/dt))
#  	for j=t_e_ind1:t_e_ind2
#  		M2[j,7]=Temp_e[i]
#  	end
#  	for j=t_i_ind1:t_i_ind2
#  		M2[j,7]=Temp_i[i]
#  	end
#  	global K2=Kmat(M2,Con,dl,dt)
#  	global x1=modconjgrad(K2,M2,dt,dl,0,x_in)
#  	global M3=[["x_coord" "y_coord"];M2[1:(size(x)[1]),3:4]]
#  	global Z0_vect=zeros(size(M3)[1]-1)
#  	global Z0_vect=["z_coord";Z0_vect]
#  	global xx=["temp";x1]
#  	writedlm("results/transient/1Dt/1D_transient$i.csv",  [M3 Z0_vect xx], ',')
#  	global x_in=copy(x1)
#  end
#  println("End of calculation!")
#  end
#
# ###Up to steady state input1###
#  println("Calculation started!")
#  @time begin
#  include("inputs/input.jl")
#  M=mesh(E,dl)
#  M2=copy(M)
#  Ns=nsort(E,M,dl)
#  Con=con1(E,M,dl,Ns)
#  K=Kmat(M,Con,dl,dt)
#  x_in=10*ones(size(K)[1])
#  #Temp_e=readdlm("inputs/Cos_1_240.csv",',')
#  #t_e_ind1=floor(Int64,Ns[4][1,1])	#global IDs of left - exterior BC
#  #t_e_ind2=floor(Int64,Ns[4][end,1])	#
#  #Temp_i=20*ones(size(K)[1])
#  #t_i_ind1=floor(Int64,Ns[6][1,1])	#global IDs of right - interior BC
#  #t_i_ind2=floor(Int64,Ns[6][end,1])	#
#  for i=1:floor(Int64,round(T/dt))
#  	# for j=t_e_ind1:t_e_ind2
#  	# 	M2[j,7]=Temp_e[i]
#  	# end
#  	# for j=t_i_ind1:t_i_ind2
#  	# 	M2[j,7]=Temp_i[i]
#  	# end
#  	global K2=Kmat(M2,Con,dl,dt)
#  	@time begin
#  	global x1=modconjgrad(K2,M2,dt,dl,0,x_in)
#  	time end
#  	global M3=[["x_coord" "y_coord"];M2[1:(size(x1)[1]),3:4]]
#  	global Z0_vect=zeros(size(M3)[1]-1)
#  	global Z0_vect=["z_coord";Z0_vect]
#  	global xx=["temp";x1]
#  	writedlm("results/transient/input_ss_to_tr/iss$i.csv",  [M3 Z0_vect xx], ',')
#  	global x_in=copy(x1)
#  end
#  println("End of calculation!")
#  end
#
# ##---2D steady-state heat transfer SparseArrays---###
#  include("inputs/input.jl")
#  @time begin
#  M=mesh(E,dl) #Mesh matrix -- function on matrix of geometry, materials and boundary conditions (E) included from input.jl
#  end
#  #M=readdlm("projects/HM_T/mesh.csv",',');	#Reading already defined mesh
#  @time begin
#  Ns=nsort(E,M,dl)
#  end
#  @time begin
#  Con=con1(E,M,dl,Ns)
#  end
#  @time begin
#  Spar=Kspar_ss(M,Con,dl)
#  end
#  K=Spar[1]
#  f=Spar[2]
#  @time begin
#  x=cg_lanczos(K,f)[1]
#  end
#
#  # @time begin
#  # K2=Kmat2(M,Con,dl)
#  # k2=K2[1];f=K2[2]
#  # kinv=inv(k2)
#  # x=kinv*f
#  M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#M2=[["x_coord" "y_coord"];M[:,3:4]]#M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#
#  Z0_vect=zeros(size(M2)[1]-1)#Z0_vect=zeros(size(M)[1])#Z0_vect=zeros(size(M2)[1]-1)#
#  Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]#
#  x1=["temp";x]#x1=["temp";x;M[size(x)[1]+1:end,7]]#x1=["temp";x]#
#  writedlm("results/paraview_input.csv",  [M2 Z0_vect x1], ',')#writedlm("results/paraview_results.csv",  [M2 Z0_vect x1], ',')
#  # end
#  xres=[x;M[size(x)[1]+1:end,7]]; #results with BC
#
# ###---2D steady-state heat transfer SparseArrays---###
#
#  include("inputs/input.jl")
#  @time begin
#  M=mesh(E,dl) #Mesh matrix -- function on matrix of geometry, materials and boundary conditions (E) included from input.jl
#  end
#  #M=readdlm("projects/HM_T/mesh.csv",',');	#Reading already defined mesh
#  @time begin
#  Ns=nsort(E,M,dl)
#  end
#  @time begin
#  Con=con1(E,M,dl,Ns)
#  end
#  @time begin
#  Spar=Kspar_ss(M,Con,dl)
#  end
#  K=Spar[1]
#  f=Spar[2]
#  @time begin
#  x=cg_lanczos(K,f)[1]
#  end
#
#  # @time begin
#  # K2=Kmat2(M,Con,dl)
#  # k2=K2[1];f=K2[2]
#  # kinv=inv(k2)
#  # x=kinv*f
#  M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#M2=[["x_coord" "y_coord"];M[:,3:4]]#M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#
#  Z0_vect=zeros(size(M2)[1]-1)#Z0_vect=zeros(size(M)[1])#Z0_vect=zeros(size(M2)[1]-1)#
#  Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]#
#  x1=["temp";x]#x1=["temp";x;M[size(x)[1]+1:end,7]]#x1=["temp";x]#
#  writedlm("results/paraview_input.csv",  [M2 Z0_vect x1], ',')#writedlm("results/paraview_results.csv",  [M2 Z0_vect x1], ',')
#  # end
#  xres=[x;M[size(x)[1]+1:end,7]]; #results with BC
#
# ###---2D transient heat transfer SparseArrays---###
#
#  println("Calculation started!")
#  @time begin
#  include("inputs/input.jl")
#  M=mesh(E,dl)
#  M2=copy(M)
#  Ns=nsort(E,M,dl)
#  Con=con1(E,M,dl,Ns)
#  Spar=Kspar_t(M,Con,dl,dt)
#  K=Spar[1]
#  f=Spar[2]
#  x_in=10*ones(size(K)[1])
#  #Temp_e=readdlm("inputs/Cos_1_240.csv",',')
#  #t_e_ind1=floor(Int64,Ns[4][1,1])	#global IDs of left - exterior BC
#  #t_e_ind2=floor(Int64,Ns[4][end,1])	#
#  #Temp_i=20*ones(size(K)[1])
#  #t_i_ind1=floor(Int64,Ns[6][1,1])	#global IDs of right - interior BC
#  #t_i_ind2=floor(Int64,Ns[6][end,1])	#
#  for i=1:floor(Int64,round(T/dt))
#  	# for j=t_e_ind1:t_e_ind2
#  	# 	M2[j,7]=Temp_e[i]
#  	# end
#  	# for j=t_i_ind1:t_i_ind2
#  	# 	M2[j,7]=Temp_i[i]
#  	# end
#  	#global Spar=Kspar_t(M2,Con,dl,dt)
#  	#global K2=Spar[1]
#  	println("Timestep $i")
#  	@time begin
#  	global ft=init(x_in,M2,dl,dt)
#  	global f2=f+ft
#  	global x1=cg_lanczos(K,f2)[1]
#  	global M3=[["x_coord" "y_coord"];M2[1:(size(x1)[1]),3:4]]
#  	global Z0_vect=zeros(size(M3)[1]-1)
#  	global Z0_vect=["z_coord";Z0_vect]
#  	global xx=["temp";x1]
#  	writedlm("results/transient/input_ss_to_tr/iss$i.csv",  [M3 Z0_vect xx], ',')
#  	global x_in=copy(x1)
#  	end
#  end
#  println("End of calculation!")
#  end
#
###---2D transient heat and mass transfer SparseArrays---###
 println("Calculation started!")
 @time begin
 include("inputs/input.jl")
 include("environmental.jl")		#Environmental hourly BC --> Tem{Float64}(2919,1), r_H{Float64}(2919,1) Ns[20]
 M=mesh(E,dl)
 Ns=nsort(E,M,dl)
 Con=con1(E,M,dl,Ns)
 Spar=Kspar_ss(M,Con,dl)
 K=Spar[1]
 f=Spar[2]
 x_t_in=cg_lanczos(K,f)[1]
 M2=[["x_coord" "y_coord"];M[1:(size(x_t_in)[1]),3:4]]
 Z0_vect=zeros(size(M2)[1]-1)
 Z0_vect=["z_coord";Z0_vect]
 xx=["temp";x_t_in]
 writedlm("results/transient/input_HM/t/t0.csv", [M2 Z0_vect xx], ',')
 x_m_in=0.1*ones(size(x_t_in)[1])
 xx=["temp";x_m_in]
 writedlm("results/transient/input_HM/m/m0.csv",  [M2 Z0_vect xx], ',')
 t_e_ind1=floor(Int64,Ns[20][1,1])
 t_e_ind2=floor(Int64,Ns[20][end,1])
 s1_1=floor(Int64,Ns[17][1,1])
 s1_2=floor(Int64,Ns[17][end,1])
 s2_1=floor(Int64,Ns[18][1,1])
 s2_2=floor(Int64,Ns[18][end,1])
 s3_1=floor(Int64,Ns[19][1,1])
 s3_2=floor(Int64,Ns[19][end,1])
 s4_1=floor(Int64,Ns[25][1,1])
 s4_2=floor(Int64,Ns[25][end,1])
 s5_1=floor(Int64,Ns[26][1,1])
 s5_2=floor(Int64,Ns[26][end,1])
 s6_1=floor(Int64,Ns[27][1,1])
 s6_2=floor(Int64,Ns[27][end,1])
 inrH_1=floor(Int64,Ns[22][1,1])
 inrH_2=floor(Int64,Ns[23][end,1])
 Mt=copy(M)
 Mm=copy(M)
 for i=inrH_1:inrH_2
 	Mm[i,7]=0.5
 end
 for i=1:744
 	@time begin
 	println("Timestep $i")
 	for j=t_e_ind1:t_e_ind2
 		Mt[j,7]=Tem[i]
 	end
 	for j=s1_1:s1_2
 		Mt[j,7]=soil3(Mt[j,4])
 	end
 	for j=s2_1:s2_2
 		Mt[j,7]=soil3(Mt[j,4])
 	end
 	for j=s3_1:s3_2
 		Mt[j,7]=soil3(Mt[j,4])
 	end
 	for j=s4_1:s4_2
 		Mt[j,7]=soil3(Mt[j,4])
 	end
 	for j=s5_1:s5_2
 		Mt[j,7]=soil3(Mt[j,4])
 	end
 	for j=s6_1:s6_2
 		Mt[j,7]=soil3(Mt[j,4])
 	end
 	global Spar=Kspar_t(Mt,Con,dl,dt)
 	global Ktem=Spar[1]
 	global fk=Spar[2]
 	global ft=init(x_t_in,Mt,dl,dt)
 	global ftem=fk+ft
 	global x_t=cg_lanczos(Ktem,ftem)[1]
 	global M2=[["x_coord" "y_coord"];Mt[1:(size(x_t)[1]),3:4]]
 	global Z0_vect=zeros(size(M2)[1]-1)
 	global Z0_vect=["z_coord";Z0_vect]
 	global xx=["temp";x_t]
 	#writedlm("results/transient/input_HM/t/t$i.csv",  [M2 Z0_vect xx], ',')
 	global x_t_in=copy(x_t)
 	global Mm=moist(Mm,x_t_in,x_m_in)
 	for j=t_e_ind1:t_e_ind2
 		Mm[j,7]=r_H[i]
 	end
 	for j=s1_1:s1_2
 		Mm[j,7]=0.95
 	end
 	for j=s2_1:s2_2
 		Mm[j,7]=0.95
 	end
 	for j=s3_1:s3_2
 		Mm[j,7]=0.80
 	end
 	for j=s4_1:s4_2
 		Mm[j,7]=0.95
 	end
 	for j=s5_1:s5_2
 		Mm[j,7]=0.95
 	end
 	for j=s6_1:s6_2
 		Mm[j,7]=0.80
 	end
 	global Spar2=Kspar_t(Mm,Con,dl,dt)
 	global Kmoi=Spar2[1]
 	global fkm=Spar2[2]
 	global ftm=init(x_m_in,Mm,dl,dt)
 	global fmoi=fkm+ftm
 	global x_m=cg_lanczos(Kmoi,fmoi)[1]
 	if i==1
 		global M2=[["x_coord" "y_coord"];Mm[1:(size(x_m)[1]),3:4]]
 		global xx=["temp";x_m]
 		writedlm("results/transient/input_HM/m/m$i.csv",  [M2 Z0_vect xx], ',')
 	end
 	if mod(i,24)==0
 		global M2=[["x_coord" "y_coord"];Mm[1:(size(x_m)[1]),3:4]]
 		global xx=["temp";x_m]
 		writedlm("results/transient/input_HM/m/m$i.csv",  [M2 Z0_vect xx], ',')
 	end
 	global x_m_in=copy(x_m)
 	end #time
 end
 for i=745:1464
 		@time begin
 		println("Timestep $i")
 		for j=t_e_ind1:t_e_ind2
 			Mt[j,7]=Tem[i]
 		end
 		for j=s1_1:s1_2
 			Mt[j,7]=soil4(Mt[j,4])
 		end
 		for j=s2_1:s2_2
 			Mt[j,7]=soil4(Mt[j,4])
 		end
 		for j=s3_1:s3_2
 			Mt[j,7]=soil4(Mt[j,4])
 		end
 		for j=s4_1:s4_2
 			Mt[j,7]=soil4(Mt[j,4])
 		end
 		for j=s5_1:s5_2
 			Mt[j,7]=soil4(Mt[j,4])
 		end
 		for j=s6_1:s6_2
 			Mt[j,7]=soil4(Mt[j,4])
 		end
 		global Spar=Kspar_t(Mt,Con,dl,dt)
 		global Ktem=Spar[1]
 		global fk=Spar[2]
 		global ft=init(x_t_in,Mt,dl,dt)
 		global ftem=fk+ft
 		global x_t=cg_lanczos(Ktem,ftem)[1]
 		global M2=[["x_coord" "y_coord"];Mt[1:(size(x_t)[1]),3:4]]
 		global Z0_vect=zeros(size(M2)[1]-1)
 		global Z0_vect=["z_coord";Z0_vect]
 		global xx=["temp";x_t]
 		#writedlm("results/transient/input_HM/t/t$i.csv",  [M2 Z0_vect xx], ',')
 		global x_t_in=copy(x_t)
 		global Mm=moist(Mm,x_t_in,x_m_in)
 		global Spar2=Kspar_t(Mm,Con,dl,dt)
 		global Kmoi=Spar2[1]
 		global fkm=Spar2[2]
 		global ftm=init(x_m_in,Mm,dl,dt)
 		global fmoi=fkm+ftm
 		global x_m=cg_lanczos(Kmoi,fmoi)[1]
 		if mod(i,24)==0
 			global M2=[["x_coord" "y_coord"];Mm[1:(size(x_m)[1]),3:4]]
 			global xx=["temp";x_m]
 			writedlm("results/transient/input_HM/m/m$i.csv",  [M2 Z0_vect xx], ',')
 		end
 		global x_m_in=copy(x_m)
 		end #time
 end
 for i=1465:2208
 		@time begin
 		println("Timestep $i")
 		for j=t_e_ind1:t_e_ind2
 			Mt[j,7]=Tem[i]
 		end
 		for j=s1_1:s1_2
 			Mt[j,7]=soil5(Mt[j,4])
 		end
 		for j=s2_1:s2_2
 			Mt[j,7]=soil5(Mt[j,4])
 		end
 		for j=s3_1:s3_2
 			Mt[j,7]=soil5(Mt[j,4])
 		end
 		for j=s4_1:s4_2
 			Mt[j,7]=soil5(Mt[j,4])
 		end
 		for j=s5_1:s5_2
 			Mt[j,7]=soil5(Mt[j,4])
 		end
 		for j=s6_1:s6_2
 			Mt[j,7]=soil5(Mt[j,4])
 		end
 		global Spar=Kspar_t(Mt,Con,dl,dt)
 		global Ktem=Spar[1]
 		global fk=Spar[2]
 		global ft=init(x_t_in,Mt,dl,dt)
 		global ftem=fk+ft
 		global x_t=cg_lanczos(Ktem,ftem)[1]
 		global M2=[["x_coord" "y_coord"];Mt[1:(size(x_t)[1]),3:4]]
 		global Z0_vect=zeros(size(M2)[1]-1)
 		global Z0_vect=["z_coord";Z0_vect]
 		global xx=["temp";x_t]
 		#writedlm("results/transient/input_HM/t/t$i.csv",  [M2 Z0_vect xx], ',')
 		global x_t_in=copy(x_t)
 		global Mm=moist(Mm,x_t_in,x_m_in)
 		global Spar2=Kspar_t(Mm,Con,dl,dt)
 		global Kmoi=Spar2[1]
 		global fkm=Spar2[2]
 		global ftm=init(x_m_in,Mm,dl,dt)
 		global fmoi=fkm+ftm
 		global x_m=cg_lanczos(Kmoi,fmoi)[1]
 		if mod(i,24)==0
 			global M2=[["x_coord" "y_coord"];Mm[1:(size(x_m)[1]),3:4]]
 			global xx=["temp";x_m]
 			writedlm("results/transient/input_HM/m/m$i.csv",  [M2 Z0_vect xx], ',')
 		end
 		global x_m_in=copy(x_m)
 		end #time
 end
 for i=2209:2919
 		@time begin
 		println("Timestep $i")
 		for j=t_e_ind1:t_e_ind2
 			Mt[j,7]=Tem[i]
 		end
 		for j=s1_1:s1_2
 			Mt[j,7]=soil6(Mt[j,4])
 		end
 		for j=s2_1:s2_2
 			Mt[j,7]=soil6(Mt[j,4])
 		end
 		for j=s3_1:s3_2
 			Mt[j,7]=soil6(Mt[j,4])
 		end
 		for j=s4_1:s4_2
 			Mt[j,7]=soil6(Mt[j,4])
 		end
 		for j=s5_1:s5_2
 			Mt[j,7]=soil6(Mt[j,4])
 		end
 		for j=s6_1:s6_2
 			Mt[j,7]=soil6(Mt[j,4])
 		end
 		global Spar=Kspar_t(Mt,Con,dl,dt)
 		global Ktem=Spar[1]
 		global fk=Spar[2]
 		global ft=init(x_t_in,Mt,dl,dt)
 		global ftem=fk+ft
 		global x_t=cg_lanczos(Ktem,ftem)[1]
 		global M2=[["x_coord" "y_coord"];Mt[1:(size(x_t)[1]),3:4]]
 		global Z0_vect=zeros(size(M2)[1]-1)
 		global Z0_vect=["z_coord";Z0_vect]
 		global xx=["temp";x_t]
 		#writedlm("results/transient/input_HM/t/t$i.csv",  [M2 Z0_vect xx], ',')
 		global x_t_in=copy(x_t)
 		global Mm=moist(Mm,x_t_in,x_m_in)
 		global Spar2=Kspar_t(Mm,Con,dl,dt)
 		global Kmoi=Spar2[1]
 		global fkm=Spar2[2]
 		global ftm=init(x_m_in,Mm,dl,dt)
 		global fmoi=fkm+ftm
 		global x_m=cg_lanczos(Kmoi,fmoi)[1]
 		if mod(i,24)==0
 			global M2=[["x_coord" "y_coord"];Mm[1:(size(x_m)[1]),3:4]]
 			global xx=["temp";x_m]
 			writedlm("results/transient/input_HM/m/m$i.csv",  [M2 Z0_vect xx], ',')
 		end
 		if i==2919
 			global M2=[["x_coord" "y_coord"];Mm[1:(size(x_m)[1]),3:4]]
 			global xx=["temp";x_m]
 			writedlm("results/transient/input_HM/m/m$i.csv",  [M2 Z0_vect xx], ',')
 		end
 		global x_m_in=copy(x_m)
 		end #time
end
println("End of calculation!")
end #time


text:
## -------------------------------------------------------------------
# model vs exp comparison
let
    METHODS = [:FBA_Z_FIX_MIN_VG_MAX_COST, :ME_MAX_POL]

    sort_iders(iders) = sort!(deepcopy(iders),
        by = (ider) -> ider == "GLC" ? "AAA" : ider
    )

    # Plot
    gc = 1
    xticks_pos = []
    xticks_iders = []
    p = plot(;xlabel = "exchange", ylabel = "flx/exp_glc_flx")

    for srcs in [
        [
            (ChK.iJR904, ChK.KayserData)],
            [(ChN.iJR904, ChN.NanchenData)],
            [(ChF.iJR904, ChF.FolsomData)
        ],
        [
            (ChH.iJR904, ChH.HeerdenData)
        ]
    ]

        # Collect
        dat_pool = Dict()
        for (iJR, Data) in srcs

            src = nameof(Data)

            DAT = ChE.load_DAT(src)

            EXCH_FLX_IDERS = sort_iders(DAT[:EXCH_FLX_IDERS])
            EXPS = DAT[:EXPS]

            @info("Doing", src)

            # Experiments
            edat = get!(dat_pool, :exp, Dict())
            for exp in EXPS
                !haskey(DAT, :exp, :flx, "GLC", exp) && continue
                glc_exp_val = DAT[:exp, :flx, "GLC", exp]

                for ider in EXCH_FLX_IDERS
                    exp_val = abs(DAT[:exp, :flx, ider, exp] / glc_exp_val)

                    @info(:exp, exp, ider, exp_val); println()

                    eidat = get!(edat, ider, [])
                    push!(eidat, exp_val)
                end
            end

            # methods
            for method in METHODS
                mdat = get!(dat_pool, method, Dict())

                for exp in EXPS
                    !haskey(DAT, :exp, :flx, "GLC", exp) && continue
                    glc_exp_val = DAT[:exp, :flx, "GLC", exp]

                    for ider in EXCH_FLX_IDERS
                        !haskey(DAT, method, :flx, ider, exp) && continue
                        model_val = abs(DAT[method, :flx, ider, exp] / glc_exp_val)

                        # @info(method, exp, ider, model_val); println()

                        midat = get!(mdat, ider, [])
                        push!(midat, model_val)

                    end
                end
            end
        end # for dsource

        # Plos
        vline!(p, [gc]; label = "", color = :black, alpha = 0.7, ls = :dash)
        for (ii, ider) in sort_iders(ALL_IDERS) |> enumerate

            # Experimental
            edat = dat_pool[:exp]
            !haskey(edat, ider) && continue
            m = mean(abs.(edat[ider]))
            e = std(abs.(edat[ider]))

            @info(ider, m, e); println()

            scatter!(p, [gc], [m]; yerr = [e],
                color = :black, label = "", ms = 6, m = :star,
                alpha = 0.8
            )
            push!(xticks_pos, gc)
            push!(xticks_iders, ider)
            gc += 5

            # modeled
            for method in METHODS
                mdat = dat_pool[method]
                color = method_colors[method]
                μ = mean(abs.(mdat[ider]))
                e = std(abs.(mdat[ider]))
                m = (6, method_markers[method])
                scatter!(p, [gc], [μ]; yerr = [e],
                    color, label = "", m, alpha = 0.8
                )
                gc += 5
            end
            gc += 15
        end
    end # for srcs

    xticks = (xticks_pos, xticks_iders)
    xrotation = 35
    size = (1100, 500)
    thickness_scaling = 1.3
    plot!(p; xticks, xrotation, size, thickness_scaling)
    mysavefig(p, "exchs_flx_comparizon")
end


text:
using LinearAlgebra
using Symbolics
using SparseArrays

f(x) = dot(x[2:3], x[2:3])
c(x) = cos.(x)
L(x, y) = f(x) + dot(y, c(x))

@variables x[1:3], y[1:3]

Lsym = L(x, y)
hess = Symbolics.hessian(Lsym, x)
hess_func = eval(Symbolics.build_function(hess, x, y)[1])
hess_sp = Symbolics.sparsehessian(Lsym, x)
hess_sp_func = eval(Symbolics.build_function(hess_sp.nzval, x, y)[2])
hess_sp.nzval
sparsity = [findnz(hess_sp)[1:2]...]
# findnz(hess_sp)
# hess_sp.nzval[1]
# hess_sp.nzval[2]
# hess_sp[1]

x0 = rand(3)
y0 = rand(3)
h0 = zeros(size(hess))
h0_sp = zeros(length(hess_sp.nzval))

hess_sp_func(h0_sp, x0, y0)
for (i, h) in enumerate(h0_sp)
    h0[sparsity[1][i], sparsity[2][i]] = h
end

h0 - hess_func(x0, y0)