text:
using JFVM
using PyPlot
include("CRT.jl")

#Calculating error as a function of timesteps

# define the domain and constants
Lx = 1.0 # [m]
ϕ  = 0.4 # porosity
u  = 1.0/(3600*24) # m/s, eq. 1 m/day
D  = [1e-9, 1e-9, 1e-9] # m^2/s diffusivity
K_eq = 10.0
k_f = 100.0
k_b = k_f/K_eq
nu_A = 1.0 # stoichiometric coefficient of component A
nu_B = 1.0 # stoichiometric coefficient of component B
nu_C = 1.0 # stoichiometric coefficient of component C
a = 1.0 #order of reaction
b = 1.0 #order of reaction
c = 1.0 #order of reaction
N_comp = 3 # number of components

# Defining time-stepping parameters
tol = 1e-10 # convergence tolerance

# boundary conditions
C_bc = [1.0, 0.5, 0.0] # boundary concentration
C_0 = zeros(N_comp) # initial concentration

NX = [11, 51, 101, 201, 501, 1001]
err_array = zeros(5, length(NX))

#Calculating reference "true solution"
Nx_ref = 10001;
m_ref = createMesh1D(Nx_ref,Lx)
dt_ref = Lx/(u/ϕ)/Nx_ref
T_ref = Nx_ref/2*dt_ref
C_ref = CRT.CRT_GIA_eps(m_ref, Nx_ref, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt_ref, T_ref, tol)

for i=1:length(NX)
  Nx = NX[i]  # number of grids
  m = createMesh1D(Nx, Lx)  # mesh

  #Time-stepping parameters
  dt = Lx/(u/ϕ)/Nx
  T = Nx/2*dt

  #Calculating alternatives
  C_SNIA        = CRT.CRT_SNIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T)
  C_SIA         = CRT.CRT_SIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  C_GIA_eps     = CRT.CRT_GIA_eps(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  C_GIA_kin_eq  = CRT.CRT_GIA_kin_eq(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, a, b, c, C_bc, C_0, N_comp, dt, T, tol)
  C_GIA_kin_fb  = CRT.CRT_GIA_kin_fb(m, Nx, ϕ, u, D, k_f, k_b, nu_A, nu_B, nu_C, a, b, c, C_bc, C_0, N_comp, dt, T, tol)

  for j = 1:N_comp
    reference = internalCells(C_ref[j])
    err_array[1, i] += sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]-internalCells(C_SNIA[j])))/sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]))
    err_array[2, i] += sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]-internalCells(C_SIA[j])))/sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]))
    err_array[3, i] += sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]-internalCells(C_GIA_eps[j])))/sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]))
    err_array[4, i] += sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]-internalCells(C_GIA_kin_eq[j])))/sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]))
    err_array[5, i] += sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]-internalCells(C_GIA_kin_fb[j])))/sum(abs(reference[Int(1):Int((Nx_ref-1)/(Nx-1)):Int(Nx_ref)]))
  end
  println(i)
end

figure()
plot(NX, err_array[1,:])
plot(NX, err_array[2,:])
plot(NX, err_array[3,:])
plot(NX, err_array[4,:])
plot(NX, err_array[5,:])

grid("on")
legend(["SNIA", "SIA", "GIA ("L"$\epsilon$"")", "GIA ("L"$K_{eq}$"")", "GIA ("L"$k_fk_b$"")"])
xlabel(L"$N_x$")
ylabel(L"||$C_{ref}-C||_1/||C_{ref}||_1$")
savefig("error_grid_plot.png")

figure()
semilogy(NX, err_array[1,:])
semilogy(NX, err_array[2,:])
semilogy(NX, err_array[3,:])
semilogy(NX, err_array[4,:])
semilogy(NX, err_array[5,:])

grid("on")
legend(["SNIA", "SIA", "GIA ("L"$\epsilon$"")", "GIA ("L"$K_{eq}$"")", "GIA ("L"$k_fk_b$"")"])
xlabel(L"$N_x$")
ylabel(L"log(||$C_{ref}-C||_1/||C_{ref}||_1)$")
savefig("error_grid_semilogy.png")

#CRT_SIA title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (SIA)"

#CRT_SNIA title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (SNIA)"

#CRT_GIA_eps title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (FIM)"

#CRT_GIA_kin_eq title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}(K_{eq}C_A^aC_B^b-C_C^c)=0,\quad K_{eq}=$""$K_eq"" (FIM)"

#CRT_GIA_kin_fb title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}(k_{f}C_A^aC_B^b-k_{b}C_C^c)=0,\quad K_{eq}=$""$K_eq"" (FIM)"

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
# parsing state codes
@enum(ParsingStateCode
    ,es_dead=1
    ,es_start_req_or_res
    ,es_res_or_resp_H
    ,es_res_first_http_major
    ,es_res_http_major
    ,es_res_first_http_minor
    ,es_res_http_minor
    ,es_res_first_status_code
    ,es_res_status_code
    ,es_res_status_start
    ,es_res_status
    ,es_res_line_almost_done
    ,es_start_req
    ,es_req_method
    ,es_req_spaces_before_target
    ,es_req_target
    ,es_req_target_wildcard
    ,es_req_schema
    ,es_req_schema_slash
    ,es_req_schema_slash_slash
    ,es_req_server_start
    ,es_req_server
    ,es_req_server_with_at
    ,es_req_path
    ,es_req_query_string_start
    ,es_req_query_string
    ,es_req_fragment_start
    ,es_req_fragment
    ,es_req_http_start
    ,es_req_http_H
    ,es_req_http_HT
    ,es_req_http_HTT
    ,es_req_http_HTTP
    ,es_req_first_http_major
    ,es_req_http_major
    ,es_req_first_http_minor
    ,es_req_http_minor
    ,es_req_line_almost_done
    ,es_trailer_start
    ,es_header_field_start
    ,es_header_field
    ,es_header_value_discard_ws
    ,es_header_value_discard_ws_almost_done
    ,es_header_value_discard_lws
    ,es_header_value_start
    ,es_header_value
    ,es_header_value_lws
    ,es_header_almost_done
    ,es_headers_almost_done
    ,es_headers_done
    ,es_body_start
    ,es_chunk_size_start
    ,es_chunk_size
    ,es_chunk_parameters
    ,es_chunk_size_almost_done
    ,es_chunk_data
    ,es_chunk_data_almost_done
    ,es_chunk_data_done
    ,es_body_identity
    ,es_body_identity_eof
    ,es_message_done
)
for i in instances(ParsingStateCode)
    @eval const $(Symbol(string(i)[2:end])) = UInt8($i)
end

const CR = '\r'
const bCR = UInt8('\r')
const LF = '\n'
const bLF = UInt8('\n')
const CRLF = "\r\n"

ext: jl


text:
__precompile__(true)
module ePPR

include("algorithm.jl")

# export all symbols
for n in names(@__MODULE__, all=true)
    if Base.isidentifier(n) && n ∉ (nameof(@__MODULE__), :eval, :include)
        @eval export $n
    end
end

end


text:
using NLS_Fit
using Documenter

DocMeta.setdocmeta!(NLS_Fit, :DocTestSetup, :(using NLS_Fit); recursive=true)

makedocs(;
    modules=[NLS_Fit],
    authors="Vincent Picaud",
    repo="https://github.com/vincent-picaud/NLS_Fit.jl/blob/{commit}{path}#{line}",
    sitename="NLS_Fit.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://vincent-picaud.github.io/NLS_Fit.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting started" => [
            "simple_fit.md",
            "recalibration.md",
            "mapped_parameters.md",
        ],
        "model2fit_bestiary.md",
        #            "getting_started.md",
        "API.md"
    ],
)

deploydocs(;
    repo="github.com/vincent-picaud/NLS_Fit.jl.git",
    devbranch="main",
)


text:
# This file is a part of SimilaritySearch.jl

using SimilaritySearch
using Test

function test_pivs(perf, pivs)
    @test probe(perf, pivs).macrorecall >= 0.99
end

@testset "indexing vectors with PivotedSearch" begin
    db = [rand(Float32, 4) for i in 1:1000]
    queries = [rand(Float32, 4) for i in 1:30]
    ksearch = 10
    dist = L2Distance()
    seq = ExhaustiveSearch(dist, db)
    perf = Performance(seq, queries, ksearch)

    # random pivots
    test_pivs(perf, PivotedSearch(dist, db, 8))

    # pivots selected with SSS criterion
    test_pivs(perf, sss(dist, db))

    # pivots selected with distant tournament criterion
    test_pivs(perf, distant_tournament(dist, db, 8))

    # Kvp
    test_pivs(perf, Kvp(dist, db))
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

module KrigingEstimators

using GeoStatsBase
using Variography

using LinearAlgebra: Factorization, lu, cholesky, issuccess, ⋅
using Distances: Euclidean
using Distributions: Normal
using Combinatorics: multiexponents
using StaticArrays: MVector

import GeoStatsBase: fit, predict, status
import GeoStatsBase: preprocess, solve, solve_single

include("estimators.jl")

include("solvers/kriging.jl")
include("solvers/sgsim.jl")

export
  # estimators
  KrigingEstimator,
  SimpleKriging,
  OrdinaryKriging,
  UniversalKriging,
  ExternalDriftKriging,
  fit, predict, status, weights,

  # solvers
  Kriging, KrigingParam,
  SeqGaussSim, SeqGaussSimParam

end


text:
using Oceananigans
using Oceananigans.Units
using GLMakie
using Printf

grid = RectilinearGrid(size=(64, 64), extent=(64, 64), halo=(3, 3), topology=(Periodic, Flat, Bounded))
buoyancy_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(1e-8))

model = NonhydrostaticModel(grid = grid,
                            advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            closure = ScalarDiffusivity(ν=1e-4, κ=1e-4),
                            coriolis = FPlane(f=1e-4),
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            boundary_conditions = (; b=buoyancy_bcs))

bᵢ(x, y, z) = 1e-5 * z
ϵ(x, y, z) = 1e-6 * randn() # noise
set!(model, b=bᵢ, u=ϵ, v=ϵ, w=ϵ)

simulation = Simulation(model, Δt=2minutes, stop_time=24hours)

wizard = TimeStepWizard(cfl=0.8, max_change=1.1, max_Δt=2minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

progress(sim) = @printf("Iter: %d, time: %s, Δt: %s, max|w|: %s\n",
                        iteration(sim), prettytime(sim), prettytime(sim.Δt), maximum(abs, model.velocities.w))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

# Turbulent kinetic energy
u, v, w = model.velocities
e = Field((u^2 + v^2 + w^2) / 2)
ξ = Field(∂x(w) + ∂z(u))

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; e, w, ξ, b = model.tracers.b),
                                                      schedule = TimeInterval(1minute),
                                                      prefix = "ocean_convection",
                                                      field_slicer = nothing,
                                                      force = true)

run!(simulation)

# Visualization

data_path = simulation.output_writers[:fields].filepath
e = FieldTimeSeries(data_path, "e")
x, y, z = nodes(e)
Nt = length(e.times)

fig = Figure()
ax = Axis(fig[1, 1])
slider = Slider(fig[2, 1], range=1:Nt, startvalue=1)
n = slider.value

eⁿ = @lift dropdims(interior(e[$n]), dims=2)
hm = heatmap!(ax, x, z, eⁿ; colormap=:solar, colorrange=(0, maximum(e) / 2))
Colorbar(fig[1, 2], hm)

title = @lift "Turbulent kinetic energy at t = " * prettytime(e.times[$n])
Label(fig[0, :], title)

display(fig)


text:
module Server # server

using Sockets

using ..Common: Error, SSConfig, Cipher, SSConnection, close, gensubkey, read, write, ioCopy, init_read, init_write, eof

@inline function gethost(buff::Array{UInt8})
    host, port = nothing, nothing
    if buff[1] == 0x01
        host = IPv4(ntoh(unsafe_load(Ptr{UInt32}(pointer(buff) + 1))))
        port = UInt16(buff[6]) * 256 + buff[7]
    elseif buff[1] == 0x03
        len = buff[2]
        host = String(unsafe_wrap(Array{UInt8}, pointer(buff) + 2, len))
        port = UInt16(buff[len+3]) * 256 + buff[len+4]
    elseif buff[1] == 0x04
        host = IPv6(ntoh(unsafe_load(Ptr{UInt128}(pointer(buff) + 1))))
        port = UInt16(buff[18]) * 256 + buff[19]
    end

    return host, port, nothing
end

@inline function connectRemote(buff::Array{UInt8})
    if !(buff[1] in [0x01; 0x03; 0x04])
        return nothing, "Not a valid CMD"
    end

    host, port, err = gethost(buff)
    if err != nothing
        return nothing, err
    end

    client = try
        connect(host, port)
    catch err
        return nothing, err
    end

    return client, nothing
end

@inline function handleConnection(ssConn::SSConnection) # Server
    remote = nothing

    while true
        buff = Array{UInt8}(undef, 272)

        if init_read(ssConn) != nothing
            break
        end

        ~, err = read(ssConn, buff)
        if err != nothing
            break
        end

        remote, err = connectRemote(buff)
        if err != nothing
            break
        end

        @async ioCopy(ssConn, remote) # First

        if init_write(ssConn) != nothing
            break
        end

        buff = nothing
        ioCopy(remote, ssConn)

        close(ssConn)
        while !eof(ssConn)
            sleep(1)
        end

        break
    end

    remote != nothing && close(remote)
    close(ssConn)
end

@inline function tcpServer(config::SSConfig, cipher::Cipher, terminate::Condition)
    server = try
        listen(config.host, config.port)
    catch err
        return
    end

    @async while isopen(server)
        conn = try
            accept(server)
        catch err
            break
        end

        @async handleConnection(
            SSConnection(conn, cipher, zeros(UInt8, cipher.ivlen), zeros(UInt8, cipher.ivlen), Array{UInt8}(undef, cipher.taglen + 2), nothing, nothing)
        )
    end

    wait(terminate)
    close(server)
end

@inline function runServer(config::SSConfig, terminate::Condition)
    cipher = Cipher(config)

    # @async udpServer()
    tcpServer(config, cipher, terminate)
end
runServer(config::SSConfig) = runServer(config, Condition())

end # module


text:
using RecipesBase
export SeqPlot
export TestPlot
export rectangle_corners

function rectangle_corners(x::Real,y::Real,w,h; anchor=:bottomright)
    if anchor == :botttomright
        [x,x+w,x+w,x], [y,y,y+h,y+h]
    elseif anchor == :center
        [x-w/2, x+w/2, x+w/2, x-w/2], [y-h/2, y-h/2, y+h/2, y+h/2]
    else
        error("Anchor not recognised.")
    end
end

function rectangle_corners(
    x_vec::Vector{T},y_vec::Vector{T},
    w,h; anchor=:bottomleft
    ) where {T<:Real}

    x_out = Vector{Real}[]
    y_out = Vector{Real}[]
    if anchor == :bottomleft
        for (x,y) in zip(x_vec, y_vec)
            push!(x_out, [x,x+w,x+w,x])
            push!(y_out, [y,y,y+h,y+h])
        end
    elseif anchor == :center
        for (x,y) in zip(x_vec, y_vec)
            push!(x_out, [x-w/2, x+w/2, x+w/2, x-w/2])
            push!(y_out, [y-h/2, y-h/2, y+h/2, y+h/2])
        end
    else
        error("Anchor not recognised.")
    end
    x_out, y_out
end

@userplot SeqPlot
@recipe function f(
    h::SeqPlot;
    entrymargin=0.1,
    entryfontsize=20
    )
    seq = h.args[1]
    w,h = (1,1)
    y = fill(0.0, length(seq))
    x = [0.0 + i*(w+entrymargin) for i in 0:(length(seq)-1)]
    annotations := [(x,y,string(i)) for (x,y,i) in zip(x,y,seq)]
    annotationfontsize := entryfontsize
    x_cords, y_cords = rectangle_corners(x,y, w, h; anchor=:center)
    showaxis --> false
    axis --> nothing
    aspect_ratio --> 1
    legend --> false
    @series begin
        seriestype := :shape
        x_cords,y_cords
    end
end

@userplot TestPlot
@recipe function f(
    h::TestPlot
    )
    aspect_ratio --> 1
    showaxis --> false
    @series begin
        seriestype := :shape
        1:4,1:4
    end
    @series begin
        seriestype := :line
        1:4, 1:4
        annotations := [[];[(i,i,string(i)) for i in 1:4]]
        annotationhalign --> :right
    end
end


text:
module SimplerTicTacToe
export GameEnv, GameSpec
include("game.jl")
end


text:
const Float64_2 = Type{Tuple{Float64, Float64}}

# +, - relation
choose(θ, ::typeof(+), ::Type{NTuple{2, <:Real}}, ::Type{A}, ::Type{ZB}, z, b) = (z - b,)
choose(θ, ::typeof(+), ::Type{<:NTuple{2, <:Real}}, ::Type{AB}, ::Type{Z}, z) =
  let θ_ = ℝ(θ) ; (z - θ_, θ_) end
choose(θ, ::typeof(+), ::Type{Tuple{Int, Int}}, ::Type{AB}, ::Type{Z}, z) =
  let θ_ = ℝ(θ) ; (z - θ_, θ_) end

# *, / relation -- z = a * b
choose(θ, ::typeof(*), ::Float64_2, ::Type{B}, ::Type{ZA}, z, a) =
  (z / a,)

function choose(θ, ::typeof(*), ::Float64_2, ::Type{AB}, ::Type{Z}, z)
  b = 𝔹(θ)
  v = ℝ(θ)
  b ? (z/v, v) : (v, z/v)
end

choose(θ, ::typeof(*), ::Float64_2, ::Type{A}, ::Type{ZB}, z, b) =
  (a = z / b,)

choose(θ, ::typeof(/), ::Float64_2, ::Type{AB}, ::Type{Z}, z) =
  let r = ℝ(θ)
    (z * r, r)
  end

choose(θ, loc::Loc, args...) = choose(project(θ, loc), args...)
### TODO:
# Handle multiple arguments like :+(%2, %3, %4, %5) (all with potentially different types)
# Handle inverse integer multiplication and division: factoring problems
# Define more primitives as needed


text:
function initRxUHD(pointerUSRP)
	# ----------------------------------------------------
	# --- Rx Streamer
	# ----------------------------------------------------
	# --- Create a pointer related to the Rx streamer
	addressStream = Ref{uhd_rx_streamer_handle}();
	# --- Cal the init
    @assert_uhd uhd_rx_streamer_make(addressStream)
	streamerPointer = addressStream[];
	# ----------------------------------------------------
	# --- Rx Metadata
	# ----------------------------------------------------
	# --- Create a pointer related to Metadata
	addressMD = Ref{uhd_rx_metadata_handle}();
	# --- Cal the init
    @assert_uhd uhd_rx_metadata_make(addressMD)
	# --- Get the usable object
	metadataPointer = addressMD[];
# --- Pointer for counting number of samples
	pointerSamples = Ref{Csize_t}(0);
	# ----------------------------------------------------
	# --- Create the USRP wrapper object
	# ----------------------------------------------------
	uhd  = UHDRxWrapper(true,pointerUSRP,streamerPointer,metadataPointer,addressStream,addressMD,pointerSamples);
	return uhd;
end

function openUHDRx(pointerUSRP,carrierFreq,samplingRate,gain;channels=[0],antennas=["RX2"],args="",nbAntennaRx=1,cpu_format="fc32",otw_format="sc16",subdev="",bypassStreamer=false);
    # ----------------------------------------------------
    # --- Parameter checks
    # ----------------------------------------------------
    # --- MIMO mode shortcut
    (nbAntennaRx == 2) ? MIMO_MODE = true  : MIMO_MODE = false
    @assert (nbAntennaRx) ≤ length(antennas) "Mismatch in antenna configuration. Number of Rx antennas should match the string vector of antenna"
    @assert (nbAntennaRx) ≤ length(channels) "Number of Rx antennas should be lower than number of channels"
    # ----------------------------------------------------
    # --- Core structure definitions
    # ----------------------------------------------------
    # --- Init Rx stage
    uhd	  = initRxUHD(pointerUSRP);
    # --- Create the structure to be updated
    rx = UHDRx(uhd,carrierFreq,samplingRate,gain,antennas,channels,0,0,nbAntennaRx);
    # ----------------------------------------------------
    # --- Radio configuration
    # ----------------------------------------------------
    # --- RF Configuration
    for (c,currChan) in enumerate(channels[1:nbAntennaRx])
        # Set the carrier frequencies
        updateCarrierFreq!(rx,carrierFreq,currChan)
        # Set the sampling rate
        updateSamplingRate!(rx,samplingRate,currChan)
        # Set the gains
        updateGain!(rx,gain,currChan)
        # Set the antennas
        uhd_usrp_set_rx_antenna(uhd.pointerUSRP,antennas[c],currChan)
    end
    # --- Custom subdev use
    if !isempty(subdev)
        # --- Custom Subdev use
        pointerSubDev = Ref{uhd_subdev_spec_handle}()
        @assert_uhd uhd_subdev_spec_make(pointerSubDev,subdev)
        @assert_uhd uhd_usrp_set_rx_subdev_spec(rx.uhd.pointerUSRP,pointerSubDev[],0)
        uhd_subdev_spec_free(pointerSubDev)
    else
        # Default configuration for subdev. It should leads to issues with multiple antenna (at least in e310)
        if nbAntennaRx > 1
            @warn "Multiple channels (antennas)  with default subdev config may leads to recv issue. Please use subdev=\"A:0 A:1\" to configure multiple subdev"
        end
    end
    # --- Internal streamer and buffer config
    if MIMO_MODE
        # --- Clock source
        @assert_uhd LibUHD.uhd_usrp_set_clock_source(rx.uhd.pointerUSRP,"internal",0)
        @assert_uhd LibUHD.uhd_usrp_set_time_now(rx.uhd.pointerUSRP,0,0,0)
    end
    # ----------------------------------------------------
    # --- Set up streamer
    # ----------------------------------------------------
    # Streamer policy depends on MIMO mode or not
    # It can be a pain to have one ssytem for every possible configurations, so one can want to define streamer policy later
    if nbAntennaRx > 0 && bypassStreamer == false
        if MIMO_MODE
            uhd_usrp_create_stream(rx;
                                   nbAntennaRx,
                                   channels,
                                   cpu_format,
                                   otw_format,
                                   args,
                                   stream_mode = UHD_STREAM_MODE_NUM_SAMPS_AND_DONE,
                                   stream_now = false,
                                   full_sec_delay = 1,
                                   frac_sec_delay = 0.5,
                                   num_samps = 10_000_000
                                  )
        else
            uhd_usrp_create_stream(rx;
                                   nbAntennaRx,
                                   channels,
                                   cpu_format,
                                   otw_format,
                                   args,
                                   stream_mode = UHD_STREAM_MODE_START_CONTINUOUS,
                                   stream_now = true,
                                   full_sec_delay = 0,
                                   frac_sec_delay = 0,
                                   num_samps = rx.packetSize
                                  )
        end
    else
        (nbAntennaRx > 0 ) &&  @warn "Streamer has not be set up as bypassStreamer has been set to true. Rx will not work without create a dedicated streamer and set up the streamer policy. See uhd_usrp_create_stream to create a streamer with custom parameters"
    end
    return rx
end

"""
Create and launch a streamer with dedicated parameters
"""
function uhd_usrp_create_stream(rx::UHDRx;nbAntennaRx = 1, channels=[0],cpu_format="fc32",otw_format="sc16",args="",stream_mode=UHD_STREAM_MODE_START_CONTINUOUS,stream_now=true,full_sec_delay=0,frac_sec_delay=0,num_samps=-1)
        # --- Streamer arguments
        a1			   =  Base.unsafe_convert(Cstring,cpu_format);
        a2			   =  Base.unsafe_convert(Cstring,otw_format);
        a3			   =  Base.unsafe_convert(Cstring,args);
        channel        = pointer(channels)
        uhdArgs_0	   = uhd_stream_args_t(a1,a2,a3,channel,nbAntennaRx);
        pointerArgs	  = Ref{uhd_stream_args_t}(uhdArgs_0);
        pointerSamples = Ref{Csize_t}(0);
        # --- Get internal buffer size
        @assert_uhd uhd_usrp_get_rx_stream(rx.uhd.pointerUSRP,pointerArgs,rx.uhd.pointerStreamer)
        @assert_uhd uhd_rx_streamer_max_num_samps(rx.uhd.pointerStreamer,pointerSamples)
        rx.packetSize= pointerSamples[]
        # --- Launch the stream
        restartStreamer(rx;stream_mode,num_samps,stream_now,full_sec_delay,frac_sec_delay)
end

"""
Restart the USRP streamer. In some cases (especially macOS) we have an issue with streamer congestion and we need to restart it.
By now, we have added the restart function in recv! method.
"""
@inline function restartStreamer(rx::UHDRx;stream_mode=UHD_STREAM_MODE_START_CONTINUOUS,num_samps=-1,stream_now=true,full_sec_delay=0,frac_sec_delay=0)
    # ----------------------------------------------------
    # --- Start the streamer
    # ----------------------------------------------------
    (num_samps == -1) && (num_samps = rx.packetSize)
    streamCmd	= UHDBindings.uhd_stream_cmd_t(stream_mode,num_samps,stream_now,full_sec_delay,frac_sec_delay);
    pointerCmd	= Ref{uhd_stream_cmd_t}(streamCmd);
    uhd_rx_streamer_issue_stream_cmd(rx.uhd.pointerStreamer,pointerCmd)

end

function updateSamplingRate!(radio::UHDRx,samplingRate,chan=0)
	# ----------------------------------------------------
	# --- Sampling rate configuration
	# ----------------------------------------------------
	# @inforx  "Try to change rate from $(radio.samplingRate/1e6) MHz to $(samplingRate/1e6) MHz";
	# --- Update the Rx sampling rate
    uhd_usrp_set_rx_rate(radio.uhd.pointerUSRP,samplingRate,chan)
	# --- Get the Rx rate from the radio
	pointerRate  = Ref{Cdouble}(0);
    uhd_usrp_get_rx_rate(radio.uhd.pointerUSRP,chan,pointerRate)
	updateRate  = pointerRate[];
	# --- Print a flag
	if updateRate != samplingRate
		@warnrx "Effective Rate is $(updateRate/1e6) MHz and not $(samplingRate/1e6) MHz\n"
	else
		# @inforx "Effective Rate is $(updateRate/1e6) MHz\n";
	end
	radio.samplingRate = updateRate;
end

function updateGain!(radio::UHDRx,gain,chan=0)
	# ----------------------------------------------------
	# --- Sampling rate configuration
	# ----------------------------------------------------
	# @inforx  "Try to change gain from $(radio.gain) dB to $(gain) dB";
	# Update the UHD sampling rate
    uhd_usrp_set_rx_gain(radio.uhd.pointerUSRP,gain,chan,"")
	# Get the updated gain from UHD
	pointerGain	  = Ref{Cdouble}(0);
    uhd_usrp_get_rx_gain(radio.uhd.pointerUSRP,chan,"",pointerGain)
	updateGain	  = pointerGain[];
	# --- Print a flag
	if updateGain != gain
		@warnrx "Effective gain is $(updateGain) dB and not $(gain) dB\n"
	else
		# @inforx "Effective gain is $(updateGain) dB\n";
	end
	radio.gain = updateGain;
	return updateGain;
end

function updateCarrierFreq!(radio::UHDRx,carrierFreq,chan=0)
	# ----------------------------------------------------
	# --- Carrier Frequency configuration
	# ----------------------------------------------------
	# @inforx  "Try to change carrier frequency from $(radio.carrierFreq/1e6) MHz to $(carrierFreq/1e6) MHz";
	a5			   =  Base.unsafe_convert(Cstring,"");
	tuneRequest	   = uhd_tune_request_t(carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,200e6,a5)
	# tuneRequest   = uhd_tune_request_t(carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,UHD_TUNE_REQUEST_POLICY_AUTO);
	tunePointer	  = Ref{uhd_tune_request_t}(tuneRequest);
	pointerTuneResult	  = Ref{uhd_tune_result_t}();
    uhd_usrp_set_rx_freq(radio.uhd.pointerUSRP,tunePointer,chan,pointerTuneResult)
	pointerCarrierFreq = Ref{Cdouble}(0);
	# sleep(0.001);
    uhd_usrp_get_rx_freq(radio.uhd.pointerUSRP,chan,pointerCarrierFreq)
	updateCarrierFreq	= pointerCarrierFreq[];
	if updateCarrierFreq != carrierFreq
		@warnrx "Effective carrier frequency is $(updateCarrierFreq/1e6) MHz and not $(carrierFreq/1e6) MHz\n"
	else
		# @inforx "Effective carrier frequency is $(updateCarrierFreq/1e6) MHz\n";
	end
	radio.carrierFreq = updateCarrierFreq;
	return updateCarrierFreq;
end

"""
Get a single buffer from the USRP device, and create all the necessary ressources

# --- Syntax

sig	  = recv(radio,nbSamples)
# --- Input parameters
- radio	  : UHD object [UHDRx]
- nbSamples : Desired number of samples [Int]
# --- Output parameters
- sig	  : baseband signal from radio [Array{Complex{CFloat}},radio.packetSize]
"""
function recv(radio::UHDRx,nbSamples);
        # --- Populate the buffer
        sig     = [zeros(Complex{Cfloat},nbSamples) for _ ∈ 1:radio.nbAntennaRx]
        nbSamples = recv!(sig,radio);
        if radio.nbAntennaRx == 1
            # In SISO mode, we return a vector
            return sig[1]
        else
            # In MIMO mode, we return a vector of vector
            return sig
        end
end

"""
Get a single buffer from the USRP device, using the Buffer structure

# --- Syntax

recv!(sig,radio,nbSamples)
# --- Input parameters
- sig	  : Complex signal to populate [Array{Complex{Cfloat}}]
- radio	  : UHD object [UHDRx]
- buffer  : Buffer object (obtained with setBuffer(radio)) [Buffer]
# --- Output parameters
-
"""
function recv!(sig::Vector{Vector{Complex{T}}},radio::UHDRx;nbSamples=0,offset=0) where T
    # restartStreamer(radio)
	# --- Defined parameters for multiple buffer reception
	filled		= false;
	# --- Fill the input buffer @ a specific offset
	if offset == 0
		posT		= Csize_t(0);
	else
		posT 		= Csize_t(offset);
	end
	# --- Managing desired size and buffer size
	if nbSamples == 0
		# --- Fill all the buffer
		# x2 as input is complex and we feed real words
        nbSamples	= Csize_t(length(sig[1]));
	else
		# ---  x2 as input is complex and we feed real words
		nbSamples 	= Csize_t(nbSamples);
		# --- Ensure that the allocation is possible
        @assert nbSamples < (length(sig[1])+posT) "Impossible to fill the buffer (number of samples > residual size";
	end
    global SIG = sig
	while !filled
		# --- Get a buffer: We should have radio.packetSize or less
		# radio.packetSize is the complex size, so x2
		(posT+radio.packetSize> nbSamples) ? n = nbSamples - posT : n = radio.packetSize;
		# --- To avoid memcopy, we direclty feed the pointer at the appropriate solution
        ptr = [pointer(sig[k],1+posT) for k ∈ 1 : radio.nbAntennaRx]
		# --- Populate buffer with radio samples
		cSamples 	= populateBuffer!(radio,ptr,n);
		# --- Update counters
		posT += cSamples;
		# @show Int(cSamples),Int(posT)
		# --- Breaking flag
		(posT == nbSamples) ? filled = true : filled = false;
	end
	return posT
end

"""
Calling UHD function wrapper to fill a buffer. It is preferable to cal this function though the use of recv or recv!

# --- Syntax

recv!(sig,radio,nbSamples)
# --- Input parameters
- radio	  	: UHD object [UHDRx]
- ptr  		: Writable memory position [Ref{Ptr{Cvoid}}]
- nbSamples : Number of samples to acquire
# --- Output parameters
- nbSamp 	: Number of samples fill in buffer [Csize_t]
"""
function populateBuffer!(radio,ptr,nbSamples::Csize_t=0)
	# --- Getting number of samples
	# If not specified, we get back to radio.packetSize
	if (nbSamples == Csize_t(0))
		nbSamples = radio.packetSize;
	end
	#@assert nbSamples <= length(buffer.x) "Number of desired samples can not be greater than buffer size";
	# --- Effectively recover data
	pointerSamples = Ref{Csize_t}(0);
    uhd_rx_streamer_recv(radio.uhd.pointerStreamer,ptr,nbSamples,radio.uhd.addressMD,0.1,false,pointerSamples)
		# --- Pointer deferencing
	return pointerSamples[];
end#

"""
Returns the Error flag of the current UHD burst

# --- Syntax

flag = getError(radio)
# --- Input parameters
- radio : UHD object [UHDRx]
# --- Output parameters
- err	: Error Flag [error_code_t]
"""
function getError(radio::UHDRx)
	ptrErr = Ref{uhd_rx_metadata_error_code_t}();
    uhd_rx_metadata_error_code(radio.uhd.pointerMD,ptrErr)
	return err = ptrErr[];
end

"""
Return the timestamp of the last UHD burst

# --- Syntax

(second,fracSecond) = getTimestamp(radio)
# --- Input parameters
- radio	  : UHD UHD object [UHDRx]
# --- Output parameters
- second  : Second value for the flag [Int]
- fracSecond : Fractional second value [Float64]
"""
function getTimestamp(radio::UHDRx)
	ptrFullSec = Ref{FORMAT_LONG}();
	ptrFracSec = Ref{Cdouble}();
    uhd_rx_metadata_time_spec(radio.uhd.pointerMD,ptrFullSec,ptrFracSec)
	return (ptrFullSec[],ptrFracSec[]);
end

function Base.close(radio::UHDRx)
	# --- Checking realease nature
	# There is one flag to avoid double close (that leads to seg fault)
	if radio.released == 0
        @assert_uhd uhd_rx_streamer_free(radio.uhd.addressStream)
        @assert_uhd uhd_rx_metadata_free(radio.uhd.addressMD)
	else
		# print a warning
		@warn "UHD ressource was already released, abort call";
	end
	# --- Force flag value
	radio.released = 1;
end

function Base.print(radio::UHDRx,chan=0)
	# Get the gain from UHD
	pointerGain	  = Ref{Cdouble}(0);
    uhd_usrp_get_rx_gain(radio.uhd.pointerUSRP,chan,"",pointerGain)
	updateGain	  = pointerGain[];
	# Get the rate from UHD
	pointerRate	  = Ref{Cdouble}(0);
    uhd_usrp_get_rx_rate(radio.uhd.pointerUSRP,chan,pointerRate)
	updateRate	  = pointerRate[];
	# Get the freq from UHD
	pointerFreq	  = Ref{Cdouble}(0);
    uhd_usrp_get_rx_freq(radio.uhd.pointerUSRP,chan,pointerFreq)
	updateFreq	  = pointerFreq[];
	# Print message
	strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n Rx Gain: %2.2f dB\n",updateFreq/1e6,updateRate/1e6,updateGain);
	@inforx "Current UHD Configuration in Rx mode\n$strF";
end

max_stars_repo_name: RGerzaguet/UHDBindings

max_issues_repo_name: RGerzaguet/UHDBindings


text:
#=
A function f is Lipschitz continous on [a, b]
    |f(x)-f(y)| < l*|x-y|    for all x, y \in [a, b], l > 0
l is as large as the largest unsigned instantaneous
rate of change the function attains on [a,b].

Given a point (x0, f(x0)), The line f(x0)-l(x-x0) for x > x0
and f(x0) + l(x-x0) for x < x0 form a lower bound of f.
=#

struct Pt
    x
    y
end

function get_sp_intersection(A, B, L)
    t = ((A.y-B.y)-L*(A.x-B.x)) / 2L;
    return Pt(A.x+t, A.y-t*L);
end

function shubert_piyavskii(f, a, b, l, ϵ, δ=0.01)
    m = (a+b)/2;
    A, M, B = Pt(a, f(a)), Pt(m, f(m)), Pt(b, f(b))
    pts = [A, get_sp_intersection(A, M, l),
    M, get_sp_intersection(M, B, l), B];
    Δ = Inf;

    while Δ > ϵ
        i = argmin([P.y for P in pts]);
        P = Pt(pts[i].x, f(pts[i].x));
        Δ = P.y = pts[i].y;

        P_prev = get_sp_intersection(pts[i-1], P, l);
        P_next = get_sp_intersection(P, pts[i+1], l);

        deleteat!(pts, i);
        insert!(pts, i, P_next);
        insert!(pts, i, P);
        insert!(pts, i, P_prev);
    end

    intervals = [];
    P_min = pts[2*(argmin([P.y for P in pts[1:2:end]])) - 1];
    y_min = P_min.y;
    for i in 2:2:length(pts)
        if pts[i].y < y_min
            dy = y_min - pts[i].y;
            x_low = max(a, pts[i].x - dy/l);
            x_high = max(b, pts[i].x - dy/l);
            if !isempty(intervals) && intervals[end][2] + δ ≥ x_low
                intervals[end] = (intervals[end][1], x_high);
            else
                push!(intervals, (x_low, x_high));
            end
        end
    end
    return (P_min, intervals);
end


text:
module LibMagick

import Base: error, size

export MagickWand,
    constituteimage,
    exportimagepixels!,
    getblob,
    getimagealphachannel,
    getimagecolorspace,
    getimagedepth,
    getnumberimages,
    importimagepixels,
    readimage,
    resetiterator,
    setimagecolorspace,
    setimagecompression,
    setimageformat,
    writeimage

# Find the library
mpath = get(ENV, "MAGICK_HOME", "") # If MAGICK_HOME is defined, add to library search path
mpaths = isempty(mpath) ? ASCIIString[] : [mpath, joinpath(mpath, "lib")]
libnames = ["libMagickWand", "CORE_RL_wand_"]
suffixes = ["", "-Q16", "-6.Q16", "-Q8"]
options = ["", "HDRI"]
const libwand = find_library(vec(libnames.*transpose(suffixes).*reshape(options,(1,1,length(options)))), mpaths)
const have_imagemagick = !isempty(libwand)

# Initialize the library
function init()
    global libwand
    if have_imagemagick
        eval(:(ccall((:MagickWandGenesis, $libwand), Void, ())))
    else
        warn("ImageMagick utilities not found. Install for more file format support.")
    end
end

# Constants
# Storage types
const CHARPIXEL = 1
const DOUBLEPIXEL = 2
const FLOATPIXEL = 3
const INTEGERPIXEL = 4
const SHORTPIXEL = 7
IMStorageTypes = Union(Uint8, Uint16, Uint32, Float32, Float64)
storagetype(::Type{Uint8}) = CHARPIXEL
storagetype(::Type{Uint16}) = SHORTPIXEL
storagetype(::Type{Uint32}) = INTEGERPIXEL
storagetype(::Type{Float32}) = FLOATPIXEL
storagetype(::Type{Float64}) = DOUBLEPIXEL

# Image type
const IMType = ["BilevelType", "GrayscaleType", "GrayscaleMatteType", "PaletteType", "PaletteMatteType", "TrueColorType", "TrueColorMatteType", "ColorSeparationType", "ColorSeparationMatteType", "OptimizeType", "PaletteBilevelMatteType"]
const IMTypedict = Dict(IMType, 1:length(IMType))

const CStoIMTypedict = ["Gray" => "GrayscaleType", "GrayAlpha" => "GrayscaleMatteType", "RGB" => "TrueColorType", "ARGB" => "TrueColorMatteType", "CMYK" => "ColorSeparationType"]

# Colorspace
const IMColorspace = ["RGB", "Gray", "Transparent", "OHTA", "Lab", "XYZ", "YCbCr", "YCC", "YIQ", "YPbPr", "YUV", "CMYK", "sRGB"]
const IMColordict = Dict(IMColorspace, 1:length(IMColorspace))
IMColordict["GrayAlpha"] = IMColordict["Gray"]
IMColordict["RGBA"] = IMColordict["RGB"]
IMColordict["ARGB"] = IMColordict["RGB"]

function nchannels(imtype::String, cs::String, havealpha = false)
    n = 3
    if beginswith(imtype, "Grayscale") || beginswith(imtype, "Bilevel")
        n = 1
        cs = havealpha ? "GrayAlpha" : "Gray"
    elseif cs == "CMYK"
        n = 4
    else
        cs = havealpha ? "ARGB" : "RGB" # only remaining variants supported by exportimagepixels
    end
    n + havealpha, cs
end

channelorder = ["Gray" => "I", "GrayAlpha" => "IA", "RGB" => "RGB", "ARGB" => "ARGB", "RGBA" => "RGBA", "CMYK" => "CMYK"]

# Compression
const NoCompression = 1

type MagickWand
    ptr::Ptr{Void}
end

function MagickWand()
    wand = MagickWand(ccall((:NewMagickWand, libwand), Ptr{Void}, ()))
    finalizer(wand, destroymagickwand)
    wand
end

destroymagickwand(wand::MagickWand) = ccall((:DestroyMagickWand, libwand), Ptr{Void}, (Ptr{Void},), wand.ptr)

type PixelWand
    ptr::Ptr{Void}
end

function PixelWand()
    wand = PixelWand(ccall((:NewPixelWand, libwand), Ptr{Void}, ()))
    finalizer(wand, destroypixelwand)
    wand
end

destroypixelwand(wand::PixelWand) = ccall((:DestroyPixelWand, libwand), Ptr{Void}, (Ptr{Void},), wand.ptr)

const IMExceptionType = Array(Cint, 1)
function error(wand::MagickWand)
    pMsg = ccall((:MagickGetException, libwand), Ptr{Uint8}, (Ptr{Void}, Ptr{Cint}), wand.ptr, IMExceptionType)
    msg = bytestring(pMsg)
    relinquishmemory(pMsg)
    error(msg)
end
function error(wand::PixelWand)
    pMsg = ccall((:PixelGetException, libwand), Ptr{Uint8}, (Ptr{Void}, Ptr{Cint}), wand.ptr, IMExceptionType)
    msg = bytestring(pMsg)
    relinquishmemory(pMsg)
    error(msg)
end

function getsize(buffer, colorspace)
    if colorspace == "Gray"
        return size(buffer, 1), size(buffer, 2), size(buffer, 3)
    else
        return size(buffer, 2), size(buffer, 3), size(buffer, 4)
    end
end

function exportimagepixels!{T}(buffer::AbstractArray{T}, wand::MagickWand,  colorspace::ASCIIString; x = 0, y = 0)
    cols, rows, nimages = getsize(buffer, colorspace)
    ncolors = colorspace == "Gray" ? 1 : size(buffer, 1)
    p = pointer(buffer)
    for i = 1:nimages
        status = ccall((:MagickExportImagePixels, libwand), Cint, (Ptr{Void}, Cssize_t, Cssize_t, Csize_t, Csize_t, Ptr{Uint8}, Cint, Ptr{Void}), wand.ptr, x, y, cols, rows, channelorder[colorspace], storagetype(T), p)
        status == 0 && error(wand)
        nextimage(wand)
        p += sizeof(T)*cols*rows*ncolors
    end
    buffer
end

# function importimagepixels{T}(buffer::AbstractArray{T}, wand::MagickWand, colorspace::ASCIIString; x = 0, y = 0)
#     cols, rows = getsize(buffer, colorspace)
#     status = ccall((:MagickImportImagePixels, libwand), Cint, (Ptr{Void}, Cssize_t, Cssize_t, Csize_t, Csize_t, Ptr{Uint8}, Cint, Ptr{Void}), wand.ptr, x, y, cols, rows, channelorder[colorspace], storagetype(T), buffer)
#     status == 0 && error(wand)
#     nothing
# end

function constituteimage{T<:Unsigned}(buffer::AbstractArray{T}, wand::MagickWand, colorspace::ASCIIString; x = 0, y = 0)
    cols, rows, nimages = getsize(buffer, colorspace)
    ncolors = colorspace == "Gray" ? 1 : size(buffer, 1)
    p = pointer(buffer)
    for i = 1:nimages
        status = ccall((:MagickConstituteImage, libwand), Cint, (Ptr{Void}, Cssize_t, Cssize_t, Ptr{Uint8}, Cint, Ptr{Void}), wand.ptr, cols, rows, channelorder[colorspace], storagetype(T), p)
        status == 0 && error(wand)
        setimagecolorspace(wand, colorspace)
        status = ccall((:MagickSetImageDepth, libwand), Cint, (Ptr{Void}, Csize_t), wand.ptr, 8*sizeof(T))
        status == 0 && error(wand)
        p += sizeof(T)*cols*rows*ncolors
    end
    nothing
end

function getblob(wand::MagickWand, format::String)
    setimageformat(wand, format)
    len = Array(Csize_t, 1)
    ptr = ccall((:MagickGetImagesBlob, libwand), Ptr{Uint8}, (Ptr{Void}, Ptr{Csize_t}), wand.ptr, len)
    blob = pointer_to_array(ptr, int(len[1]))
    finalizer(blob, relinquishmemory)
    blob
end

function pingimage(wand::MagickWand, filename::String)
    status = ccall((:MagickPingImage, libwand), Cint, (Ptr{Void}, Ptr{Uint8}), wand.ptr, filename)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, filename::String)
    status = ccall((:MagickReadImage, libwand), Cint, (Ptr{Void}, Ptr{Uint8}), wand.ptr, filename)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, stream::IO)
    status = ccall((:MagickReadImageFile, libwand), Cint, (Ptr{Void}, Ptr{Void}), wand.ptr, fdopen(stream))
    status == 0 && error(wand)
    nothing
end

function writeimage(wand::MagickWand, filename::String)
    status = ccall((:MagickWriteImages, libwand), Cint, (Ptr{Void}, Ptr{Uint8}, Cint), wand.ptr, filename, true)
    status == 0 && error(wand)
    nothing
end

function size(wand::MagickWand)
    height = ccall((:MagickGetImageHeight, libwand), Csize_t, (Ptr{Void},), wand.ptr)
    width = ccall((:MagickGetImageWidth, libwand), Csize_t, (Ptr{Void},), wand.ptr)
    return int(width), int(height)
end

getnumberimages(wand::MagickWand) = int(ccall((:MagickGetNumberImages, libwand), Csize_t, (Ptr{Void},), wand.ptr))

nextimage(wand::MagickWand) = ccall((:MagickNextImage, libwand), Cint, (Ptr{Void},), wand.ptr) == 1

resetiterator(wand::MagickWand) = ccall((:MagickResetIterator, libwand), Void, (Ptr{Void},), wand.ptr)

newimage(wand::MagickWand, cols::Integer, rows::Integer, pw::PixelWand) = ccall((:MagickNewImage, libwand), Cint, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}), wand.ptr, cols, rows, pw.ptr) == 0 && error(wand)

# test whether image has an alpha channel
getimagealphachannel(wand::MagickWand) = ccall((:MagickGetImageAlphaChannel, libwand), Cint, (Ptr{Void},), wand.ptr) == 1

# # get number of colors in the image
# magickgetimagecolors(wand::MagickWand) = ccall((:MagickGetImageColors, libwand), Csize_t, (Ptr{Void},), wand.ptr)

# get the type
function getimagetype(wand::MagickWand)
    t = ccall((:MagickGetImageType, libwand), Cint, (Ptr{Void},), wand.ptr)
    # Apparently the following is necessary, because the type is "potential"
    ccall((:MagickSetImageType, libwand), Void, (Ptr{Void}, Cint), wand.ptr, t)
    1 <= t <= length(IMType) || error("Image type ", t, " not recognized")
    IMType[t]
end

# get the colorspace
function getimagecolorspace(wand::MagickWand)
    cs = ccall((:MagickGetImageColorspace, libwand), Cint, (Ptr{Void},), wand.ptr)
    1 <= cs <= length(IMColorspace) || error("Colorspace ", cs, " not recognized")
    IMColorspace[cs]
end

# set the colorspace
function setimagecolorspace(wand::MagickWand, cs::ASCIIString)
    status = ccall((:MagickSetImageColorspace, libwand), Cint, (Ptr{Void},Cint), wand.ptr, IMColordict[cs])
    status == 0 && error(wand)
    nothing
end

# set the compression
function setimagecompression(wand::MagickWand, compression::Integer)
    status = ccall((:MagickSetImageCompression, libwand), Cint, (Ptr{Void},Cint), wand.ptr, int32(compression))
    status == 0 && error(wand)
    nothing
end

# set the image format
function setimageformat(wand::MagickWand, format::ASCIIString)
    status = ccall((:MagickSetImageFormat, libwand), Cint, (Ptr{Void}, Ptr{Uint8}), wand.ptr, format)
    status == 0 && error(wand)
    nothing
end

# get the pixel depth
getimagedepth(wand::MagickWand) = int(ccall((:MagickGetImageDepth, libwand), Csize_t, (Ptr{Void},), wand.ptr))

pixelsetcolor(wand::PixelWand, colorstr::ByteString) = ccall((:PixelSetColor, libwand), Csize_t, (Ptr{Void},Ptr{Uint8}), wand.ptr, colorstr) == 0 && error(wand)

relinquishmemory(p) = ccall((:MagickRelinquishMemory, libwand), Ptr{Uint8}, (Ptr{Uint8},), p)

# Implementation of fdopen
# Delete this once 0.3 or greater is common, and use the CFILE mechanism
function fdopen(s::IO)
    @unix_only FILEp = ccall(:fdopen, Ptr{Void}, (Cint, Ptr{Uint8}), convert(Cint, fd(s)), "r")
    @windows_only FILEp = ccall(:_fdopen, Ptr{Void}, (Cint, Ptr{Uint8}), convert(Cint, fd(s)), "r")
    systemerror("fdopen failed", FILEp == 0)
    status = ccall(:fseek, Cint, (Ptr{Void}, Clong, Cint), FILEp, convert(Clong, position(s)), int32(0))
    systemerror("fseek failed", status == 0)
    FILEp
end

end

LibMagick.have_imagemagick


text:
#FCL.jl - an implementation of the "fully connected layer"

module FCL

import ..@import_interface
@import_interface

import ..TF.sigmoid

immutable FullyConnectedLayer{F, inputlength, outputlength, transferfunction} <: Layer{F}
  bias::Vector{F}
  transition::Matrix{F}
end

(fcl::FullyConnectedLayer{F,i,o,tf}){F,i,o,tf}(v::Array) = ml_call(fcl, v)

function (::Type{FullyConnectedLayer{F,i,o,tf}}){F,i,o,tf}(randomizationfunction::Function = zeros)
  biasvector = randomizationfunction(F, o)
  transitionmatrix = randomizationfunction(F, o, i)
  FullyConnectedLayer{F,i,o,tf}(biasvector, transitionmatrix)
end

#implementation of the interfaces
inputs{F, i, o, tf}(::Type{FullyConnectedLayer{F, i, o, tf}}) = i
outputs{F, i, o, tf}(::Type{FullyConnectedLayer{F, i, o, tf}}) = o
parameters{F, i, o, tf}(::Type{FullyConnectedLayer{F, i, o, tf}})  = (i + 1) * o

include("./serialization.jl")
include("./evaluation.jl")
include("./backpropagation.jl")

end #module

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
module TestClassWithManyMethods

using CxxCall
using Test
using Libdl

dir = mktempdir()
lib = joinpath(dir, "libClassWithManyMethods")
filepath = joinpath(dir, "ClassWithManyMethods.cxx")
eval(cxxsetup())
eval(cxxnewfile(filepath,
"""
template <typename T>
struct ClassWithManyMethods {
    T number;
    T zero_more()      {return number+0;}
    T one_more()   {return number+1;}
    T two_more()   {return number+2;}
    T three_more() {return number+3;}
    T four_more()  {return number+4;}
    T five_more()  {return number+5;}
};
"""))
struct ClassWithManyMethods{T}
    _data::T
end
function CxxCall.tocxx(::Type{ClassWithManyMethods{T}}) where {T}
    "ClassWithManyMethods<$(tocxx(T))>"
end
for T in [Int, Float64]
    @cxx lib zero_more(self::ClassWithManyMethods{T})::T = "return self.zero_more();"
    # lets define the other methods in a loop
    for f in [
        :one_more,
        :two_more,
        :three_more,
        :four_more,
        :five_more,
        ]
        CT = ClassWithManyMethods{T}
        @cxx lib function $f(self::CT)::T
            """
            // lets check we are ABI compatible
            static_assert(sizeof($(tocxx(CT))) == $(sizeof(CT)));
            return self.$f();
            """
        end
    end
end

@testset "ClassWithManyMethods" begin
    @test !ispath(filepath)
    cxx_write_code!()
    @test isfile(filepath)
    libpath = lib * ".so"
    run(`g++ -shared -fPIC $filepath -o $libpath`)
    @test isfile(libpath)
    dlopen(libpath)

    obj = ClassWithManyMethods(42)
    @test zero_more(obj)  === 42+0
    @test one_more(obj)   === 42+1
    @test two_more(obj)   === 42+2
    @test three_more(obj) === 42+3
    @test four_more(obj)  === 42+4
    @test five_more(obj)  === 42+5

    obj = ClassWithManyMethods(42.0)
    @test zero_more(obj)  === 42.0+0
    @test one_more(obj)   === 42.0+1
    @test two_more(obj)   === 42.0+2
    @test three_more(obj) === 42.0+3
    @test four_more(obj)  === 42.0+4
    @test five_more(obj)  === 42.0+5
end

end#module


text:
### truncation
Base.trunc(dt::Date,p::Type{Year}) = Date(UTD(totaldays(year(dt),1,1)))
Base.trunc(dt::Date,p::Type{Month}) = firstdayofmonth(dt)
Base.trunc(dt::Date,p::Type{Day}) = dt

Base.trunc(dt::DateTime,p::Type{Year}) = DateTime(trunc(Date(dt),Year))
Base.trunc(dt::DateTime,p::Type{Month}) = DateTime(trunc(Date(dt),Month))
Base.trunc(dt::DateTime,p::Type{Day}) = DateTime(Date(dt))
Base.trunc(dt::DateTime,p::Type{Hour}) = dt - Minute(dt) - Second(dt) - Millisecond(dt)
Base.trunc(dt::DateTime,p::Type{Minute}) = dt - Second(dt) - Millisecond(dt)
Base.trunc(dt::DateTime,p::Type{Second}) = dt - Millisecond(dt)
Base.trunc(dt::DateTime,p::Type{Millisecond}) = dt

# Adjusters
firstdayofweek(dt::Date) = Date(UTD(value(dt) - dayofweek(dt) + 1))
firstdayofweek(dt::DateTime) = DateTime(firstdayofweek(Date(dt)))
lastdayofweek(dt::Date) = Date(UTD(value(dt) + (7-dayofweek(dt))))
lastdayofweek(dt::DateTime) = DateTime(lastdayofweek(Date(dt)))

@vectorize_1arg TimeType firstdayofweek
@vectorize_1arg TimeType lastdayofweek

firstdayofmonth(dt::Date) = Date(UTD(value(dt)-day(dt)+1))
firstdayofmonth(dt::DateTime) = DateTime(firstdayofmonth(Date(dt)))
function lastdayofmonth(dt::Date)
    y,m,d = yearmonthday(dt)
    return Date(UTD(value(dt)+daysinmonth(y,m)-d))
end
lastdayofmonth(dt::DateTime) = DateTime(lastdayofmonth(Date(dt)))

@vectorize_1arg TimeType firstdayofmonth
@vectorize_1arg TimeType lastdayofmonth

firstdayofyear(dt::Date) = Date(UTD(value(dt)-dayofyear(dt)+1))
firstdayofyear(dt::DateTime) = DateTime(firstdayofyear(Date(dt)))
function lastdayofyear(dt::Date)
    y,m,d = yearmonthday(dt)
    return Date(UTD(value(dt)+daysinyear(y)-dayofyear(y,m,d)))
end

@vectorize_1arg TimeType firstdayofyear
@vectorize_1arg TimeType lastdayofyear

function firstdayofquarter(dt::Date)
    y,m = yearmonth(dt)
    mm = m < 4 ? 1 : m < 7 ? 4 : m < 10 ? 7 : 10
    return Date(y,mm,1)
end
function lastdayofquarter(dt::Date)
    y,m = yearmonth(dt)
    mm,d = m < 4 ? (3,31) : m < 7 ? (6,30) : m < 10 ? (9,30) : (12,31)
    return Date(y,mm,d)
end
firstdayofquarter(dt::DateTime) = DateTime(firstdayofquarter(Date(dt)))
lastdayofquarter(dt::DateTime) = DateTime(lastdayofquarter(Date(dt)))
@vectorize_1arg TimeType firstdayofquarter
@vectorize_1arg TimeType lastdayofquarter

# Temporal Adjusters
immutable DateFunction
    f::Function
    # validate boolean, single-arg inner constructor
    function DateFunction(f::Function,negate::Bool,dt::TimeType)
        try
            f(dt) in (true,false) || throw(ArgumentError("Provided function must take a single TimeType argument and return true or false"))
        catch e
            throw(ArgumentError("Provided function must take a single TimeType argument"))
        end
        n = negate ? (!) : identity
        return new(@eval x->$n($f(x)))
    end
end
Base.show(io::IO,df::DateFunction) = println(df.f.code)

# Core adjuster
function adjust(df::DateFunction,start,step,limit)
    for i = 1:limit
        df.f(start) && return start
        start += step
    end
    throw(ArgumentError("Adjustment limit reached: $limit iterations"))
end

function adjust(func::Function,start;step::Period=Day(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,start),start,step,limit)
end

# Constructors using DateFunctions
function Date(func::Function,y,m=1,d=1;step::Period=Day(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,Date(y,m,d)),Date(y,m,d),step,limit)
end

function DateTime(func::Function,y,m=1;step::Period=Day(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,DateTime(y,m)),DateTime(y,m),step,limit)
end
function DateTime(func::Function,y,m,d;step::Period=Hour(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,DateTime(y)),DateTime(y,m,d),step,limit)
end
function DateTime(func::Function,y,m,d,h;step::Period=Minute(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,DateTime(y)),DateTime(y,m,d,h),step,limit)
end
function DateTime(func::Function,y,m,d,h,mi;step::Period=Second(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,DateTime(y)),DateTime(y,m,d,h,mi),step,limit)
end
function DateTime(func::Function,y,m,d,h,mi,s;step::Period=Millisecond(1),negate::Bool=false,limit::Int=10000)
    return adjust(DateFunction(func,negate,DateTime(y)),DateTime(y,m,d,h,mi,s),step,limit)
end

# Return the next TimeType that falls on dow
ISDAYOFWEEK = [Mon=>DateFunction(ismonday,false,Date(0)),
               Tue=>DateFunction(istuesday,false,Date(0)),
               Wed=>DateFunction(iswednesday,false,Date(0)),
               Thu=>DateFunction(isthursday,false,Date(0)),
               Fri=>DateFunction(isfriday,false,Date(0)),
               Sat=>DateFunction(issaturday,false,Date(0)),
               Sun=>DateFunction(issunday,false,Date(0))]

# "same" indicates whether the current date can be considered or not
tonext(dt::TimeType,dow::Int;same::Bool=false) = adjust(ISDAYOFWEEK[dow],same ? dt : dt+Day(1),Day(1),7)
# Return the next TimeType where func evals true using step in incrementing
function tonext(func::Function,dt::TimeType;step::Period=Day(1),negate::Bool=false,limit::Int=10000,same::Bool=false)
    return adjust(DateFunction(func,negate,dt),same ? dt : dt+step,step,limit)
end

toprev(dt::TimeType,dow::Int;same::Bool=false) = adjust(ISDAYOFWEEK[dow],same ? dt : dt+Day(-1),Day(-1),7)
function toprev(func::Function,dt::TimeType;step::Period=Day(-1),negate::Bool=false,limit::Int=10000,same::Bool=false)
    return adjust(DateFunction(func,negate,dt),same ? dt : dt+step,step,limit)
end

# Return the first TimeType that falls on dow in the Month or Year
function tofirst(dt::TimeType,dow::Int;of::Union(Type{Year},Type{Month})=Month)
    dt = of <: Month ? firstdayofmonth(dt) : firstdayofyear(dt)
    return adjust(ISDAYOFWEEK[dow],dt,Day(1),366)
end

# Return the last TimeType that falls on dow in the Month or Year
function tolast(dt::TimeType,dow::Int;of::Union(Type{Year},Type{Month})=Month)
    dt = of <: Month ? lastdayofmonth(dt) : lastdayofyear(dt)
    return adjust(ISDAYOFWEEK[dow],dt,Day(-1),366)
end

function recur{T<:TimeType}(fun::Function,start::T,stop::T;step::Period=Day(1),negate::Bool=false,limit::Int=10000)
    ((start != stop) & ((step > zero(step)) != (stop > start))) && return T[]
    a = T[]
    check = start <= stop ? 1 : -1
    df = Dates.DateFunction(fun,negate,start)
    while true
        next = Dates.adjust(df,start,step,limit)
        cmp(next,stop) == check && break
        push!(a,next)
        start = next + step
    end
    return a
end
function recur{T<:TimeType}(fun::Function,dr::StepRange{T};negate::Bool=false,limit::Int=10000)
    return recur(fun,first(dr),last(dr);step=step(dr),negate=negate,limit=limit)
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMIImport
using Test
import Random
using FMIZoo

using FMIImport.FMICore: fmi2Integer, fmi2Boolean, fmi2Real, fmi2String

exportingToolsWindows = [("Dymola", "2022x")]
exportingToolsLinux = [("Dymola", "2022x")]

function runtests(exportingTool)
    ENV["EXPORTINGTOOL"] = exportingTool[1]
    ENV["EXPORTINGVERSION"] = exportingTool[2]

    @testset "Testing FMUs exported from $exportingTool" begin
        @testset "Functions for fmi2Component" begin
            @testset "Variable Getters / Setters" begin
                include("getter_setter.jl")
            end
            @testset "State Manipulation" begin
                include("state.jl")
            end
            @testset "Directional derivatives" begin
                include("dir_ders.jl")
            end
        end

        @testset "Model Description Parsing" begin
            include("model_description.jl")
        end

        @testset "Logging" begin
            include("logging.jl")
        end
    end
end

@testset "FMIImport.jl" begin
    if Sys.iswindows()
        @info "Automated testing is supported on Windows."
        for exportingTool in exportingToolsWindows
            runtests(exportingTool)
        end
    elseif Sys.islinux()
        @info "Automated testing is supported on Linux."
        for exportingTool in exportingToolsLinux
            runtests(exportingTool)
        end
    elseif Sys.isapple()
        @warn "Test-sets are currrently using Windows- and Linux-FMUs, automated testing for macOS is currently not supported."
    end
end


text:
function fwrite_note_quake!(C::Union{SeisHdr, SeisSrc}, method::String, fname::String, opts::String)
  wstr = string(timestamp(), " ¦ write ¦ ", method, "(H", opts, ") ¦ wrote to file ", fname)
  push!(C.notes, wstr)
  return nothing
end


text:
function initialize!(integrator,cache::AitkenNevilleCache)
  integrator.kshortsize = 2
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1

  cache.step_no = 1
  cache.cur_order = max(integrator.alg.init_order, integrator.alg.min_order)
end

function perform_step!(integrator,cache::AitkenNevilleCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,T,utilde,atmp,dtpropose,cur_order,A = cache

  @muladd @. u = uprev + dt*fsalfirst

  T[1,1] = copy(u)

  halfdt = dt/2
  for i in 2:size(T)[1]
    # Solve using Euler method
    @muladd @. u = uprev + halfdt*fsalfirst
    f(k, u, p, t+halfdt)
    integrator.destats.nf += 1
    for j in 2:2^(i-1)
      @muladd @. u = u + halfdt*k
      f(k, u, p, t+j*halfdt)
      integrator.destats.nf += 1
    end
    T[i,1] = copy(u)
    # Richardson Extrapolation
    for j in 2:i
      T[i,j] = ((2^(j-1))*T[i,j-1] - T[i-1,j-1])/((2^(j-1)) - 1)
    end

    halfdt = halfdt/2
  end

  if integrator.opts.adaptive
      minimum_work = Inf
      range_start = max(2,cur_order - 1)
      if cache.step_no == one(cache.step_no)
          range_start = 2
      end

      for i = range_start:min(size(T)[1], cur_order + 1)

          A = 2^(i-1)
          @. utilde = T[i,i] - T[i,i-1]
          atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
          EEst = integrator.opts.internalnorm(atmp,t)

          beta1 = integrator.opts.beta1
          e = integrator.EEst
          qold = integrator.qold

          integrator.opts.beta1 = 1/(i+1)
          integrator.EEst = EEst
          dtpropose = step_accept_controller!(integrator,integrator.alg,stepsize_controller!(integrator,integrator.alg))
          integrator.EEst = e
          integrator.opts.beta1 = beta1
          integrator.qold = qold

          work = A/dtpropose

          if work < minimum_work
              integrator.opts.beta1 = 1/(i+1)
              cache.dtpropose = dtpropose
              cache.cur_order = i
              minimum_work = work
              integrator.EEst = EEst
          end
      end
  end

  # using extrapolated value of u
  @. u = T[cache.cur_order, cache.cur_order]
  cache.step_no = cache.step_no + 1
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

function initialize!(integrator,cache::AitkenNevilleConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  cache.step_no = 1
  cache.cur_order = max(integrator.alg.init_order, integrator.alg.min_order)
end

function perform_step!(integrator,cache::AitkenNevilleConstantCache,repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  @unpack dtpropose, T, cur_order, work, A = cache
  @muladd u = @. uprev + dt*integrator.fsalfirst
  T[1,1] = u
  halfdt = dt/2
  for i in 2:min(size(T)[1], cur_order+1)
    # Solve using Euler method
    @muladd u = @. uprev + halfdt*integrator.fsalfirst
    k = f(u, p, t+halfdt)
    integrator.destats.nf += 1

    for j in 2:2^(i-1)
      @muladd u = @. u + halfdt*k
      k = f(u, p, t+j*halfdt)
      integrator.destats.nf += 1
    end
    T[i,1] = u

    # Richardson Extrapolation
    for j in 2:i
      T[i,j] = ((2^(j-1))*T[i,j-1] - T[i-1,j-1])/((2^(j-1)) - 1)
    end
    halfdt = halfdt/2
  end

  if integrator.opts.adaptive
      minimum_work = Inf
      range_start = max(2,cur_order - 1)
      if cache.step_no == one(cache.step_no)
          range_start = 2
      end

      for i = range_start:min(size(T)[1], cur_order + 1)

          A = 2^(i-1)
          utilde = T[i,i] - T[i,i-1]
          atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
          EEst = integrator.opts.internalnorm(atmp,t)

          beta1 = integrator.opts.beta1
          e = integrator.EEst
          qold = integrator.qold

          integrator.opts.beta1 = 1/(i+1)
          integrator.EEst = EEst
          dtpropose = step_accept_controller!(integrator,integrator.alg,stepsize_controller!(integrator,integrator.alg))
          integrator.EEst = e
          integrator.opts.beta1 = beta1
          integrator.qold = qold

          work = A/dtpropose

          if work < minimum_work
              integrator.opts.beta1 = 1/(i+1)
              cache.dtpropose = dtpropose
              cache.cur_order = i
              minimum_work = work
              integrator.EEst = EEst
          end
      end
  end

  cache.step_no = cache.step_no + 1
  integrator.u = T[cache.cur_order,cache.cur_order]

  k = f(integrator.u, p, t+dt)
  integrator.destats.nf += 1
  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast

  # Use extrapolated value of u
end


text:
using CoordinateTransformations: Translation, LinearMap
using StaticArrays
using LinearAlgebra: BLAS
using RigidBodyDynamics
using RigidBodyDynamics.Spatial: vector_to_skew_symmetric
using MeshCat
using ColorTypes: RGBA, RGB
using MeshCatMechanisms
using Rotations
using PressureFieldContact
using PressureFieldContact.Radau
using LinearAlgebra

set_zero_subnormals(true)
BLAS.set_num_threads(1)  # NOTE: comment out this line if using IntelMKL

is_bristle = true  # set to false to use REGURALIZED friction between the pencil and the pads

# TUNABLE BRISTLE FRICTION PARAMATERS
k̄ = 8.0e4
τ = 0.01
magic = 1.0e-2

# TUNABLE REGURALIZED FRICTION PARAMATERS
v_tol = 1.0e-5

# Geometric Constants
pad_rad = 0.0035
penci_length = 0.16
penci_rad = 0.0035

mutable struct GripperIns
	z::Float64  # height of the arm
	slip_allowed::Float64  # see comments below GripperIns
	θ_arm::Float64  # desired angle of arm
	Δt::Float64  # when should it look for the next instruction
	t0::Float64  # is set to the current time the first time the instruction is used
	rot_0::Float64  # is set to the current ***pencil*** rotation angle the first time the instruction is used
	function GripperIns(; z::Float64=z_high, slip_allowed::Float64=0.0, θ_arm::Float64=0.0, Δt::Float64=2.0)
		# This object contains information about how the gripper should behave.
		# Much of this behavior is encoded in the paramater called slip_allowed as follows:
		# -Inf    --> no grip (fingers are separated )
		#  0.0    --> no slip allowed (fingers grip tightly)
		#  finite --> finite slip allowed (fingers allow slip until a certain slip is reached)
		#  Inf    --> infinite slip allowed (the fingers grip loosely and allow the pencil to rotate)
		return new(z, slip_allowed, θ_arm, Δt, Inf, Inf)
	end
end

# setpoints for the arm height and rotation anngle
z_high = 1.2 * penci_length
z_low = penci_rad
θ_arm_1 =  pi / 2
θ_arm_2 = -pi / 2

v_gripper_ins = [
	### Pick up Pencil
	# lower gripper to pencil
	GripperIns(Δt=1.2, slip_allowed=-Inf, z=z_low,),
	# grip pencil while lowered
	GripperIns(Δt=0.25, slip_allowed=+Inf, z=z_low),
	# let pencil slip as gripper lifts up
	GripperIns(Δt=1.3, slip_allowed=+Inf),

	### Rotate Arm 90 degrees CCW
	# grip pencil hard to prevent future slip
	GripperIns(Δt=0.5, slip_allowed=0.0),
	# rotate pencil ccw 90 degrees
	GripperIns(Δt=2.0, slip_allowed=0.0, θ_arm=θ_arm_1),
	# loosen grip and let pencil slip
	GripperIns(Δt=1.5, slip_allowed=+Inf, θ_arm=θ_arm_1),

	### Rotate Arm 180 degrees CW
	# grip pencil hard to prevent future slip
	GripperIns(Δt=0.5, slip_allowed=0.0, θ_arm=θ_arm_1),
	# rotate pencil cw 180 degrees
	GripperIns(Δt=2.5, slip_allowed=0.0, θ_arm=θ_arm_2),

	### Control Slipping
	# let pencil slip a bit then stop
	GripperIns(Δt=2.5, slip_allowed=0.2, θ_arm=θ_arm_2),
	# let pencil slip a bit then stop
	GripperIns(Δt=2.5, slip_allowed=0.2, θ_arm=θ_arm_2),
	# let pencil rotate until it falls out of grasp
	GripperIns(Δt=Inf, slip_allowed=Inf, θ_arm=θ_arm_2),
]

function grip_control!(x::Vector{Float64}, m::MechanismScenario, t::Float64)
	joint_to_id(joint::Joint) = mech_scen.float.state.v[joint].indices[1]

	function calc_q̈(j::Joint, q_des::Float64, q̇_des=0.0; q̈_max::Float64=Inf)
		# Calculates the desired acceleration for a PD control law
		q = mech_scen.float.state.q[j][1]
		q̇ = mech_scen.float.state.v[j][1]
		q_err = q - q_des
		q̇_err = q̇ - 0.0
		q̈_des = - gain_d * q̇_err - gain_p * q_err
		q̈_des = clamp(q̈_des, -q̈_max, +q̈_max)
		return q̈_des, q_err, q̇_err
	end

	function get_grip_ins(t::Float64, rot_0::Float64)
		# get the correct item in v_gripper_ins based on how much time has passed
		grip_ins = v_gripper_ins[1]
		(grip_ins.t0 == Inf) && (grip_ins.t0 = t)
		if grip_ins.Δt < (t - grip_ins.t0)  # instruction expired
			popfirst!(v_gripper_ins)
			grip_ins = v_gripper_ins[1]
			grip_ins.t0 = t
			grip_ins.rot_0 = rot_0
		end
		return grip_ins
	end

	gain_p, gain_d = make_pd_gains(0.75, 1.0)  # 0.75s settling time, 1.0 damping ratio

	# the joint is in the named tuple associated with each joint
	j_tra_z = nt_tra_z.joint
	j_rev_y = nt_rev_y.joint
	j_pad_n = nt_pad_n.joint
	j_pad_p = nt_pad_p.joint
	j_penci = nt_penci.joint

	# Each joint has an ID
	id_tra_z = joint_to_id(j_tra_z)
	id_rev_y = joint_to_id(j_rev_y)
	id_pad_n = joint_to_id(j_pad_n)
	id_pad_p = joint_to_id(j_pad_p)

	id_actuate = vcat(id_tra_z, id_rev_y, id_pad_n, id_pad_p)  # only the DoF of the gripper are actuated

	# update dynamics quantaties
	tm = mech_scen.float
	state = tm.state
    copyto!(tm, x)
    H = tm.result.massmatrix
    mass_matrix!(H, state)
    dynamics_bias!(tm.result, state)
    configuration_derivative!(tm.result.q̇, state)
	H = mech_scen.float.result.massmatrix
	H4 = H[id_actuate, id_actuate]
	C4 = mech_scen.float.result.dynamicsbias[id_actuate]

	angle_penci = AngleAxis(rotation(transform_to_root(state, nt_penci.body))).theta
	grip_ins = get_grip_ins(t, angle_penci)

	q̈_tra_z, q_err_tra_z 	= calc_q̈(j_tra_z, grip_ins.z, q̈_max=1.0)
	q̈_rev_y, q_rev_y		= calc_q̈(j_rev_y, grip_ins.θ_arm)
	q̈_pad_n, _ 			= calc_q̈(j_pad_n, 0.0)
	q̈_pad_p, _ 			= calc_q̈(j_pad_p, 0.0)

	q̈ = [q̈_tra_z,	q̈_rev_y, q̈_pad_n, q̈_pad_p]
	m.τ_ext[id_actuate] .= H4 * q̈ + C4  # See Featherstone

	# use the slip_allowed field on the GripperIns to figure out how hard to pinch
	τ_soft = 0.3
	τ_hard = 2.5
	slip_allowed = grip_ins.slip_allowed
	if slip_allowed == -Inf
		τ_grip = 0.0
	elseif slip_allowed == 0.0
		τ_grip = τ_hard
	elseif slip_allowed == Inf
		τ_grip = τ_soft
	elseif isfinite(slip_allowed)
		Δ_ang_penci = abs(angle_penci - grip_ins.rot_0)
		τ_grip = ifelse(Δ_ang_penci < grip_ins.slip_allowed, τ_soft, τ_hard)
	end
	m.τ_ext[id_pad_n] .+= τ_grip  # just add grip forces onto the contant-ignorant forces
	m.τ_ext[id_pad_p] .+= τ_grip
end

### Create mechanism and temporary structure
mech_scen = MechanismScenario(discrete_controller=DiscreteControl(grip_control!, 0.01))
my_mechanism = mech_scen.float.state.mechanism

c_prop_compliant = ContactProperties(Ē=1.0e6)
i_prop_compliant = InertiaProperties(400.0)
i_prop_pad       = InertiaProperties(16000.0)
i_prop_rigid     = InertiaProperties(400.0, d=penci_rad)

eM_plane = eMesh_half_plane()
transform!(eM_plane, 0.6)

eM_pad_n = as_tet_eMesh(eMesh_sphere(pad_rad, 4))
transform!(eM_pad_n, SVector(0.0, -(pad_rad + penci_rad), 0.0))
transform!(eM_pad_n, SMatrix{3,3,Float64,9}(2.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0))

name_rev_y = "rev_y"
eM_pad_p = as_tet_eMesh(eMesh_sphere(pad_rad, 4))
transform!(eM_pad_p, SVector(0.0, +(pad_rad + penci_rad), 0.0))
transform!(eM_pad_p, SMatrix{3,3,Float64,9}(2.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0))
eM_rev_y = as_tri_eMesh(eMesh_box(pad_rad * SVector(1, 7, 1), pad_rad * SVector(0, -4, 8)))

tip_turn = 0.013
eM_penci = as_tri_eMesh(create_swept_mesh(f_swept_triv, [0.0, tip_turn, penci_length], [0.0, penci_rad, penci_rad], 12, true, rot_half=true))
transform!(eM_penci, SVector(0.0, 0.0, penci_rad))

joint_tra_z = Prismatic(SVector(0.0,  0.0, 1.0))
joint_rev_y = Revolute( SVector(0.0,  1.0, 0.0))
joint_n     = Prismatic(SVector(0.0, +1.0, 0.0))
joint_p     = Prismatic(SVector(0.0, -1.0, 0.0))

nt_plane = add_contact!(mech_scen, "plane", as_tet_eMesh(eM_plane), c_prop=c_prop_compliant)

mii_base = MeshInertiaInfo(1 .* ones(SMatrix{3,3,Float64,9}), zeros(SVector{3,Float64}), 1.0, NaN)
nt_tra_z = add_body_from_inertia!(my_mechanism, "tra_z", mii_base, joint=joint_tra_z)
nt_rev_y = add_body_from_inertia!(my_mechanism, name_rev_y, mii_base, joint=joint_rev_y, body=nt_tra_z.body)
nt_rev_y_id = add_contact!(mech_scen, name_rev_y, eM_rev_y, body=nt_rev_y.body)

nt_pad_n = add_body_contact!(mech_scen, "pad_n", eM_pad_n, c_prop=c_prop_compliant, i_prop=i_prop_pad, joint=joint_n, body=nt_rev_y.body)
nt_pad_p = add_body_contact!(mech_scen, "pad_p", eM_pad_p, c_prop=c_prop_compliant, i_prop=i_prop_pad, joint=joint_p, body=nt_rev_y.body)
nt_penci = add_body_contact!(mech_scen, "name", eM_penci, i_prop=i_prop_rigid)

if is_bristle
	c_ins_pad_n = add_friction_bristle!(mech_scen, nt_penci.id, nt_pad_n.id, μd=0.5, χ=0.6, k̄=k̄, magic=magic, τ=τ)
	c_ins_pad_p = add_friction_bristle!(mech_scen, nt_penci.id, nt_pad_p.id, μd=0.5, χ=0.6, k̄=k̄, magic=magic, τ=τ)
else
	c_ins_pad_n = add_friction_regularize!(mech_scen, nt_penci.id, nt_pad_n.id, μd=0.5, χ=0.6, v_tol=v_tol)
	c_ins_pad_p = add_friction_regularize!(mech_scen, nt_penci.id, nt_pad_p.id, μd=0.5, χ=0.6, v_tol=v_tol)
end
add_friction_regularize!(mech_scen, nt_penci.id, nt_plane.id, μd=0.5, χ=0.6, v_tol=v_tol)
add_friction_regularize!(mech_scen, nt_pad_n.id, nt_pad_p.id, μd=0.0, χ=0.6, v_tol=v_tol)

finalize!(mech_scen)

###############

q_tra_z_0 = 0.10
q_rev_y_0 = pi/4
set_configuration!(mech_scen, nt_tra_z.joint, [q_tra_z_0])
set_configuration!(mech_scen, nt_rev_y.joint, [q_rev_y_0])
set_state_spq!(mech_scen, nt_penci.joint, trans=SVector(penci_length * -0.8, 0.0, 0.0), rot=RotZ(-pi/2))

###############

x = get_state(mech_scen)

### Add meshcat visualizer
if !@isdefined vis
    vis = Visualizer()
    open(vis)
end

if true  ### visualization
	color_gray  = RGBA{Float32}(0.5, 0.5, 0.5, 1.0)
	color_red   = RGBA{Float32}(1.0, 0.0, 0.0, 1.0)
	color_green = RGBA{Float32}(0.0, 1.0, 0.0, 1.0)
	color_blue  = RGBA{Float32}(0.0, 0.0, 1.0, 0.5)
	mvis = MechanismVisualizer(my_mechanism, vis)

	### Add meshes to visualizer
	set_body_mesh_visual!(mvis, mech_scen, nt_pad_n.id,    color_gray)
	set_body_mesh_visual!(mvis, mech_scen, nt_pad_p.id,    color_gray)
	set_body_mesh_visual!(mvis, mech_scen, nt_rev_y_id.id, color_gray)
	set_body_mesh_visual!(mvis, mech_scen, nt_penci.id,    color_blue)

	set_mesh_visual!(     mvis, mech_scen, nt_plane.id,    color_green)
	set_configuration!(mech_scen, mvis)
end

### Run forward dynamics
rr = Radau_for_MechanismScenario(mech_scen)
rr.step.h_min = 1.0e-14
rr.step.h_max = 0.01
rr.step.h = rr.step.h_max

t_final   = 18.0
max_steps = 8200
println()

@time data_time, data_state = integrate_scenario_radau(rr, t_final=t_final, max_steps=max_steps)
println("t_sim: ", data_time[end])
println("simulated steps: ", length(data_time))
println("last time steps: ", (data_time[end] - data_time[end-10]) / 10)

### Move camera
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 7)
settransform!(vis["/Cameras/default"], Translation(0.0, 0.0, 0.30) ∘ LinearMap(RotZ(-π * 0.35)))

### Playback data
play_recorded_data(mvis, mech_scen, data_time, data_state, slowdown=1.0)


text:
using Lenny
using Test


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d

include("../LF_common/sparse_LF3d.jl")

using PyPlot
const plt = PyPlot

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  lin2xyz = LF.lin2xyz
  D2jl_x = LF.LFx.D2jl
  D2jl_y = LF.LFy.D2jl
  D2jl_z = LF.LFz.D2jl

  nnzc = Nx + Ny + Nz - 2
  NNZ = nnzc*Npoints

  rowval = zeros( Int64, NNZ )
  nzval = zeros( Float64, NNZ )
  colptr = zeros( Int64, Npoints + 1 )

  rowGbl_x_orig = zeros( Int64, Nx )
  rowGbl_y_orig = zeros( Int64, Ny )
  rowGbl_z_orig = zeros( Int64, Nz )

  rowGbl_x = zeros( Int64, Nx )
  rowGbl_y = zeros( Int64, Ny )
  rowGbl_z = zeros( Int64, Nz )

  for iz = 1:Nz
    rowGbl_z_orig[iz] = iz
  end

  rowGbl_y_orig[1] = 1
  for iy = 2:Ny
    rowGbl_y_orig[iy] = rowGbl_y_orig[iy-1] + Nz
  end

  rowGbl_x_orig[1] = 1
  for ix = 2:Nx
    rowGbl_x_orig[ix] = rowGbl_x_orig[ix-1] + Ny*Nz
  end

  # use @inbounds ?
  colGbl = 1
  idx_diag = 1

  ip = 0
  # rowLoc = 1:Nx, 1:Ny, 1:Nz for each gblCol
  for colGbl = 1:Npoints

    colLoc_x = ceil(Int,colGbl/(Ny*Nz))
    #println( rowGbl_x_orig + colGbl - 1 - (colLoc_x-1)*Nz*Ny )

    izz = ceil(Int,colGbl/Nz)
    #println( rowGbl_y_orig + colGbl - 1 - (izz-1)*Nz + (colLoc_x-1)*Ny*Nz )

    #println( rowGbl_z_orig + (izz-1)*Nz )

    yy = colGbl - (colLoc_x-1)*Ny*Nz
    colLoc_y = ceil( Int, yy/Nz )

    colLoc_z = colGbl - (izz-1)*Nz

    #println("coLLoc x, y, z:", [colLoc_x, colLoc_y, colLoc_z])
    #println("")

    # diagonal element
    ip = ip + 1
    rowval[ip] = colGbl
    nzval[ip] = D2jl_x[colLoc_x,colLoc_x] + D2jl_y[colLoc_y,colLoc_y] +  D2jl_z[colLoc_z,colLoc_z]
    #println("diagonal: ", rowval[ip])

    # non
    for ir = 1:Nx
      if ir != colLoc_x
        ip = ip + 1
        rowval[ip] = rowGbl_x_orig[ir] + colGbl - 1 - (colLoc_x-1)*Nz*Ny
        nzval[ip] = D2jl_x[ir,colLoc_x]
        #println("x: ", rowval[ip])
      end
    end

    for ir = 1:Ny
      if ir != colLoc_y
        ip = ip + 1
        rowval[ip] = rowGbl_y_orig[ir] + colGbl - 1 - (izz-1)*Nz + (colLoc_x-1)*Ny*Nz
        nzval[ip] = D2jl_y[ir,colLoc_y]
        #println("y: ", rowval[ip])
      end
    end

    for ir = 1:Nz
      if ir != colLoc_z
        ip = ip + 1
        rowval[ip] = rowGbl_z_orig[ir] + (izz-1)*Nz
        nzval[ip] = D2jl_z[ir,colLoc_z]
        #println("z: ", rowval[ip])
      end
    end

  end

  # colptr
  colptr[1] = 1
  for i = 2:Npoints+1
    colptr[i] = colptr[i-1] + nnzc
  end

  #println("ip = ", ip)
  #println("NNZ = ", NNZ)

  #println(colptr)
  #println(rowval)

  nabla2 = SparseMatrixCSC( Npoints, Npoints, colptr, rowval, nzval )

  nabla2_v2 = get_Laplacian3d_kron( LF )

  #println(nabla2)
  #println(nabla2_v2)

  #println( full(nabla2) - full(nabla2_v2) )

  #println( nabla2 - nabla2_v2 )

  srand(123)
  x = rand(Npoints)

  Lx = nabla2 * x
  Lx_v2 = nabla2_v2 * x

  println( sum(Lx-Lx_v2) )

  #plt.clf()
  #plt.spy(full(nabla2))
  #plt.title("nabla")
  #plt.savefig("nabla2.pdf")

  #plt.clf()
  #plt.spy(full(nabla2_v2))
  #plt.title("nabla_v2")
  #plt.savefig("nabla2_v2.pdf")

end

test_main( [30, 30, 30] )


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

# prototypes for space transforms and cutoffs
function transform end
function transform_d end
function fcut end
function fcut_d end

"""
`VecOrTup = Union{AbstractVector, Tuple}`
"""
const VecOrTup = Union{AbstractVector, Tuple}

abstract type ScalarBasis{T} <: IPBasis end

abstract type OneParticleBasis{T} <: IPBasis end

abstract type OnepBasisFcn end

# ------------------------------------------------------------
#  Abstract polynomial degree business

"""
`AbstractDegree` : object specifying a degree can be called via
`degree(D, arg)` or via `D(arg)`
"""
abstract type AbstractDegree end

(D::AbstractDegree)(args...) = degree(D, args...)

"""
`function degree(D::AbstractDegree, arg)` : compute some notion of degree of
the `arg` argument.
"""
function degree end

"""
interface functions for `OneParticleBasis`
"""
function add_into_A! end

"""
interface functions for `OneParticleBasis`
"""
function add_into_A_dA! end

"""
`function scaling(b, p)`:

a scaling factor for a basis functions ϕ, which gives a rought estimate on
   the magnitude of ∇ᵖϕ e.g.,
```
ϕ = r^n Ylm
```
has scaling factor `n^p + l^p`, though sharper estimates are also possible.
"""
function scaling end

using LinearAlgebra: Diagonal

diagonal_regulariser(basis; diff = 0) = Diagonal(scaling(basis, diff))

"""
every scalar basis must implement this
"""
function rand_radial end


text:
using PyPlot

function extract_planet(data,i)
eobs = zeros(0)
tobs = zeros(0)
sobs = zeros(0)
nobs = 0
for j=1:size(data,1)
  if data[j,1] == i
# Extract data:
    if data[j,4] != 0.0
      push!(eobs,data[j,2])
      push!(tobs,data[j,3])
      push!(sobs,data[j,4])
      nobs += 1
    end
  end
end
return eobs,tobs,sobs,nobs
end

function plot_ttv(data,elements,tt1)

# Now make some plots:
fig,axes = subplots(3,2,sharex="col",figsize=(12,8))
plabel = ["T1b","T1c","T1d","T1e","T1f","T1g","T1h"]
range = [1.0,1.0,1.0,2.0,10.0,8.0,3.0]
# Now make some plots:
#data = readdlm(fname,',')
for ip=1:nplanet
  if ip == 1
    ax = axes[ip]
  else
    ax = axes[ip-1]
  end
  coeff = zeros(Float64,2)
  # Plot each "Mode"
  eobs,tobs,sobs,nobs = extract_planet(data,ip)
  fn = zeros(2,nobs)
  fn[1,:] .= 1.0
  fn[2,:] .= eobs
  coeff,cov = regress(fn,tobs,sobs)
  println(ip," ",coeff)
  # Now plot TTVFast:
  fn = zeros(Float64,2,count1[ip+1])
  sig = ones(count1[ip+1])
  tti1 = tt1[ip+1,1:count1[ip+1]]
  tt_ref1 = zeros(count1[ip+1])
  for j=1:count1[ip+1]
    fn[1,j] = 1.0
    fn[2,j] = round(Int64,(tti1[j]-elements[ip+1,3])/elements[ip+1,2])
    epoch = round((tti1[j]-coeff[1])/coeff[2])
    tt_ref1[j] = coeff[1]+coeff[2]*epoch
  end
  coeff,cov = regress(fn,tti1,sig)
  tt_ref1 = coeff[1] .+coeff[2]*fn[2,:]
  ttv1 = (tti1 .-tt_ref1)*24*60.
  ax.plot(tti1,ttv1,label=plabel[ip],zorder=-32)
  ttv_obs = (tobs .- coeff[1] .- coeff[2]*eobs)*24*60.
#  ax[:plot](tobs,ttv_obs)
  ax.errorbar(tobs, ttv_obs, sobs*24*60.,fmt=".")
  if ip == 4 || ip == 7
    ax.set_xlabel(L"BJD$_\mathrm{TDB}$-2,450,000 [d]")
  end
  ax.set_ylabel("TTV [min]")
#  ax.axis([7260,8860,-range[ip],range[ip]])
#  if ip != 2 && ip != 3 && ip != 6
#    ax.legend(loc="upper right")
#  else
#    ax.legend(loc="lower right")
    ax.legend()
#  end
end
tight_layout()
return
end


text:
tofortran(val::Bool) =  val ? ".TRUE."  :  ".FALSE."
tofortran(val::Number) = val
tofortran(val::AbstractString) = "'$val'"

function dict2namelist(io::IO, dict::AbstractDict)
    println(io, "# this file has been generated automatically")
    for (groupname, group) in dict
        println(io, "&$groupname")
        if !(group === nothing)
            for (key, value) in group
                println(io, "  $key = $(tofortran(value)),")
            end
        end
        println(io, "/")
    end
end


text:
struct InputError <: Exception
    s::String
end
Base.showerror(io::IO, ex::InputError) = print(io, ex.s)


text:
#=
# Some things based on HOMG Tensor class
=#

# 2D routines
function tensor_IAX(A, x)
    N = size(A)[1];
    y = A * reshape(x, N, N);
    y = y[:];
    return y;
end

function tensor_AIX(A, x)
    N = size(A)[1];
    y = A * reshape(x, N, N)';
    y = y';
    y = y[:];
    return y;
end

# 3D routines
function tensor_IIAX(A, x)
    N = size(A)[1];
    y = A * reshape(x, N, N*N);
    y = y[:];
    return y;
end

function tensor_IAIX(A, x)
    N = size(A)[1];
    q = reshape(x, N, N, N);
    y = zeros(N,N,N);
    for i=1:N
        y[i,:,:] = A * q[i,:,:];
    end
    y = y[:];
    return y;
end

function tensor_AIIX(A, x)
    N = size(A)[1];
    y = reshape(x, N*N, N) * A';
    y = y[:];
    return y;
end

function tensor_grad2(A, x)
    dx = tensor_IAX(A, x);
    dy = tensor_AIX(A, x);
    return (dx, dy);
end

function tensor_grad3(A, x)
    dx = tensor_IIAX(A, x);
    dy = tensor_IAIX(A, x);
    dz = tensor_AIIX(A, x);
    return (dx, dy, dz);
end


text:
module Registrator

using UUIDs, LibGit2

import Base: PkgId

include("slack.jl")
include("regedit/RegEdit.jl")
include("pull_request.jl")
include("Messaging.jl")
include("RegService.jl")
include("commentbot/CommentBot.jl")
include("webui/WebUI.jl")

end # module