text:
abstract type ShaderStage end
struct Vertex <: ShaderStage end
struct Geometry <: ShaderStage end
struct Fragment <: ShaderStage end

struct Program
    context::AbstractContext
    vertexarray::AbstractArray
    uniforms::Dict{Symbol, Any}
    vertex_source::String
    fragment_source::String
end

struct InstancedProgram
    program::Program
    per_instance::AbstractArray
end

function getter_function(io, T, t_str, name)
    println(io, t_str, " get_$(name)(){return $name;}")
end

function getter_function(io, ::AbstractSampler, t_str, name)
end

function input_element(context::AbstractContext, stage::Vertex, io::IO, element::AbstractVector{T}, name::Symbol, uniforms) where {T}
    t_str = type_string(context, T)
    println(io, "in ", t_str, " $name;")
    getkey = Symbol(string(name, "_", "getter"))
    if haskey(uniforms, getkey)
        println(io, uniforms[getkey])
    else
        getter_function(io, T, t_str, name)
    end
end

function input_block(context::AbstractContext, io, input_array, uniforms)
    for name in propertynames(input_array)
        element = getproperty(input_array, name)
        input_element(context, Vertex(), io, element, name, uniforms)
    end
end

function InstancedProgram(
        context::AbstractContext,
        vertshader, fragshader,
        instance::AbstractVector,
        per_instance::AbstractVector;
        uniforms...
    )
    instance_attributes = sprint() do io
        println(io, "\n// Per instance attributes: ")
        input_block(context, io, per_instance, uniforms)
        println(io)
    end

    p = Program(
        context,
        instance_attributes * vertshader,
        fragshader,
        instance; uniforms...
    )
    return InstancedProgram(p, per_instance)
end

function vertex_header(context::AbstractContext)
    return """
    #version 300 es
    precision mediump int;
    precision mediump float;
    precision mediump sampler2D;
    precision mediump sampler3D;
    """
end

function fragment_header(context::AbstractContext)
    return """
    #version 300 es
    precision mediump int;
    precision mediump float;
    precision mediump sampler2D;
    precision mediump sampler3D;

    out vec4 fragment_color;
    """
end

function Program(
        context::AbstractContext,
        vertshader, fragshader,
        instance::AbstractVector;
        uniforms...
    )
    c_uniforms = Dict{Symbol, Any}()
    uniform_block = sprint() do io
        println(io, "\n// Uniforms: ")
        for (name, v) in uniforms
            endswith(string(name), "_getter") && continue
            vc = convert_uniform(context, v)
            t_str = try
               type_string(context, vc)
            catch e
                @error("Type $(typeof(vc)) isn't supported for uniform: $(name)")
                rethrow(e)
            end
            println(io, "uniform ", t_str, " $name;")
            getkey = Symbol(string(name, "_", "getter"))
            if haskey(uniforms, getkey)
                println(io, uniforms[getkey])
            else
                getter_function(io, vc, t_str, name)
            end
            c_uniforms[name] = vc
        end
        println(io)
    end
    src = sprint() do io
        println(io, "// Instance inputs: ")
        input_block(context, io, instance, uniforms)
        println(io, uniform_block)
        println(io)
        println(io, vertshader)
    end
    Program(
        context, instance, c_uniforms,
        vertex_header(context) * src,
        fragment_header(context) * uniform_block * fragshader
    )
end


text:
@testset "Testing IndexZero" begin
    IndexZero = CounterTools.IndexZero
    x = CounterTools.IndexZero(0)
    @test CounterTools.value(x) == 0

    # `Integers` should have 1 subtracted from them.
    # `IndexZero`s should just pass through
    y = CounterTools.indexzero(1)
    z = CounterTools.indexzero(x)
    @test CounterTools.value(y) == 0
    @test y == z

    A = [1,2,3]
    B = (1,2,3)
    @test A[y] == A[1]
    @test B[y] == B[1]

    # Iterator interface
    @test collect(y) == [y]
end

ext: jl


text:
mutable struct Bus
  bus_i::Int
  bustype::Int
  Pd::Float64
  Qd::Float64
  Gs::Float64
  Bs::Float64
  area::Int
  Vm::Float64
  Va::Float64
  baseKV::Float64
  zone::Int
  Vmax::Float64
  Vmin::Float64
end

mutable struct Line
  from::Int
  to::Int
  r::Float64
  x::Float64
  b::Float64
  rateA::Float64
  rateB::Float64
  rateC::Float64
  ratio::Float64 #TAP
  angle::Float64 #SHIFT
  status::Int
  angmin::Float64
  angmax::Float64
end
Line() = Line(0,0,0.,0.,0.,0.,0.,0.,0.,0.,0,0.,0.)

mutable struct Gener
  # .gen fields
  bus::Int
  Pg::Float64
  Qg::Float64
  Qmax::Float64
  Qmin::Float64
  Vg::Float64
  mBase::Float64
  status::Int
  Pmax::Float64
  Pmin::Float64
  Pc1::Float64
  Pc2::Float64
  Qc1min::Float64
  Qc1max::Float64
  Qc2min::Float64
  Qc2max::Float64
  ramp_agc::Float64
  # .gencost fields
  gentype::Int
  startup::Float64
  shutdown::Float64
  n::Int
  coeff::Array
end

mutable struct OPFData
  buses::Array{Bus}
  lines::Array{Line}
  generators::Array{Gener}
  bus_ref::Int
  baseMVA::Float64
  BusIdx::Dict{Int,Int}    #map from bus ID to bus index
  FromLines::Array         #From lines for each bus (Array of Array)
  ToLines::Array           #To lines for each bus (Array of Array)
  BusGenerators::Array     #list of generators for each bus (Array of Array)
end

function opf_loaddata(case_name, lineOff=Line())
  #
  # load buses
  #
  # bus_arr = readdlm("data/" * case_name * ".bus")
  bus_arr = readdlm(case_name * ".bus")
  num_buses = size(bus_arr,1)
  buses = Array{Bus}(undef, num_buses)
  bus_ref=-1
  for i in 1:num_buses
    @assert bus_arr[i,1]>0  #don't support nonpositive bus ids
    buses[i] = Bus(bus_arr[i,1:13]...)
    buses[i].Va *= pi/180
    if buses[i].bustype==3
      if bus_ref>0
        error("More than one reference bus present in the data")
      else
         bus_ref=i
      end
    end
    #println("bus ", i, " ", buses[i].Vmin, "      ", buses[i].Vmax)
  end

  #
  # load branches/lines
  #
  # branch_arr = readdlm("data/" * case_name * ".branch")
  branch_arr = readdlm(case_name * ".branch")
  num_lines = size(branch_arr,1)
  lines_on = findall((branch_arr[:,11].>0) .& ((branch_arr[:,1].!=lineOff.from) .| (branch_arr[:,2].!=lineOff.to)) )
  num_on   = length(lines_on)

  if lineOff.from>0 && lineOff.to>0
    println("opf_loaddata: was asked to remove line from,to=", lineOff.from, ",", lineOff.to)
    #println(lines_on, branch_arr[:,1].!=lineOff.from, branch_arr[:,2].!=lineOff.to)
  end
  if length(findall(branch_arr[:,11].==0))>0
    println("opf_loaddata: ", num_lines-length(findall(branch_arr[:,11].>0)), " lines are off and will be discarded (out of ", num_lines, ")")
  end

  lines = Array{Line}(undef, num_on)

  lit=0
  for i in lines_on
    @assert branch_arr[i,11] == 1  #should be on since we discarded all other
    lit += 1
    lines[lit] = Line(branch_arr[i, 1:13]...)
    if lines[lit].angmin>-360 || lines[lit].angmax<360
      error("Bounds of voltage angles are still to be implemented.")
    end

  end
  @assert lit == num_on

  #
  # load generators
  #
  # gen_arr = readdlm("data/" * case_name * ".gen")
  gen_arr = readdlm(case_name * ".gen")
  # costgen_arr = readdlm("data/" * case_name * ".gencost")
  costgen_arr = readdlm(case_name * ".gencost")
  num_gens = size(gen_arr,1)

  baseMVA=100

  @assert num_gens == size(costgen_arr,1)

  gens_on=findall(x->x!=0, gen_arr[:,8]); num_on=length(gens_on)
  if num_gens-num_on>0
    println("loaddata: ", num_gens-num_on, " generators are off and will be discarded (out of ", num_gens, ")")
  end

  generators = Array{Gener}(undef, num_on)
  i=0
  for git in gens_on
    i += 1
    generators[i] = Gener(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, Array{Int}(undef, 0)) #gen_arr[i,1:end]...)

    generators[i].bus      = gen_arr[git,1]
    generators[i].Pg       = gen_arr[git,2] / baseMVA
    generators[i].Qg       = gen_arr[git,3] / baseMVA
    generators[i].Qmax     = gen_arr[git,4] / baseMVA
    generators[i].Qmin     = gen_arr[git,5] / baseMVA
    generators[i].Vg       = gen_arr[git,6]
    generators[i].mBase    = gen_arr[git,7]
    generators[i].status   = gen_arr[git,8]
    @assert generators[i].status==1
    generators[i].Pmax     = gen_arr[git,9]  / baseMVA
    generators[i].Pmin     = gen_arr[git,10] / baseMVA
    generators[i].Pc1      = gen_arr[git,11]
    generators[i].Pc2      = gen_arr[git,12]
    generators[i].Qc1min   = gen_arr[git,13]
    generators[i].Qc1max   = gen_arr[git,14]
    generators[i].Qc2min   = gen_arr[git,15]
    generators[i].Qc2max   = gen_arr[git,16]
    generators[i].gentype  = costgen_arr[git,1]
    generators[i].startup  = costgen_arr[git,2]
    generators[i].shutdown = costgen_arr[git,3]
    generators[i].n        = costgen_arr[git,4]
    if generators[i].gentype == 1
      generators[i].coeff = costgen_arr[git,5:end]
      error("Piecewise linear costs remains to be implemented.")
    else
      if generators[i].gentype == 2
        generators[i].coeff = costgen_arr[git,5:end]
        #println(generators[i].coeff, " ", length(generators[i].coeff), " ", generators[i].coeff[2])
      else
        error("Invalid generator cost model in the data.")
      end
    end
  end

  # build a dictionary between buses ids and their indexes
  busIdx = mapBusIdToIdx(buses)

  # set up the FromLines and ToLines for each bus
  FromLines,ToLines = mapLinesToBuses(buses, lines, busIdx)

  # generators at each bus
  BusGeners = mapGenersToBuses(buses, generators, busIdx)

  #println(generators)
  #println(bus_ref)
  return OPFData(buses, lines, generators, bus_ref, baseMVA, busIdx, FromLines, ToLines, BusGeners)
end

function  computeAdmitances(lines, buses, baseMVA)
  nlines = length(lines)
  YffR=Array{Float64}(undef, nlines)
  YffI=Array{Float64}(undef, nlines)
  YttR=Array{Float64}(undef, nlines)
  YttI=Array{Float64}(undef, nlines)
  YftR=Array{Float64}(undef, nlines)
  YftI=Array{Float64}(undef, nlines)
  YtfR=Array{Float64}(undef, nlines)
  YtfI=Array{Float64}(undef, nlines)

  for i in 1:nlines
    @assert lines[i].status == 1
    Ys = 1/(lines[i].r + lines[i].x*im)
    #assign nonzero tap ratio
    tap = lines[i].ratio==0 ? 1.0 : lines[i].ratio

    #add phase shifters
    tap *= exp(lines[i].angle * pi/180 * im)

    Ytt = Ys + lines[i].b/2*im
    Yff = Ytt / (tap*conj(tap))
    Yft = -Ys / conj(tap)
    Ytf = -Ys / tap

    #split into real and imag parts
    YffR[i] = real(Yff); YffI[i] = imag(Yff)
    YttR[i] = real(Ytt); YttI[i] = imag(Ytt)
    YtfR[i] = real(Ytf); YtfI[i] = imag(Ytf)
    YftR[i] = real(Yft); YftI[i] = imag(Yft)
    #@printf("[%4d]  tap=%12.9f   %12.9f\n", i, real(tap), imag(tap));
  end

  nbuses = length(buses)
  YshR = Array{Float64}(undef, nbuses)
  YshI = Array{Float64}(undef, nbuses)
  for i in 1:nbuses
    YshR[i] = buses[i].Gs / baseMVA
    YshI[i] = buses[i].Bs / baseMVA
    #@printf("[%4d]   Ysh  %15.12f + %15.12f i \n", i, YshR[i], YshI[i])
  end

  @assert 0==length(findall(isnan.(YffR)))+length(findall(isinf.(YffR)))
  @assert 0==length(findall(isnan.(YffI)))+length(findall(isinf.(YffI)))
  @assert 0==length(findall(isnan.(YttR)))+length(findall(isinf.(YttR)))
  @assert 0==length(findall(isnan.(YttI)))+length(findall(isinf.(YttI)))
  @assert 0==length(findall(isnan.(YftR)))+length(findall(isinf.(YftR)))
  @assert 0==length(findall(isnan.(YftI)))+length(findall(isinf.(YftI)))
  @assert 0==length(findall(isnan.(YtfR)))+length(findall(isinf.(YtfR)))
  @assert 0==length(findall(isnan.(YtfI)))+length(findall(isinf.(YtfI)))
  @assert 0==length(findall(isnan.(YshR)))+length(findall(isinf.(YshR)))
  @assert 0==length(findall(isnan.(YshI)))+length(findall(isinf.(YshI)))

  return YffR, YffI, YttR, YttI, YftR, YftI, YtfR, YtfI, YshR, YshI
end

# Builds a map from lines to buses.
# For each line we store an array with zero or one element containing
# the  'From' and 'To'  bus number.
function mapLinesToBuses(buses, lines, busDict)
  nbus = length(buses)
  FromLines = [Int[] for b in 1:nbus]
  ToLines   = [Int[] for b in 1:nbus]
  for i in 1:length(lines)
    busID = busDict[lines[i].from]
    @assert 1<= busID <= nbus
    push!(FromLines[busID], i)

    busID = busDict[lines[i].to]
    @assert 1<= busID  <= nbus
    push!(ToLines[busID], i)
  end
  return FromLines,ToLines
end

# Builds a mapping between bus ids and bus indexes
#
# Returns a dictionary with bus ids as keys and bus indexes as values
function mapBusIdToIdx(buses)
  dict = Dict{Int,Int}()
  for b in 1:length(buses)
    @assert !haskey(dict,buses[b].bus_i)
    dict[buses[b].bus_i] = b
  end
  return dict
end

# Builds a map between buses and generators.
# For each bus we keep an array of corresponding generators number (as array).
#
# (Can be more than one generator per bus)
function mapGenersToBuses(buses, generators,busDict)
  gen2bus = [Int[] for b in 1:length(buses)]
  for g in 1:length(generators)
    busID = busDict[ generators[g].bus ]
    #@assert(0==length(gen2bus[busID])) #at most one generator per bus
    push!(gen2bus[busID], g)
  end
  return gen2bus
end


text:
struct pyconvert_rule_ctypessimplevalue{R,S} <: Function end

function (::pyconvert_rule_ctypessimplevalue{R,SAFE})(::Type{T}, x::Py) where {R,SAFE,T}
    ptr = GC.@preserve x C.PySimpleObject_GetValue(Ptr{R}, getptr(x))
    ans = unsafe_load(ptr)
    if SAFE
        pyconvert_return(convert(T, ans))
    else
        pyconvert_tryconvert(T, ans)
    end
end

const CTYPES_SIMPLE_TYPES = [
    ("char", Cchar),
    ("wchar", Cwchar_t),
    ("byte", Cchar),
    ("ubyte", Cuchar),
    ("short", Cshort),
    ("ushort", Cushort),
    ("int", Cint),
    ("uint", Cuint),
    ("long", Clong),
    ("ulong", Culong),
    ("longlong", Clonglong),
    ("ulonglong", Culonglong),
    ("size_t", Csize_t),
    ("ssize_t", Cssize_t),
    ("float", Cfloat),
    ("double", Cdouble),
    ("char_p", Ptr{Cchar}),
    ("wchar_p", Ptr{Cwchar_t}),
    ("void_p", Ptr{Cvoid}),
]

function init_ctypes()
    for (t,T) in CTYPES_SIMPLE_TYPES
        isptr = endswith(t, "_p")
        isreal = !isptr
        isnumber = isreal
        isfloat = t in ("float", "double")
        isint = isreal && !isfloat
        isuint = isint && (startswith(t, "u") || t == "size_t")

        name = "ctypes/c_$t"
        rule = pyconvert_rule_ctypessimplevalue{T, false}()
        saferule = pyconvert_rule_ctypessimplevalue{T, true}()

        t == "char_p" && pyconvert_add_rule(name, Cstring, saferule)
        t == "wchar_p" && pyconvert_add_rule(name, Cwstring, saferule)
        pyconvert_add_rule(name, T, saferule)
        isuint && pyconvert_add_rule(name, UInt, sizeof(T) ≤ sizeof(UInt) ? saferule : rule)
        isuint && pyconvert_add_rule(name, Int, sizeof(T) < sizeof(Int) ? saferule : rule)
        isint && !isuint && pyconvert_add_rule(name, Int, sizeof(T) ≤ sizeof(Int) ? saferule : rule)
        isint && pyconvert_add_rule(name, Integer, rule)
        isfloat && pyconvert_add_rule(name, Float64, saferule)
        isreal && pyconvert_add_rule(name, Real, rule)
        isnumber && pyconvert_add_rule(name, Number, rule)
        isptr && pyconvert_add_rule(name, Ptr, saferule)
    end
end


text:
function ungetfield(dy, s::Tuple, f::Int)
    T = typeof(s)
    return Tangent{T}([i == f ? dy : ZeroTangent() for i=1:length(s)]...)
end

function ungetindex!(dx::AbstractArray, x::AbstractArray, dy::AbstractArray, I...)
    idx = CartesianIndices(size(x))[I...]
    idx = idx isa CartesianIndex ? [idx] : idx
    return NNlib.scatter!(+, dx, dy, idx)
end

_array_type_only(A::AT) where AT <: AbstractArray{T, N} where {T, N} = AT

function ungetindex!(dx::AbstractArray, x::AbstractArray, dy::Number, I...)
    d_dy = similar(dx, (1,))
    fill!(d_dy, dy)
    ungetindex!(dx, x, d_dy, I...)
end

function ungetindex(x::AbstractArray, dy, I...)
    dx = zero(x)
    return ungetindex!(dx, x, dy, I...)
end

function ungetindex(x::Tuple, dy, I...)
    dx = map(1:length(x)) do i
        i in I ? dy : zero(x[i])
    end
    return dx
end

"""
    _getfield(value, fld)
This function can be used instead of getfield() to bypass Yota rules
during backpropagation.
"""
_getfield(value, fld) = getfield(value, fld)

# function ∇sum(x::AbstractArray, dy)
#     dx = similar(x)
#     dx .= dy
#     return dx
# end

# function ∇mean(x::AbstractArray, dy, dims=1:ndims(x))
#     dx = similar(x)
#     dx .= dy ./ prod(size(x, d) for d in dims)
#     return dx
# end

# function sum_dropdims(x::AbstractArray, dims)
#     return dropdims(sum(x; dims=dims); dims=dims)
# end

# unbroadcast from Flux
# in in-place version we can consider sum!(similar(x), ds),
# but we need to carefully measure performance in each case

# reshape Δ to be consistent with x
trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

function unbroadcast(x::AbstractArray, Δ)
    if size(x) == size(Δ)
        return Δ
    elseif length(x) == length(Δ)
        return trim(x, Δ)
    else
        sum_dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))
        return trim(x, sum(Δ, dims=sum_dims))
    end
end

unbroadcast(::Number, Δ) = sum(Δ)

function unbroadcast_prod_x(x::AbstractArray, y::AbstractArray, Δ)
    if size(x) == size(Δ)
        return Δ .* y
    elseif length(x) == length(Δ)
        return trim(x, Δ .* y)
    else
        sum_dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))
        return trim(x, sum(Δ.* y, dims=sum_dims))
    end
end
unbroadcast_prod_y(x::AbstractArray, y::AbstractArray, Δ) = unbroadcast_prod_x(y, x, Δ)

# device_like(example, a) = (device = guess_device([example]); device(a))

# unbroadcast_prod_x(x::Number, y::AbstractArray, Δ) = unbroadcast_prod_x(device_like(y, [x]), y, Δ)[1]
unbroadcast_prod_x(x::Number, y::AbstractArray, Δ) = unbroadcast_prod_x(to_same_device([x], y), y, Δ)[1]
unbroadcast_prod_x(x::AbstractArray, y::Number, Δ) = unbroadcast_prod_x(x, to_same_device([y], x), Δ)
unbroadcast_prod_y(x::AbstractArray, y::Number, Δ) = unbroadcast_prod_y(x, to_same_device([y], x), Δ)[1]
unbroadcast_prod_y(x::Number, y::AbstractArray, Δ) = unbroadcast_prod_y(to_same_device([x], y), y, Δ)

untranspose_vec(ds::Transpose{T, <:AbstractVector{T}}) where T = transpose(ds)
untranspose_vec(ds::Adjoint{T, <:AbstractVector{T}}) where T = adjoint(ds)
untranspose_vec(ds::AbstractMatrix) = dropdims(transpose(ds); dims=2)

# function unvcat(dy::AbstractArray, n::Int, arrs::AbstractArray...)
#     a = arrs[n]
#     from = n == 1 ? 1 : sum(size(arr, 1) for arr in arrs[1:n-1]) + 1
#     to = from + size(a, 1) - 1
#     return dy[from:to, [(:) for i=1:length(size(dy)) - 1]...]
# end

# function unhcat(dy::AbstractArray, n::Int, arrs::AbstractArray...)
#     a = arrs[n]
#     from = n == 1 ? 1 : sum(size(arr, 2) for arr in arrs[1:n-1]) + 1
#     to = from + size(a, 2) - 1
#     return dy[:, from:to, [(:) for i=1:length(size(dy)) - 2]...]
# end

function uncat(dy::AbstractArray, n::Int, arrs::AbstractArray...; dims)
    @assert(dims isa Integer, "Can only undo cat() over a single dimension, " *
            "but dimensions $dims were provided")
    dim = dims
    a = arrs[n]
    from = n == 1 ? 1 : sum(size(arr, dim) for arr in arrs[1:n-1]) + 1
    to = from + size(a, dim) - 1
    return dy[[(:) for i=1:dim - 1]..., from:to, [(:) for i=1:length(size(dy)) - dim]...]
end

namedtuple(names, values) = NamedTuple{names}(values)
namedtuple(d::Dict) = NamedTuple{tuple(keys(d)...)}(values(d))

function rev_perm(perm::NTuple{N, Int}) where N
    rperm = Vector{Int}(undef, length(perm))
    for (i, j) in enumerate(perm)
        rperm[j] = i
    end
    return tuple(rperm...)
end

function ∇permutedims(dy, perm)
    rperm = rev_perm(perm)
    return permutedims(dy, rperm)
end


text:
@testset "batched_mul" begin
    using NNlib: batched_mul, batched_mul!, batched_vec,
                 batched_adjoint, batched_transpose

    A = randn(Float32, 3,3,2);
    B = randn(Float32, 3,3,2);

    C = batched_mul(A, B)
    @test CuArray(C) ≈ batched_mul(CuArray(A), CuArray(B))

    Ct = batched_mul(batched_transpose(A), B)
    @test CuArray(Ct) ≈ batched_mul(batched_transpose(CuArray(A)), CuArray(B))

    Ca = batched_mul(A, batched_adjoint(B))
    @test CuArray(Ca) ≈ batched_mul(CuArray(A), batched_adjoint(CuArray(B)))

    # 5-arg batched_mul!
    C .= pi
    batched_mul!(C, A, B, 2f0, 3f0)
    cuCpi = CuArray(similar(C)) .= pi
    @test CuArray(C) ≈ batched_mul!(cuCpi, CuArray(A), CuArray(B), 2f0, 3f0)

    # PermutedDimsArray
    @test CuArray(Ct) ≈ batched_mul(PermutedDimsArray(CuArray(A), (2,1,3)), CuArray(B))

    D = permutedims(B, (1,3,2))
    Cp = batched_mul(batched_adjoint(A), B)
    @test CuArray(Cp) ≈ batched_mul(batched_adjoint(CuArray(A)), PermutedDimsArray(CuArray(D), (1,3,2)))

    # Methods which reshape
    M = randn(Float32, 3,3)

    Cm = batched_mul(A, M)
    @test CuArray(Cm) ≈ batched_mul(CuArray(A), CuArray(M))

    Cv = batched_vec(permutedims(A,(3,1,2)), M)
    @test CuArray(Cv) ≈ batched_vec(PermutedDimsArray(CuArray(A),(3,1,2)), CuArray(M))
end

@testset "NNlib storage_type etc." begin
    using LinearAlgebra
    using NNlib: is_strided, are_strided, storage_type

    M = cu(ones(10,10))

    @test is_strided(M)
    @test is_strided(view(M, 1:2:5,:))
    @test is_strided(PermutedDimsArray(M, (2,1)))

    @test !is_strided(reshape(view(M, 1:2:10,:), 10,:))
    @test !is_strided((M .+ im)')
    @test !is_strided(Diagonal(cu(ones(3))))

    @test storage_type(M) <: CuArray{Float32,2}
    @test storage_type(reshape(view(M, 1:2:10,:), 10,:)) <: CuArray{Float32,2}
end


text:
#include("../get_home_dir.jl")

using DataFrames
using CSV
using PyPlot
using Random
using StatsBase

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams");
rcParams["lines.linewidth"] = 2;
rcParams["font.size"] = 24;
rcParams["mathtext.default"] = "regular";

#leafRMSE =  [0.04030680922286562, 0.1114426954163669, 0.19222915929366224]
#etRMSE = [5.11686931714865e-05, 3.992463510637098e-05, 5.691494713989705e-05]*18.02/1000*60*60*24;
#rzRMSE = [0.0030338018951135853, 0.003058330701815462, 0.003075149421206762];

leafRMSE = [0.04030680922286563, 0.11144269541636687, 0.19222915929366224];
etRMSE = [5.11686931714865e-05, 3.992463510637098e-05, 5.691494713989705e-05];
rzRMSE = [0.0012556391476988972, 0.003711247616340684, 0.0016033746735104248];

mylabels = ["All", "1 AM/PM", "6 AM/PM"];

#=
figure()
subplot(131)
bar([1,2,3],leafRMSE)
xticks([1,2,3], mylabels)
title("Leaf water potential (MPa)")

subplot(132)
bar([1,2,3],rzRMSE)
xticks([1,2,3], mylabels)
title("Column-average soil moisture")

subplot(133)
bar([1,2,3],etRMSE)
xticks([1,2,3], mylabels)
title(L"Daily ET $(mm/m^2/s)$")
ticklabel_format(style="plain",axis="y")

#tight_layout()
=#

gsRMSE = [0.0009707089849704303, 0.000809590695874667, 0.0012329717941626317]
leafscaleRMSE = [0.02781028047973449, 0.01783460829728346, 0.03495551737372624]
betaRMSE = [0.01874038306283515, 0.013533028689382927, 0.023987804499910826]

figure()
subplot(231)
bar([1,2,3],leafRMSE)
xticks([1,2,3], mylabels)
title("Leaf water potential (MPa)")

subplot(232)
bar([1,2,3],leafscaleRMSE)
xticks([1,2,3], mylabels)
title("LWP divided by P63")
ticklabel_format(style="plain",axis="y")

subplot(233)
bar([1,2,3],betaRMSE)
xticks([1,2,3], mylabels)
title("Beta factor")
ticklabel_format(style="plain",axis="y")

subplot(234)
bar([1,2,3],gsRMSE)
xticks([1,2,3], mylabels)
title(L"Stomatal conductance $(mol/m^2/s)$")

subplot(235)
bar([1,2,3],etRMSE)
xticks([1,2,3], mylabels)
title(L"Daily ET $(mm/m^2/s)$")
ticklabel_format(style="plain",axis="y")

leafmean = convert(Array, CSV.read("leaf_mean_tab.csv",DataFrame,header=false));
leafscale = convert(Array, CSV.read("leaf_scale_tab.csv",DataFrame,header=false));

figure()
for i in 1:3
	plot(leafmean[i,:], 1 ./ leafscale[i,:],"o",label=mylabels[i])
end
plot(leafmean[1,end], 1 ./ leafscale[1,end],"k+",label="True model",markersize=25,markeredgewidth=3)
xlabel("Mean LWP (MPa)")
ylabel("Stomatal P63 (MPa)")
legend()


text:
using BenchmarkTools, DataFrames, MixedModels, RData, Tables

const SUITE = BenchmarkGroup()

const dat = Dict(Symbol(k) => v for (k, v) in load(joinpath(
    dirname(pathof(MixedModels)),
    "..",
    "test",
    "dat.rda",
)));

categorical!(dat[:ml1m], [:G,:H]);  # forgot to convert these grouping factors

const mods = Dict{Symbol,Vector{Expr}}(
    :Alfalfa => [:(1 + A * B + (1 | G)), :(1 + A + B + (1 | G))],
    :Animal => [:(1 + (1 | G) + (1 | H))],
    :Arabidopsis => [],              # glmm and rename variables
    :Assay => [:(1+A+B*C+(1|G)+(1|H))],
    :AvgDailyGain => [:(1 + A * U + (1 | G)), :(1 + A + U + (1 | G))],
    :BIB => [:(1 + A * U + (1 | G)), :(1 + A + U + (1 | G))],
    :Bond => [:(1 + A + (1 | G))],
    :Chem97 => [:(1 + (1 | G) + (1 | H)), :(1 + U + (1 | G) + (1 | H))],
    :Contraception => [],            # glmm and rename variables
    :Cultivation => [:(1 + A * B + (1 | G)), :(1 + A + B + (1 | G)), :(1 + A + (1 | G))],
    :Demand => [:(1 + U + V + W + X + (1 | G) + (1 | H))],
    :Dyestuff => [:(1 + (1 | G))],
    :Dyestuff2 => [:(1 + (1 | G))],
    :Early => [:(1 + U + U & A + (1 + U | G))],
    :Exam => [:(1 + A * U + B + (1 | G)), :(1 + A + B + U + (1 | G))],
    :Gasoline => [:(1 + U + (1 | G))],
    :Gcsemv => [:(1 + A + (1 | G))], # variables must be renamed
    :Genetics => [:(1 + A + (1 | G) + (1 | H))],
    :HR => [:(1 + A * U + V + (1 + U | G))],
    :Hsb82 => [:(1 + A + B + C + U + (1 | G))],
    :IncBlk => [:(1 + A + U +  + W + Z + (1 | G))],
    :InstEval => [:(1 + A + (1 | G) + (1 | H) + (1 | I)), :(1 + A * I + (1 | G) + (1 | H))],
    :KKL => [],                      # variables must be renamed
    :KWDYZ => [],                    # variables must be renamed
    :Mississippi => [:(1 + A + (1 | G))],
    :Mmmec => [],                    # glmm (and offset) and variables renamed
    :Multilocation => [:(1 + A + (0 + A | G) + (1 | H))],
    :Oxboys => [:(1 + U + (1 + U | G))],
    :PBIB => [:(1 + A + (1 | G))],
    :Pastes => [:(1 + (1 | G) + (1 | H))],
    :Penicillin => [:(1 + (1 | G) + (1 | H))],
    :Pixel => [:(1 + U + V + (1 + U | G) + (1 | H))],  # variables must be renamed
    :Poems => [:(1 + U + V + W + (1 | G) + (1 | H) + (1 | I))],
    :Rail => [:(1 + (1 | G))],
    :SIMS => [:(1 + U + (1 + U | G))],
    :ScotsSec => [:(1 + A + U + V + (1 | G) + (1 | H))],
    :Semi2 => [:(1 + A + (1 | G) + (1 | H))],
    :Semiconductor => [:(1 + A * B + (1 | G))],
    :Socatt => [],                   # variables must be renamed - binomial glmm?
    :TeachingII => [:(1 + A + T + U + V + W + X + Z + (1 | G))],
    :VerbAgg => [:(1 + A + B + C + U + (1 | G) + (1 | H))], # Bernoulli glmm and rename variables
    :Weights => [:(1 + A * U + (1 + U | G))],
    :WWheat => [:(1 + U + (1 + U | G))],
    :bdf => [],                      # rename variables and look up model
    :bs10 => [:(1 + U + V + W + ((1 + U + V + W) | G) + ((1 + U + V + W) | H))],
    :cake => [:(1 + A * B + (1 | G))],
    :cbpp => [:(1 + A + (1 | G))],   # Binomial glmm, create and rename variables
    :d3 => [
        :(1 + U + (1 | G) + (1 | H) + (1 | I)),
        :(1 + U + (1 + U | G) + (1 + U | H) + (1 + U | I)),
    ],
    :dialectNL => [:(1 + A + T + U + V + W + X + (1 | G) + (1 | H) + (1 | I))],
    :egsingle => [:(1 + A + U + V + (1 | G) + (1 | H))],
    :epilepsy => [],                 # unknown origin
    :ergoStool => [:(1 + A + (1 | G))],
    :gb12 => [:(1 + S + T + U + V + W + X + Z + ((1 + S + U + W) | G) +
                ((1 + S + T + V) | H))],
    :grouseticks => [],              # rename variables, glmm needs formula
    :guImmun => [],                  # rename variables, glmm needs formula
    :guPrenat => [],                 # rename variables, glmm needs formula
    :kb07 => [
        :(1 + S + T + U + V + W + X + Z + ((1 + S + T + U + V + W + X + Z) | G) +
          ((1 + S + T + U + V + W + X + Z) | H)),
        :(1 + S + T + U + V + W + X + Z +
          zerocorr((1 + S + T + U + V + W + X + Z) | G) +
          zerocorr((1 + S + T + U + V + W + X + Z) | H)),
    ],
    :ml1m => [:(1 + (1 | G) + (1 | H))],
    :paulsim => [:(1 + S + T + U + (1 | H) + (1 | G))],  # names of H and G should be reversed
    :sleepstudy => [:(1 + U + (1 + U | G)), :(1 + U + zerocorr(1 + U | G))],
    :s3bbx => [],                    # probably drop this one
    :star => [],                     # not sure it is worthwhile working with these data
);

fitbobyqa(rhs::Expr, dsname::Symbol) =
    fit(MixedModel, @eval(@formula(Y ~ $rhs)), dat[dsname])
compactstr(ds, rhs) = replace(string(ds, ':', rhs), ' ' => "")

SUITE["simplescalar"] = BenchmarkGroup(["single", "simple", "scalar"])
for ds in [
    :Alfalfa,
    :AvgDailyGain,
    :BIB,
    :Bond,
    :cake,
    :Cultivation,
    :Dyestuff,
    :Dyestuff2,
    :ergoStool,
    :Exam,
    :Gasoline,
    :Hsb82,
    :IncBlk,
    :Mississippi,
    :PBIB,
    :Rail,
    :Semiconductor,
    :TeachingII,
]
    for rhs in mods[ds]
        SUITE["simplescalar"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
            $(QuoteNode(rhs)),
            $(QuoteNode(ds)),
        )
    end
end

SUITE["singlevector"] = BenchmarkGroup(["single", "vector"])
for ds in [:Early, :HR, :Oxboys, :SIMS, :sleepstudy, :Weights, :WWheat]
    for rhs in mods[ds]
        SUITE["singlevector"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
            $(QuoteNode(rhs)),
            $(QuoteNode(ds)),
        )
    end
end

SUITE["nested"] = BenchmarkGroup(["multiple", "nested", "scalar"])
for ds in [:Animal, :Chem97, :Genetics, :Pastes, :Semi2]
    for rhs in mods[ds]
        SUITE["nested"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
            $(QuoteNode(rhs)),
            $(QuoteNode(ds)),
        )
    end
end

SUITE["crossed"] = BenchmarkGroup(["multiple", "crossed", "scalar"])

for ds in [
    :Assay,
    :Demand,
    :InstEval,
    :Penicillin,
    :ScotsSec,
    :dialectNL,
    :egsingle,
    :ml1m,
    :paulsim,
]
    for rhs in mods[ds]
        SUITE["crossed"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
            $(QuoteNode(rhs)),
            $(QuoteNode(ds)),
        )
    end
end

SUITE["crossedvector"] = BenchmarkGroup(["multiple", "crossed", "vector"])
for ds in [:bs10, :d3, :gb12, :kb07]
    for rhs in mods[ds]
        SUITE["crossedvector"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
            $(QuoteNode(rhs)),
            $(QuoteNode(ds)),
        )
    end
end


text:
#=
Main functions to build MPS-circuit.
=#
export MPSpar, MPSC, MPSDCpar

"""
    MPSpar(nBitA::Int64, vBit::Int64, rBit::Int64)
Construct parameters that MPSC needs.
\n
\nFields:
\n`nBitA::Int64`:  Number of lines(bits) in the QMPS Extended circuit.
\n`vBit::Int64`:   Number of virtual bits in the QMPS circuit.
\n`rBit::Int64`:   Number of reusable bits in the QMPS circuit.
\n`nBit::Int64`:   Number of lines(bits) in the QMPS circuit.
\n`nBlock::Int64`: Number of  MPS blocks in the QMPS circuit.
"""
struct MPSpar
    nBitA::Int64  # Number of lines(bits) in the QMPS Extended circuit.
    vBit::Int64   # Number of virtual bits in the QMPS circuit.
    rBit::Int64   # Number of reusable bits in the QMPS circuit.
    nBit::Int64   # Number of lines(bits) in the QMPS circuit.
    nBlock::Int64 # Number of  MPS blocks in the QMPS circuit.

    function MPSpar(nBitA::Int64, vBit::Int64, rBit::Int64)
        if (nBitA - vBit -rBit) % rBit != 0
            error("Error: nBlock is not an integer.")
        end
        nBlock = Int((nBitA - vBit) / rBit)
        nBit = rBit + vBit
        new(nBitA, vBit, rBit, nBit, nBlock)
    end
end

"""
    MPSDCpar(circuit::ChainBlock)
Get the circuit parameters of a differentiable QMPS circuit (QMPS-DC) or of a QMPS-DC extended circuit.
\n
\nFields:
\n`nBitA::Int64`:  Number of qubits(lines) in the QMPS-DC extended circuit.
\n`vBit::Int64`:   Number of virtual qubits in the circuit.
\n`rBit::Int64`:   Number of reusable qubits in the circuit.
\n`nBit::Int64`:   Number of qubits(lines) in the circuit.
\n`nBlock::Int64`: Number of  MPS blocks in the circuit.
\n`depth::Int64`:  Depth(number of layers) of each DC block in the circuit.
"""
struct MPSDCpar
    nBitA::Int64  # Number of qubits(lines) in the MPS-DC extended circuit.
    vBit::Int64   # Number of virtual qubits in the MPS-DC circuit.
    rBit::Int64   # Number of reusable qubits in the MPS-DC circuit.
    nBit::Int64   # Number of qubits(lines) in the MPS-DC circuit.
    nBlock::Int64 # Number of  MPS blocks in the MPS-DC circuit or MPS-DC extended circuit.
    depth::Int64  # Depth(number of layers) of each DC block in the MPS-DC circuit.

    function MPSDCpar(circuit::ChainBlock)
        rBit = 0
        nBit = 0
        vBit = 0
        depth = 0
        nBitA = 0
        nBlock = 0
        try
            if typeof(circuit[1]) <: AbstractContainer
                # MPSC().cExtend[1]::AbstractContainer
                rBit = content(circuit[1])[2][1].locs[1]
                nBit = nqubits(content(circuit[1]))
                vBit = nBit - rBit
                depth = length(content(circuit[1])[1])
                nBitA = nqubits(circuit)
            elseif typeof(circuit[1]) <: CompositeBlock
                # MPSC().circuit[1]::ChainBlock
                rBit = circuit[1][1][2][1].locs[1]
                nBit = nqubits(circuit)
                vBit = nBit - rBit
                depth = length(circuit[1][1][1])
                nBitA = vBit + rBit*length(circuit)
            end
            nBlock = Int((nBitA - vBit) / rBit)
        catch err
            if isa(err, DomainError) || (0 in (nBitA, vBit, rBit, nBit, nBlock, depth))
                error("ERROR: Input circuit is not supported by the function!")
            end
        end
        new(nBitA, vBit, rBit, nBit, nBlock, depth)
    end
end

"""
    MPSC(blockT::Union{String, Tuple{String, Int64}}, nBitA::Int64, vBit::Int64, rBit::Int64=1; dBlocksPar::Array{Float64,1}=[0.0])
Structure of elements related to a QMPS circuit.
\n`blockT::Union{String, Tuple{String, Int64}}`
\n1) `blockT = ("DC", depth)`
\n   "DC" stands for "differentiable circuit".
\n2) `blockT = "CS"`
\n   "CS" stands for "cluster state".
\n`dBlocksPar::Array{Float64,1}`
\n   The array of the parameters (for differentiable blocks) in the circuit.
\n
\nFields:
\n`circuit::ChainBlock`:                QMPS circuit.
\n`mpsBlocks::Array{CompositeBlock,1}`: Array of all the MPS blocks in the QMPS circuit.
\n`cExtend::ChainBlock`:                The circuit QMPS circuit is extended back to that doesn't reuse any qubit.
\n`cEBlocks::Array{CompositeBlock,1}`:  Array of all the MPS blocks in the Extended circuit.
\n`dGates::Array{QMPS.QDiff,1}`:        Differentiable gates of the QMPS circuit if applicable.
\n`nBit::Int64`:                        Number of lines (bits) of the QMPS circuit.
\n`nBlock::Int64`:                      Number of blocks in the QMPS circuit.
"""
struct MPSC
    circuit::ChainBlock                # QMPS circuit.
    mpsBlocks::Array{CompositeBlock,1} # Array of all the MPS blocks in the QMPS circuit.
    cExtend::ChainBlock                # The circuit QMPS circuit is extended back to that doesn't reuse any qubit.
    cEBlocks::Array{CompositeBlock,1}  # Array of all the MPS blocks in the Extended circuit.
    dGates::Array{QMPS.QDiff,1}        # Differentiable gates of the QMPS circuit if applicable.
    nBit::Int64                        # Number of lines (bits) of the QMPS circuit.
    nBlock::Int64                      # Number of blocks in the QMPS circuit.

    function MPSC(blockT::Union{String, Tuple{String, Int64}}, nBitA::Int64, vBit::Int64, rBit::Int64; dBlocksPar::Array{Float64,1}=[0.0])
        par2nd = MPSpar(nBitA, vBit, rBit)
        nBlock = par2nd.nBlock
        nBit = par2nd.nBit
        MPS = MPSbuilder(nBitA, vBit, rBit, blockT)
        circuit = MPS.circuit
        cExtend = MPS.cExtend
        dispatch!(circuit, :random)
        dGates = collect_blocks(QMPS.QDiff, circuit)
        if dBlocksPar != [0.0]
            if length(dGates) == length(dBlocksPar)
                pars = [parameters(dGates[i])[1] for i=1:length(dGates)]
                pars .= dBlocksPar
                dispatch!.(dGates, pars)
            elseif blockT == "CS"
                error("Cluster states don't support optional argument `dBlocksPar`.")
            else
                error("The number of elements in `dBlocksPar` is not correct.")
            end
        end
        mpsBlocks = CompositeBlock[]
        for i = 1:length(circuit)-1
            push!(mpsBlocks, circuit[i][1])
        end
        push!(mpsBlocks, circuit[end])
        new(circuit, mpsBlocks, cExtend, cExtend.blocks, dGates, nBit, nBlock)
    end
end


text:
using TrussMorph
tm = TrussMorph
using Optim
using Printf

function compute_morph_path(t0::tm.Truss, t1::tm.Truss, load::Matrix{Float64},
    node_dof::Int, full_node_dof::Int, design_var_ids::Vector{Int}; path_disc::Int=5, parm_smooth::Float64=100.0, parm_weight::Float64=50.0, do_opt::Bool=true, kwargs...)

    @assert(size(t0.X, 1) == size(t1.X, 1))
    @assert(t0.T == t1.T)

    n_v = size(t0.X,1)
    n_e = size(t0.T,1)
    # initial cross secs
    r = ones(size(t0.T,1)) * 0.01
    F = tm.assemble_load_vector(load, n_v, node_dof)
    perm_vec, perm, n_dof_free = tm.dof_permutation(t0.S, n_v, node_dof)

    F_perm = perm * F
    F_m = Array(F_perm[1:n_dof_free])

    var_chuck = length(design_var_ids)::Int
    X0_var = reshape(t0.X', prod(size(t0.X)))[design_var_ids]
    X1_var = reshape(t1.X', prod(size(t1.X)))[design_var_ids]

    # TODO: symmetric variable
    weight_fn = tm.get_weight_calculation_fn(t0.X, r, t0.T, F_m, perm, n_dof_free,
    node_dof, full_node_dof, t0.mp, design_var_ids)

    # set up Xpath0
    Xpath0 = Vector{Float64}(undef, path_disc * var_chuck)
    for i=1:path_disc
        Xpath0[(i-1)*var_chuck + 1 : i*var_chuck] =
            Float64.((path_disc - i + 1) / (path_disc + 1)) * X0_var +
            Float64.(i / (path_disc + 1)) * X1_var
    end

    ### for |d^2ϕ/dt^2|
    # L_I = collect(1:path_disc+1)
    # L_J = collect(2:path_disc+2)
    # L_V = ones(Float64, path_disc+1)
    # L = sparse(L_I, L_J, L_V, path_disc+2, path_disc+2)
    # L = L + L'
    # rowsums = vec(sum(L,dims=2))
    # L = L - spdiagm(0 => rowsums)

    function path_energy(Xpath::Vector{Float64})
        Xmat = vcat(X0_var', reshape(Xpath, (var_chuck, path_disc))', X1_var')
        dXpath_dt = Xmat[2:end, :] - Xmat[1:end-1, :]
        path_weight_vec = zeros(Float64, size(Xmat, 1))
        for i=1:size(Xmat, 1)
            path_weight_vec[i] = weight_fn(Xmat[i,:])
        end
        # ddt = (L * path_weight_vec).^2

        # ∑ |dϕ/dt|^2 + |WF(X(t))|
        # return parm_smooth * sum(dXpath_dt.^2) + parm_weight * sum(path_weight_vec)

        # ∑ |dϕ/dt|^2 + |dWF/dt|^2
        # return parm_smooth * sum(dXpath_dt.^2) + parm_weight * sum((path_weight_vec[2:end] - path_weight_vec[1:end-1]).^2)

        # ∑ |dWF/dt|^2
        return sum((path_weight_vec[2:end] - path_weight_vec[1:end-1]).^2)

        # ∑ |d^2WF/dt^2|^2
        # return sum(ddt[2:end-1])
        # return sum(ddt)
    end

    @show init_sumE = path_energy(Xpath0)
    opt_sumE = init_sumE

    # run opt
    if do_opt
        if :box_constraint in keys(kwargs) && kwargs[:box_constraint] != undef
            bc_var = kwargs[:box_constraint]
            @assert(size(bc_var, 1) == length(design_var_ids))
            @assert(size(bc_var, 2) == 2)
            println("using box constraint...")

            Xpath_bc_lower = reshape(bc_var[:,1] * ones(path_disc)', path_disc * var_chuck)
            Xpath_bc_upper = reshape(bc_var[:,2] * ones(path_disc)', path_disc * var_chuck)
            # @assert(length(Xpath_bc_lower) == length(Xpath0))
            # @assert(length(Xpath_bc_upper) == length(Xpath0))

            # inner_optimizer = GradientDescent()
            # res = Optim.optimize(path_energy, box_constraint[:,1], box_constraint[:,2], Xpath0, Fminbox(inner_optimizer))
            res = Optim.optimize(path_energy, Xpath_bc_lower, Xpath_bc_upper, Xpath0)
        else
            println("unconstrained...")
            # Optim constructs an approximate gradient for us using central finite differencing
            res = Optim.optimize(path_energy, Xpath0, LBFGS())
        end

        @show Optim.summary(res)
        @show Optim.iterations(res)
        @show Optim.iteration_limit_reached(res)
        @show Optim.f_calls(res)
        @show Optim.converged(res)

        X_var_star_vec = Optim.minimizer(res)
        opt_sumE = Optim.minimum(res)
        @show opt_sumE
    else
        X_var_star_vec = Xpath0
    end

    X_var_star = reshape(X_var_star_vec, (var_chuck, path_disc))'
    X_var_init = reshape(Xpath0, (var_chuck, path_disc))'

    X_template = copy(t0.X)
    dim = size(t0.X,2)
    design_var_id_map = zeros(Int, var_chuck, 2)
    for i=1:length(design_var_ids)
        design_var_id_map[i,1] = Int.((design_var_ids[i] + design_var_ids[i]%dim) / dim)
        design_var_id_map[i,2] = dim - design_var_ids[i]%dim
        X_template[design_var_id_map[i,1], design_var_id_map[i,2]] = 0
    end

    morph_path = Array{Matrix{Float64}}(undef, path_disc + 2)
    morph_path[1] = t0.X
    morph_path[end] = t1.X

    for i=2:path_disc+1
        morph_path[i] = copy(X_template)
        for j=1:var_chuck
            morph_path[i][design_var_id_map[j,1], design_var_id_map[j,2]] = X_var_star[i-1, j]
        end
    end

    init_morph_path = Array{Matrix{Float64}}(undef, path_disc + 2)
    init_morph_path[1] = t0.X
    init_morph_path[end] = t1.X

    for i=2:path_disc+1
        init_morph_path[i] = copy(X_template)
        for j=1:var_chuck
            init_morph_path[i][design_var_id_map[j,1], design_var_id_map[j,2]] = X_var_init[i-1, j]
        end
    end

    # function map_full_X_to_design_vars(origX::Matrix{Float64})::Vector{Float64}
    #     return reshape(origX', prod(size(origX)))[design_var_ids]
    # end
    function cal_ptwise_energies(X_var_vec::Vector{Float64})
        X_var = reshape(X_var_vec, (var_chuck, path_disc))'
        ptwise_smoothness_energy = Float64[]
        ptwise_weight_energy = Float64[]
        ptwise_total_energy = Float64[]
        push!(ptwise_smoothness_energy, 0.5 * sum((X0_var - X_var[1,:]).^2))
        push!(ptwise_weight_energy, weight_fn(X0_var))
        push!(ptwise_total_energy, parm_smooth * ptwise_smoothness_energy[end] +
                                   parm_weight * ptwise_weight_energy[end])

        for i=1:path_disc
            push!(ptwise_weight_energy, weight_fn(X_var[i, :]))
            if 1 == i
                prevX_var = X0_var
            else
                prevX_var = X_var[i-1,:]
            end
            if path_disc == i
                nextX_var = X1_var
            else
                nextX_var = X_var[i+1,:]
            end
            push!(ptwise_smoothness_energy, 0.5 * (sum((X_var[i, :] - prevX_var).^2) + sum((X_var[i, :] - nextX_var).^2)))
            push!(ptwise_total_energy, parm_smooth * ptwise_smoothness_energy[end] + parm_weight * ptwise_weight_energy[end])
        end

        push!(ptwise_smoothness_energy, 0.5 * sum((X1_var - X_var[end,:]).^2))
        push!(ptwise_weight_energy, weight_fn(X1_var))
        push!(ptwise_total_energy, parm_smooth * ptwise_smoothness_energy[end] +
                                   parm_weight * ptwise_weight_energy[end])
        return ptwise_smoothness_energy, ptwise_weight_energy, ptwise_total_energy
    end

    opt_sE, opt_wE, opt_totE = cal_ptwise_energies(X_var_star_vec)
    init_sE, init_wE, init_totE = cal_ptwise_energies(Xpath0)

    # return intial morph path and energy as a comparision
    return morph_path, init_morph_path, opt_sE, opt_wE, opt_totE, init_sE, init_wE, init_totE, opt_sumE, init_sumE, [parm_smooth, parm_weight]
end


text:
module Symmetry

using CoordinateTransformations: AffineMap, Translation, LinearMap
using LinearAlgebra: I, diagm
using StaticArrays: SVector, SMatrix, SDiagonal

using Crystallography

import LinearAlgebra: det, tr

export SeitzOperator,
    symmetrytype,
    isidentity,
    istranslation,
    ispointsymmetry,
    genpath,
    pearsonsymbol,
    arithmeticclass

pearsonsymbol(::Triclinic) = "a"
pearsonsymbol(::Monoclinic) = "m"
pearsonsymbol(::Orthorhombic) = "o"
pearsonsymbol(::Tetragonal) = "t"
pearsonsymbol(::Cubic) = "c"
pearsonsymbol(::Hexagonal) = "h"
pearsonsymbol(::Trigonal) = "h"
pearsonsymbol(::Primitive) = "P"
pearsonsymbol(::BaseCentering{T}) where {T} = string(T)
pearsonsymbol(::BodyCentering) = "I"
pearsonsymbol(::FaceCentering) = "F"
pearsonsymbol(::RhombohedralCentering) = "R"
pearsonsymbol(b::Bravais) = pearsonsymbol(crystalsystem(b)) * pearsonsymbol(centering(b))

arithmeticclass(::Triclinic) = "-1"
arithmeticclass(::Monoclinic) = "2/m"
arithmeticclass(::Orthorhombic) = "mmm"
arithmeticclass(::Tetragonal) = "4/mmm"
arithmeticclass(::Hexagonal) = "6/mmm"
arithmeticclass(::Cubic) = "m-3m"
arithmeticclass(::Primitive) = "P"
arithmeticclass(::BaseCentering) = "S"
arithmeticclass(::BodyCentering) = "I"
arithmeticclass(::FaceCentering) = "F"
arithmeticclass(::RhombohedralCentering) = "R"
arithmeticclass(b::Bravais) =
    arithmeticclass(crystalsystem(b)) * arithmeticclass(centering(b))
arithmeticclass(::RCenteredHexagonal) = "-3mR"

abstract type PointSymmetry end
struct Identity <: PointSymmetry end
struct RotationAxis{N} <: PointSymmetry end
struct Inversion <: PointSymmetry end
struct RotoInversion{N} <: PointSymmetry end
const Mirror = RotoInversion{2}
RotationAxis(N::Int) =
    N ∈ (2, 3, 4, 6) ? RotationAxis{N}() :
    throw(ArgumentError("rotation axis must be either 2, 3, 4 or 6!"))
RotoInversion(N::Int) =
    N ∈ (2, 3, 4, 6) ? RotoInversion{N}() :
    throw(ArgumentError("rotoinversion axis must be either 2, 3, 4 or 6!"))

struct PointSymmetryPower{T<:PointSymmetry,N} end
PointSymmetryPower(X::PointSymmetry, N::Int) = PointSymmetryPower{typeof(X),N}()
function PointSymmetryPower(::RotationAxis{N}, M::Int) where {N}
    if M >= N
        if M == N
            Identity()
        else
            PointSymmetryPower(RotationAxis(N), M - N)  # Recursive call
        end
    else  # Until M < N
        PointSymmetryPower{RotationAxis{N},M}()
    end
end

function symmetrytype(trace, determinant)
    return Dict(
        (3, 1) => Identity(),
        (-1, 1) => RotationAxis(2),
        (0, 1) => RotationAxis(3),
        (1, 1) => RotationAxis(4),
        (2, 1) => RotationAxis(6),
        (-3, -1) => Inversion(),
        (1, -1) => Mirror(),
        (0, -1) => RotoInversion(3),
        (-1, -1) => RotoInversion(4),
        (-2, -1) => RotoInversion(6),
    )[(trace, determinant)]
end
symmetrytype(op::PointSymmetry) = symmetrytype(tr(op), det(op))

"""
    SeitzOperator(m::AbstractMatrix)

Construct a `SeitzOperator` from a 4×4 matrix.
"""
struct SeitzOperator{T}
    data::SMatrix{4,4,T,16}
end
SeitzOperator(m::AbstractMatrix) = SeitzOperator(SMatrix{4,4}(m))
function SeitzOperator(l::LinearMap)
    m = l.linear
    @assert size(m) == (3, 3)
    x = diagm(ones(eltype(m), 4))
    x[1:3, 1:3] = m
    return SeitzOperator(x)
end # function PointSymmetryOperator
function SeitzOperator(t::Translation)
    v = t.translation
    @assert length(v) == 3
    x = diagm(ones(eltype(v), 4))
    x[1:3, 4] = v
    return SeitzOperator(x)
end
"""
    SeitzOperator(l::LinearMap, t::Translation)
    SeitzOperator(a::AffineMap)
    SeitzOperator(l::LinearMap)
    SeitzOperator(t::Translation)

Construct a `SeitzOperator` from `LinearMap`s, `Translation`s or `AffineMap`s by the
following equation:

```math
S = \\left(
\\begin{array}{ccc|c}
& & & \\\\
& R & & {\\boldsymbol \\tau} \\\\
& & & \\\\
\\hline 0 & 0 & 0 & 1
\\end{array}
\\right),
```

where ``R`` is a `LinearMap` and ``{\\boldsymbol \\tau}`` is a `Translation`.
"""
function SeitzOperator(l::LinearMap, t::Translation)
    m, v = l.linear, t.translation
    return SeitzOperator(vcat(hcat(m, v), [zeros(eltype(m), 3)... one(eltype(v))]))
end # function SeitzOperator
SeitzOperator(a::AffineMap) = SeitzOperator(LinearMap(a.linear), Translation(a.translation))
"""
    SeitzOperator(s::SeitzOperator, pos::AbstractVector)

Construct a `SeitzOperator` that locates at `pos` from a `SeitzOperator` passing through the
origin.
"""
function SeitzOperator(s::SeitzOperator, pos::AbstractVector)
    @assert length(pos) == 3
    t = SeitzOperator(Translation(pos))
    return t * s * inv(t)
end # function SeitzOperator
(op::SeitzOperator)(v::AbstractVector) = (op.data*[v; 1])[1:3]

isidentity(op::SeitzOperator) = op.data == I

function istranslation(op::SeitzOperator)
    m = op.data
    if m[1:3, 1:3] != I || !(iszero(m[4, 1:3]) && isone(m[4, 4]))
        return false
    end
    return true
end

function ispointsymmetry(op::SeitzOperator)
    m = op.data
    if !(
        iszero(m[4, 1:3]) &&
        iszero(m[1:3, 4]) &&
        isone(m[4, 4]) &&
        abs(det(m[1:3, 1:3])) == 1
    )
        return false
    end
    return true
end

"""
    genpath(nodes, densities)
    genpath(nodes, density::Integer, iscircular = false)

Generate a reciprocal space path from each node.

# Arguments
- `nodes::AbstractVector`: a vector of 3-element k-points.
- `densities::AbstractVector{<:Integer}`: number of segments between current node and the next node.
- `density::Integer`: assuming constant density between nodes.

# Examples
```jldoctest
julia> nodes = [
    [0.0, 0.0, 0.5],
    [0.0, 0.0, 0.0],
    [0.3333333333, 0.3333333333, 0.5],
    [0.3333333333, 0.3333333333, -0.5],
    [0.3333333333, 0.3333333333, 0.0],
    [0.5, 0.0, 0.5],
    [0.5, 0.0, 0.0]
];

julia> genpath(nodes, 100, true)  # Generate a circular path
700-element Array{Array{Float64,1},1}:
...

julia> genpath(nodes, 100, false)  # Generate a noncircular path
600-element Array{Array{Float64,1},1}:
...

julia> genpath(nodes, [10 * i for i in 1:6])  # Generate a noncircular path
210-element Array{Array{Float64,1},1}:
...
```
"""
function genpath(nodes, densities)
    if length(densities) == length(nodes) || length(densities) == length(nodes) - 1
        path = similar(nodes, sum(densities .- 1) + length(densities))
        s = 0
        for (i, (thisnode, nextnode, density)) in
            enumerate(zip(nodes, circshift(nodes, -1), densities))
            step = @. (nextnode - thisnode) / density
            for j in 1:density
                path[s+j] = @. thisnode + j * step
            end
            s += density
        end
        return path
    else
        throw(
            DimensionMismatch(
                "the length of `densities` is either `length(nodes)` or `length(nodes) - 1`!",
            ),
        )
    end
end
genpath(nodes, density::Integer, iscircular::Bool = false) =
    genpath(nodes, (density for _ in 1:(length(nodes)-(iscircular ? 0 : 1))))

Base.getindex(A::SeitzOperator, I::Vararg{Int}) = getindex(A.data, I...)

Base.one(::Type{SeitzOperator{T}}) where {T} =
    SeitzOperator(SDiagonal(SVector{4}(ones(T, 4))))
Base.one(A::SeitzOperator) = one(typeof(A))

Base.inv(op::SeitzOperator) = SeitzOperator(Base.inv(op.data))

Base.:*(::Identity, ::Identity) = Identity()
Base.:*(::Identity, x::Union{PointSymmetryPower,PointSymmetry}) = x
Base.:*(x::Union{PointSymmetryPower,PointSymmetry}, ::Identity) = x
Base.:*(::Inversion, ::Inversion) = Identity()
Base.:*(::RotationAxis{N}, ::Inversion) where {N} = RotoInversion(N)
Base.:*(i::Inversion, r::RotationAxis) = r * i
Base.:*(::RotationAxis{N}, ::RotationAxis{N}) where {N} =
    PointSymmetryPower(RotationAxis(N), 2)
Base.:*(::RotationAxis{2}, ::RotationAxis{2}) = Identity()
Base.:*(::PointSymmetryPower{RotationAxis{N},M}, ::RotationAxis{N}) where {N,M} =
    PointSymmetryPower(RotationAxis(N), M + 1)
Base.:∘(a::SeitzOperator, b::SeitzOperator) = SeitzOperator(a.data * b.data)

function Base.convert(::Type{Translation}, op::SeitzOperator)
    @assert(istranslation(op), "operator is not a translation!")
    return Translation(collect(op.data[1:3, 4]))
end
function Base.convert(::Type{LinearMap}, op::SeitzOperator)
    @assert(ispointsymmetry(op), "operator is not a point symmetry!")
    return LinearMap(collect(op.data[1:3, 1:3]))
end

tr(::Identity) = 3
tr(::RotationAxis) where {N} = N - 3  # 2, 3, 4
tr(::RotationAxis{6}) = 2
tr(::Inversion) = -3
tr(::RotoInversion{N}) where {N} = -tr(RotationAxis(N))
det(::Union{Identity,RotationAxis}) = 1
det(::Inversion) = -1
det(::RotoInversion) = -1

end


text:
using Multisets
using Test

A = Multiset(1:5)
B = A+A
@test length(B)==2*length(A)
@test A<B
@test A<=B
@test B>A
@test B>=A
B[1]=0
@test Set(B)==Set(2:5)
@test length(A|B) == 9
@test length(A&B) == 4
@test length(A*B) == 40
@test length(A+B) == 13
@test collect(A) == collect(1:5)
@test hash(A) != hash(B)


text:
using KitBase, Plots

set = Setup(
    "gas", # matter
    "sod", # case
    "1d2f1v", # space
    "kfvs", # flux
    "bgk", # collision
    1, # species
    2, # interpolation order
    "vanleer", # limiter
    "fix", # boundary
    0.5, # cfl
    0.15, # simulation time
)
ps = PSpace1D(0.0, 1.0, 100, 1)
vs = VSpace1D(-5.0, 5.0, 72)
gas = Gas(1e-4, 0.0, 1.0, 2.0)
fw, ff, bc = config_ib(set, ps, vs, gas)
ib = IB1F(fw, ff, bc)

ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks, ks.ps, structarray = true)

t = 0.0
dt = KitBase.timestep(ks, ctr, t)
nt = ks.set.maxTime ÷ dt |> Int
for iter = 1:nt
    reconstruct!(ks, ctr)
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, zeros(3))
end

plot(ks, ctr)
plot_line(ks, ctr)


text:
# Functions
import Polynomials: ImmutablePolynomial
# Macros
import MLStyle: @match
# Modules
import MPI

mutable struct PayloadRange
    a::Int32
    b::Int32
end

# PayloadRange(r::NTuple{2,<:Integer}) = PayloadRange(r[1], r[2])

function Base.copy(r::PayloadRange)
    return PayloadRange(r.a, r.b)
end

function Base.:(==)(r::PayloadRange, t::NTuple{2,<:Integer})
    return r.a == t[1] && r.b == t[2]
end

struct PayloadGroup
    range::PayloadRange
    ratio::Float64
end

PayloadGroup(range::Vector{<:Integer}, ratio::T) where {T} =
    PayloadGroup(PayloadRange(range[1], range[2]), get_ratio(ratio))

abstract type EvolutionFunction end

struct DataStreamConfig
    name::String
    evolve::EvolutionFunction
    payload_groups::Vector{PayloadGroup}
end

function get_payload_group_id(
    rank::Integer,
    nranks::Integer,
    cfg::DataStreamConfig,
)
    percentile = (rank + 1) / nranks
    current = 0.0
    for id = 1:length(cfg.payload_groups)
        grp = cfg.payload_groups[id]
        current += grp.ratio
        if percentile <= current
            return id
        end
    end
    # TODO: throw an error here ?
end

function get_payload_group_id(cfg::DataStreamConfig)
    comm = MPI.COMM_WORLD
    return get_payload_group_id(MPI.Comm_rank(comm), MPI.Comm_size(comm), cfg)
end

abstract type DataObject end

mutable struct DataVector <: DataObject
    vec::Vector{Float64}
end

DataVector() = DataVector(Float64[])

mutable struct DataStream
    initial_range::PayloadRange
    range::PayloadRange
    evolve::EvolutionFunction
    data::DataObject
end

DataStream(range::PayloadRange, evolve::EvolutionFunction) =
    DataStream(range, copy(range), evolve, DataVector())

DataStream(cfg::DataStreamConfig) =
    DataStream(cfg.payload_groups[get_payload_group_id(cfg)].range, cfg.evolve)

struct GrowthFactorEvFn <: EvolutionFunction
    factor::Float64
end

GrowthFactorEvFn(params::Vector{<:Real}) = GrowthFactorEvFn(params[1])

function evolve_payload_range!(
    stream::DataStream,
    step::Integer,
    fn::GrowthFactorEvFn,
)
    growth = fn.factor^step
    stream.range.a = round(stream.initial_range.a * growth)
    stream.range.b = round(stream.initial_range.b * growth)
end

struct PolynomialEvFn <: EvolutionFunction
    poly::ImmutablePolynomial
end

PolynomialEvFn(params::Vector{<:Real}) =
    PolynomialEvFn(ImmutablePolynomial(vcat(0, params)))

function evolve_payload_range!(
    stream::DataStream,
    step::Integer,
    fn::PolynomialEvFn,
)
    growth = fn.poly(step)
    stream.range.a = stream.initial_range.a + growth
    stream.range.b = stream.initial_range.b + growth
end

function evolve_payload_range!(stream::DataStream, step::Integer)
    evolve_payload_range!(stream, step, stream.evolve)
end

function check_length(expected::Integer, params::Vector{<:Real})
    if length(params) != expected
        @error "Wrong number of parameters"
    end
    return params
end

function get_evolution_function(evcfg::Config)
    params = evcfg[:params]
    @match evcfg[:function] begin
        "GrowthFactor" => return GrowthFactorEvFn(check_length(1, params))
        "Polynomial" => return PolynomialEvFn(params)
        "Linear" => return PolynomialEvFn(check_length(1, params))
        _ => @error "Unsupported stream size evolution function"
    end
end

get_evolution_function(nothing) = GrowthFactorEvFn(1.0)


text:
function get_data(num_classes::Int = 10)

    if dataset == "3dmnist"

        class_file = Dict()
        shapeids = Dict()
        class = Dict()
        datapath = Dict()

        classes = Dict{String, Int}()
        for i in 0:9
            classes["$(i)"] = i
        end

        shapeids["train"] = ["$(i)" for i in 1:5000]
        shapeids["test"] = ["$(i)" for i in 1:1000]

        for split_ in ["train","test"]
            class_file[split_] = joinpath(root, "labels_$(split_).txt")
            class[split_] = [line::String for line in readlines(class_file[split_])]
            datapath[split_] = [(class[split_][i], joinpath(root, split_, (shapeids[split_][i])*".txt")) for i in 1:length(shapeids[split_])]
        end

        #shuffling trainset Point
        rng = MersenneTwister(1234);
        shuffle!(rng, datapath["train"])

        # Fetching the train and validation data and getting them into proper shape
        datapath, classes = get_data(num_classes)
        X = [datapoint(p, classes) for p in datapath["train"]]
        data = [X[i][1] for i in 1:length(X)]
        labels = onehotbatch(cat([X[i][2][1] for i in 1:length(X)]..., dims=1), 0:9)
        train = [(cat(data[i]..., dims = 3), labels[:,i]) for i in partition(1:length(data), batch_size)]

        VAL = [datapoint(p, classes) for p in datapath["test"]]
        valX = cat([VAL[i][1] for i in 1:length(VAL)]..., dims=3)
        valY = onehotbatch(cat([VAL[i][2][1] for i in 1:length(VAL)]..., dims=1), 0:9)

    elseif dataset == "modelnet"
        catFile = joinpath(root, "modelnet$(num_classes)_shape_names.txt")
        categories = [line::String for line in readlines(catFile)]

        classes = Dict{String, Int}()
        for i in 1:length(categories)
            classes[categories[i]] = i
        end

        shapeids = Dict()
        shape_names = Dict()
        datapath = Dict()

        shapeids["train"] = [line::String for line in readlines(joinpath(root, "modelnet$(num_classes)_train.txt"))]
        shapeids["test"] = [line::String for line in readlines(joinpath(root, "modelnet$(num_classes)_test.txt"))]

        for split_ in ["train","test"]
            shape_names[split_] = [join(split(shapeids[split_][i], "_")[1:end-1], "_") for i in 1:length(shapeids[split_])]
            datapath[split_] = [(shape_names[split_][i], joinpath(root, shape_names[split_][i], (shapeids[split_][i])*".txt")) for i in 1:length(shapeids[split_])]
        end

        #shuffling trainset Point
        rng = MersenneTwister(1234);
        shuffle!(rng, datapath["train"])

        # Fetching the train and validation data and getting them into proper shape
        datapath, classes = get_data(num_classes)
        X = [datapoint(p, classes) for p in datapath["train"]]
        data = [X[i][1] for i in 1:length(X)]
        labels = onehotbatch(cat([X[i][2] for i in 1:length(X)]..., dims=1), 1:num_classes)
        train = [(cat(data[i]..., dims = 3), labels[:,i]) for i in partition(1:length(data), batch_size)]

        VAL = [datapoint(p, classes) for p in datapath["test"]]
        valX = cat([VAL[i][1] for i in 1:length(VAL)]..., dims=3)
        valY = onehotbatch(cat([VAL[i][2] for i in 1:length(VAL)]..., dims=1), 1:num_classes)
    else
        print("invalid dataset name choose between modelnet and 3dmnist")
        exit()
    end

    return train, (valX, valY)
end

function pcNormalize(pc)
    centroid = mean((x->x), pc, dims = 1)
    pc = pc .- centroid
    m = max((sum(pc.^2, dims=2).^.5)...)
    pc = pc ./ convert(Float32, m)
    return pc
end

function datapoint(fn, classes)
    #fn = datapath[index]
    cls = [classes[fn[1]]]
    Pointset = Array{Float32}(undef, npoints, 3)
    stream = open(fn[2], "r")
    for i in 1:npoints
        tmp = readline(stream, keep=false)
        Pointset[i, :] = map((x->parse(Float32, x)), split(tmp, " ")[1:3])
    end
    Pointset = pcNormalize(Pointset)
    return (Pointset,cls)
end

max_stars_repo_name: nirmal-suthar/DGCNN

max_issues_repo_name: nirmal-suthar/DGCNN

max_forks_repo_name: nirmal-suthar/DGCNN


text:
"""
    HyperComplex{T <: Real} = Nilplexification{DualComplex{T}}

A hyper-complex number is a nilplexification of dual-complex numbers.
"""
const HyperComplex{T <: Real} = Nilplexification{DualComplex{T}}

function HyperComplex(z::DualComplex{T}) where T <: Real
    HyperComplex{T}(z, zero(DualComplex{T}))
end

function HyperComplex(z::Cmplex{T}) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(z, zero(Cmplex{T})),
        zero(DualComplex{T})
    )
end

function HyperComplex(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        DualComplex{T}(
            Cmplex{T}(f, g),
            Cmplex{T}(h, j)
        )
    )
end

function HyperComplex(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        DualComplex{T}(
            Cmplex{T}(f, g),
            Cmplex{T}(h, zero(T))
        )
    )
end

function HyperComplex(a::T, b::T, c::T, d::T, f::T, g::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        DualComplex{T}(
            Cmplex{T}(f, g),
            zero(Cmplex{T})
        )
    )
end

function HyperComplex(a::T, b::T, c::T, d::T, f::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        DualComplex{T}(
            Cmplex{T}(f, zero(T)),
            zero(Cmplex{T})
        )
    )
end

function HyperComplex(a::T, b::T, c::T, d::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        zero(DualComplex{T})
    )
end

function HyperComplex(a::T, b::T, c::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, zero(T))
        ),
        zero(DualComplex{T})
    )
end

function HyperComplex(a::T, b::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, b),
            zero(Cmplex{T})
        ),
        zero(DualComplex{T})
    )
end

function HyperComplex(a::T) where T <: Real
    HyperComplex{T}(
        DualComplex{T}(
            Cmplex{T}(a, zero(T)),
            zero(Cmplex{T})
        ),
        zero(DualComplex{T})
    )
end

function HyperComplex(a::Real, b::Real)
    HyperComplex(promote(a, b)...)
end

function HyperComplex(a::Real, b::Real, c::Real)
    HyperComplex(promote(a, b, c)...)
end

function HyperComplex(a::Real, b::Real, c::Real, d::Real)
    HyperComplex(promote(a, b, c, d)...)
end

function HyperComplex(a::Real, b::Real, c::Real, d::Real, f::Real)
    HyperComplex(promote(a, b, c, d, f)...)
end

function HyperComplex(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    HyperComplex(promote(a, b, c, d, f, g)...)
end

function HyperComplex(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    HyperComplex(promote(a, b, c, d, f, g, h)...)
end

function HyperComplex(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    HyperComplex(promote(a, b, c, d, f, g, h, j)...)
end

"""
    unreal(z::HyperComplex)

The unreal part of a hyper-complex number is a 7-vector.
"""
function unreal(z::HyperComplex)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::HyperComplex)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", H: ")
    print(io, z.l.l.r)
    print(io, ", A: ")
    print(io, z.l.r.l)
    print(io, ", HA: ")
    print(io, z.l.r.r)
    print(io, ", B: ")
    print(io, z.r.l.l)
    print(io, ", HB: ")
    print(io, z.r.l.r)
    print(io, ", AB: ")
    print(io, z.r.r.l)
    print(io, ", HAB: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function random(::Type{HyperComplex{T}}) where T <: Real
    HyperComplex{T}(random(DualComplex{T}), random(DualComplex{T}))
end

function barH(z::HyperComplex{T}) where T <: Real
    HyperComplex{T}(barH(z.l), barH(z.r))
end

function barA(z::HyperComplex{T}) where T <: Real
    HyperComplex{T}(barA(z.l), barA(z.r))
end

function barB(z::HyperComplex{T}) where T <: Real
    HyperComplex{T}(z.l, -(z.r))
end

function quadrance(z::HyperComplex)
    (z * barB(z)).l
end

function norm(z::HyperComplex)
    abs2(quadrance(quadrance(z)))
end

function iszerodivisor(z::HyperComplex)
    iszerodivisor(quadrance(z))
end

function inv(z::HyperComplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    q = quadrance(z)
    qq = quadrance(q)
    barB(z) * (barA(q) * barH(qq)) / abs2(qq)
end


text:
@testset "EUtils" begin

    @testset "einfo" begin
        res = einfo(db="pubmed")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"
    end

    @testset "esearch" begin
        res = esearch(db="pubmed", term="""(Asthma[MeSH Major Topic]) AND
                                        ("1/1/2018"[Date - Publication] :
                                        "3000"[Date - Publication])""")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"
    end

    @testset "epost" begin
        ctx = Dict()
        res = epost(ctx, db="protein", id="NP_005537")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test haskey(ctx, :WebEnv)
        @test haskey(ctx, :query_key)
    end

    @testset "esummary" begin
        # esummary doesn't seem to support accession numbers
        res = esummary(db="protein", id="15718680,157427902,119703751")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"

        res = esummary(db="protein", id=["15718680", "157427902", "119703751"])
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"

        # esearch then esummary
        query = "asthma[mesh] AND leukotrienes[mesh] AND 2009[pdat]"
        ctx = Dict()
        res = esearch(ctx, db="pubmed", term=query, usehistory=true, retmode="xml")
        @test res.status == 200
        res = esummary(ctx, db="pubmed")
        @test res.status == 200

        ctx = Dict()
        res = esearch(ctx, db="pubmed", term=query, usehistory=true, retmode="json")
        @test res.status == 200
        res = esummary(ctx, db="pubmed")
        @test res.status == 200
    end

    @testset "efetch" begin
        res = efetch(db="nuccore", id="NM_001178.5", retmode="xml", idtype="acc")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        @test isa(parse_xml(String(res.body)), XMLDict.XMLDictElement)

        # epost then efetch
        ctx = Dict()
        res = epost(ctx, db="protein", id="NP_005537")
        @test res.status == 200
        res = efetch(ctx, db="protein", retmode="xml")
        @test res.status == 200

        # esearch then efeth for large number of ids
        retmax = 1000
        search_term = """(Asthma[MeSH Major Topic])
                        AND ("1/1/2018"[Date - Publication] :
                        "3000"[Date - Publication])"""
        res = esearch(db = "pubmed", term = search_term,
        retstart = 0, retmax = retmax, tool = "BioJulia")

        #convert xml to dictionary
        esearch_dict = parse_xml(String(res.body))

        #get the list of ids and perfom a fetch
        ids = [parse(Int64, id_node) for id_node in esearch_dict["IdList"]["Id"]]

        res = efetch(db = "pubmed", tool = "BioJulia", retmode = "xml", rettype = "null", id = ids)
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")

        body = parse_xml(String(res.body))

        @test haskey(body, "PubmedArticle")
        @test length(body["PubmedArticle"]) == retmax

    end

    @testset "elink" begin
        res = elink(dbfrom="protein", db="gene", id="NM_001178.5")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        @test isa(parse_xml(String(res.body)), XMLDict.XMLDictElement)
    end

    @testset "egquery" begin
        res = egquery(term="""(Asthma[MeSH Major Topic]) AND
                              ("1/1/2018"[Date - Publication] :
                              "3000"[Date - Publication])""")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        @test isa(parse_xml(String(res.body)), XMLDict.XMLDictElement)
    end

    @testset "espell" begin
        res = espell(db="pmc", term="fiberblast cell grwth")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        replaced_spell = body["SpelledQuery"]["Replaced"][2]
        @test replaced_spell == "growth"
    end

    @testset "ecitmatch" begin
        res = ecitmatch(
            db="pubmed",
            retmode="xml",
            bdata="proc natl acad sci u s a|1991|88|3248|mann bj|Art1|")
        @test res.status == 200
        # ECitMatch is mysterious: it returns a plain text data even though it
        # requires retmode="xml".
        # @show res.headers["Content-Type"]
    end
end


text:
using FMMTrees
using Test

N = FMMTrees.SimpleTrees.TreeNode
nodes = N[
    N(8,"a"),
        N(2,"b"),
            N(0,"c"),
            N(0,"d"),
        N(4,"e"),
            N(0,"f"),
            N(2,"g"),
                N(0,"h"),
                N(0,"i")]

tree = FMMTrees.SimpleTrees.SimpleTree(nodes)
FMMTrees.print_tree(tree)

chdit = FMMTrees.children(tree)

@test data(tree) == "a"
@test collect(data(ch) for ch in FMMTrees.children(tree)) == ["b","e"]
@test collect(data(lv) for lv in FMMTrees.leaves(tree)) == ["c","d","f","h","i"]

const N2 = FMMTrees.PointerBasedTrees.Node
nodes2 = N2[
    N2("a",2,-1,-1,2),
        N2("b",2,5,1,3),
            N2("c",0,4,2,-1),
            N2("d",0,-1,2,-1),
        N2("e",2,-1,1,6),
            N2("f",0,7,5,-1),
            N2("g",2,-1,5,8),
                N2("h",0,9,7,-1),
                N2("i",0,-1,7,-1)]
tree2 = FMMTrees.PointerBasedTrees.PointerBasedTree(nodes2,1)

@test data(tree2) == "a"
@test collect(data(tree2,ch) for ch in FMMTrees.children(tree2)) == ["b","e"]
@test collect(data(tree2,lv) for lv in FMMTrees.leaves(tree2)) == ["c","d","f","h","i"]

FMMTrees.print_tree(tree2)


text:
using Test

include("test_utils.jl")
include("test_converter.jl")
include("test_resnet.jl")
include("test_run_resnet.jl")
include("test_mnist.jl")


text:
#############################################################
## NOTE: This file is written in the standard GTK convention;
## documentation is not necessary in this particular case.
#############################################################

include("blurrFunctions.jl")
include("colorFunctions.jl")
include("transformationFunctions.jl")

using Gtk, Images, ImageView, FileIO

# GLOBAL VARIABLES AND GENERAL FUNCTIONS

bld = GtkBuilder(filename="projekt/GUILayout.glade")
saving_path = ""
save_flag = false
undo_counter = -1
transits_counter = 0
transit_given_reg = []

current_image = Array{RGB{Normed{UInt8,8}},2}
current_image_u = Array{RGB{Normed{UInt8,8}},2}
current_image_uu = Array{RGB{Normed{UInt8,8}},2}
function new_current_image(new_image, canvas)
    global current_image_uu = current_image_u
    global current_image_u = current_image
    global current_image = new_image
    global undo_counter += 1
    imshow(canvas, new_image)
end

function undo_image(canvas)
    if undo_counter > 0
        global current_image = current_image_u
        global current_image_u = current_image_uu
        global undo_counter -= 1
        imshow(canvas, current_image)
    else
        println("ERROR - cant undo")
    end
end

rgb_choice = ""
gray_choice = ""
blur_choice = ""
orig_choice = ""
range_right_val = 0
range_left_val = 0
range_up_val = 0
range_down_val = 0
selection_ru = (Nothing, Nothing)
selection_ld = (Nothing, Nothing)

# ELEMENTS BUILDING

# main window
mainW = bld["mainW"]
frame, cnv = ImageView.frame_canvas(:auto)
image_grid = bld["image_grid"]
image_grid[2,2] = frame

b_open = bld["b_open"]
b_save = bld["b_save"]
b_saveas = bld["b_saveas"]
b_quit = bld["b_quit"]
b_undo = bld["b_undo"]

b_brctr = bld["b_brctr"]
b_hsl = bld["b_hsl"]
b_rgb = bld["b_rgb"]
b_gray = bld["b_gray"]
b_negative = bld["b_negative"]
b_blur = bld["b_blur"]
b_sharp = bld["b_sharp"]
b_transit = bld["b_transit"]
b_xmirr = bld["b_xmirr"]
b_ymirr = bld["b_ymirr"]

scale_left = bld["scale_left"]
scale_right = bld["scale_right"]
scale_up = bld["scale_up"]
scale_down = bld["scale_down"]

a_range_up = bld["a_range_up"]
a_range_down = bld["a_range_down"]
a_range_left = bld["a_range_left"]
a_range_right = bld["a_range_right"]
b_range_ok = bld["b_range_ok"]
b_range_cancel = bld["b_range_cancel"]

# brightness and conterast window
brctrW = bld["brctrW"]
b_brctr_cancel = bld["b_brctr_cancel"]
a_brightness = bld["a_brightness"]
a_contrast = bld["a_contrast"]
b_brctr_ok = bld["b_brctr_ok"]
# hue satur light window
hslW = bld["hslW"]
b_hsl_cancel = bld["b_hsl_cancel"]
a_hue = bld["a_hue"]
a_satur = bld["a_satur"]
a_light = bld["a_light"]
b_hsl_ok = bld["b_hsl_ok"]
# rgb components
rgbW = bld["rgbW"]
b_rgb_cancel = bld["b_rgb_cancel"]
rgb_sr = bld["rgb_sr"]
rgb_sg = bld["rgb_sg"]
rgb_sb = bld["rgb_sb"]
rgb_hr = bld["rgb_hr"]
rgb_hg = bld["rgb_hg"]
rgb_hb = bld["rgb_hb"]
b_rgb_ok = bld["b_rgb_ok"]
# grayscale
grayW = bld["grayW"]
b_gray_cancel = bld["b_gray_cancel"]
gray_r = bld["gray_r"]
gray_g = bld["gray_g"]
gray_b = bld["gray_b"]
gray_rg = bld["gray_rg"]
gray_gb = bld["gray_gb"]
gray_br = bld["gray_br"]
gray_rgb = bld["gray_rgb"]
b_gray_ok = bld["b_gray_ok"]
# blur
blurW = bld["blurW"]
b_blur_cancel = bld["b_blur_cancel"]
a_blur_radius_s = bld["a_blur_radius_s"]
a_blur_intens_s = bld["a_blur_intens_s"]
blur_mask_aver = bld["blur_mask_aver"]
blur_mask_circ = bld["blur_mask_circ"]
blur_mask_lp3 = bld["blur_mask_lp3"]
b_blur_ok = bld["b_blur_ok"]
# sharp
sharpW = bld["sharpW"]
b_sharp_cancel = bld["b_sharp_cancel"]
a_sharp_radius_s = bld["a_sharp_radius_s"]
a_sharp_intens_s = bld["a_sharp_intens_s"]
b_sharp_ok = bld["b_sharp_ok"]
# transform
transitW = bld["transitW"]
a_transl_vect_x = bld["a_transl_vect_x"]
a_transl_vect_y = bld["a_transl_vect_y"]
a_rotat_angle = bld["a_rotat_angle"]
a_scale_ratio_x = bld["a_scale_ratio_x"]
a_scale_ratio_y = bld["a_scale_ratio_y"]
b_add_transl = bld["b_add_transl"]
b_add_rotat = bld["b_add_rotat"]
b_add_scale = bld["b_add_scale"]

b_transit_cancel = bld["b_transit_cancel"]
b_transit_ok = bld["b_transit_ok"]
list_transit = bld["list_transit"]
item_transit_1 = bld["item_transit_1"]
item_transit_2 = bld["item_transit_2"]
item_transit_3 = bld["item_transit_3"]
item_transit_4 = bld["item_transit_4"]
item_transit_5 = bld["item_transit_5"]
label_transit_1 = bld["label_transit_1"]
label_transit_2 = bld["label_transit_2"]
label_transit_3 = bld["label_transit_3"]
label_transit_4 = bld["label_transit_4"]
label_transit_5 = bld["label_transit_5"]
b_minus_transit_1 = bld["b_minus_transit_1"]
b_minus_transit_2 = bld["b_minus_transit_2"]
b_minus_transit_3 = bld["b_minus_transit_3"]
b_minus_transit_4 = bld["b_minus_transit_4"]
b_minus_transit_5 = bld["b_minus_transit_5"]
transit_list_reg = [[item_transit_1, label_transit_1, b_minus_transit_1],
                    [item_transit_2, label_transit_2, b_minus_transit_2],
                    [item_transit_3, label_transit_3, b_minus_transit_3],
                    [item_transit_4, label_transit_4, b_minus_transit_4],
                    [item_transit_5, label_transit_5, b_minus_transit_5],]
orig_c = bld["orig_c"]
orig_l = bld["orig_l"]
orig_r = bld["orig_r"]
orig_u = bld["orig_u"]
orig_d = bld["orig_d"]
orig_ru = bld["orig_ru"]
orig_rd = bld["orig_rd"]
orig_lu = bld["orig_lu"]
orig_ld = bld["orig_ld"]

transitLimitW = bld["transitLimitW"]
b_transit_limit_cancel = bld["b_transit_limit_cancel"]

# CALLBACK FUNCTIONS
function open_fileopen(w)
    path = open_dialog("Pick an image file", GtkNullContainer(), (GtkFileFilter("*.jpg", name="All supported formats"), "*.jpg"))
    if path != ""
        img = load(path)
        new_current_image(img, cnv)
    end
end
function save_filesaveas(w)
    path_raw = save_dialog("Save as...", mainW, ("*.jpg", GtkFileFilter("*.jpg", name="All supported formats")))

    if path_raw[end - 3:end] != ".jpg"
        global saving_path = path_raw * ".jpg"
    else
        global saving_path = path_raw
    end
    Images.save(saving_path, current_image)
    global save_flag = true
end
function save_filesave(w)
    if save_flag == true
        Images.save(saving_path, current_image)
    else
        save_filesaveas(w)
    end
end
function quit_app(w)
    destroy(mainW)
end
function undo_clb(w)
    undo_image(cnv)
end

function brctr_open(w)
    show(brctrW)
    set_gtk_property!(a_brightness, :value, 0)
    set_gtk_property!(a_contrast, :value, 0)
end
brctr_close(w) = hide(brctrW)
function hsl_open(w)
    show(hslW)
    set_gtk_property!(a_hue, :value, 0)
    set_gtk_property!(a_satur, :value, 0)
    set_gtk_property!(a_light, :value, 0)
end
hsl_close(w) = hide(hslW)

function rgb_open(w)
    show(rgbW)
end
rgb_close(w) = hide(rgbW)
function rgb_update(w)
    if get_gtk_property(w, :active, Bool)
        global rgb_choice = w
    end
end

function gray_open(w)
    show(grayW)
end
function gray_update(w)
    if get_gtk_property(w, :active, Bool)
        global gray_choice = w
    end
end
gray_close(w) = hide(grayW)

function blur_open(w)
    set_gtk_property!(a_blur_intens_s, :sensitive, false)
    show(blurW)
    if get_gtk_property(blur_mask_circ, :active, Bool)
        set_gtk_property!(a_blur_intens_s, :sensitive, false)
    end
end
blur_close(w) = hide(blurW)
function blur_update(w)
    if get_gtk_property(w, :active, Bool)
        global blur_choice = w
    end
    if w == blur_mask_circ || w == blur_mask_aver
        set_gtk_property!(a_blur_intens_s, :sensitive, false)
    else
        set_gtk_property!(a_blur_intens_s, :sensitive, true)
    end
end

function sharp_open(w)
    show(sharpW)
end
sharp_close(w) = hide(sharpW)

function show_selection(image, range_left, range_right, range_up, range_down)

    color_matrices = generateMatricesRGB(image)

    img_height = size(color_matrices[1])[1]
    img_width = size(color_matrices[1])[2]
    ru = (max(Int(floor(range_up * img_height / 100)), 1), max(Int(floor((100 - range_right) * img_width / 100)), 1))
    ld = (max(Int(floor((100 - range_down) * img_height / 100)), 1), max(Int(floor(range_left * img_width / 100)), 1))

    for i in 1:3
        original_matrix = copy(color_matrices[i])

        color_matrices[i] = original_matrix .* 0.3
        color_matrices[i][ru[1]:ld[1], ld[2]:ru[2]] = original_matrix[ru[1]:ld[1], ld[2]:ru[2]]
    end

    global selection_ru = ru
    global selection_ld = ld

    selection_image = matriceRGB(color_matrices...)
    imshow(cnv, selection_image)
end # RGB

function range_open(w)

    global range_left_val = get_gtk_property(a_range_left, :value, Int)
    global range_right_val = get_gtk_property(a_range_right, :value, Int)
    global range_up_val = get_gtk_property(a_range_up, :value, Int)
    global range_down_val = get_gtk_property(a_range_down, :value, Int)
    show_selection(current_image, range_left_val, range_right_val, range_up_val, range_down_val)

    show(scale_left)
    show(scale_right)
    show(scale_down)
    show(scale_up)
    show(b_range_cancel)
    show(b_range_ok)
end

function range_accept(w)
    show(transitW)
end

function range_close(w)
    hide(scale_left)
    hide(scale_right)
    hide(scale_down)
    hide(scale_up)
    hide(b_range_cancel)
    hide(b_range_ok)
    imshow(cnv, current_image)

end

function update_range_left(w)
    global range_left_val = get_gtk_property(w, :value, Int)
    if range_left_val + range_right_val >= 99
        set_gtk_property!(w, :value, 99 - range_right_val)
    end
    show_selection(current_image, range_left_val, range_right_val, range_up_val, range_down_val)
end
function update_range_right(w)
    global range_right_val = get_gtk_property(w, :value, Int)
    if range_right_val + range_left_val >= 99
        set_gtk_property!(w, :value, 99 - range_left_val)
    end
    show_selection(current_image, range_left_val, range_right_val, range_up_val, range_down_val)
end
function update_range_up(w)
    global range_up_val = get_gtk_property(w, :value, Int)
    if range_up_val + range_down_val >= 99
        set_gtk_property!(w, :value, 99 - range_down_val)
    end
    show_selection(current_image, range_left_val, range_right_val, range_up_val, range_down_val)
end
function update_range_down(w)
    global range_down_val = get_gtk_property(w, :value, Int)
    if range_down_val + range_up_val >= 99
        set_gtk_property!(w, :value, 99 - range_up_val)
    end
    show_selection(current_image, range_left_val, range_right_val, range_up_val, range_down_val)
end

function transit_open(w)
    global transits_counter = 0
    for i in 1:5
        hide(transit_list_reg[i][1])
        global transit_given_reg = []
    end
    show(transitW)
end
transit_close(w) = hide(transitW)

function draw_transit_list()
    for i in 1:length(transit_given_reg)
        if transit_given_reg[i][1] == 't'
            set_gtk_property!(transit_list_reg[i][2], :label,
                string("Translation: (", string(transit_given_reg[i][2]), ", ", string(transit_given_reg[i][3]), ")"))
        elseif transit_given_reg[i][1] == 's'
            set_gtk_property!(transit_list_reg[i][2], :label,
             string("Scaling: (", string(transit_given_reg[i][2]), ", ", string(transit_given_reg[i][3]), ")"))
        elseif transit_given_reg[i][1] == 'r'
            set_gtk_property!(transit_list_reg[i][2], :label,
             string("Rotation: ", string(transit_given_reg[i][2]), "°"))
        end

        show(transit_list_reg[i][1])
    end
    for i in length(transit_given_reg) + 1:5
        hide(transit_list_reg[i][1])
    end
end

function transl_add(w)
    if transits_counter >= 5
        show(transitLimitW)
    else
        vect_x = get_gtk_property(a_transl_vect_x, :value,  Float32)
        vect_y = get_gtk_property(a_transl_vect_y, :value, Float32)
        if !(vect_x == 0.0 && vect_y == 0.0)
            push!(transit_given_reg, ('t', vect_x, vect_y))
            global transits_counter += 1
        end
    end
    draw_transit_list()
end

function rotat_add(w)
    if transits_counter >= 5
        show(transitLimitW)
    else
        angle = get_gtk_property(a_rotat_angle, :value, Float32)
        if angle != 0
            push!(transit_given_reg, ('r', angle))
            global transits_counter += 1
        end
    end
    draw_transit_list()
end

function scale_add(w)
    if transits_counter >= 5
        show(transitLimitW)
    else
        ratio_x = get_gtk_property(a_scale_ratio_x, :value, Float32)
        ratio_y = get_gtk_property(a_scale_ratio_y, :value,  Float32)
        if !(ratio_x == 1.0 && ratio_y == 1.0)
            push!(transit_given_reg, ('s', ratio_x, ratio_y))
            global transits_counter += 1
        end
    end
    draw_transit_list()
end

function transit_del_elem_1(w)
    global transits_counter -= 1
    hide(transit_list_reg[1][1])
    splice!(transit_given_reg, 1)
    draw_transit_list()
end
function transit_del_elem_2(w)
    global transits_counter -= 1
    hide(transit_list_reg[2][1])
    splice!(transit_given_reg, 2)
    draw_transit_list()
end
function transit_del_elem_3(w)
    global transits_counter -= 1
    hide(transit_list_reg[3][1])
    splice!(transit_given_reg, 3)
    draw_transit_list()
end
function transit_del_elem_4(w)
    global transits_counter -= 1
    hide(transit_list_reg[4][1])
    splice!(transit_given_reg, 4)
    draw_transit_list()
end
function transit_del_elem_5(w)
    global transits_counter -= 1
    hide(transit_list_reg[5][1])
    splice!(transit_given_reg, 5)
    draw_transit_list()
end
function orig_update(w)
    if get_gtk_property(w, :active, Bool)
        global orig_choice = w
    end
end

function transit_limit_dialog_open(w)
    show(transitLimitW)
end
transit_limit_dialog_close(w) = hide(transitLimitW)

# BACKEND FUNCTIONS CALL
function call_brctr(w)
    rgb = generateMatricesRGB(current_image)
    brightness_fact = get_gtk_property(a_brightness, :value, Float64)
    contrast_fact = get_gtk_property(a_contrast, :value, Float64)

    if brightness_fact != 0
        rgb = changeBrightness(rgb, brightness_fact)
    end
    if contrast_fact != 0
        rgb = changeContrast(rgb, contrast_fact)
    end

    new_current_image(matriceRGB(rgb...), cnv)
    hide(brctrW)
end

function call_hsl(w)
    rgb = generateMatricesRGB(current_image)
    hue_fact = get_gtk_property(a_hue, :value, Float64)
    satur_fact = get_gtk_property(a_satur, :value, Float64)
    light_fact =  get_gtk_property(a_light, :value, Float64)
    if hue_fact < 0
        hue_fact += 360
    end
    if hue_fact != 0
        rgb = changeColors(rgb, hue_fact)
    end
    if satur_fact != 0
        rgb = changeSaturation(rgb, satur_fact)
    end
    if light_fact != 0
        rgb = changeLightness(rgb, light_fact)
    end

    new_current_image(matriceRGB(rgb...), cnv)
    hide(hslW)
end

function call_rgb(w)
    rgb = generateMatricesRGB(current_image)
    if rgb_choice == rgb_sr
        rgb = onlyRed(rgb)
    elseif rgb_choice == rgb_sg
        rgb = onlyGreen(rgb)
    elseif rgb_choice == rgb_sb
        rgb = onlyBlue(rgb)
    elseif rgb_choice == rgb_hr
        rgb = withoutRed(rgb)
    elseif rgb_choice == rgb_hg
        rgb = withoutGreen(rgb)
    elseif rgb_choice == rgb_hb
        rgb = withoutBlue(rgb)
    end
    new_current_image(matriceRGB(rgb...), cnv)
    hide(rgbW)
end

function call_gray(w)
    rgb = generateMatricesRGB(current_image)
    if gray_choice == gray_r
        rgb = redAsAGrayscale(rgb)
    elseif gray_choice == gray_g
        rgb = greenAsAGrayscale(rgb)
    elseif gray_choice == gray_b
        rgb = blueAsAGrayscale(rgb)
    elseif gray_choice == gray_rg
        rgb = onlyBlueAndGrayscale(rgb)
    elseif gray_choice == gray_gb
        rgb = onlyRedAndGrayscale(rgb)
    elseif gray_choice == gray_br
        rgb = withoutGreen(rgb)
    elseif gray_choice == gray_rgb
        rgb = grayscaleLuminosity(rgb)
    end
    new_current_image(matriceRGB(rgb...), cnv)
    hide(grayW)
end

function call_negative(w)
    rgb = generateMatricesRGB(current_image)
    rgb = negative(rgb)
    new_current_image(matriceRGB(rgb...), cnv)
end

function call_blur(w)
    rgb = generateMatricesRGB(current_image)
    blur_radius = get_gtk_property(a_blur_radius_s, :value, Int)
    blur_intens = get_gtk_property(a_blur_intens_s, :value, Int)
    if blur_radius % 2 == 0
        blur_radius += 1
    end
    if blur_choice == blur_mask_aver
        rgb = converting(rgb, maskAverage(blur_radius))
    elseif blur_choice == blur_mask_circ
        rgb = converting(rgb, circle(blur_radius))
    elseif blur_choice == blur_mask_lp3
        rgb = converting(rgb, LP3(blur_radius, blur_intens))
    end
    new_current_image(matriceRGB(rgb...), cnv)
    hide(blurW)
end

function call_sharp(w)
    rgb = generateMatricesRGB(current_image)
    sharp_radius = get_gtk_property(a_sharp_radius_s, :value, Int)
    sharp_intens = get_gtk_property(a_sharp_intens_s, :value, Int)
    if sharp_radius % 2 == 0
        sharp_radius += 1
    end
    rgb = converting(rgb, meanRemoval(sharp_radius, sharp_intens))

    new_current_image(matriceRGB(rgb...), cnv)
    hide(sharpW)
end

function call_affin(w)
    rgb = generateMatricesRGB(current_image)
    if length(transit_given_reg) > 0
        origin = (0, 0)
        s_r_list = []
        df_origin = defaultOrigin(selection_ld, selection_ru)
        if orig_choice == orig_c
            origin = df_origin
        elseif orig_choice == orig_l
            origin = (df_origin[1], selection_ld[2])
        elseif orig_choice == orig_r
            origin = (df_origin[1], selection_ru[2])
        elseif orig_choice == orig_u
            origin = (selection_ru[1], df_origin[2])
        elseif orig_choice == orig_d
            origin = (selection_ld[1], df_origin[2])
        elseif orig_choice == orig_ru
            origin = selection_ru
        elseif orig_choice == orig_ld
            origin = selection_ld
        elseif orig_choice == orig_rd
            origin = (selection_ld[1], selection_ru[2])
        elseif orig_choice == orig_lu
            origin = (selection_ru[1], selection_ld[2])
        end
        transl_x = 0
        transl_y = 0
        for i in 1:length(transit_given_reg)
            if transit_given_reg[i][1] == 't'
                transl_x += transit_given_reg[i][2]
                transl_y -= transit_given_reg[i][3]
            elseif transit_given_reg[i][1] == 's'
                push!(s_r_list, ('s', (transit_given_reg[i][3],  transit_given_reg[i][2])))
            elseif transit_given_reg[i][1] == 'r'
                push!(s_r_list, ('r', transit_given_reg[i][2]))
            end
        end
        transl_vect = (transl_y, transl_x)
        rgb = selectionTransform(rgb, selection_ld, selection_ru, [(origin, s_r_list, transl_vect)])
        new_current_image(matriceRGB(rgb...), cnv)
    else
        new_current_image(current_image, cnv)
    end
    hide(transitW)

end

function call_xmirr(w)
    rgb = generateMatricesRGB(current_image)
    rgb = mirror(rgb, :x)
    new_current_image(matriceRGB(rgb...), cnv)
end

function call_ymirr(w)
    rgb = generateMatricesRGB(current_image)
    rgb = mirror(rgb, :y)
    new_current_image(matriceRGB(rgb...), cnv)
end

# SIGNAL CONNECTING
signal_connect(open_fileopen, b_open, "activate")
signal_connect(save_filesave, b_save, "activate")
signal_connect(save_filesaveas, b_saveas, "activate")
signal_connect(quit_app, b_quit, "activate")
signal_connect(undo_clb, b_undo, "activate")

signal_connect(brctr_open, b_brctr, "clicked")
signal_connect(brctr_close, b_brctr_cancel, "clicked")
signal_connect(call_brctr, b_brctr_ok, "clicked")

signal_connect(hsl_open, b_hsl, "clicked")
signal_connect(hsl_close, b_hsl_cancel, "clicked")
signal_connect(call_hsl, b_hsl_ok, "clicked")

signal_connect(rgb_open, b_rgb, "clicked")
signal_connect(rgb_close, b_rgb_cancel, "clicked")
signal_connect(rgb_update, rgb_sr, "toggled")
signal_connect(rgb_update, rgb_sg, "toggled")
signal_connect(rgb_update, rgb_sb, "toggled")
signal_connect(rgb_update, rgb_hr, "toggled")
signal_connect(rgb_update, rgb_hg, "toggled")
signal_connect(rgb_update, rgb_hb, "toggled")
signal_connect(call_rgb, b_rgb_ok, "clicked")

signal_connect(gray_open, b_gray, "clicked")
signal_connect(gray_close, b_gray_cancel, "clicked")
signal_connect(gray_update, gray_r, "toggled")
signal_connect(gray_update, gray_g, "toggled")
signal_connect(gray_update, gray_b, "toggled")
signal_connect(gray_update, gray_rg, "toggled")
signal_connect(gray_update, gray_gb, "toggled")
signal_connect(gray_update, gray_br, "toggled")
signal_connect(gray_update, gray_rgb, "toggled")
signal_connect(call_gray, b_gray_ok, "clicked")

signal_connect(call_negative, b_negative, "clicked")

signal_connect(blur_open, b_blur, "clicked")
signal_connect(blur_close, b_blur_cancel, "clicked")
signal_connect(blur_update, blur_mask_aver, "toggled")
signal_connect(blur_update, blur_mask_circ, "toggled")
signal_connect(blur_update, blur_mask_lp3, "toggled")
signal_connect(call_blur, b_blur_ok, "clicked")

signal_connect(sharp_open, b_sharp, "clicked")
signal_connect(sharp_close, b_sharp_cancel, "clicked")
signal_connect(call_sharp, b_sharp_ok, "clicked")

signal_connect(range_open, b_transit, "clicked")
signal_connect(range_accept, b_range_ok, "clicked")
signal_connect(range_close, b_range_cancel, "clicked")

signal_connect(update_range_left, a_range_left, "value-changed")
signal_connect(update_range_right, a_range_right, "value-changed")
signal_connect(update_range_up, a_range_up, "value-changed")
signal_connect(update_range_down, a_range_down, "value-changed")
signal_connect(transit_open, b_range_ok, "clicked")
signal_connect(transit_close, b_transit_cancel, "clicked")

signal_connect(transl_add, b_add_transl, "clicked")
signal_connect(rotat_add, b_add_rotat, "clicked")
signal_connect(scale_add, b_add_scale, "clicked")
signal_connect(orig_update, orig_c, "toggled")
signal_connect(orig_update, orig_l, "toggled")
signal_connect(orig_update, orig_r, "toggled")
signal_connect(orig_update, orig_d, "toggled")
signal_connect(orig_update, orig_u, "toggled")
signal_connect(orig_update, orig_ld, "toggled")
signal_connect(orig_update, orig_ru, "toggled")
signal_connect(orig_update, orig_rd, "toggled")
signal_connect(orig_update, orig_lu, "toggled")
signal_connect(transit_del_elem_1, transit_list_reg[1][3], "clicked")
signal_connect(transit_del_elem_2, transit_list_reg[2][3], "clicked")
signal_connect(transit_del_elem_3, transit_list_reg[3][3], "clicked")
signal_connect(transit_del_elem_4, transit_list_reg[4][3], "clicked")
signal_connect(transit_del_elem_5, transit_list_reg[5][3], "clicked")
signal_connect(call_affin, b_transit_ok, "clicked")

signal_connect(transit_limit_dialog_close, b_transit_limit_cancel, "clicked")

signal_connect(call_xmirr, b_xmirr, "clicked")
signal_connect(call_ymirr, b_ymirr, "clicked")

# INITIAL OPERATIONS EXECUTE

global rgb_choice = rgb_sr
global gray_choice = gray_rgb
global blur_choice = blur_mask_aver
global orig_choice = orig_c
showall(mainW)
hide(scale_left)
hide(scale_right)
hide(scale_down)
hide(scale_up)
hide(b_range_cancel)
hide(b_range_ok)
open_fileopen(mainW)

max_stars_repo_name: o-Mateo-o/ImageEditorApp

max_issues_repo_name: o-Mateo-o/ImageEditorApp

max_forks_repo_name: o-Mateo-o/ImageEditorApp